(new pb) : invert a large matrix with lapack

Open discussion regarding features, bugs, issues, vendors, etc.

(new pb) : invert a large matrix with lapack

Postby cecile » Thu Oct 14, 2010 6:48 pm

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.
Last edited by cecile on Mon Oct 25, 2010 10:43 pm, edited 1 time in total.
cecile
 
Posts: 4
Joined: Thu Oct 14, 2010 6:31 pm

Re: invert a matrix with lapack

Postby cottrell » Fri Oct 15, 2010 12:06 pm

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;
}
cottrell
 
Posts: 68
Joined: Thu Jan 15, 2009 1:40 pm

Re: invert a matrix with lapack

Postby cecile » Fri Oct 15, 2010 1:27 pm

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
 
Posts: 4
Joined: Thu Oct 14, 2010 6:31 pm

Re: invert a matrix with lapack

Postby cecile » Fri Oct 15, 2010 8:13 pm

Ok problem solved !
cecile
 
Posts: 4
Joined: Thu Oct 14, 2010 6:31 pm

Re: (new pb) : invert a large matrix with lapack

Postby cecile » Mon Oct 25, 2010 11:20 pm

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.
cecile
 
Posts: 4
Joined: Thu Oct 14, 2010 6:31 pm

Re: (new pb) : invert a large matrix with lapack

Postby cottrell » Sun Nov 07, 2010 8:52 pm

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!
cottrell
 
Posts: 68
Joined: Thu Jan 15, 2009 1:40 pm

Re: (new pb) : invert a large matrix with lapack

Postby tmatematikas » Wed May 22, 2013 6:58 am

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.
tmatematikas
 
Posts: 1
Joined: Wed May 22, 2013 6:48 am


Return to User Discussion

Who is online

Users browsing this forum: Yahoo [Bot] and 4 guests