Newton619
2021-09-19 22:36
采纳率: 0%
浏览 98

fortran中mkl库pardiso计算错误

大家好,我想利用pardiso计算Ax=B问题求x,想计算A是4*4矩阵,16个元素都是6.25, x=(1,1,1,1)T,然后用mkl_dcsrsymv计算得到B=(25,25,25,25)T,但是为什么将A和B带入pardiso,得到的x是(0,0,0,4)T呢?代码如下谢谢大家!!

program mkl
    implicit none
    integer, parameter        :: n = 4 
    real(kind=8)              ::  x(n), y(n),A(n)
     integer(kind=4)           ::i(5),j(10)
     real(kind=8)           ::csr(10)     
      character(len=1)          :: uplo = 'u'
     integer(kind=8)                 :: ipt(64)
    integer(kind=8)                 :: idum(n)
    integer(kind=4)                 :: maxfct, mnum, mtype, phase, nrhs, error, msglvl
    integer(kind=4)                 :: iparm(64)
  
     x=(/1,1,1,1/)
        i=(/1,5,8,10,11/)
     j=(/1,2,3,4,2,3,4,3,4,4/)
     csr=(/6.25,6.25,6.25,6.25,6.25,6.25,6.25,6.25,6.25,6.25/)
     call mkl_dcsrsymv( uplo, n, csr, i, j, x, y )
 
    ipt    = 0
    maxfct = 1
    mnum   = 1
    mtype  = -2 
    phase  = 13  
    nrhs   = 1
    msglvl = 0 
    error  = 0    
    iparm  = 0          
    idum=0             
   call pardiso(ipt, maxfct, mnum, mtype, phase, n, csr, i, j, idum, nrhs, iparm, msglvl,y, A, error)
  end program mkl

  • 收藏

3条回答 默认 最新

  • 技术专家团-Joel 2021-09-20 17:13

    没有错误哦,同学,
    你看你用的是全6.25的矩阵,这种4x4的矩阵就是秩亏矩阵,秩为1,所以按理来说是无数个解的,之所以出现某个特解,是跟求解器所使用的求解方法有关。

    打赏 评论
  • 冰岛少年 2021-09-20 08:11

    是正常的,x有无限个解

    img

    打赏 评论
  • 孙叫兽 2021-09-22 11:21

    答案是对的

    打赏 评论

相关推荐 更多相似问题