## Inversing a Dynamic Matrix using dgetri_ and dgetrf_

Post here if you have a question about LAPACK or ScaLAPACK algorithm or data format

### Inversing a Dynamic Matrix using dgetri_ and dgetrf_

I am currently trying to inverse a dynamic square (n by n) matrix. The code works for a 2*2 matrix if i specify the elements. However when i tried to make the code work for a dynamic matrix which prompts the user to enter the elements, it crashes. I do not know how dgetri_really works so its proving difficult to solve this problem. Help please!!! Here is my code so far.

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

void dgetri_ ( int* n, double** A, int* lda,int* ipiv, double* work, int* lwork, int* info );

void dgetrf_ ( int* m, int* n, double** a, int* lda,int* ipiv, int* info );

// this allocates memory to the dynamic matrix//
double** allomat (int row, int col)
{

double **matrix;
int i;
matrix = (double **) malloc (row*sizeof(double *));
if (!matrix) return NULL;
matrix[0] = (double *) malloc (row*col*sizeof(double));
if (!matrix[0])
{
free(matrix);
return NULL;
}
for (i = 1; i < row; i++)
matrix[i] = matrix[i-1] + col;
return matrix;
}

// this part of the code carries out matrix inversion//
void inverse (double** A, int n) {
int i,n,j;
int lda;
int *ipiv;
int lwork;
int info;
double* work;

lda=n;
ipiv = (int*)malloc(20*n*sizeof(int));
work = (double*)malloc(100*n*(sizeof(double)));

dgetrf_ ( &n, &n, A, &lda, ipiv, &info );
printf ("Matrix A after dgetrf\n");
for (i=0;i<n;i++) {
for (j=0;j<n;j++)
{
printf ("A[%d][%d]:%f\t", i,j,A[i][j]);
}
}
printf ("\ninfo for dgetrf_:%d\n\n", info);

lda=n;
lwork=100*n; //magic number - dimension of work

dgetri_ (&n, A, &lda, ipiv, work, &lwork, &info);
printf ("Matrix A after dgetri\n");
for (i=0;i<n;i++) {
for (j=0;j<n;j++){
printf ("A[%d][%d]:%f\t", i,j,A[i][j]);
}
}
printf ("\ninfo for dgetri_:%d\n\n", info);
free (ipiv);
free (work);

}

int main ()
{
int i,j,n;
double **A;
double **Ainv;
printf
scanf("%d",&n);
A = allomat(n,n);
Ainv = allomat(n,n);
if (!A||!Ainv)
{printf("error in malloc\n");}
for (i=0;i < n; i++){
for(j=0;j<n;j++)
{ scanf("%lf",&A[i][j]);

}
}
for (i=0;i < 2; i++){
for(j=0;j<2;j++)
{ Ainv[i][j]=A[i][j];
}
}
printf("crap\n");
for (i=0;i < 2; i++){
for(j=0;j<2;j++)
{ printf("%lf\n", Ainv[i][j]);
}
}

inverse(A,2);
free (A);
free(Ainv);
return 0;

}
David12

Posts: 3
Joined: Mon Mar 07, 2011 12:33 pm

### Re: Inversing a Dynamic Matrix using dgetri_ and dgetrf_

[ I just had a quick look at your code, there might be more problems.

1) the main problem is your matrix A, your declare it as being a ( double ** ), and access the element (i,j) with A[i][j], that's not the way to do it,

2) you need to declare A as ( double * ), allocate it of size n*n, and access element (i,j) with A[ i + j*n ] (this is called column major format )

You can also consider using lapacke to have a native C interface to lapack (this can wait)

for the workspace, 100*n, sounds good for the moment, but please consider using a workspace query:
1) allocate work ( double *) of size 1,
2) call lapack dgetri with the lwork argument set to -1
3) read the value in work[0], then cast it in lwork: lwork = (int) work[0];
4) free work
5) reallocate work of size lwork
now do the "real" call to dgetri, work has the good size, lwork the good value

j
Julien Langou

Posts: 824
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

### Re: Inversing a Dynamic Matrix using dgetri_ and dgetrf_

Following Julien's advie to use LAPACKE , here is a quick code with LAPACKE (the new Standard C language APIs for LAPACK)
LAPACKE is available at http://netlib.org/lapack/#_standard_c_l ... for_lapack

Code: Select all
`#include<stdio.h>#include<stdlib.h>#include "lapacke.h"#include "cblas.h"int main (int argc, char **argv){ lapack_int n,lda,info;int i;lapack_int *ipiv;double alpha,beta,norm;double *A,*Acopy,*MatId;n=1000;lda=1000;// Allocating Memory   A = (double *)malloc(lda*n*sizeof(double)) ;if (A==NULL){ printf("error of memory allocation\n"); exit(0); }Acopy = (double *)malloc(lda*n*sizeof(double)) ;if (Acopy==NULL){ printf("error of memory allocation\n"); exit(0); }MatId = (double *)malloc(lda*n*sizeof(double)) ;if (MatId==NULL){ printf("error of memory allocation\n"); exit(0); }ipiv = (lapack_int *)malloc(n*sizeof(lapack_int)) ;if (ipiv==NULL){ printf("error of memory allocation\n"); exit(0); }// Generating A as a Random Matrix and MatId as the Identity Matrixfor(i=0;i<lda*n;i++)   {     A[i] = ((double) rand()) / ((double) RAND_MAX) - 0.5;     Acopy[i]=A[i];     MatId[i]=0.0;    }for(i=0;i<lda;i++)     MatId[i*n]=1.0;     // Call to dgetrfinfo=LAPACKE_dgetrf(LAPACK_COL_MAJOR, n, n, A, lda, ipiv );printf("\tDGETRF [n=%d] INFO=%d \n",n,info);info=LAPACKE_dgetri(LAPACK_COL_MAJOR, n, A, lda, ipiv); // Call to dgetriprintf("\tDGETRI [n=%d] INFO=%d \n",n,info);// Checkprintf("CHECK: Computing the norm OF A*Ainv-MatID (Should be zero)\n");alpha = 1.0; beta = -1.0;cblas_dgemm ( CblasColMajor, CblasNoTrans, CblasNoTrans, n, n, n, alpha, Acopy, lda, A, lda, beta, MatId, lda);norm = cblas_dnrm2(n,MatId,1);printf("\tNORM=%f\n",norm);free(A);free(Acopy);free(MatId);free(ipiv);return 0;}`

This code is also using the cblas for checking purpose, just to make sure that the inversion is correct.
CBLAS is available at http://www.netlib.org/blas/blast-forum/cblas.tgz
gcc -c prog.c
gfortran prog.o -o prog.exe lapacke.a cblas.a lapack.a blas.a

Posts: 609
Joined: Wed Dec 08, 2004 7:07 pm

### Re: Inversing a Dynamic Matrix using dgetri_ and dgetrf_

Julien Langou wrote:[ I just had a quick look at your code, there might be more problems.

1) the main problem is your matrix A, your declare it as being a ( double ** ), and access the element (i,j) with A[i][j], that's not the way to do it,

2) you need to declare A as ( double * ), allocate it of size n*n, and access element (i,j) with A[ i + j*n ] (this is called column major format )

You can also consider using lapacke to have a native C interface to lapack (this can wait)

for the workspace, 100*n, sounds good for the moment, but please consider using a workspace query:
1) allocate work ( double *) of size 1,
2) call lapack dgetri with the lwork argument set to -1
3) read the value in work[0], then cast it in lwork: lwork = (int) work[0];
4) free work
5) reallocate work of size lwork
now do the "real" call to dgetri, work has the good size, lwork the good value

j

I have made the changes you asked me to do (making the matrix a column major array and changing the prototype to a fortran based one). And this has solved the problem. Thank you very much. I am also gonna look into lapacke. Thanks again.
The new code is below if any one ever needs it, its simple but does the job.
Code: Select all
`#include<stdio.h>#include<stdlib.h>#include "f2c.h"#include "clapack.h"// void sgetri_ ( int* n, const void* A, int* lda,  //  int* ipiv, const void* work, int* lwork, int* info );  //void sgetrf_ ( int* m, int* n, const void* a, int* lda,    //int* ipiv, int* info ); // double A[2][2] = {{1.0,3.0},{2.0,4.0}};  double** allomat (int row, int col){double **matrix;int i;matrix = (double **) malloc (row*sizeof(double *));if (!matrix) return NULL;matrix[0] = (double *) malloc (row*col*sizeof(double));if (!matrix[0]){free(matrix); return NULL; }for (i = 1; i < row; i++)matrix[i] = matrix[i-1] + col;return matrix;}  void printmat (double** mat, int row, int col){   int i,j;   for (i=0; i< row; i++)   {      for (j=0; j< col; j++)         printf ("%9.3lf", mat [i][j]);         printf("\n");   }  } void inverse (doublereal *A, integer N) {    integer i,n,j;    integer lda;    integer  *ipiv;    integer lwork;    integer  info;    doublereal* work;        lda=N;    ipiv = (integer *)malloc(20*N*sizeof(integer));   work = (doublereal*)malloc(10*N*(sizeof(doublereal)));        dgetrf_ ( &N, &N, A, &lda, ipiv, &info );    printf ("Matrix A after dgetrf\n");    for (i=0;i<N*N;i++) {      //for (j=0;j<2;j++)      //{        printf ("A[%d]:%f\t", i,A[i]);    //}   }    printf ("\ninfo for dgetrf_:%d\n\n", info);    //n=2;    //lda=2;    lwork=100*N; //magic number - dimension of work       dgetri_ (&N, A, &lda, ipiv, work, &lwork, &info);    printf ("Matrix A after dgetri\n");    for (i=0;i<N*N;i++) {      //for (j=0;j<2;j++){            printf ("A[%d]:%f\t", i,A[i]);    //}   }    printf ("\ninfo for dgetri_:%d\n\n", info);   free (ipiv);   free (work);  } int main () {      int i,j,N;   // order of matrix (in this case 2*2)    //double A[2][2] = {{3.0,3.0},{3.0,3.0}};    double **A;    double *Ainv;    N = 2;     A = allomat(N,N);    Ainv = (double*)malloc(N*N*sizeof(double));    if (!A||!Ainv)    {printf("error in malloc\n");}    for (i=0;i < N; i++){      for(j=0;j<N;j++)      { scanf("%lf",&A[i][j]);         //A[i][j]=3;      }   }    // thiss part of the code transform the matrix to a column major array    for (i=0;i < N; i++){      for(j=0;j<N;j++)      {   Ainv[i+(N*j)]=A[i][j];      }   }    printf("crap\n");for (i=0;i < N*N; i++){               printf("Ainv[%d]%lf",i, Ainv[i]);         printf("\n");   }   // the matrix order N = 2    inverse(Ainv,N);    for (i=0;i < N; i++){      for(j=0;j<N;j++)      {   A[i][j]=Ainv[i+(N*j)];      }   }    printf("this is the inverse of A\n");    printmat(A,N,N);    free (A);    free(Ainv);    return 0; }`
David12

Posts: 3
Joined: Mon Mar 07, 2011 12:33 pm

### Re: Inversing a Dynamic Matrix using dgetri_ and dgetrf_

Yes, please check out lapacke and look at Julie's example. If you main application code is written in C is definitely worth switching to lapacke. There is plenty of benefit: NaN checks, automatic mem allocation, prototypes already existing, and native C interface. So it's really worth it. Julien.
Julien Langou

Posts: 824
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA