main
|
Last change
on this file was ea777aa, checked in by Alex Wilton <awilton@…>, 3 years ago |
|
Moved examples, include, build_default.properties, common.xml, and README out from dev.civl.com into the root of the repo.
git-svn-id: svn://vsl.cis.udel.edu/civl/trunk@5704 fb995dde-84ed-4084-dfe6-e5aef3e2452c
|
-
Property mode
set to
100755
|
|
File size:
1.9 KB
|
| Line | |
|---|
| 1 | SUBROUTINE LU(A,N,ldim,IR,IC)
|
|---|
| 2 | C IT IS THE FIRST SUBROUTINE TO COMPUTE THE MX. INV.
|
|---|
| 3 | DIMENSION A(ldim,1),IR(1),IC(1)
|
|---|
| 4 | DO 10I=1,N
|
|---|
| 5 | IR(I)=I
|
|---|
| 6 | IC(I)=I
|
|---|
| 7 | 10 CONTINUE
|
|---|
| 8 | K=1
|
|---|
| 9 | L=K
|
|---|
| 10 | M=K
|
|---|
| 11 | XMAX=ABS(A(K,K))
|
|---|
| 12 | DO 100I=K,N
|
|---|
| 13 | DO 100J=K,N
|
|---|
| 14 | Y=ABS(A(I,J))
|
|---|
| 15 | IF(XMAX.GE.Y) GOTO 100
|
|---|
| 16 | XMAX=Y
|
|---|
| 17 | L=I
|
|---|
| 18 | M=J
|
|---|
| 19 | 100 CONTINUE
|
|---|
| 20 | DO 1000K=1,N
|
|---|
| 21 | IRL=IR(L)
|
|---|
| 22 | IR(L)=IR(K)
|
|---|
| 23 | IR(K)=IRL
|
|---|
| 24 | ICM=IC(M)
|
|---|
| 25 | IC(M)=IC(K)
|
|---|
| 26 | IC(K)=ICM
|
|---|
| 27 | IF(L.EQ.K) GOTO 300
|
|---|
| 28 | DO 200J=1,N
|
|---|
| 29 | B=A(K,J)
|
|---|
| 30 | A(K,J)=A(L,J)
|
|---|
| 31 | A(L,J)=B
|
|---|
| 32 | 200 CONTINUE
|
|---|
| 33 | 300 IF(M.EQ.K) GOTO 500
|
|---|
| 34 | DO 400I=1,N
|
|---|
| 35 | B=A(I,K)
|
|---|
| 36 | A(I,K)=A(I,M)
|
|---|
| 37 | A(I,M)=B
|
|---|
| 38 | 400 CONTINUE
|
|---|
| 39 | 500 C=1./A(K,K)
|
|---|
| 40 | A(K,K)=C
|
|---|
| 41 | IF(K.EQ.N) GOTO 1000
|
|---|
| 42 | K1=K+1
|
|---|
| 43 | XMAX=ABS(A(K1,K1))
|
|---|
| 44 | L=K1
|
|---|
| 45 | M=K1
|
|---|
| 46 | DO 600I=K1,N
|
|---|
| 47 | A(I,K)=C*A(I,K)
|
|---|
| 48 | 600 CONTINUE
|
|---|
| 49 | DO 800I=K1,N
|
|---|
| 50 | B=A(I,K)
|
|---|
| 51 | DO 800J=K1,N
|
|---|
| 52 | A(I,J)=A(I,J)-B*A(K,J)
|
|---|
| 53 | Y=ABS(A(I,J))
|
|---|
| 54 | IF(XMAX.GE.Y) GOTO 800
|
|---|
| 55 | XMAX=Y
|
|---|
| 56 | L=I
|
|---|
| 57 | M=J
|
|---|
| 58 | 800 CONTINUE
|
|---|
| 59 | 1000 CONTINUE
|
|---|
| 60 | RETURN
|
|---|
| 61 | END
|
|---|
| 62 | SUBROUTINE SOLVE(F,A,K,N,ldim,IR,IC)
|
|---|
| 63 | C IT IS THE SECOND PART OF THE MATRIX INVERSION
|
|---|
| 64 | DIMENSION A(ldim,1),F(ldim,1),IR(1),IC(1)
|
|---|
| 65 | COMMON /CTMPG/ G(2000)
|
|---|
| 66 | C
|
|---|
| 67 | C
|
|---|
| 68 | IF (N.GT.2000) THEN
|
|---|
| 69 | write(6,*) 'Abort IN Subrtouine SOLVE: N>2000, N=',N
|
|---|
| 70 | call exitt
|
|---|
| 71 | ENDIF
|
|---|
| 72 | C
|
|---|
| 73 | N1=N+1
|
|---|
| 74 | DO 1000KK=1,K
|
|---|
| 75 | DO 100I=1,N
|
|---|
| 76 | IRI=IR(I)
|
|---|
| 77 | G(I)=F(IRI,KK)
|
|---|
| 78 | 100 CONTINUE
|
|---|
| 79 | DO 400I=2,N
|
|---|
| 80 | I1=I-1
|
|---|
| 81 | B=G(I)
|
|---|
| 82 | DO 300J=1,I1
|
|---|
| 83 | B=B-A(I,J)*G(J)
|
|---|
| 84 | 300 CONTINUE
|
|---|
| 85 | G(I)=B
|
|---|
| 86 | 400 CONTINUE
|
|---|
| 87 | DO 700IT=1,N
|
|---|
| 88 | I=N1-IT
|
|---|
| 89 | I1=I+1
|
|---|
| 90 | B=G(I)
|
|---|
| 91 | IF(I.EQ.N) GOTO 701
|
|---|
| 92 | DO 600J=I1,N
|
|---|
| 93 | B=B-A(I,J)*G(J)
|
|---|
| 94 | 600 CONTINUE
|
|---|
| 95 | 701 G(I)=B*A(I,I)
|
|---|
| 96 | 700 CONTINUE
|
|---|
| 97 | DO 900I=1,N
|
|---|
| 98 | ICI=IC(I)
|
|---|
| 99 | F(ICI,KK)=G(I)
|
|---|
| 100 | 900 CONTINUE
|
|---|
| 101 | 1000 CONTINUE
|
|---|
| 102 | RETURN
|
|---|
| 103 | END
|
|---|
Note:
See
TracBrowser
for help on using the repository browser.
