Using sparse solver like SuperLU etc. to find the matrix inverse in fortran

Question

Ok i want to compare the result using inverse of matrix for dense rectangular matrix non symmetric. Usually using DGETRF and DGETRI Blas for get the matrix inverse. Let say [2000x2000] double precision matrix A, i want to solve to find the matrix - inverse.

And then i got the SuperLU solver. The idea using sparse solver to get inverse matrix, is to solve A*X=I, where I is the identity matrix. If there is a solution, X will be the inverse matrix. So if the A[2000x2000], Ainvers aka X[2000x2000], identity I[2000x2000]

While the input for calling superlu lib in fortran

call c_fortran_dgssv ( iopt, n, ncc, nrhs, acc, icc,ccc, b, ldb, factors, info )

How to use it? Since the output of the superlu only b matrices, 1 dimension — realbabilu, Mar 21 '22 at 00:22
It seems that you have a computational science question and not a programming question. I suggest to ask at https://scicomp.stackexchange.com/ instead. — Vladimir F Героям слава, Mar 21 '22 at 07:06
You could call your linear system solver N-times (2000x, for every column of X and I) but it would ve very inefficient. Think twice whether you actually need the inverse matrix. Often, it is not very useful for any actual computation. Think about what the A*X=I actually means. — Vladimir F Героям слава, Mar 21 '22 at 07:08
2000x2000 is not very big. It's almost certain that the dense solver will out perform the sparse one. Is there a particular reason you want to use the sparse solver? — Ian Bush, Mar 21 '22 at 07:24
Basically i just want to compare head to head, direct solver vs iterative solver for each size and matrix condition (dense, sparse, ill conditioned, etc.) — realbabilu, Mar 21 '22 at 11:24
Ok Vladimir it worked one by one column as manual does. How can it done in paralel calling superlu ? — realbabilu, Mar 22 '22 at 03:05
How to paralell call superlu, it is nonpure if use do concurrent — realbabilu, Mar 22 '22 at 16:13

realbabilu · Answer 1 · 2022-03-23T02:49:09.883

B in superlu can be (n x n) matrices, A * B = x if A (n x n) input, then inverse A is B (n x n) with x (n x n) as identity matrix, then all B, and X should be (n x n).

  ! a matrix(n,n) bmatrix(n,n) identity matrix(n,n)
  ! find nonzero val in a, get ncc
  ! convert a matrix to ccc matrix first
  allocate (icc(ncc),acc(ncc))
  ....
 
  nrhs = n ! THIS B matrix size
  ldb = n
  allocate (b(n,n))
  b = 0.d0
  !insert identity matrix as B in
  do i = 1, n
    b(i,i) = 1.d0
  end do

  !Factor the matrix in superlu.
  iopt = 1
  call c_fortran_dgssv ( iopt, n, ncc, nrhs, acc, icc,ccc, b, ldb, factors, info )

  !Solve the factored system in superlu.
  iopt = 2
  call c_fortran_dgssv ( iopt, n, ncc, nrhs, acc, icc,ccc, b, ldb, factors, info )

  ! Free memory in superlu.
  iopt = 3
  call c_fortran_dgssv ( iopt, n, ncc, nrhs, acc, icc, ccc, b, ldb, factors, info )
  
  !b out is the inverse
  a = b
  deallocate (b,icc,acc)

Alternatively using B (n) as 1 dimensional matrix for saving memory. as Vladimir pointed out. Do per column of A matrix

do k=1,n
nrhs = 1
ldb = n

do i = 1, n
 if (k.eq.i) b(i) = 1.d0
 if (k.ne.i) b(i) = 0.d0   
end do

!  Factor the matrix.
  iopt = 1
  call c_fortran_dgssv ( iopt, n, ncc, nrhs, acc, icc,ccc, b, ldb, factors, info )
!  Solve the factored system.
  iopt = 2
  call c_fortran_dgssv ( iopt, n, ncc, nrhs, acc, icc,ccc, b, ldb, factors, info )
  do i = 1, n
    a(i,k)= b(i)
  end do
!  Free memory.
  iopt = 3
  call c_fortran_dgssv ( iopt, n, ncc, nrhs, acc, icc, ccc, b, ldb, factors, info )
end do !k

  deallocate (icc,acc)

Using sparse solver like SuperLU etc. to find the matrix inverse in fortran

1 Answers1