(new pb) : invert a large matrix with lapack

[cecile]

Hi,

I am a new user of lapack, and I tried as a simple first step to invert a matrix. The functions seems to work correctly, but when I check the result with the result that matlab gives me, it doesn't match.

Here is what I tried :
this is my initial matrix (the numbers are floats)

Code: Select all

mat =

    0.1000    2.2000    2.0000    4.2000    1.1000
    0.3000    3.1000         0    8.9000    0.2000
    7.1000    2.0000    1.0000    8.1000    3.5400
    0.5000    4.7000    0.8000    4.1000    1.2000
    2.2000    7.1000    6.4000    2.3000    1.3000

then i call the lapack routines :

Code: Select all

integer N1=5;
integer INFO = 0; 
float *A;
integer  *IPIV;
float *work;
sgetrf_(&N1,&N1,A,&N1,IPIV,&INFO);
sgetri_(&N1,A,&N1,IPIV,work,&lwork,&INFO);

where A is my previous matrix written as a vector column by column.

Here is the result I get from lapack :

Code: Select all

inv =

         0   -0.0845   -0.0322   -0.0756    0.3057
         0    0.6369   -0.0986   -0.1324   -0.1767
         0    0.0342    0.1292    0.0192   -0.1093
         0   -0.2103   -0.2495    0.3188    0.1274
         0         0         0         0    1.3000

And here is the correct answer, given by matlab :

Code: Select all

ans =

   -0.3986    0.1077    0.1124   -0.1158    0.1215
   -0.2446    0.0111   -0.0304    0.2592    0.0488
    0.2405    0.0085   -0.0206   -0.2919    0.1208
    0.0834    0.1159    0.0053   -0.0950   -0.0150
    0.6789   -0.4899    0.0679    0.3858   -0.2710

I don't know if this forum is a good place to post my problem, but so far it's the only forum about Lapack that I found.
Could you please help me finding where my error comes from ?

Thank you.

[cottrell]

You don't give the details of your code so we can't tell exactly
what you have done wrong, but the following C program prints
the solution as given by matlab. (For convenience, I'm transposing
on both input and output.)

Code: Select all

#include <stdio.h>

int main (void)
{
    int N1 = 5;
    int INFO = 0;
    int IPIV[5];
    int lwork = 5;
    float work[5];
    float A[25] = {
        0.1, 2.2, 2.0, 4.2, 1.1,
        0.3, 3.1, 0.0, 8.9, 0.2,
        7.1, 2.0, 1.0, 8.1, 3.54,
        0.5, 4.7, 0.8, 4.1, 1.2,
        2.2, 7.1, 6.4, 2.3, 1.3
    };
    int i, j, k;

    sgetrf_(&N1, &N1, A, &N1, IPIV, &INFO);
    sgetri_(&N1, A, &N1, IPIV, work, &lwork, &INFO);

    k = 0;
    for (i=0; i<N1; i++) {
        for (j=0; j<N1; j++) {
            printf("%10.4f", A[k++]);
        }
        putchar('\n');
    }
   
    return 0;
}


[cecile]

Thank you for your help, this also works for me.
I think the problem is in the way I define A, since I define it as a float* and then read every value from a file. (I do this because that's what I will have to do later).
Here is how I define A :

Code: Select all

   float *A;
   if ((A=(float*)malloc( N*N*sizeof(float) )) == NULL)
   {
      printf("memory");
      exit(1);
   }
   LapackRead(A,N,N,file1);

And the LapackRead function is like this :

Code: Select all

void LapackRead(float *data,int N, int M, char *filename)
{
  int x,y;
  FILE *fp;
  char *err_msg=NULL;

  if ((fp=fopen(filename,"r"))== NULL)
    {
      err_msg = strerror(errno);
      fprintf(stderr, "Error opening file in LapackRead: %s\n", err_msg);
      exit(1);
    }

  for(y=1;y<=M*N;y++){
     fscanf(fp,"%f", &data[y]);
  }

  fclose(fp);
}

Since when I write to a file, I write column by column, I am also reading column by column here. When I print the output of the LapackRead function, I get the correct matrix.

Is this where my problem comes from ? How can I solve it ?
Thank you.

[cecile]

Ok problem solved !

[cecile]

Hi again !

As I said before, inverting a small matrix works, but now I am trying to invert a bigger matrix (1024*1024) and I get results that are not close at all to what I should get.
Here is how I do it :

Code: Select all

float *Kg;
   Kg = LapackCov(g,1,N*N);
   
   integer N1=N*N;
   integer INFO = 0;
   integer  *IPIV;
   float *work;
   char *memString = "\nBuy more memory!\n";
   
   if ((IPIV=(integer*)malloc( N1*sizeof(int) )) == NULL)
   {
      printf(memString);
      exit(1);
   }

   if ((work=(float*)malloc( N1*sizeof(float) )) == NULL)
   {
      printf(memString);
      exit(1);
   }
   
   printf("Calling Lapack\n");

   sgetrf_(&N1,&N1,Kg,&N1,IPIV,&INFO);

   integer lwork = -1;
   sgetri_(&N1,Kg,&N1,IPIV,work,&lwork,&INFO);
   if (INFO != 0 || work[0] <= 0.0) {
        fprintf(stderr, "dgetri (workspace): info = %d\n", INFO);
        exit(1);
        }

   lwork = (int) work[0];

    realloc(work, lwork*sizeof(float));
    if (work == NULL) {
        fputs("out of memory\n", stderr);
        exit(1);
    } 

   sgetri_(&N1,Kg,&N1,IPIV,work,&lwork,&INFO);

   printf("Lapack Done\n");

   /* Check the output flag. */
   if (INFO==0)
        printf("\nComputation successful\n");
   else
   {
      printf("\nThe tatterdemalions have run amok on our code.\n");
      //free(KgLapack);
      free(IPIV);
      exit(1);
   }

   free(IPIV);
   free(work);

LapackCov is a home function that computes a covariance matrix and put it in the lapack format.
This function returns values that match the ones I get with Matlab.

When I run my code and set breakpoints to track different values here is what I see :
- before calling sgetrf, Kg has a "good" value : 0.21770304 float
- before calling sgetri for the first time, Kg still has a "good" value : 4.1318846 float
- before calling sgetri for the second time, Kg still has a "good" value : 4.1318846 float, nothing has changed, which is what we expected since we called sgetri for the first time with a lwork = -1.
- after calling sgetri for the second time, Kg has a total nonsense value : 18053338. float

To me this looks like a memory error : some of my variables are getting overwritten, and I don't know how I can prevent this from happening !
How can I allocate enough memory space for my computation so that none of my variables get corrupted ?

I would greatly appreciate any help on this problem.
Thanks.

[cottrell]

Size doesn't matter ;-)

Seriously, Lapack is well tested on huge matrices. You don't post enough of
your own code for a third party to detect the error it contains. This forum
is primarily for discussing questions relating to Lapack itself. In general, we
cannot diagnose bugs in programs that happen to call Lapack functions. We
may be able to do so if you compose and post a minimal but complete test
case -- otherwise, your crystal ball is as good as ours!

[tmatematikas]

Hello,
This might be an old thread, but I have one idea and a question. The results of LAPACK and Mathcad may be different, because Mathcad uses double precision arithmetic and LAPACK uses single (eg. sgetrf like in this case) or double (dgetrf).
And here is my question. My sample program works with dgetrf, but do not with sgetrf (compiler does not find it). Why is that? installation problem, link issues,... ? And how to use qgetrf (which I assume would be quadruple precision)? I use -llapack, -lblas in DevC++ compiler options. Sorry for slightly ot.