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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 程序语言 » Fortran » 逆矩阵运算
211/3返回列表123下一页
查看: 2718  |  回复: 20
只看楼主@他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

kimileegdut

捐助贵宾 (小有名气)

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

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

» 猜你喜欢
1楼2015-05-08 22:34:28
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
回帖支持 ( 显示支持度最高的前 50 名 )

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的回帖
普通回帖

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的回帖

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), ...

你好，我试了一下，还是不太对，我把数据传上来，麻烦你帮我看一下，谢谢！
赞一下 回复此楼

» 本帖附件资源列表

  • 欢迎监督和反馈：小木虫仅提供交流平台，不对该内容负责。
    本内容由用户自主发布，如果其内容涉及到知识产权问题，其责任在于用户本人，如对版权有异议，请联系邮箱：xiaomuchong@tal.com
  • 附件 1 : 系数矩阵G
    • 2015-05-13 15:20:05, 150.39 K
5楼2015-05-13 15:21:11
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

pippi6

铁杆木虫 (著名写手)

工程和科学数值计算咨询
引用回帖:
5楼: Originally posted by kimileegdut at 2015-05-13 15:21:11
你好，我试了一下，还是不太对，我把数据传上来，麻烦你帮我看一下，谢谢！...

什么格式？100*100？
赞一下 回复此楼
6楼2015-05-13 16:34:58
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

kimileegdut

捐助贵宾 (小有名气)
引用回帖:
6楼: Originally posted by pippi6 at 2015-05-13 16:34:58
什么格式？100*100？...

单精度，是另外的程序计算之后按标准格式输出的，100*100
赞一下 回复此楼
7楼2015-05-13 20:19:10
 已阅   回复此楼   关注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 ...

好的，很感谢你啊！可能是我调用子程序的时候形参和实参没有对应上！gaussr_inverse 的形参表各个变量是代表什么？
赞一下 回复此楼
9楼2015-05-14 09:48:38
 已阅   回复此楼   关注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=a*c的时候，为什么不直接用矩阵相乘语句b=matmul(a,c),而使用了循环？
赞一下 回复此楼
10楼2015-05-14 10:47:47
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 kimileegdut 的主题更新
211/3返回列表123下一页
普通表情 龙 兔 虎 猫
信息提示
关闭
请填处理意见
关闭