I am trying to write a code using lapacke libraries to invert a complex matrix in C. However I am stuck with a segmentation fault which seems to depend on the size N of the matrix. What's more is that the size for which the segmentation fault happens varies each time I compile the program or I touch anything. This makes me think that somewhere the code is trying to access badly allocated or forbidden memory. Unfortunately, I dont' understand how this happens since it seems to be related with the LAPACKE functions themeselves. In fact, when the function /*MatrixComplexInv(invA,A,N);*/ ( in which the LAPACKE functions are called for the inversion) is commented the segmentation fault doesn't happen.
Below there is the working code that can be compiled and run on its own.
#include <stdio.h>
#include <lapacke.h>
#include <complex.h>
#include <stdlib.h>
#include <math.h>
void Ctranspose( double complex *, double complex * ,int );
void MatrixComplexInv(double complex *, double complex *, int );
int main(int argc, const char * * argv) {
int i,j,k,N = 4;/*if N> bigger than a small number 4,6,7.. it gives segmentation fault*/
double complex *A = calloc(N*N,sizeof(double complex)),
*b = calloc(N*N,sizeof(double complex)),
*Ap =calloc(N*N,sizeof(double complex));
double complex *invA =calloc(N*N,sizeof(double complex));
for(i=0;i<N;i++){
for(j=0;j<N;j++){
A[i*N+j] = 1+sin(i*j)*i+I*j;
Ap[i*N+j] = 1+sin(i*j)*i+I*j;
}
}
/*Segmentation fault in this function, due to
*
LAPACKE_zgetrf(LAPACK_ROW_MAJOR, n, n, tempA , n,&n);
LAPACKE_zgetri(LAPACK_ROW_MAJOR, n, tempA , n, &n );
*
* both.
*/
MatrixComplexInv(invA,A,N);
for(i=0;i<N;i++){
for(j=0;j<N;j++){
for(k = 0;k<N;k++){
b[i*N+j]+=invA[i*N + k]*Ap[k*N + j];
}
printf("(%lf,%lf)\t", creal(b[i*N + j]),cimag(b[i*N + j]));/*tests that the result produces the inverse matrix A^{-1}A = 1*/
}
printf("\n");
}
return 0;
}
void Ctranspose( 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 MatrixComplexInv(double complex *invA, double complex *A, int n)
{
double complex *tempA = (double complex*) malloc( n*n*sizeof(double complex) );
Ctranspose(tempA,A,n);
/*SEGMENTATION HAPPEN IN THESE TWO FUNCTIONS*/
LAPACKE_zgetrf(LAPACK_ROW_MAJOR, n, n, tempA , n,&n);
LAPACKE_zgetri(LAPACK_ROW_MAJOR, n, tempA , n, &n );
Ctranspose(invA,tempA,n);
free(tempA);
}
mallocbut freeing it withLAPACKE_free?