implicit real*8 (a-h,o-z)
character*12 fname,filename(20)
common /aa/ ia(15750000)
common /bb/ ib(7870000)
common /cc/ ic(3930000)
open(1,file=' ',form='unformatted',status='old')
read(1) knode,kdgof
close(1)
MAXT=15750000/2

C.......OPEN SYS File
OPEN (2,FILE=' ',FORM='UNFORMATTED',STATUS='OLD')
read(2) NUMEL,NEQ,NUMTYP,MAXA
CLOSE (2)

IF (MAXA.GT.MAXT) THEN
WRITE(*,*) 'MATRIX A EXCEED CORE MEMERY .... ',MAXA
WRITE(*,*) 'REQUIRED CORE MEMERY ........... ',MAXT
STOP 0000
ENDIF
KVAR=KNODE*KDGOF
KCOOR=3
KELEM=500000
WRITE(*,*) 'KNODE,KDGOF,KVAR,KCOOR,KELEM ='
WRITE(*,'(1X,6I7)') KNODE,KDGOF,KVAR,KCOOR,KELEM

kna1=kdgof*knode*1
if (kna1/2*2 .lt. kna1) kna1=kna1+1
knc3=kdgof*knode*2
knc1=kcoor*knode*2
knc6=kdgof*knode*2
knc2=neq*2
knb1=maxa*2
kna2=neq*1
if (kna2/2*2 .lt. kna2) kna2=kna2+1
kna3=kelem*1
if (kna3/2*2 .lt. kna3) kna3=kna3+1
knc7=100000*2
knc5=neq*2
knc4=kdgof*knode*2
kna0=1
kna1=kna1+kna0
kna2=kna2+kna1
kna3=kna3+kna2
knb0=1
knb1=knb1+knb0
knc0=1
knc1=knc1+knc0
knc2=knc2+knc1
knc3=knc3+knc2
knc4=knc4+knc3
knc5=knc5+knc4
knc6=knc6+knc5
knc7=knc7+knc6
call ebara(knode,kdgof,kvar,kcoor,numtyp,numel,
*neq,kelem,maxa,maxt,neq1,ib(kna0),ib(kna1),
*ib(kna2),ia(knb0),ic(knc0),ic(knc1),ic(knc2),
*ic(knc3),ic(knc4),ic(knc5),ic(knc6),
*filename)
end
subroutine ebara(knode,kdgof,kvar,kcoor,numtyp,numel,
*neq,kelem,maxa,maxt,neq1,nodvar,jdiag,node,
*a,coor,f,u,ubf,u1,eu,sml,
*filename)
implicit real*8 (a-h,o-z)
character*12 filename(20)
DIMENSION NODVAR(KDGOF,KNODE),U(KDGOF,KNODE),COOR(KCOOR,KNODE),
*eu(kdgof,knode),
& F(NEQ),A(MAXA),JDIAG(NEQ),EMATE(3000),
& NODE(KELEM),SML(100000),u1(neq),UBF(KDGOF,KNODE)

6 FORMAT (1X,15I5)
7 FORMAT (1X,5e15.5)
1001 FORMAT(1X,9I7)

C.......OPEN TIME File
OPEN(1,FILE=' ',FORM='UNFORMATTED',STATUS='OLD')
READ(1) TMAX,DT,TIME,IT
WRITE(*,*) ' TMAX,DT,TIME,IT =',TMAX,DT,TIME,IT
CLOSE(1)

C.......OPEN NODVAR file
OPEN (1,FILE=' ',FORM='UNFORMATTED',STATUS='OLD')
READ (1) ((NODVAR(I,J),I=1,KDGOF),J=1,KNODE)
CLOSE (1)
cc WRITE(*,*) 'KDGOF =',KDGOF,' KNODE =',KNODE
cc WRITE (*,*) 'NODVAR ='
cc WRITE (*,6) ((NODVAR(I,J),I=1,KDGOF),J=1,KNODE)

C.......OPEN COOR file
OPEN (1,FILE=' ',FORM='UNFORMATTED',STATUS='OLD')
READ (1) NUMNOD,NCOOR,((COOR(I,J),I=1,NCOOR),J=1,NUMNOD)
CLOSE(1)
cc WRITE(*,*) 'NUMNOD,NCOOR=',NUMNOD,NCOOR

C.......OPEN BF file
OPEN (1,FILE=' ',FORM='UNFORMATTED',STATUS='OLD')
READ (1) ((UBF(J,I),J=1,KDGOF),I=1,KNODE)
CLOSE (1)
cc WRITE (*,*) 'BF ='
cc WRITE(*,7) ((U(J,I),J=1,KDGOF),I=1,KNODE)

numtyp = 1
IF (IT.EQ.0) THEN
ELSE
ENDIF

C.......OPEN DIAG file
OPEN (2,FILE=' ',FORM='UNFORMATTED',STATUS='OLD')
READ(2) (JDIAG(I),I=1,NEQ)
CLOSE(2)

C.......OPEN ELEM0 file
OPEN (3,FILE=' ',FORM='UNFORMATTED',STATUS='OLD')

itime=0
1 continue
itime=itime+1
if (itime.gt.1) then
write(*,*) 'Nonlinear Iteration Times ========',itime
rewind(3)
endif

DO 111 I=1,KNODE
DO 111 J=1,KDGOF
U(J,I) = UBF(J,I)
111 CONTINUE
cc WRITE (*,*) 'BF ='
cc WRITE(*,7) ((U(J,I),J=1,KDGOF),I=1,KNODE)
DO 112 I=1,MAXA
A(I) = 0.0
112 CONTINUE

DO 2300 I=1,NEQ
2300 CONTINUE
NUMEL=0
C.......OPEN EMATE+ENODE+ELOAD file
C OPEN (3,FILE=' ',FORM='UNFORMATTED',STATUS='OLD')
DO 2000 ITYP=1,NUMTYP
C.......INPUT ENODE
READ (3) NUM,NNODE,
* ((NODE((I-1)*NNODE+J),J=1,NNODE),I=1,NUM)
cc WRITE(*,*) 'NUM =',NUM,' NNODE =',NNODE
cc WRITE(*,*) 'NODE ='
cc WRITE(*,6) ((NODE((I-1)*NNODE+J),J=1,NNODE),I=1,NUM)
NNE = NNODE
nne = nne-1
K=0
DO 115 J=1,NNE
JNOD = NODE(J)
DO 115 L=1,KDGOF
IF (NODVAR(L,JNOD).NE.0) K=K+1
115 CONTINUE
WRITE(*,*) 'K =',K
kk=k*k
k0=1
k1=k0+k*k
k2=k1+k
k3=k2+k
k4=k3+k*k
k5=k4+k*k
CALL ETSUB(KNODE,KDGOF,IT,KCOOR,KELEM,K,KK,NNODE,NNE,
* ITYP,NCOOR,NUM,TIME,DT,neq,maxa,NODVAR,COOR,NODE,EMATE,
& A,JDIAG,
*sml(k0),sml(k1),sml(k2),sml(k3),sml(k4),
&eu,
*U)
2000 CONTINUE

C CLOSE(1)
C CLOSE(2)
c CLOSE(3)
C NEQ = 0
DO 2050 IJ=1,NEQ
if (itime.le.1) u1(IJ) = 0.0
F(IJ)=0.0D0
2050 CONTINUE
DO 2200 I=1,KNODE
DO 2100 J=1,KDGOF
IJ=NODVAR(J,I)
IF (IJ.LE.0) GOTO 2100
C IF (IJ.GT.NEQ) NEQ = IJ
F(IJ)=F(IJ)+U(J,I)
U1(IJ)=F(IJ)
2100 CONTINUE
2200 CONTINUE

CC IF (IT.GT.0) THEN
cc WRITE (*,*) 'U ='
cc WRITE (*,7) ((U(J,I),J=1,KDGOF),I=1,KNODE)
cc WRITE (*,*) 'NEQ =',NEQ,' F ='
cc WRITE(*,7) (F(I),I=1,NEQ)

if (itime.le.1) then
C.......OPEN LMATRIX FILE
OPEN (2,FILE=' ',FORM='UNFORMATTED',STATUS='unknown')
CLOSE (2)
endif

WRITE(*,*) 'SIN_SOLVER .... MAXA =',MAXA
c IF (MAXA.GT.MAXT) THEN
c WRITE(*,*) 'MATRIX A EXCEED CORE MEMORY'
c WRITE(*,*) 'CORE MEMORY = ',MAXT
c STOP 0000
c ENDIF
CALL REDU(A,U1,JDIAG,NEQ,MAXA,1)
C SUBROUTINE REDU(A,B,U,JDIAG,NEQ,MAXA,KKK)
C WRITE(*,*) ' U1 = '
C WRITE(*,7) (A(I),I,MAXA)
C WRITE(*,7) (F(I),I=1,NEQ)

NOUT = 20
OPEN(NOUT,FILE=' ',FORM='FORMATTED',STATUS='unknown')
DO 3200 INOD=1,KNODE
DO 3100 IDFG=1,KDGOF
N=NODVAR(IDFG,INOD)
C WRITE (*,*) 'N =',N
cc IF (N.eq.0) GOTO 3100
if(n.le.0) then
eu(IDFG,INOD)=u(IDFG,INOD)
else
eu(IDFG,INOD)=u1(N)
endif
C WRITE (*,*) 'F(N) =',F(N)
3100 CONTINUE
3200 CONTINUE
DO 3400 N=1,KNODE
WRITE (NOUT,3600) N,(eu(I,N),I=1,KDGOF)
3400 CONTINUE
3600 FORMAT (1X,I5,1X,6E11.4,9(/6X,6E11.4))
CLOSE (NOUT)

open(10,file='unod',form='unformatted',status='unknown')
write(10) ((eu(j,i),i=1,knode),j=1,kdgof)
close(10)
CLOSE(3)

RETURN
END

SUBROUTINE ETSUB(KNODE,KDGOF,IT,KCOOR,KELEM,K,KK,NNODE,NNE,
*ITYP,NCOOR,NUM,TIME,DT,neq,maxa,NODVAR,COOR,NODE,EMATE,
&A,JDIAG,
*es,em,ef,Estifn,Estifv,eu,
*U)
implicit real*8 (a-h,o-z)
DIMENSION NODVAR(KDGOF,KNODE),COOR(KCOOR,KNODE),NODE(KELEM),
*U(KDGOF,KNODE),EMATE(300),
&A(MAXa),JDIAG(neq),
*es(k,k),em(k),ef(k),eu(kdgof,knode),
*Estifn(k,k),Estifv(kk),
*R(500),PRMT(500),COEF(500),LM(500)
17 FORMAT (1X,15I5)
18 FORMAT (1X,8e9.2)

READ (3) MMATE,NMATE,((EMATE((I-1)*NMATE+J),J=1,NMATE),
* I=1,MMATE)
WRITE(*,*) 'MMATE =',MMATE,' NMATE =',NMATE
WRITE (*,*) 'EMATE ='
WRITE (*,18) ((EMATE((I-1)*NMATE+J),J=1,NMATE),
* I=1,MMATE)
DO 1000 NE=1,NUM
NR=0
DO 130 J=1,NNE
JNOD = NODE((NE-1)*NNODE+J)
IF (JNOD.LT.0) JNOD = -JNOD
DO 120 I=1,NCOOR
NR=NR+1
120 R(NR) = COOR(I,JNOD)
130 CONTINUE
IMATE = NODE(NNODE*NE)
DO 140 J=1,NMATE
140 PRMT(J) = EMATE((IMATE-1)*NMATE+J)
PRMT(NMATE+1)=TIME
PRMT(NMATE+2)=DT
prmt(nmate+3)=imate
prmt(nmate+4)=num
goto 1
1 call ael2(r,coef,prmt,es,em,ec,ef,ne)
goto 2
2 continue

C WRITE(*,*) 'ES EM EF ='
C DO 555 I=1,K
C555 WRITE(*,18) (ES(I,J),J=1,K)
C WRITE(*,18) (EM(I),I=1,K)
C WRITE(*,18) (EF(I),I=1,K)

CC IF (IT.GT.0) THEN
do 201 i=1,k
do 201 j=1,k
Estifn(i,j)=0.0
201 continue
do 202 i=1,k
Estifn(i,i)=Estifn(i,i)
do 202 j=1,k
Estifn(i,j)=Estifn(i,j)+es(i,j)
202 continue

L=0
M=0
I=0
DO 700 INOD=1,NNE
NODI=NODE((NE-1)*NNODE+INOD)
DO 600 IDGF=1,KDGOF
INV=NODVAR(IDGF,NODI)
IF (INV.EQ.0) GOTO 600
I=I+1
IF (INV.LT.0) GOTO 305
L=L+1
LM(L)=INV
U(IDGF,NODI)=U(IDGF,NODI)
*+ef(i)
305 J=0
DO 500 JNOD=1,NNE
NODJ=NODE((NE-1)*NNODE+JNOD)
DO 400 JDGF=1,KDGOF
JNV=NODVAR(JDGF,NODJ)
IF (JNV.EQ.0) GOTO 400
J=J+1

IF (JNV.LT.0) GOTO 400
IF (INV.LT.0) GOTO 310
M=M+1
Estifv(m)=Estifn(i,j)
310 CONTINUE


IF (INV.LT.0)
* U(JDGF,NODJ)=U(JDGF,NODJ)-ESTIFN(J,I)*U(IDGF,NODI)
400 CONTINUE
500 CONTINUE
600 CONTINUE
700 CONTINUE
C WRITE (*,*) 'U ='
C WRITE (*,18) ((U(J,I),J=1,KDGOF),I=1,KNODE)

LRD=M
NER=NUMEL+NE
CCCC WRITE (1) L,(LM(I),I=1,L)
CCCC#WDIST.sub
C WRITE(*,*) '**************************'
C WRITE(*,*) (ESTIFV(I),I=1,LRD)
C WRITE (*,*) 'Einform ............'
C WRITE (*,'(1X,15I5)') L,LRD,(LM(I),I=1,L)
CC ENDIF
DO 800 I=1,L
J=LM(I)
800 CONTINUE
call ADDA(A,JDIAG,L,LM,ESTIFV,NEQ,MAXA)
1000 CONTINUE
RETURN
END

SUBROUTINE ADDA(A,JDIAG,ND,LM,ESTIF,NEQ,MAXA)
implicit real*8 (a-h,o-z)
DIMENSION A(MAXA),JDIAG(NEQ),LM(ND),ESTIF(ND,ND)
C READ(NT1) ND,(LM(I),I=1,ND)
C READ(NT2) ((ESTIF(I,J),J=1,ND),I=1,ND)
C WRITE (*,*) ND, (LM(I),I=1,ND)
C WRITE (*,*) ((ESTIF(I,J),J=1,ND),I=1,ND)
DO 300 I=1,ND
II = LM(I)
DO 200 J=1,I
JJ = LM(J)
IF (II.LT.JJ) THEN
K = JDIAG(JJ) - JJ + II
ELSE
K = JDIAG(II) - II + JJ
ENDIF
A(K) = A(K) + ESTIF(J,I)
200 CONTINUE
300 CONTINUE
RETURN
END

SUBROUTINE REDU(A,U,JDIAG,NEQ,MAXA,KKK)
implicit real*8 (a-h,o-z)
DIMENSION A(MAXA),JDIAG(NEQ),U(NEQ)
DOUBLE PRECISION C,AA
C.......COMPUTE L D & U, L*D*L(T) = A, A*U = F

DO 500 I=2,NEQ
I1 = I-1
NI = JDIAG(I)
LI = I-NI+JDIAG(I-1)+1
IF (KKK.EQ.2) GOTO 333
DO 200 J=LI,I1
NJ = JDIAG(J)
LJ = J-NJ+1
IF (J.GT.1) LJ = LJ+JDIAG(J-1)
LIJ = MAX0(LI,LJ)
J1 = J-1
AA = 0.0
DO 100 L=LIJ,J1
100 AA = AA + A(NI-I+L)*A(NJ-J+L)
200 A(NI-I+J) = A(NI-I+J) - AA
C = 0.0
DO 300 J=LI,I1
NJ = JDIAG(J)
K = NI-I+J
T = A(K)
A(K) = T/A(NJ)
300 C = C+T*A(K)
A(NI) = A(NI) - C

IF (KKK.EQ.3) GOTO 500
333 C = 0.0
DO 400 L=LI,I1
400 C = C + A(NI-I+L)*U(L)
U(I) = U(I)-C
500 CONTINUE
IF (KKK.EQ.3) GOTO 999

DO 600 I=1,NEQ
K = JDIAG(I)
600 U(I) = U(I)/A(K)

J = MAXA + 1
N = NEQ + 1
700 N = N - 1
J = J - 1
IF (N.EQ.1) GOTO 999
M = J-JDIAG(N-1)-1
DO 800 I=1,M
J = J - 1
800 U(N-I) = U(N-I) - A(J)*U(N)
GOTO 700
999 RETURN
END

【Close】