Solving a system of linear equations using Lapack's dgesv

Question

I want to solve a linear equation system using the Lapack package in C++. I plan to implement it like this using the routines from here, namely dgesv.

This is my code:

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;
}

I have the following problem:

How can it be that dgesv calculates a vector x with the elements {3, 3}? The solution should be {0.6, 0} (checked it with matlab).

Greetings

Edit: dgesv might work for square matrices. My solution below shows how to solve overdetermined systems with dgels.

Argh, those interfaces you're using are just horrible. Please do yourself a favor and use the official C bindings which come with recent lapack versions: http://www.netlib.org/lapack/lapacke.html — janneb, Oct 09 '12 at 11:21
The corresponding function seems to be LAPACKE_dgels, but that does neither modify a passed-by-reference vector x nor return a vector x, do I see that correctly? So, how to get the solution x? — bogus, Oct 09 '12 at 11:25
I'd say LAPACKE_dgesv is the wrapper for the Fortran routine DGESV, what makes you think it should be LAPACKE_dgels? Per the documentation, on output the solution x is stored into B. — janneb, Oct 09 '12 at 11:38
@janneb I adopted all the suggestions, but what I don't really get is how to achieve the solution vector x whcih contains the "unknowns". My A has 19 rows/2 cols, b has 19 rows/1 col, so x should have 2 rows/1 col. The values don't seem to be in b! — bogus, Oct 09 '12 at 13:46
With those dimensions your system is overdetermined, and you cannot use dgesv (that is, dgesv requires A to be square). Use dgels instead, if you want an approximate solution in the least-squares sense (that is IIRC what the matlab '\' operator does for non-square A). — janneb, Oct 09 '12 at 14:04

score 1 · Answer 1 · answered Oct 09 '12 at 10:28

1

Lapack assumes the matrix is stored in column major format. It will not work with the pointer-to-pointer format of your matrix A. All elements of the matrix need to be stored continuously in memory.

Try declaring A as such:

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

And the for-loop:

for(int i=0; i<n; i++)
{
    A[i+0*n]=5;
    A[i+1*n]=1;
    b[i]=3;
}

answered Oct 09 '12 at 10:28

digitalvision

2,552
19
23

Ok, I wasn't aware of Lapack expecting column major format. Having adapted your remark, I get the following message: Error: argument of type "double *" is incompatible with parameter of type "double **". Does dgesv in fact require a pointer-to-pointer (i.e. a twodimensional input)? – bogus Oct 09 '12 at 11:04
I assumed you were using a regular lapack C wrapper. It seems that the code you are using is in fact doing the conversion for you. See the dgesv_ctof() function in your dgesv.h – digitalvision Oct 09 '12 at 11:30
The link provided by @janneb states that Lapack, depending on a parameter, can accept both row and column major format. EDIT: This refers to my interface, not the official lapack c binding! – bogus Oct 09 '12 at 11:31
There are so many different Lapack bindings. Regarding row major, the official ones state: "Note that using row-major ordering may require more memory and time than column-major ordering, because the routine must transpose the row-major order to the column-major order required by the underlying LAPACK routine." – digitalvision Oct 09 '12 at 11:34

bogus · Accepted Answer · 2012-10-10T11:32:02.963

The problem was for once the mixup of row-major/column major but also the fact that apparently dgesv is not suitable for non-square matrices.

An overdetermined system of linear equations can be solved in the least-squares style, using dgels (#include "lapacke.h").

Noteworthy, to me not apparent at first sight: The solution vector (usually denoted by x) is then stored in b.

My example system consists of a matrix A holding the values 1-10 in its first column and a vector b with the value i^2 for {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;

siimurik · Answer 3 · 2022-12-10T15:18:09.793

This what the code would look like in C:

#include <stdio.h>
#include <stdlib.h>
#include <lapacke.h>
#include <time.h>

#define N 3
#define NRHS 1
#define LDA N
#define LDB NRHS

int main (){
    int n = N, nrhs = NRHS, lda = LDA, ldb = LDB, info;
    int ipiv[N];

    // Define a system of linear equations
    // 3x + 2y - z = 2
    // 2x - 2y + 4z = -3
    // -x + y/2 - z = z/2
    double M[LDA * N] = {
        3.0,  2.0, -1.0, 
        2.0, -2.0,  4.0, 
        -1.0,  0.5, -1.0
    };

    double V[LDB * N] = {
        2.0, 
        -3.0, 
        0.5
    };

    // Record the start time
    clock_t start = clock();

    // Call LAPACKE_dgesv
    info = LAPACKE_dgesv(LAPACK_ROW_MAJOR, n, nrhs, M, lda, ipiv, V, ldb);

    // Record the end time
    clock_t end = clock();

    // Calculate the elapsed time in seconds
    double elapsed_time = (end - start) / (double)CLOCKS_PER_SEC;

    printf("Solution to the system of linear equations:");
    for (int i = 0; i < N; i++) {
        printf("\nx[%d] = %f", i, V[i]);
    }
    printf("\nElapsed time: %e s\n", elapsed_time);
    printf("\n");
    return 0;
}

If you haven't installed the LAPACK library, this is how you can do it on Ubuntu:

sudo apt-get install liblapack-dev

Then you can compile your code with this command:

gcc laplin.c -o lin -llapacke

Solving a system of linear equations using Lapack's dgesv

3 Answers3