使用 Fortran 中的 Strumpack 求解稀疏线性系统

Question

我必须解决 CSR 格式的稀疏线性系统。我的代码使用库 Strumpack。 我不明白为什么它不返回解决方案。 CSR 格式的矩阵由 col_index、row_index 和 valeurs（存储在矩阵中的实数值）构成。代码从文件中读取参数，组装矩阵，然后进入“do while”循环。

program main
use, intrinsic :: ISO_C_BINDING
use strumpack
use hermite
use linalg
implicit none




integer(c_int) :: ierr
integer(c_int) :: k      ! grid dimension
integer, target :: matrice_dimension

integer  ::   Nx, i, troncature , nbblocs,j, Nv,r
real(rp) ::   Tfin ,dx,dxdemi,dt,date, tzero, tauzero, eps, dv,xj,vr , delta

integer, dimension(:), allocatable, target ::  row_index, col_index
integer  :: rc  
real(rp), dimension(:), allocatable, target :: D, Dnew, psi_v , diag , valeurs
real(rp), dimension(:,:), allocatable :: Mat  , D_matrice , solution_finale


! sparse solver object
type(STRUMPACK_SparseSolver) :: S



!!! lecture des parametres
open(unit = 10, file = "parametres.txt")
read(10,*)  Nx, Nv
read(10,'(E15.8)')  Tfin
read(10,*)  troncature
read(10,'(E15.8)') tzero
read(10,'(E15.8)') delta
read(10,'(E15.8)') eps
read(10,'(E15.8)') tauzero



!!!!! il y a nbblocs*nbblocs dans la matrice 
nbblocs = troncature + 1


matrice_dimension = (troncature+1)*Nx

!!!!! allocation 
!allocate(Mat((troncature+1)*Nx,(troncature+1)*Nx))
!allocate(diag(Nx*nbblocs))

allocate(D_matrice(Nx,nbblocs))
allocate(psi_v(0:troncature))
allocate(D(matrice_dimension))
allocate(Dnew(matrice_dimension))
allocate(solution_finale(Nx,Nv))



allocate(valeurs(8*Nx+7*(troncature-1)*Nx))
allocate(col_index(8*Nx+7*(troncature-1)*Nx))
allocate(row_index((troncature+1)*Nx+1))


!!!! pas de temps et d'espace et de vitesse
dx = (xdroite-xgauche)/real(Nx,rp)
dv = (vdroite-vgauche)/real(Nv-1,rp)


dt = 0.001_rp 


print*,""
print*," !!!!!!!!!!!! PARAMETRES  !!!!!!!!!!!!!!"
print*,"Tfin = ",Tfin
print*,"Intervalle [",xgauche,",",xdroite,"]"
print*,"dx =",dx,",dt =",dt
print*,"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"
print*," " 





!!!!!! initialisation 

!!!!! D(t=0,x)
!call vect_init_cas1(nbblocs,Nx,dx,delta,c0,tzero,D)


!call vect_init(nbblocs,troncature,Nx,xgauche,dx,tzero,c0,D)
!print*,D

!call vect_init_trapeze(nbblocs,troncature,Nx,xgauche,dx,tzero,c0,D)



D = 1.0_rp
Dnew=5.0_rp



print*,D

call matricecreuse(dx,dt,tauzero,eps,tzero,Nx,troncature,valeurs,col_index,row_index)

!call STRUMPACK_init_mt(S, STRUMPACK_DOUBLE, STRUMPACK_MT, 0, c_null_ptr, 1)

call STRUMPACK_set_csr_matrix(S, c_loc(matrice_dimension), c_loc(row_index), c_loc(col_index), c_loc(valeurs), 0);

date = 0.0_rp


open(unit=50,file="erreur1.dat")
open(unit=51,file="erreur2.dat")

write(50,"(E15.8,x,E15.8)") date , sqrt(dx/(xdroite-xgauche)) * norm2(D(Nx+1:))
write(51,"(E15.8,x,E15.8)") date , erreur2(troncature,Nx,dx,tzero,D)



!!!!!!!!!!!!!!! RESOLUTION MATRICE PLEINE !!!!!!!!!!!!!!!!!!!!!

do while( date<Tfin) 
    dt = min(dt,Tfin-date)
    date = date  + dt
    print*,"date =",date
    ierr = STRUMPACK_solve(S, c_loc(D), c_loc(Dnew), 0);
    print*,"ierr=",ierr
    print*,Dnew
    D = Dnew
    write(50,"(E15.8,x,E15.8)") date , sqrt(dx/(xdroite-xgauche)) * norm2(D(Nx+1:))
    write(51,"(E15.8,x,E15.8)") date , erreur2(troncature,Nx,dx,tzero,D)
    exit
end do
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

close(50)
close(51)

!!!!!!!!! EVALUATION DE LA SOLUTION SUR UNE GRILLE xj, j=1,...,Nx et vr, r=1,...,Nv


D_matrice = reshape(D, (/Nx,nbblocs/))

do j=1,Nx
    xj = (xgauche+0.5_rp*dx) + real(j-1,rp)*dx
    D_matrice(j,:) = sqrt(rho_infini(tzero,c0,xj)) * D_matrice(j,:)
end do

do r=0,Nv-1
    vr = vgauche + r*dv
    call hermite_polynome(troncature,vr,psi_v)
    psi_v = maxwellien(tzero,vr)*psi_v
    solution_finale(:,r+1) = matmul(D_matrice,psi_v) 
end do


print*,"Minimum de la solution : ",minval(solution_finale)


!!!!!! ECRITURE DES FICHIERS 
open(unit=50,file="approximation.dat")
open(unit=51,file="stationnaire.dat")
open(unit=52,file="erreur.dat")
open(unit=53,file="initial.dat")

do r=1,Nv
    do j=1,Nx
        xj = (xgauche+0.5_rp*dx) + real(j-1,rp)*dx
        vr = vgauche + real(r-1,rp)*dv
        write(50,"(E15.8,x,E15.8,x,E15.8)") xj,vr,solution_finale(j,r)
        write(51,"(E15.8,x,E15.8,x,E15.8)") xj,vr, stationnaire(c0,tzero,xj,vr)
        write(52,"(E15.8,x,E15.8,x,E15.8)") xj,vr, abs(solution_finale(j,r)-stationnaire(c0,tzero,xj,vr))

    end do
    write(50,*) ""
    write(51,*) ""
    write(52,*) ""
end do


do r=1,Nv
    vr = vgauche + real(r-1,rp)*dv
    write(53,"(E15.8,x,E15.8,x,E15.8)") vr,solution_finale(1,r),f_init(c0,tzero,1.0_rp,vr)
end do

close(50)
close(51)
close(52)
close(53)

call execute_command_line ("gnuplot trace.p")





!!!!!!!! Désallocation mémoire, fin du code
deallocate(D_matrice,D,Dnew,psi_v)



end program main

我确定矩阵是正确的（我在上面使用了 Scipy，效果很好）。我试图调试它，方法是在循环之前在 D 和 Dnew 中设置值，并打印 vector 和 ierr。一轮迭代后，Dnew还是5，ierr是0。我认为这意味着Strumpack没有解决这个问题。因此，我很可能在进入循环之前误用了 Strumpack。我很难在 the strumpack 的网站 上找到解释。

有人有答案吗？

致以诚挚的问候。