用IMSL中DLINRG程序求矩阵逆的问题

GWinston

我用IMSL库的子程序DLINRG来求矩阵A的逆（放在AINV中），然后演算结果是否正确，发现：
    AINV*A=I（单位阵），但A*AINV≠I
不知是怎么回事？计算过程中并没有警告出现。

程序如下：

        REAL*8 A(80,80),AINV(80,80),SUM
        OPEN(2,FILE='G1.TXT')
        N=80
        READ(2,*) ((A(I,J),I=1,N),J=1,N)
        CALL DLINRG (N, A, N, AINV, N)
        AINV=MATMUL(A,AINV)
        SUM=0.0d0
        DO I=1,N
              SUM=SUM+DABS(A(I,I)-1.0d0)
        ENDDO
        WRITE(*,*) "SUM |A(i,i)-1| =",SUM
        SUM=0.0d0
        DO 2,I=1,N
        DO 2,J=1,N
               IF(I.EQ.J) GOTO 2
               SUM=SUM+DABS(A(I,J))
2       CONTINUE
        WRITE(*,*) "SUM |A(i,j)| =",SUM
        CLOSE(2)
        END

运行结果：
SUM |A(i,i)-1| =  1.973305759124577E-003
SUM |A(i,j)| =    2562.604705789780000
Press any key to continue

SUM |A(i,j)| 不为零。输入文件G1.TXT

jinguangyang

楼主的程序我看了半天，感觉是不是有一个笔误，AINV=MATMUL(A,AINV)
应改为A=MATMUL(A,AINV)，但楼主所提出的问题依然存在，我的计算结果是
SUM |A(i,i)-1| =  4.798695419008236E-004
SUM |A(i,j)| =   347.968499911844
Press any key to continue

GWinston

感谢jinguangyang！
是我写错了，应该是A=MATMUL(A,AINV)
我用的是POWER STATION4.0。
可能矩阵太病态了吧，不过颠倒过来：A=MATMUL(AINV,A)，运算结果是

SUM |A(i,i)-1| =  3.274622591531440E-006
SUM |A(i,j)| =  7.778592744525364E-006
Press any key to continue

jinguangyang

是的，我的结果和你的运行结果是一样的，我用的是cvf6.6，我还自己编了一下矩阵相乘的子程序，计算结果与MATMUL的计算结果完全相同，这就说明了不是MATMUL精度不够的问题。这个问题挺有趣的，大家应该多探讨一下。

jacobi

将你的输入矩阵输出生成G2.TXT，导入matlab，并求其条件数得：
cond(G2)＝4.447557796477102e+009
求其逆的条件数
cond(inv(G2))＝ 4.447557609789049e+009
次条件数均忒大！！！计算结果一般不可信
两种乘法得：
>> cond(G2*inv(G2))
ans =
    1.170700575612592e+002
>> cond(inv(G2)*G2)
ans =
   1.00000040544797
可见:inv(G2)*G2=1.0
但是，G2*inv(G2)不对。

ps：matlab输入
>> format long
>> load G2.txt;
>> a=G2;
>> cond(a)
ans =
    4.447557796477102e+009
>> cond(inv(a))
ans =
    4.447557609789049e+009
>> cond(inv(a)*a)
ans =
   1.00000040544797
>> cond(a*inv(a))
ans =
    1.170700575612592e+002
>>

偶自己做了一个矩阵，求其逆后，分别进行AINV*A与A*AINV计算，结果无误。
platform VF6.6C and coded in free format
ps：源文件如下。

jacobi

G2：

jacobi

G2：

GWinston

谢谢jacobi和jinguangyang的讨论！
其实矩阵G1中含有惩罚因子，该因子理论上应趋于无穷大，我一般取1.0d5左右，所以G1肯定呈病态，无法避免。
DLINRG是通过LU分解来求逆的，现在我程序中用的是从一本书中查到的求逆程序，用的是全主元高斯-约当消去法，计算结果还较满意。

雷雨9

看了你们的讨论，很受益，谢谢拉！

BG

原帖由 GWinston 于 2005-5-8 13:32 发表
谢谢jacobi和jinguangyang的讨论！
其实矩阵G1中含有惩罚因子，该因子理论上应趋于无穷大，我一般取1.0d5左右，所以G1肯定呈病态，无法避免。
DLINRG是通过LU分解来求逆的，现在我程序中用的是从一本书中查到的求 ...

我感觉对于一些相对病态的矩阵，采用高斯约当消元法，甚至传统的迭代法都比一些所谓高等的算法的适用性要强，我曾经就遇到过一个方程求解问题采用MATLAB中的共轭梯度法无法求解，但是自己编的迭代法程序可以计算，所以有时候选用商业数学库中的子程序最好在实际应用前作一些数值实验测试一下适用性，如果盲目采用的化很可能造成最后结论错误而且自己根本就不知道的问题。