| 查看: 2893 | 回复: 20 | ||
| 当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖 | ||
[求助]
逆矩阵运算 已有2人参与
|
||
| 求助各位大神，本人最近在编写一个程序，涉及到逆矩阵运算，而且矩阵维数很大（100*100），采用了西尔曼-莫里生公式求逆矩阵，但是后续计算结果有误，于是验算逆矩阵的正误，发现原来的矩阵和逆矩阵相乘不为单位阵，想请教一下各位前辈，这种情况求逆矩阵应该怎么处理？ |
回复此楼
5楼2015-05-13 15:21:11
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
