Решение системы линейных уравнений с использованием дэгсва Лапака
Я хочу решить систему линейных уравнений, используя пакет Лапака в C ++.
Я планирую реализовать это как этот используя процедуры из Вот, а именно дгесв.
Это мой код:
unsigned short int n=profentries.size();
double **A, *b, *x;
A= new double*[n];
b=new double[n];
x=new double[2];
// Generate the matrices for my equation
for(int i = 0; i < n; ++i)
{
A[i] = new double[2];
}
for(int i=0; i<n; i++)
{
A[i][0]=5;
A[i][1]=1;
b[i]=3;
}
std::cout<< "printing matrix A:"<<std::endl;
for(int i=0; i<n; i++)
{
std::cout<< A[i][0]<<" " <<A[i][1]<<std::endl;
}
// Call the LAPACK solver
x = new double[n];//probably not necessary
dgesv(A, b, 2, x); //wrong result for x!
std::cout<< "printing vector x:"<<std::endl;
/*prints
3
3
but that's wrong, the solution is (0.6, 0)!
*/
for(int i=0; i<2; i++)
{
std::cout<< x[i]<<std::endl;
}
У меня есть следующая проблема:
Как так получается, что dgesv вычисляет вектор x с элементами {3, 3}? Решение должно быть {0.6, 0} (проверил его с помощью matlab).
Привет
Изменить: dgesv может работать для квадратных матриц. Мое решение ниже показывает, как решать переопределенные системы с помощью dgels.
Решение
Лапак предполагает, что матрица хранится в основном формате столбца. Он не будет работать с форматом указатель-указатель вашей матрицы A. Все элементы матрицы должны постоянно храниться в памяти.
Попробуйте объявить A так:
double *A;
A = new double[2*n];
И для цикла:
for(int i=0; i<n; i++)
{
A[i+0*n]=5;
A[i+1*n]=1;
b[i]=3;
}
Другие решения
Проблема заключалась в том, что на этот раз произошло смешение мажор-строка-мажор, а также тот факт, что dgesv явно не подходит для неквадратных матриц.
Переопределенная система линейных уравнений может быть решена в стиле наименьших квадратов с использованием dgels (#include «lapacke.h»).
Примечательно, что для меня это неочевидно на первый взгляд: вектор решения (обычно обозначаемый x) затем сохраняется в b.
Мой пример системы состоит из матрицы A, содержащей значения 1-10 в первом столбце, и вектора b со значением i ^ 2 для {i | 1<= я<= 10}:
double a[10][2]={1,1,2,1,3,1,4,1,5,1,6,1,7,1,8,1,9,1,10,1};
double b[10][1]={1,4,9,16,25,36,49,64,81,100};
lapack_int info,m,n,lda,ldb,nrhs;
int i,j;
// description here: http://www.netlib.org/lapack/double/dgesv.f
m = 10;
n = 2;
nrhs = 1;
lda = 2;
ldb = 1;
for(int i=0; i<m; i++)
{
for(int k=0; k<lda; k++)
{
std::cout << a[i][k]<<" ";
}
std::cout <<"" << std::endl;
}
std::cout<< "printing vector b:"<<std::endl;
for(int i=0; i<m; i++)
{
std::cout<< *b[i]<<std::endl;
}
std::cout<< "\nStarting calculation..."<<std::endl;
info = LAPACKE_dgels(LAPACK_ROW_MAJOR,'N',m,n,nrhs,*a,lda,*b,ldb);
std::cout<< "Done."<<std::endl;
std::cout<< " "<<std::endl;
std::cout<< "printing vector b:"<<std::endl;
for(int i=0; i<2; i++)
{
std::cout<< *b[i]<<std::endl;
}
std::cout<< "Used values: "<< m << ", " << n << ", " << nrhs << ", " << lda << ", " << ldb << std::endl;
detector