Re: [cython-users] Using gsl to invert a matrix

[Hanno Klemm]

On 24. Nov 2017, at 12:35, Zahra Sheikh <sheikh...@gmail.com> wrote:

Hi,

What is the fastest and most efficient way to invert a huge 2-d matrix in cython? I tried to use gsl library and wrote the following code but it crashes (Segmentation fault (core dumped)) when I tried to set gsl matrix from a 2d pointer array

import numpy as np
cimport numpy as np
cimport cython
from libc.stdlib cimport malloc, free
from cython_gsl cimport *
from libc.string cimport memset
cdef void freeMemory(double** arr, int rows):
     cdef int i
     if (arr!=NULL):
        for i from 0 <= i < rows:
            if (arr[i]!=NULL):
               free(arr[i])
        free(arr)
               
cdef double** gsl_matrix_set_inverse(double** Nm, int num):
     cdef:
         gsl_matrix *m1
         gsl_matrix *inverse
         int i,j
         gsl_permutation *p
         double** NmInv
     p = gsl_permutation_alloc(num)
   
     for i from 0 <= i < num:
         for j from 0 <= j < num:
             gsl_matrix_set(m1, i, j, Nm[i][j])
    
     gsl_linalg_LU_invert(m1, p, inverse)
     NmInv= <double**> malloc(sizeof(double*) * num)
     for i from 0 <= i < num:
         NmInv[i] = <double*> malloc(sizeof(double) * num)
     #convert gsl matrix back to 2d pointer array
     for i from 0 <= i < num:
        for j from 0 <= j < num:
            NmInv[i][j] = gsl_matrix_get(inverse, i, j)
     gsl_matrix_free(m1)      
     return NmInv
    
def run( int n):
    cdef int i,j
    cdef double** matrix = <double**> malloc(sizeof(double*) * n)

    for i from 0 <= i < n:
        matrix[i] = <double*> malloc(sizeof(double) * n)
        memset(matrix[i], 0, sizeof(double) * n)

    for i from 0 <= i < n:
         for j from 0 <= j < n:
             if ((i!=0) and (j!=0)):
                matrix[i][j]=matrix[i-1][j-1]+1.
             else:
                matrix[i][j]=1.
    print "inverse of a matrix:"
    cdef double** matrixInverse=gsl_matrix_set_inverse(matrix, n)
    for i from 0 <= i < n:
         for j from 0 <= j < n:
             print matrixInverse[i][j]
    freeMemory(matrixInverse,n)
   
    freeMemory(matrix,n)
Any suggestion why this happens? Is it the best way to invert a matrix?

Cheers

--

Depending on the mathematical library NumPy is build with in your Python installation, just using np.linalg.inv ( or scipy.linalg.inv) might be quite efficient already. Usually, they hand off these tasks to the underlying Blas/Lapack implementation. Any particular reason why you want to use gsl?

Hanno