Решение системы линейных уравнений с помощью dgesv Лапака

Я хочу решить систему линейных уравнений, используя пакет Lapack на C++. Я планирую реализовать его, как этим, используя подпрограммы из здесь, а именно dgesv.

Это мой код:

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.


c++ linear-algebra equation-solving lapack
person bogus    schedule 09.10.2012    source источник
comment
Аргх, те интерфейсы, которые ты используешь, просто ужасны. Пожалуйста, сделайте себе одолжение и используйте официальные привязки C, которые поставляются с последними версиями lapack: netlib.org/ лапак/лапак.html   -  person janneb    schedule 09.10.2012
comment
Соответствующая функция выглядит как LAPACKE_dgels, но она не изменяет переданный по ссылке вектор x и не возвращает вектор x, правильно ли я понимаю? Итак, как получить решение x?   -  person bogus    schedule 09.10.2012
comment
Я бы сказал, что LAPACKE_dgesv — это оболочка для подпрограммы Фортрана DGESV, почему вы думаете, что это должен быть LAPACKE_dgels? Согласно документации, на выходе решение x сохраняется в B.   -  person janneb    schedule 09.10.2012
comment
@janneb Я принял все предложения, но чего я действительно не понимаю, так это того, как получить вектор решения x, который содержит неизвестные. Мой A имеет 19 строк/2 столбца, b имеет 19 строк/1 столбец, поэтому x должен иметь 2 строки/1 столбец. Значения, кажется, не в b!   -  person bogus    schedule 09.10.2012
comment
С этими размерностями ваша система переопределена, и вы не можете использовать dgesv (то есть dgesv требует, чтобы A был квадратным). Вместо этого используйте dgels, если вам нужно приблизительное решение в смысле наименьших квадратов (это IIRC, что делает оператор matlab '\' для неквадрата A).   -  person janneb    schedule 09.10.2012

Ответы (2)


arrow_upward
1
arrow_downward

Лапак предполагает, что матрица хранится в основном формате столбца. Он не будет работать с форматом указателя на указатель вашей матрицы A. Все элементы матрицы должны постоянно храниться в памяти.

Попробуйте объявить A как таковой:

double *A;
A = new double[2*n];

И цикл for:

for(int i=0; i<n; i++)
{
    A[i+0*n]=5;
    A[i+1*n]=1;
    b[i]=3;
}
person digitalvision    schedule 09.10.2012
comment
Хорошо, я не знал, что Лапак ожидает основной формат столбца. Адаптировав ваше замечание, получаю следующее сообщение: Ошибка: аргумент типа double * несовместим с параметром типа double **. Действительно ли для dgesv требуется указатель на указатель (т.е. двумерный ввод)? - person bogus; 09.10.2012
comment
Я предположил, что вы используете обычную оболочку lapack C. Похоже, код, который вы используете, на самом деле выполняет преобразование за вас. См. функцию dgesv_ctof() в файле dgesv.h. - person digitalvision; 09.10.2012
comment
В ссылке, предоставленной @janneb, говорится, что Lapack, в зависимости от параметра, может принимать как основной формат строки, так и столбца. РЕДАКТИРОВАТЬ: Это относится к моему интерфейсу, а не к официальной привязке lapack c! - person bogus; 09.10.2012
comment
Существует так много разных креплений Lapack. Что касается основного порядка строк, официальные утверждают: Обратите внимание, что использование упорядоченного упорядочения строк может потребовать больше памяти и времени, чем упорядочение столбцов, потому что подпрограмма должна транспонировать основной порядок строк в порядок столбцов, требуемый основной процедурой LAPACK. . - person digitalvision; 09.10.2012

arrow_upward
1
arrow_downward

На этот раз проблема заключалась в смешении основных строк и основных столбцов, а также в том, что, по-видимому, dgesv не подходит для неквадратных матриц.

Переопределенная система линейных уравнений может быть решена методом наименьших квадратов с использованием dgels (#include "lapacke.h").

Примечательно, что для меня это не очевидно с первого взгляда: вектор решения (обычно обозначаемый x) затем сохраняется в b.

Моя примерная система состоит из матрицы A, содержащей значения 1-10 в первом столбце, и вектора b со значением i^2 для {i|1‹=i‹=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;
person bogus    schedule 10.10.2012