应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 程序语言 » Fortran » 逆矩阵运算
北京石油化工学院2026年研究生招生接收调剂公告
52/1返回列表
查看: 2887  |  回复: 20
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看
当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

kimileegdut

捐助贵宾 (小有名气)

[求助] 逆矩阵运算 已有2人参与

求助各位大神，本人最近在编写一个程序，涉及到逆矩阵运算，而且矩阵维数很大（100*100），采用了西尔曼-莫里生公式求逆矩阵，但是后续计算结果有误，于是验算逆矩阵的正误，发现原来的矩阵和逆矩阵相乘不为单位阵，想请教一下各位前辈，这种情况求逆矩阵应该怎么处理？
回复此楼

» 猜你喜欢
1楼 2015-05-08 22:34:28
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

kimileegdut

捐助贵宾 (小有名气)
引用回帖:
8楼: Originally posted by pippi6 at 2015-05-13 21:27:12
我试了一下，gaussr_inverse 应该工作。 逆矩阵和原矩阵的乘积与单位矩阵的最大差值为  2.442490654175344E-015。


subroutine  InverseMatrix
  implicit none
  integer, parameter :: nn=100,m=20
  inte ...

你好，我还有一个不理解的地方，就是在调用子程序计算逆矩阵前，为什么要先把原来的矩阵元素值赋给b?如果没有b=a，这一句，调用子程序时实参用a，计算出来的结果有误，这什么原因导致的？
赞一下 回复此楼
19楼2015-05-15 10:54:50
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 21 个回答

pippi6

铁杆木虫 (著名写手)

工程和科学数值计算咨询

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
感谢参与，应助指数 +1
kimileegdut: 金币+20, ★有帮助 2015-05-11 11:49:44
试试这个 全面选主元的gauss消去。不行的话。把你的数据传上来


!flag
subroutine gaussr_inverse(a,n,np,b,m,mp)
  implicit none
  integer :: n,np,m,mp,i,j,k,irow,icol,l
  real*8 a(np,np),b(np,mp),pivinv,cc
  real*8, allocatable, dimension( :: dum
  real*8 big
  integer, allocatable, dimension( :: ipiv,indxr,indxc
  real*8, allocatable, dimension( :: c1

  allocate(ipiv(1:n),indxr(1:n),indxc(1:n),dum(1:max(n,m)),c1(1:n))
  ipiv=0

  !     row normalization (implicite pivoting)
  do i=1,n
     c1(i)=1/sqrt(sum(abs(a(i,1:n))**2))
     a(i,1:n)=a(i,1:n)*c1(i)
     b(i,1:m)=b(i,1:m)*c1(i)
  end do

  do i=1,n

     !     Full pivoting, however, not necessarily starting from 1st col
     big=0.d0
     do j=1,n
        if(ipiv(j)==1)cycle
        do k=1,n
           if (ipiv(k)==1) cycle
           if (abs(a(j,k)) <= big) cycle
           big=abs(a(j,k))
           irow=j
           icol=k
        end do
     end do

     !     ipiv(icol)=1 means the column icol has been selected once
     ipiv(icol)=ipiv(icol)+1
     if (irow.ne.icol) then

        !     exchange a(irow, with a(icol,
        !     so a(irow,irow) is the selected pivot
        dum(1:n)= a(irow,1:n)  
        a(irow,1:n)= a(icol,1:n)  
        a(icol,1:n)=dum(1:n)

        dum(1:m)= b(irow,1:m)  
        b(irow,1:m)= b(icol,1:m)  
        b(icol,1:m)=dum(1:m)
     endif

     !     bookkeeping row- and column-indices
     indxr(i)=irow
     indxc(i)=icol
     if (abs(a(icol,icol))==0.d0) pause 'singular matrix.'
     pivinv=1.d0/a(icol,icol)
     a(icol,icol)=1.d0
     a(icol,1:n)=a(icol,1:n)*pivinv
     b(icol,1:m)=b(icol,1:m)*pivinv

     do l=1,n
        if(l==icol)cycle  ! only for row l not equal to icol
        cc=a(l,icol)
        a(l,icol)=0.d0    ! clever design
        a(l,1:n)=a(l,1:n)-a(icol,1:n)*cc
        b(l,1:m)=b(l,1:m)-b(icol,1:m)*cc
     end do
  end do

  do l=1,n
     if(indxr(l) == indxc(l)) cycle
     dum(1:n)=a(1:n,indxr(l))
     a(1:n,indxr(l))=a(1:n,indxc(l))
     a(1:n,indxc(l))=dum(1:n)
  end do

  do i=1,n
     a(1:n,i)=a(1:n,i)*c1(i)
  end do

  deallocate(ipiv,indxr,indxc,dum,c1)
end subroutine gaussr_inverse
赞一下 回复此楼
2楼2015-05-10 08:27:12
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

kimileegdut

捐助贵宾 (小有名气)
引用回帖:
2楼: Originally posted by pippi6 at 2015-05-10 08:27:12
试试这个 全面选主元的gauss消去。不行的话。把你的数据传上来


!flag
subroutine gaussr_inverse(a,n,np,b,m,mp)
  implicit none
  integer :: n,np,m,mp,i,j,k,irow,icol,l
  real*8 a(np,np),b(np,mp), ...

ok,我先试试，谢谢！
赞一下 回复此楼
3楼2015-05-11 11:47:39
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

kimileegdut

捐助贵宾 (小有名气)
分享一个求逆矩阵的代码


SUBROUTINE INVERMAT(AA,N,AV)
!
!程序说明：利用西尔曼-莫里生(Sherman-Morrison)
!          公式求方阵[A]的逆矩阵，能够克服矩阵元素由于阶差较大引起的困难.
!          [A]存放在数组AA(N,N),逆矩阵存放在数组AV(N,N)中.
!
!
!定义数据类型：注意采用双精度浮点数
REAL(kind=8),DIMENSION(N,N)::AA,TA,AV
REAL(kind=8),DIMENSION(N)::VT,WK,DD
!
!
    DO i=1,N
      DO j=1,N
        TA(i,j)=AA(i,j)
      END DO
    END DO
!
!
    DO i=1,N
      DD(i)=1.D0
      IF(TA(i,i).NE.0.) DD(i)=SQRT(ABS(TA(i,i)))
      DO j=1,N
        TA(i,j)=TA(i,j)/DD(i)
        TA(j,i)=TA(j,i)/DD(i)
        AV(i,j)=0.D0
      END DO
      AV(i,i)=1.D0
      TA(i,i)=TA(i,i)-1.D0
      WK(i)=1.D0
    END DO
!
    LJK=0
30  ISW=0
    DO i=1,N
      IF(WK(i).LE.1.E-40) GOTO 40
      DO j=1,N
        VT(j)=0.D0
        DO k=1,N
          VT(j)=VT(j)+TA(i,k)*AV(k,j)
        END DO
      END DO
      VTA=1.D0+VT(i)
      IF(ABS(VTA).LE.1.E-40) GO TO 70
      WK(i)=0.
      DO j=1,N
        HH=AV(j,i)
        DO k=1,N
          AV(j,k)=AV(j,k)-HH*VT(k)/VTA
        END DO
      END DO
      GO TO 40
70    ISW=ISW+1
    END DO
40  CONTINUE
!
    IF(ISW.EQ.LJK) GO TO 80
    LJK=ISW
    GO TO 30
80  CONTINUE
    DO i=1,N
      DO j=1,N
        AV(i,j)=AV(i,j)/DD(i)/DD(j)
      END DO
    END DO
!
!
RETURN
END SUBROUTINE INVERMAT
赞一下 回复此楼
4楼2015-05-11 18:03:28
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 21 个回答
普通表情 龙 兔 虎 猫
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 336求调剂 +4 kiyy 2026-04-01 4/200 2026-04-02 06:51 by ilovexiaobin
[考研] 求生物学调剂 +8 15172915737 2026-04-01 8/400 2026-04-02 06:49 by ilovexiaobin
[考研] 271求调剂 +14 勒布朗@ 2026-03-31 19/950 2026-04-02 00:01 by 勒布朗@
[考研] 求调剂推荐 材料 304 +21 荷包蛋hyj 2026-03-26 21/1050 2026-04-01 21:09 by lijunpoly
[考研] 材料调剂 +14 一样YWY 2026-04-01 14/700 2026-04-01 21:07 by lijunpoly
[考研] 342求调剂 +12 Mary Keen 2026-03-28 13/650 2026-04-01 21:02 by 流情牧豪
[考研] 【求调剂】085601材料工程专硕 | 总分272 | +10 脚滑的守法公民 2026-03-27 10/500 2026-04-01 17:23 by pies112
[考研] 26考研调剂 +3 Wnz.20030617 2026-04-01 3/150 2026-04-01 17:22 by salamander`
[考研] 085600，材料与化工321分，求调剂 +11 大馋小子 2026-03-27 11/550 2026-04-01 16:10 by cymywx
[考研] 330分求调剂 +11 qzenlc 2026-03-29 11/550 2026-04-01 14:32 by chenqifeng666
[考研] 一志愿北交材料工程总分358 +5 cs0106 2026-04-01 7/350 2026-04-01 11:45 by wangjy2002
[硕博家园] 博一被送出联培感觉不适应怎么办 +3 全村的狗 2026-03-31 3/150 2026-04-01 10:44 by 328838485
[考研] 考研调剂 +9 小蜡新笔 2026-03-29 10/500 2026-03-31 19:52 by Dyhoer
[考研] 0856求调剂 +9 楒桉 2026-03-28 9/450 2026-03-31 19:06 by 暮泽12
[考研] 江苏苏北高校诚邀调剂同学 +3 zzll406 2026-03-31 3/150 2026-03-31 16:54 by 及时行乐fan
[考研] 哈尔滨工业大学材料与化工专硕378求调剂 +3 塔比乌斯 2026-03-30 3/150 2026-03-30 22:55 by 无际的草原
[考研] 297 地理学070500 复试求调剂 +3 小圆圈圈ooo 2026-03-30 3/150 2026-03-30 21:05 by 余震yz
[考研] 312，生物学求调剂 +3 小译同学abc 2026-03-28 3/150 2026-03-28 15:32 by 落睿可思
[考研] 292求调剂 +4 求求了收下我吧� 2026-03-26 4/200 2026-03-27 10:37 by zhshch
[考研] 341求调剂 +7 青柠檬1 2026-03-26 7/350 2026-03-27 00:19 by wxiongid
信息提示
关闭
请填处理意见
关闭