Complex matrix inverse: C/Lapack

Here a simple example on how to use Lapack within C to invert matrices using the LU decomposition.

The C++ code is different and I will post that another day

[The html code for this post was generated by http://www.bedaux.net/cpp2html/]

#include<stdio.h>
#include<math.h>

#include<complex.h>
#include <stdlib.h>

//........................................................................................
void zgeTranspose( double complex *Transposed, double complex *M ,int n)
{

int i,j;
for(i=0;i<n;i++)

for(j=0;j<n;j++) Transposed[i+n*j] = M[i*n+j];
}

//.........................................................................................
void MatrixComplexInverse(double complex *invA, double complex *A, int n)
{

int LWORK=10*n;

int *permutations;

double complex *WORK, *tempA;

tempA = (double complex*) malloc( n*n*sizeof(double complex) );

permutations = (int*) malloc( 2*n*sizeof(int) );

WORK = (double complex *)malloc(LWORK*sizeof(double complex));


int INFO;

zgeTranspose(tempA,A,n);


zgetrf_( &n, &n, tempA , &n, permutations , &INFO );

if (INFO != 0) {
 printf("ComplexMatrixInverse: Error at zgetrf  \n"); exit(0);
 }



zgetri_( &n, tempA , &n, permutations , WORK, &LWORK, &INFO );

if (INFO != 0) {
 printf("ComplexMatrixInverse: Error at zgetri  \n"); exit(0);
 }

zgeTranspose(invA,tempA,n);

free(WORK);

free(tempA);
free(permutations);

}


/////////////////////////////////////////////////////////////////////////////////////////


int main()
{
int i,j;
const int N = 3;


double complex A[] = { 1.+I , 2. ,  3 , 4. , 5.+I , 6. , 7., 8., 9. + I};

double complex invA[N*N];

MatrixComplexInverse(invA,A,N);


for(i=0;i<N;i++){
 for(j=0;j<N;j++) printf(" (%f,%f) \t", invA[i*N + j]);

 printf("\n");
}

printf("------------------------------------------------------------\n");
return 0;

}


The make file is

#
CC = gcc

#edit LAPACK_PATH if necessary
LAPACK_PATH = /usr/lib64/atlas

a,out: main.c
$(CC)  main.c  -L$(LAPACK_PATH) -llapack -lblas  -lgfortran -lm