使用 SuperLU 等稀疏求解器在 fortran 中求逆矩阵

Question

好的，我想使用非对称密集矩形矩阵的矩阵逆来比较结果。 通常使用 DGETRF 和 DGETRI Blas 来获取矩阵逆。 假设 [2000x2000] 双精度矩阵 A，我想求解找到矩阵 - 逆矩阵。

然后我得到了 SuperLU 解算器。 使用稀疏求解器获得逆矩阵的想法是求解 A*X=I，其中 I 是单位矩阵。如果有解，X就是逆矩阵。因此，如果 A[2000x2000]，Ainvers 又名 X[2000x2000]，恒等式 I[2000x2000]

而在 fortran 中调用 superlu lib 的输入

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

Answer 1

superlu 中的 B 可以是 (nxn) 矩阵，
A * B = x
如果输入A(nxn)，则A的逆矩阵是B(nxn)，以x(nxn)为单位矩阵，则所有B和X应该是(nxn)。

  ! 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)

或者使用 B(n) 作为一维矩阵以节省内存。
正如弗拉基米尔指出的那样。 A 矩阵每列的执行次数

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)