1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
| PROGRAM MULTISCALE_1
IMPLICIT NONE
INTEGER A,N,INDX,X
C
C n is the number of row minus 1 in the data file (the first line contains n)
C
OPEN (22,file='Cijkl.dat',form='formatted',status='old')
READ (22,*) A
CALL MIGS(A,N,X,INDX)
PRINT *,X
END
SUBROUTINE MIGS(A,N,X,INDX)
C
C Subroutine to invert matrix A(N,N) with the inverse stored
C in X(N,N) in the output.
C
DIMENSION A(N,N),X(N,N),INDX(N),B(N,N)
C
DO 20 I = 1, N
DO 10 J = 1, N
B(I,J) = 0.0
10 CONTINUE
20 CONTINUE
DO 30 I = 1, N
B(I,I) = 1.0
30 CONTINUE
C
CALL ELGS(A,N,INDX)
C
DO 100 I = 1, N-1
DO 90 J = I+1, N
DO 80 K = 1, N
B(INDX(J),K) = B(INDX(J),K)
* -A(INDX(J),I)*B(INDX(I),K)
80 CONTINUE
90 CONTINUE
100 CONTINUE
C
DO 200 I = 1, N
X(N,I) = B(INDX(N),I)/A(INDX(N),N)
DO 190 J = N-1, 1, -1
X(J,I) = B(INDX(J),I)
DO 180 K = J+1, N
X(J,I) = X(J,I)-A(INDX(J),K)*X(K,I)
180 CONTINUE
X(J,I) = X(J,I)/A(INDX(J),J)
190 CONTINUE
200 CONTINUE
C
RETURN
END
C
SUBROUTINE ELGS(A,N,INDX)
C
C Subroutine to perform the partial-pivoting Gaussian elimination.
C A(N,N) is the original matrix in the input and transformed
C matrix plus the pivoting element ratios below the diagonal in
C the output. INDX(N) records the pivoting order.
C
DIMENSION A(N,N),INDX(N),C(N)
C
C Initialize the index
C
DO 50 I = 1, N
INDX(I) = I
50 CONTINUE
C
C Find the rescaling factors, one from each row
C
DO 100 I = 1, N
C1= 0.0
DO 90 J = 1, N
C1 = AMAX1(C1,ABS(A(I,J)))
90 CONTINUE
C(I) = C1
100 CONTINUE
C
C Search the pivoting (largest) element from each column
C
DO 200 J = 1, N-1
PI1 = 0.0
DO 150 I = J, N
PI = ABS(A(INDX(I),J))/C(INDX(I))
IF (PI.GT.PI1) THEN
PI1 = PI
K = I
ELSE
ENDIF
150 CONTINUE
C
C Interchange the rows via INDX(N) to record pivoting order
C
ITMP = INDX(J)
INDX(J) = INDX(K)
INDX(K) = ITMP
DO 170 I = J+1, N
PJ = A(INDX(I),J)/A(INDX(J),J)
C
C Record pivoting ratios below the diagonal
C
A(INDX(I),J) = PJ
C
C Modify other elements accordingly
C
DO 160 K = J+1, N
A(INDX(I),K) = A(INDX(I),K)-PJ*A(INDX(J),K)
160 CONTINUE
170 CONTINUE
200 CONTINUE
C
RETURN
END |
