Inversing a Dynamic Matrix using dgetri_ and dgetrf_

[David12]

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;

}

[Julien Langou]

[ 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 )

3) please change your prototype

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

[admin]

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 Matrix
for(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 dgetrf
info=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 dgetri
printf("\tDGETRI [n=%d] INFO=%d \n",n,info);

// Check
printf("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
The compile and link is

gcc -c prog.c
gfortran prog.o -o prog.exe lapacke.a cblas.a lapack.a blas.a

[David12]

[ 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 )

3) please change your prototype

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;

 }

[Julien Langou]

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.