[Fortran] 纯文本查看 复制代码
program main22
IMPLICIT REAL*8(A-H,O-Z)
IMPLICIT INTEGER*4(I-N)
dimension a(7,7)
INTEGER, DIMENSION(7,1)::MAXA(7,1)=RESHAPE((/1,2,4,7,9,10,14/),(/7,1/))
DIMENSION V(7,1)
open (1,file='E\1.txt')
DO i=1,7
read(1,*)(a(i,j),j=1,7)
END DO
write(*,100)((a(i,j),j=1,7),i=1,7)
100 FORMAT (2x,7f6.2/6(2x,7f6.2/))
WRITE(*,*)'INPUT THE VALUE OF V'
DO i=1,7
READ(*,*)V(7,1)
END DO
write(*,100)(V(i,1),i=1,7)
call LDLT(a,MAXA,7,0,1,NWK,8)
call RESOLVE(a,V,MAXA,7,NWK,8)
END
SUBROUTINE LDLT(A,MAXA,NN,ISH,IOUT,NWK,NNM)
IMPLICIT REAL*8(A-H,O-Z)
IMPLICIT INTEGER*4(I-N)
DIMENSION A(NWK), MAXA(NNM)
IF(NN.EQ.1)RETURN
DO 200 N=1,NN
KN=MAXA(N)
KL=KL+1
KU=MAXA(N+1)-1
KH=KU-KL
IF(KH)304,240,210
210 K=N-KH
IC=0
KLT=KU
DO 260 J=1,KH
KLT=KLT-1
IC=IC+1
KI=MAXA(K)
ND=MAXA(K+1)-KI-1
IF(ND)260,260,270
270 KK=MINO(IC,ND)
C=0.0
DO 280 L=1,KK
280 C=C+A(KI+L)*A(KLT+L)
A(KLT)=A(KLT)-C
260 K=K+1
240 K=N
B=0.0
DO 300 KK=KL,KU
K=K-1
KI=MAXA(K)
C=A(KK)/A(KI)
IF(ABS(C).LT.1.0E+07) GOTO 290
WRITE(IOUT,2010) N,C
STOP
290 B=B+C*A(KK)
300 A(KK)=C
A(KN)=A(KN)-B
304 IF(A(KN)) 310,310,200
310 IF(ISH.EQ.0) GOTO 320
IF(A(KN).EQ.0.0) A(KN)=-1.0E-16
GOTO 200
320 WRITE(IOUT,2000) N,A(KN)
STOP
200 CONTINUE
RETURN
2000 FORMAT(//'Stop-stiffness matrix not positive definite'//'nonpositive pivot for equation', I4//'pivot= ',E20.10)
2010 FORMAT(//'Stop-sturm sequence check failed because of multiplier growth for column number', I4,//'Multiplier= ',E20.8)
END
SUBROUTINE RESOLVE(A,V,MAXA,NN,NWK,NNM)
IMPLICIT REAL*8(A-H,O-Z)
IMPLICIT INTEGER*4(I-N)
DIMENSION A(NWK),V(NN,1),MAXA(NNM)
NIP=1
DO IP=1,NIP
DO 400 N=1,NN
KL=MAXA(N)+1
KU=MAXA(N+1)-1
IF(KU-KL)400,410,410
410 K=N
C=0.0
DO 420 KK=KL,KU
K=K-1
420 C=C+A(KK)*V(I,IP)
V(N,IP)=V(N,IP)-C
400 CONTINUE
DO 480 N=1,NN
K=MAXA(N)
480 V(N,IP)=V(N,IP)-C
IF(NN.EQ.1)RETURN
N=NN
DO 500 L=2,NN
KL=MAXA(N)+1
KU=MAXA(N+1)-1
IF(KU-KL)500,510,510
510 K=N
DO 520 KK=KL,KU
K=K-1
520 V(K,IP)=V(K,IP)-A(KK)*V(N,IP)
500 N=N-1
ENDDO
WRITE(*,*)(V(I,1),I=1,7)
END