关于mkl中lapack的线性代数问题

frankgty

想解决的是一个Ax=b的问题，第一次使用mkl，选用的函数是SYSV(不确定是否选择正确)。A(Vollmatrix)是一个2400x2400的对称稀疏矩阵，x则是待验证向量Results，右侧Right hand site为RHS。下面是自己试着编的代码，想先试着输出x，然后再跟Results比较。在红色行触发断点，错误0警告0。因为对函数不太了解，输入参数也不知道是否正确，想请教大神帮看一眼，万分感谢。

  Program kontrolle

  use f95_precision
  use lapack95

  Implicit None

  Integer :: n , nrhs , info
  Integer,allocatable :: ipiv(:)
  Real,allocatable :: A(:,:) , B(:,:)
  n = 2400
  nrhs = 1
  Allocate( A(n,n) , B(n,nrhs) , ipiv(n) )

  Open(12,File="Vollmatix.txt")
  Open(13,File="RHS.txt")

  Read(12,*) A(:,:)
  Read(13,*) B(:,:)

  Call SYSV( A(:,:) , B(:,:) , 'L' , ipiv , info )
  Write(*,*) B(:,:)
  Write(*,*) 'IPIV=', ipiv
  Write(*,*) 'INFO=', info

  End Program kontrolle

li913

[Fortran] 纯文本查看 复制代码

! MKL库函数求解线性方程组
include 'lapack.f90'
program test 
use lapack95
integer,parameter::n=3
integer i
real(4)::a(n,n),b(n)
integer ipiv(n)
a=reshape([1,0,1,2,1,0,3,3,0],[n,n])
b(:)=[14,11,1]
call sgetrf(n,n,a,n,ipiv,i) !LU分解
call sgetrs( 'N', n, n, a, n, ipiv, b,n, i )
print*,b !1,2,3
pause
end program 

frankgty

[Fortran] 纯文本查看 复制代码

include 'lapack.f90'
Program kontrolle

use lapack95

Implicit None

Integer,parameter :: n=2400 , nrhs=1
Integer :: i
Integer,allocatable :: ipiv(:)
Real,allocatable :: A(:,:) , B(:,:)

Allocate( A(n,n) , B(n,nrhs) , ipiv(n) )

Open(12,File="Vm.txt")
Open(13,File="RHS.txt")

Read(12,*) A
Read(13,*) B
Close(12)
Close(13)

Call sgetrf(n,n,A,n,ipiv,i)
Call sgetrs('N',n,n,A,n,ipiv,B,n,i)
Print*,B
Pause
End Program kontrolle

您好，十分感谢您的指导，我参照您的例子改写了代码，出现了错误（0x00007FF93CC2B811 (mkl_core.dll)处(位于 Kontrolle.exe 中)引发的异常: 0xC0000005: 读取位置 0x000001D26396D000 时发生访问冲突。）

li913

给你的例子，你能得到正确结果不？我自己测试了一下，新版本ivf的结果不对。所以，你需要试试新的代码。

[Fortran] 纯文本查看 复制代码

! MKL库函数求解线性方程组
! 依赖 mkl_lapack95.lib
program test 
use lapack95
integer,parameter::n=3
integer i
real(4)::a(n,n), b(n)
integer ipiv(n)
a=reshape([1,0,1,2,1,0,3,3,0],[n,n])
b(:)=[14,11,1]
call getrf(a,ipiv,i)   !LU分解
call getrs(a,ipiv, b)  !a,b均被覆盖
print*,b !1,2,3
pause
end program 

frankgty

li913 发表于 2019-9-28 16:32
给你的例子，你能得到正确结果不？我自己测试了一下，新版本ivf的结果不对。所以，你需要试试新的代码。[mw ...

非常感谢您的耐心指导，已得到正确结果