Line data Source code
1 4013 : subroutine leq(a,b,nn,mm,ia,ib,err)
2 : c
3 : c this routine will find the inverse of a matrix by the method of
4 : c partial pivoting and gaussian elimination if b is set equal to
5 : c the identity matrix
6 : c
7 : c nn - dimension of segment of a to be used
8 : c mm - number of right hand columns of b to be used
9 : c ia - the total number of rows in large array a
10 : c ib - the total number of rows in large array b
11 : c the matrix equation is ax=b
12 : c err = det(a)
13 : c if mm = 0 leq calculates err = det(a)
14 : c
15 : c note on modification on 23/1 1985:
16 : c
17 : c previously, the routine contained the statements
18 : c
19 : c do 14 k=i2,m1,ib
20 : c b(k)=b(k)+b(i1)*r
21 : c 14 i1=i1+ib
22 : c
23 : c this caused problems, on some computers, for m = ib = 1.
24 : c then m1 = 1 and the loop was skipped when i2 .gt. 1.
25 : c this has been changed.
26 : c
27 : c Double precision version.
28 : c ++++++++++++++++++++++++
29 : c
30 : c Dated 10/3/90.
31 : c
32 : c Note: this double precision version of the routine has same name
33 : c as single precision version, but is distinguished by its file name
34 : c
35 : implicit double precision(a-h,o-z)
36 : dimension a(100),b(100)
37 : c
38 : c common defining standard input and output
39 : c
40 : common/cstdio/ istdin, istdou, istdpr, istder
41 4013 : 10000 n=nn
42 4013 : m=mm
43 4013 : err=0.0d0
44 4013 : detsc=1.0d0
45 : c
46 : c treat n .le. 1 separately
47 4013 : if(n-1) 280,281,285
48 : c no equation
49 0 : 280 if(istdpr.gt.0) write(istdpr,295) n
50 0 : return
51 : c n = 1
52 0 : 281 err=a(1)
53 0 : if(m.le.0) return
54 0 : ai=1.d0/err
55 0 : m1=ib*m
56 0 : do 282 j=1,m1,ib
57 0 : 282 b(j)=b(j)*ai
58 4013 : return
59 : c
60 : c
61 : c find maximum element in each row and divide the row by it
62 4013 : 285 n1=ia*n
63 4013 : ia1=ia+1
64 4013 : m1=ib*m
65 21487 : do 1 i=1,n
66 17474 : r= abs(a(i))
67 17474 : do 2 j= 1,n1,ia
68 77098 : ij=j+i-1
69 77098 : 2 r=max(r, abs(a(ij)))
70 17474 : if(r)31,30,31
71 0 : 30 if(istdpr.gt.0) write(istdpr,298)i
72 17474 : return
73 17474 : 31 do 3 j=1,n1,ia
74 77098 : ij=j+i-1
75 77098 : 3 a(ij)=a(ij)/r
76 17474 : if(m.eq.0) go to 1
77 15522 : do 4 j=1,m1,ib
78 15522 : ij=j+i-1
79 15522 : 4 b(ij)=b(ij)/r
80 21487 : 1 detsc=detsc*r
81 : c
82 : c
83 : c find maximum element in the i'th column
84 4013 : n2=n-1
85 17474 : do 5 i=1,n2
86 13461 : ialow=(i-1)*ia+i
87 13461 : iaup=(i-1)*ia+n
88 13461 : r= abs(a(ialow))
89 13461 : ialow1=ialow+1
90 13461 : imax=ialow
91 43273 : do 6 j=ialow1,iaup
92 29812 : if(r- abs(a(j)))7,6,6
93 7525 : 7 imax=j
94 29812 : r= abs(a(j))
95 13461 : 6 continue
96 13461 : if(imax-ialow)8,8,9
97 : c replace the i'th row with the row that has the maximum element in
98 : c the respective column and put the i'th row in its place
99 5036 : 9 im=imax
100 13479 : 72 if(im-ia)70,70,71
101 8443 : 71 im=im-ia
102 13479 : go to 72
103 5036 : 70 do 10 j=1,n1,ia
104 23628 : jj=i+j-1
105 23628 : ji=im+j-1
106 23628 : r=a(jj)
107 23628 : a(jj)=a(ji)
108 23628 : 10 a(ji)=r
109 : c change sign of determinant
110 5036 : detsc=-detsc
111 : c
112 5036 : if(m.eq.0) go to 8
113 4366 : do 11 j=1,m1,ib
114 4366 : jj=i+j-1
115 4366 : ji=im+j-1
116 4366 : r=b(jj)
117 4366 : b(jj)=b(ji)
118 12791 : 11 b(ji)=r
119 : c multiply the i'th row by (the negative of each i'th column element
120 : c below the diagonal divided by the diagonal element) and add the
121 : c resulting row to the respective row of the element used
122 13461 : 8 iaup1=iaup-1
123 : c
124 : c
125 43273 : do 12 j=ialow,iaup1
126 29812 : if(a(ialow))32,33,32
127 0 : 33 joy=i
128 0 : if(a(ialow1))81,82,81
129 0 : 82 if(istdpr.gt.0) write(istdpr,299)joy,joy
130 0 : return
131 0 : 81 if(istdpr.gt.0) write(istdpr,297)joy,joy
132 0 : do 34 k=1,n1,ia
133 0 : jj=joy+k-1
134 0 : ji=joy+k
135 0 : if(joy+1-n)35,36,36
136 0 : 35 if(istdpr.gt.0) write(istdpr,296)
137 0 : return
138 0 : 36 r=a(jj)
139 0 : a(jj)=a(ji)
140 0 : 34 a(ji)=r
141 : c change sign of determinant
142 0 : detsc=-detsc
143 : c
144 0 : if(m.eq.0) go to 8
145 0 : do 37 k=1,m1,ib
146 0 : jj=joy+k-1
147 0 : ji=joy+k
148 0 : r=b(jj)
149 0 : b(jj)=b(ji)
150 0 : 37 b(ji)=r
151 29812 : go to 8
152 29812 : 32 j1=j+1
153 29812 : r=-a(j1)/a(ialow)
154 29812 : i1=ialow
155 29812 : do 13 k=j1,n1,ia
156 109068 : a(k)=a(k)+a(i1)*r
157 109068 : 13 i1=i1+ia
158 : c
159 : c loop to reset b has been modified, 25/1/1985.
160 : c
161 29812 : if(m.eq.0) go to 12
162 26884 : i1=i
163 26884 : i2=j-ialow+i+1
164 26884 : do 14 k=1,m1,ib
165 26884 : b(i2)=b(i2)+b(i1)*r
166 26884 : i1=i1+ib
167 26884 : 14 i2=i2+ib
168 13461 : 12 continue
169 : c
170 : c
171 4013 : 5 continue
172 : c
173 : c
174 : c the matrix is now in triangular form
175 : c first set err=1.0d0
176 4013 : err=1.0d0
177 : c if(any diagonal element of a is zero x cannot be solved for
178 21487 : do 15 i=1,n
179 17474 : idiag=(i-1)*ia+i
180 17474 : err=err*a(idiag)
181 17474 : if(err) 15,16,15
182 0 : 16 if(istdpr.gt.0) write(istdpr,299)i,i
183 17474 : return
184 4013 : 15 continue
185 : c scale determinant
186 4013 : err=err*detsc
187 : c
188 4013 : if(m.eq.0) return
189 : c find solution to ax=b
190 7050 : do 18 k=1,m
191 3525 : ka=(n-1)*ia+n
192 3525 : kb=(k-1)*ib+n
193 3525 : b(kb)=b(kb)/a(ka)
194 15522 : do 19 l=1,n2
195 11997 : i=n-l
196 11997 : r=0.0d0
197 11997 : imax=i+1
198 38881 : do 20 j=imax,n
199 26884 : jj=i+n+1-j
200 26884 : ja=(jj-1)*ia+i
201 26884 : jb=(k-1)*ib+jj
202 38881 : 20 r=r+a(ja)*b(jb)
203 11997 : la=(i-1)*ia+i
204 11997 : lb=(k-1)*ib+i
205 15522 : 19 b(lb)=(b(lb)-r)/a(la)
206 3525 : 18 continue
207 : c
208 : c
209 3525 : return
210 : c
211 : 295 format(///20h leq called with n =,i4)
212 : 296 format(///48h the row cannot be changed with the row below it ,
213 : .40h because it it the last row in the array /
214 : .55h the solution, matrix x, cannot be found by this method )
215 : 297 format(5h1 a(,i4,1h,,i4,13h) equals zero /
216 : .47h try switching this row with the row below it ,
217 : .48h and go back to statement number 8 and try again )
218 : 298 format(26h1 all the elements in row,i5,20h are zero therefore,
219 : .55h the solution, matrix x, cannot be found by this method )
220 : 299 format(50h1 the solution, matrix x, cannot be found by this ,
221 : .57h method because there is a zero array element in the main ,
222 : .9h diagonal / 30x,2ha(,i4,1h,,i4,8h) = zero )
223 : end
224 0 : subroutine leqsd(a,b,n,m,ia,ib,det,isa,isb)
225 : implicit double precision (a-h,o-z)
226 : dimension a(ia,1),b(ib,1)
227 : c
228 0 : call leq(a,b,n,m,ia,ib,det)
229 0 : return
230 : end
|
