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);
}
In LAPACKE_zgetrf(LAPACK_ROW_MAJOR, n, n, tempA , n,&n);
, the last argument of LAPACKE_zgetrf points to n, a single integer. On the contrary, the argument ipiv
should be a pointer to an array of integer of dimension max(m,n), to store pivot indices This could explain a segmentation fault.
The ipiv
computed by LAPACKE_zgetrf()
must also be provided to LAPACKE_zgetri()
as input, to get the correct inverse matrix.