LAPACKE functioning oddly

[pkassebaum]

I'm new to LAPACKE, so I wanted to run a simple test program that would compute the inverse of a 2 by 2 matrix. My compiler makes no complaints, but the results are obviously wrong. In trying to invert the matrix

a11 = 1.0, a12 = -1.0,
a21 = 1.0, a22 = 1.0,

my program returns a successful run (info = 0) but spits out the supposed inverse

a11 = 1382417208487871530309079876829361519920890827655236399202446373087004409281425975828759961223649754730531017408099857515546836392400417046590332715399634071080623695707883976652716068383393396038903725080647755086048184917192852687494574532544186355366261864363207133093848253625865070671519163486482464768.0
a12 = 1.00
a21 = 0.00
a22 = 1.00.

Any ideas on what I'm doing so horribly wrong? Thanks in advance! Here's my code:

Code: Select all

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

void printMatrix( int lda, int n, double *a );

int main( void )
{
   double *a;
   int n, lda, info, *ipiv;

   lda  = 2;
   n    = 2;
   a    = calloc( lda*n, sizeof(double) );
   ipiv = calloc( n,     sizeof(int) );

   a[0] = 1.0;   a[1] =-1.0;
   a[2] = 1.0;   a[3] = 1.0;

   printMatrix( lda, n, a );
   info = LAPACKE_dgetri( LAPACK_ROW_MAJOR, n, a, lda, ipiv);
   printf("info = %i\n", info);
   printMatrix( lda, n, a );

   free(a);
   free(ipiv);

   return 0;
}
//==========================================================

void printMatrix( int lda, int n, double *a )
{
   int i, j;
   for(i = 0; i < lda; i++) {
      for(j = 0; j < n; j++) {
         printf("a[%i][%i] = %.2f",
               i, j, a[i*n + j]
         );
         putchar('\t');
      }
      putchar('\n');
   }
}


[admin]

Hi,
you just need to call LAPACK_dgetrf before calling LAPACK_dgetri.
Try adding this line before the call to LAPACK_dgetri

Code: Select all

   info = LAPACKE_dgetrf( LAPACK_COL_MAJOR, n, n, a, lda, ipiv);


LAPACK_dgetri needs the factorization given by LAPACK_dgetrf and the resulting IPIV vector.

See http://www.netlib.org/lapack/explore-ht ... ource.html for more info
Julie

[pkassebaum]

It works! Thank you so much!