| 查看: 2887 | 回复: 20 | |||
[求助]
逆矩阵运算 已有2人参与
|
|||
| 求助各位大神，本人最近在编写一个程序，涉及到逆矩阵运算，而且矩阵维数很大（100*100），采用了西尔曼-莫里生公式求逆矩阵，但是后续计算结果有误，于是验算逆矩阵的正误，发现原来的矩阵和逆矩阵相乘不为单位阵，想请教一下各位前辈，这种情况求逆矩阵应该怎么处理？ |
回复此楼» 猜你喜欢
336求调剂
已经有4人回复
-
求生物学调剂
已经有8人回复
-
08生物与医药专硕初试346找调剂
已经有6人回复
-
生物学327，求调剂
已经有6人回复
-
0710生物学求调剂
已经有7人回复
-
385分 生物学（071000）求调剂
已经有4人回复
-
材料考研调剂
已经有4人回复
-
366求调剂一志愿东北大学
已经有3人回复
-
08工科275分求调剂
已经有13人回复
-
材料工程322分
已经有6人回复
pippi6
铁杆木虫 (著名写手)
工程和科学数值计算咨询
- 应助: 413 (硕士)
- 贵宾: 0.002
- 金币: 7116.5
- 散金: 15
- 红花: 63
- 帖子: 1639
- 在线: 798.9小时
- 虫号: 2469437
- 注册: 2013-05-14
- 专业: 计算数学与科学工程计算
|
我试了一下，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 |
» 本帖已获得的红花（最新10朵）
kimileegdut8楼2015-05-13 21:27:12
pippi6
铁杆木虫 (著名写手)
工程和科学数值计算咨询
- 应助: 413 (硕士)
- 贵宾: 0.002
- 金币: 7116.5
- 散金: 15
- 红花: 63
- 帖子: 1639
- 在线: 798.9小时
- 虫号: 2469437
- 注册: 2013-05-14
- 专业: 计算数学与科学工程计算
【答案】应助回帖
★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
感谢参与，应助指数 +1
kimileegdut: 金币+20, ★有帮助 2015-05-11 11:49:44
感谢参与，应助指数 +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(  :: dumreal*8 big integer, allocatable, dimension( :: ipiv,indxr,indxcreal*8, allocatable, dimension( :: c1allocate(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
3楼2015-05-11 11:47:39
|
分享一个求逆矩阵的代码 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
5楼2015-05-13 15:21:11
pippi6
铁杆木虫 (著名写手)
工程和科学数值计算咨询
- 应助: 413 (硕士)
- 贵宾: 0.002
- 金币: 7116.5
- 散金: 15
- 红花: 63
- 帖子: 1639
- 在线: 798.9小时
- 虫号: 2469437
- 注册: 2013-05-14
- 专业: 计算数学与科学工程计算
6楼2015-05-13 16:34:58
7楼2015-05-13 20:19:10
9楼2015-05-14 09:48:38
10楼2015-05-14 10:47:47
