应《网络安全法》要求，自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

捐助贵宾 (小有名气)
引用回帖:
15楼: Originally posted by pippi6 at 2015-05-14 21:58:42
无所谓，应该也可以的...

ok!还有一个问题想请教您，就是我的主程序里面的数据都是单精度，但是子程序里面的计算涉及到双精度，我在调用子程序前需要把数据先变为双精度吗？然后计算出来的逆矩阵，在放回到主程序被应用的话，是不是要再变回单精度？这样会有问题吗？
赞一下 回复此楼
17楼2015-05-14 23:08:17
已阅   回复此楼   关注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 个回答
普通表情 龙 兔 虎 猫
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 311求调剂 +13 zchqwer 2026-04-01 15/750 2026-04-01 22:30 by 2026材料调剂
[考研] 298求调剂 +4 什么是胖头鱼 2026-03-30 6/300 2026-04-01 22:06 by 客尔美德
[考研] 重庆大学材料与化工085600，初试370+，求求调剂建议 +3 shzhou_ 2026-04-01 3/150 2026-04-01 20:41 by 无懈可击111
[考研] 一志愿085600中科院宁波所276分求调剂 +22 材料学257求调剂 2026-03-28 23/1150 2026-04-01 19:03 by 客尔美德
[考研] 0703化学/290求调剂/本科经历丰富/工科也可 +14 丹青奶盖 2026-03-26 16/800 2026-04-01 15:58 by oooqiao
[考研] 0856初试324分求调剂 +5 想上学求调 2026-04-01 5/250 2026-04-01 15:28 by 考研上岸呀一定�
[考研] 085600，材料与化工321分调剂 +5 大馋小子 2026-03-27 7/350 2026-04-01 15:13 by oooqiao
[考研] 环境工程调剂 +9 hyzzzzzzz. 2026-04-01 9/450 2026-04-01 14:20 by salamander`
[考研] 材料调剂 +10 Eujd1 2026-03-31 11/550 2026-04-01 11:23 by ivanqyq
[考研] 08工科，295，接受跨专业调剂 +6 lmnlzy 2026-03-31 6/300 2026-04-01 11:02 by 逆水乘风
[考研] 考研生物与医药调剂 +7 铁憨憨123425 2026-03-31 7/350 2026-04-01 08:45 by JourneyLucky
[考研] 0856材料与化工调剂，339 +10 10213207 2026-03-31 10/500 2026-04-01 08:28 by Linda Hu
[考研] 一志愿华东理工大学，080500学硕，317分，求调剂 +9 s1145 2026-03-31 9/450 2026-03-31 18:47 by JourneyLucky
[考研] 322求调剂：一志愿湖南大学 材料与化工（085600），已过六级。 +10 XX小邓 2026-03-29 10/500 2026-03-31 16:46 by 不吃魚的貓
[考研] 一志愿中海洋材料357 +4 麦恩莉. 2026-03-30 4/200 2026-03-31 14:35 by 记事本2026
[考博] 材料专业申博 +5 杜雨婷dyt 2026-03-29 5/250 2026-03-31 11:19 by oooqiao
[考研] 287求调剂 +17 land xuxu 2026-03-26 17/850 2026-03-31 11:16 by Zzxxxs
[考研] 08工科，295，接受跨专业调剂 +6 lmnlzy 2026-03-30 6/300 2026-03-31 10:04 by cal0306
[考研] 085701求调剂初试286分 +5 secret0328 2026-03-28 5/250 2026-03-30 12:54 by fangnagu
[考研] 086000生物与医药调剂 +5 Feisty。 2026-03-28 9/450 2026-03-29 12:02 by longlotian
信息提示
关闭
请填处理意见
关闭