Fortran Coder

标题: 关于mkl中lapack的线性代数问题 [打印本页]

作者: frankgty 时间: 2019-5-9 03:03
标题: 关于mkl中lapack的线性代数问题
本帖最后由 frankgty 于 2019-9-28 00:10 编辑

想解决的是一个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

RHS.txt

32.81 KB, 下载次数: 0

右手侧b

Results.txt

5.81 KB, 下载次数: 0

Vollmatrix.zip

169.98 KB, 下载次数: 0

矩阵A

Vm.zip

145.51 KB, 下载次数: 0

矩阵A

作者: li913 时间: 2019-5-9 10:04

[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 时间: 2019-9-28 00:06
本帖最后由 frankgty 于 2019-9-28 00:19 编辑

[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 时间: 2019-9-28 16:32
本帖最后由 li913 于 2019-9-28 16:57 编辑

给你的例子，你能得到正确结果不？我自己测试了一下，新版本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

QQ截图20190928123041.png (114.18 KB, 下载次数: 194)

作者: frankgty 时间: 2019-9-28 19:09

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

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

欢迎光临 Fortran Coder (http://bbs.fcode.cn/)