求解类型的矩阵方程，其中是以对称矩阵。
通过 Bill 给出的解决方案和对因式分解的更多研究，我使用了SH Lee提供的原型代码以及相同的测试代码：
parameter (neqs = 5)
real*8 au(7), ad(neqs), b(neqs), energy
integer jp(neqs)
data au/-2.,-2.,-3.,-1.,0.,0.,4./
data ad/2.,3.,5.,10.,10./
data b/0.,1.,0.,0.,0./
data jp/0,1,2,3,7/
CALL REDUCl(neqs,jp,ad,au)
CALL SOLVE1(neqs,jp,ad,au,b)
结果令人满意：
解 636. 619. 292. 74. 34.
子程序 REDUC1 和 SOLVE1 如下：
SUBROUTINE REDUCl (N,IP,D,A)
IMPLICIT REAL(8) (A-H,O-Z)
DIMENSION IP(*),D(*),A(*)
!FUNCTION : CROUT DECOMPOSITION A = LDU
!REDUCTION COLUMN BY COLUMN
DO 11 K=2,N
K1=K-1
LK=IP(K)-K1
KH=IP(K1)-LK+1
S=D(K)-A(LK+KH)*A(LK+KH)/D(KH)
DO 22 J=KH+1,K1
J1=J-1
LJ=IP(J)-J1
JH=IP(J1)-LJ+1
IF (KH.GT.JH) JH=KH
T=A(LK+J)
DO 33 M=JH,J1
33 T=T-A(LJ+M)*A(LK+M)
A(LK+J)=T
22 S=S-T*T/D(J)
DO 44 J=KH,K1
L=LK+J
44 A(L)=A(L)/D(J)
11 D(K)=S
RETURN
END
SUBROUTINE SOLVE1 (N,IP,D,A,B)
IMPLICIT REAL(8) (A-H,O-Z)
DIMENSION IP(*),D(*),A(*),B(*)
!FUNCTION : SOLVE FOR X, LDUx = b where U = transpose of L
!1. FORWARD SUBSTITUTION : Lz = b
DO 11 J=2,N
J1=J-1
LJ=IP(J)-J1
JH=IP(J1)-LJ+1
T=B(J)
DO 22 M=JH,J1
22 T=T-A(LJ+M)*B(M)
11 B(J)=T
!2. DIVIDING BY DIAGONAL ELEMENTS : Dy = 2
DO 33 K=1,N
33 B(K)=B(K)/D(K)
!BACKWARD SUBSTITUTION : UX = y
DO 44 K=N,2,-1
K1=K-1
LK=IP(K)-K1
KH=IP(K1)-LK+1
T=B(K)
DO 44 J=KH,K1
44 B(J)=B(J)-T*A(LK+J)
RETURN
END
但是，现在我遇到了一个新问题！当我的矩阵与上面类似但在以下两种情况下失败（我得到错误的答案）时，这可以正常工作：
案例 1：
data au/0.6801,0.3846,0.4820,1.4074,1.2731,0.7047,0.1078,0.6171,0.9063,0.5895,0.,0.5260,0.3439,0.8178/
data ad/0.3818,0.4524,0.,1.4605,1.7593,1.1887/
data b/1.,1.,1.,1.,1.,1./
data jp/0,0,2,5,9,14/
CALL REDUCl(neqs,jp,ad,au)
CALL SOLVE1(neqs,jp,ad,au,b)
案例二：
!!C Solve A*U = b
!!C [1 0 0 0 0 0 ]
!!C [0 1 0 0 0 0 ]
!!C A = [0 0 1 0 0 0 ]
!! [0 0 0 1 0 0 ]
!!C [0 0 0 1.6667E10 -0.4167E10 ]
!!C [0 0 0 0 -0.4167E10 0.1041e10 ]
data au/-0.4166667E10/
data ad/1.,1.,1.,1.,1.6667E10,0.1041667E10/
data b/1.,1.,1.,1.,1.,1./
data jp/0,0,0,0,0,1/
Print*,au
Print*,ad
CALL REDUCl(neqs,jp,ad,au)
CALL SOLVE1(neqs,jp,ad,au,b)
我哪里错了？是由于某些分配或使用精度造成的吗？