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

24小时热门版块排行榜    

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

kimileegdut

捐助贵宾 (小有名气)

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

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

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

pippi6

铁杆木虫 (著名写手)

工程和科学数值计算咨询
我试了一下，gaussr_inverse 应该工作。 逆矩阵和原矩阵的乘积与单位矩阵的最大差值为  2.442490654175344E-015。


subroutine  InverseMatrix
  implicit none
  integer, parameter :: nn=100,m=20
  integer :: i,j,k,imax,jmax
  real*8,dimension(nn,nn) :: a,b,c,e
  real*8 :: s,dm,d

  open(20,file='coeff.dat')
  read(20,*) a
  close(20)
  b=a

  ! e us unit matrix
  e=0
  do i=1,nn
     e(i,i)=1
  end do
  c=e

  ! c is inversed matrix
  call gaussr_inverse(b,nn,nn,c,nn,nn)

  ! b=a*c, and check if maximal derivation from unit matrix e

  dm=0
  do i=1,nn
     do j=1,nn
        s=0
        do k=1,nn
           s=s+a(i,k)*c(k,j)
        end do
        b(i,j)=s

        d=b(i,j)-e(i,j)
        if(abs(d) > abs(dm) ) then
           dm=d
           imax=i
           jmax=j
        end if
     end do
  end do

  ! print first 20x20 sub-block of b just to get a sense
  write(6,1) b(1:m,1:m)
1 format('b='/(<m>e11.3))

  print*,' imax=',imax,' jmax=',jmax,' dm=',dm
  print*,' b(imax,jmax)=',b(imax,jmax)
!  print*,' a(nn,nn)=',a(nn,nn)
!  print*,' a(nn-1,nn)=',a(nn-1,nn)

end subroutine InverseMatrix
赞一下(1人) 回复此楼

» 本帖已获得的红花（最新10朵）

kimileegdut kimileegdut kimileegdut
8楼2015-05-13 21:27:12
已阅   回复此楼   关注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 个回答
普通表情 龙 兔 虎 猫
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 085601一志愿中山大学深圳材料工程330求调剂 +7 pipiver 2026-03-30 7/350 2026-04-02 07:58 by 尚水阁主
[考研] 化学工程专硕324分，一志愿中国矿业大学求调剂 +7 耿耿1314 2026-04-01 7/350 2026-04-02 07:40 by 尚水阁主
[考研] 生物学308分求调剂（一志愿华东师大） +4 相信必会光芒万� 2026-03-31 4/200 2026-04-02 06:52 by ilovexiaobin
[考研] 319求调剂 +9 太容易1018 2026-04-01 9/450 2026-04-01 21:58 by 我的船我的海
[考研] 070300化学279求调剂 +15 哈哈哈^_^ 2026-03-31 17/850 2026-04-01 21:37 by 给你你注意休息
[考研] 362求调剂 +13 西南交材料专硕3 2026-03-31 13/650 2026-04-01 17:38 by JYD2011
[考研] 350求调剂 +7 阿佳～ 2026-03-31 7/350 2026-04-01 16:12 by yanflower7133
[考研] 0703求调剂 +4 zizimo 2026-03-31 4/200 2026-04-01 16:04 by yanflower7133
[考研] 生物学296求调剂 +10 汤圆包 2026-03-29 14/700 2026-04-01 10:44 by 求调剂zz
[考研] 282求调剂 +6 呼吸都是减肥 2026-04-01 6/300 2026-04-01 08:58 by laoshidan
[考研] 张芳铭-中国农业大学-环境工程专硕-298 +9 手机用户 2026-03-26 9/450 2026-03-31 18:09 by 544594351
[考研] 材料求调剂 一志愿哈工大总分298分，前三科223分 +11 dongfang59 2026-03-27 11/550 2026-03-31 16:51 by Wang200018
[考研] 286求调剂 +6 Faune 2026-03-30 6/300 2026-03-31 14:37 by jp9609
[考研] 272求调剂，接受跨专业调剂！ +3 闲鱼卢 2026-03-31 3/150 2026-03-31 13:00 by 替代品000
[考研] 08工科，295，接受跨专业调剂 +6 lmnlzy 2026-03-30 6/300 2026-03-31 10:04 by cal0306
[考研] 22408 359分调剂 +4 Qshers 2026-03-27 8/400 2026-03-31 08:53 by Qshers
[考研] 327求调剂 +5 小卡不卡. 2026-03-29 5/250 2026-03-30 19:30 by Wang200018
[考研] 材料化工340求调剂 +3 jhx777 2026-03-30 3/150 2026-03-30 17:54 by JourneyLucky
[考研] 085602 化学工程专硕 340分求调剂 +4 qianbai11 2026-03-29 4/200 2026-03-30 11:34 by 唐沐儿
[考研] 药学105500求调剂 +3 Ssun。。 2026-03-28 3/150 2026-03-28 11:24 by lxf170613
信息提示
关闭
请填处理意见
关闭