complex_matrix_inversion.c

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <complex.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_linalg.h>

int main() {
    int Ndrive = 3;
    complex double **overlap_matrix = (complex double**)malloc(Ndrive * sizeof(complex double*));
    complex double **overlap_matrix_inverse = (complex double**)malloc(Ndrive * sizeof(complex double*));

    for(int i = 0; i < Ndrive; i++) {
        overlap_matrix[i] = (complex double*)malloc(Ndrive * sizeof(complex double));
        overlap_matrix_inverse[i] = (complex double*)malloc(Ndrive * sizeof(complex double));
    }

    overlap_matrix[0][0] = 1.0 + I * 0.0;
    overlap_matrix[0][1] = 0.0 + I * 0.0;
    overlap_matrix[0][2] = 2.0 + I * 3.0;
    overlap_matrix[1][0] = 0.0 + I * 0.0;
    overlap_matrix[1][1] = 3.0 + I * 0.0;
    overlap_matrix[1][2] = 4.0 + I * 0.0;
    overlap_matrix[2][0] = 2.0 - I * 3.0;
    overlap_matrix[2][1] = 4.0 + I * 0.0;
    overlap_matrix[2][2] = 5.0 + I * 0.0;

    /*Inverse is:
        {{ 0.025        , -0.2 -i0.3 ,  0.15 + i0.225},
         {-0.2 + i0.3   ,  0.2       ,  0.1          },
         { 0.15 - i0.225,  0.1       , -0.075        }};
    */


    //GSL Complex Matrix Inversion to find inverse of overlap matrix
    size_t matSize = Ndrive;
    gsl_permutation *perm = gsl_permutation_calloc(matSize);
    gsl_matrix_complex *mat = gsl_matrix_complex_calloc(matSize, matSize);
    gsl_matrix_complex *inv = gsl_matrix_complex_calloc(matSize, matSize);

    gsl_complex temp;
    for(int i = 0; i < matSize; i++) {
        for(int j = 0;j < matSize; j++) {
            GSL_SET_COMPLEX(&temp, creal(overlap_matrix[i][j]), cimag(overlap_matrix[i][j]));
            gsl_matrix_complex_set(mat,i,j,temp);
        }
    }      

    int s;
    gsl_linalg_complex_LU_decomp(mat, perm, &s);
    gsl_linalg_complex_LU_invert(mat, perm, inv);

    for(int i = 0; i < matSize; i++) {
        for(int j = 0; j < matSize; j++) {
            gsl_complex temp_result = gsl_matrix_complex_get(inv, i, j);
            overlap_matrix_inverse[i][j] = GSL_REAL(temp_result) + I * GSL_IMAG(temp_result);
            printf("%f + %fi \n", creal(overlap_matrix_inverse[i][j]), cimag(overlap_matrix_inverse[i][j]));
        }
    }
    gsl_matrix_complex_free(inv);
    gsl_matrix_complex_free(mat);
    gsl_permutation_free(perm);

    return 0;
}