1 2!***** Double precision linear algebra ******************************* 3 4 subroutine loupd(x,n,nm,iperm,nperm) 5! Lower-upper decomposition of matrix x, logical size n x n, 6! stored as size nm x nm. Returns matrix x(i,j),i<=j as upper 7! triangular matrix and x(i,j),i>j as lower diagonal matrix with 1 on 8! the diagonal. The product of these two gives the original matrix x. 9! Nperm is the order of the permutations of rows, either +1 or -1, 10! unless the matrix is singular, in which case nperm is returned as 0 11! and evaluation of the subroutine is stopped. 12 implicit none 13 integer n,nm,nperm 14 integer iperm(nm) 15 integer i,ii,j,jj,k,kk,im 16 double precision x(nm,nm),temp 17 double precision cmax,scl(nm),check 18! scl(i) is scaling factor for row i, to normalize largest element to 1 19 nperm=1 20 do i=1,n 21 cmax=0.d0 22 do j=1,n 23 if(abs(x(i,j)).gt.cmax) cmax=abs(x(i,j)) 24 enddo 25 if (cmax.eq.0.d0) then 26 nperm=0 27 return 28 endif 29 scl(i)=1.d0/cmax 30 enddo 31 do j=1,n 32 do i=1,j-1 33 temp=x(i,j) 34 do k=1,i-1 35 temp=temp-x(i,k)*x(k,j) 36 enddo 37 x(i,j)=temp 38 enddo 39 cmax=0.d0 40 do i=j,n 41 temp=x(i,j) 42 do k=1,j-1 43 temp=temp-x(i,k)*x(k,j) 44 enddo 45 x(i,j)=temp 46 check=scl(i)*abs(temp) 47 if (check.ge.cmax) then 48 im=i 49 cmax=check 50 endif 51 enddo 52 if (j.ne.im) then 53 do k=1,n 54 temp=x(im,k) 55 x(im,k)=x(j,k) 56 x(j,k)=temp 57 enddo 58 nperm=-nperm 59 check=scl(j) 60 scl(j)=scl(im) 61 scl(im)=check 62 endif 63 iperm(j)=im 64 if(j.ne.n) then 65 if (x(j,j).ne.0) then 66 temp=(1.d0,0.d0)/x(j,j) 67 do i=j+1,n 68 x(i,j)=x(i,j)*temp 69 enddo 70 else 71 nperm=0 72 return 73 endif 74 endif 75 enddo 76 return 77 end 78 79 subroutine baksublu(x,v,n,nm,iperm) 80! Solves a linear equation x'*w=v for w. X is the lower-upper 81! decomposition of x' from subroutine loupd. Vector v is replaced by 82! the solutions w as output. 83 implicit none 84 integer n,nm,iperm(nm),i,j,ii,jj 85 double precision x(nm,nm),v(nm),temp 86 ii=0 87 do i=1,n 88 temp=v(iperm(i)) 89 v(iperm(i))=v(i) 90 if (ii.ne.0) then 91 do j=ii,i-1 92 temp=temp-x(i,j)*v(j) 93 enddo 94 elseif (temp.ne.0.d0) then 95 ii=i 96 endif 97 v(i)=temp 98 enddo 99 do i=n,1,-1 100 temp=v(i) 101 do j=i+1,n 102 temp=temp-x(i,j)*v(j) 103 enddo 104 v(i)=temp/x(i,i) 105 enddo 106 return 107 end 108 109 110 subroutine inverse(x,y,n,nm) 111! Returns matrix y as the inverse of matrix x, both 112! of dimension n x n. 113 implicit none 114 integer n,nm 115 double precision x(nm,nm),y(nm,nm),v(nm) 116 integer i,j,k,ii,jj,kk 117 integer iperm(nm),nperm 118 do i=1,n 119 do j=1,n 120 if (i.eq.j) then 121 y(i,j)=1.d0 122 else 123 y(i,j)=0.d0 124 endif 125 enddo 126 enddo 127 call loupd(x,n,nm,iperm,nperm) 128 do j=1,n 129 call baksublu(x,y(1,j),n,nm,iperm) 130 enddo 131 return 132 end 133 134!***** Complex linear algebra ********************************************* 135 136 subroutine loupdc(x,n,nm,iperm,nperm) 137! Lower-upper decomposition of matrix x, logical size n x n, 138! stored as size nm x nm. Returns matrix x(i,j),i<=j as upper 139! triangular matrix and x(i,j),i>j as lower diagonal matrix with 1 on 140! the diagonal. The product of these two gives the original matrix x. 141! Nperm is the order of the permutations of rows, either +1 or -1, 142! unless the matrix is singular, in which case nperm is returned as 0 143! and evaluation of the subroutine is stopped. 144 implicit none 145 integer n,nm,nperm 146 integer iperm(nm) 147 integer i,ii,j,jj,k,kk,im 148 double complex x(nm,nm),temp 149 double precision cmax,scl(nm),check 150! scl(i) is scaling factor for row i, to normalize largest element to 1 151 nperm=1 152 do i=1,n 153 cmax=0.d0 154 do j=1,n 155 if(abs(x(i,j)).gt.cmax) cmax=abs(x(i,j)) 156 enddo 157 if (cmax.eq.0.d0) then 158 nperm=0 159 return 160 endif 161 scl(i)=1.d0/cmax 162 enddo 163 do j=1,n 164 do i=1,j-1 165 temp=x(i,j) 166 do k=1,i-1 167 temp=temp-x(i,k)*x(k,j) 168 enddo 169 x(i,j)=temp 170 enddo 171 cmax=0.d0 172 do i=j,n 173 temp=x(i,j) 174 do k=1,j-1 175 temp=temp-x(i,k)*x(k,j) 176 enddo 177 x(i,j)=temp 178 check=scl(i)*abs(temp) 179 if (check.ge.cmax) then 180 im=i 181 cmax=check 182 endif 183 enddo 184 if (j.ne.im) then 185 do k=1,n 186 temp=x(im,k) 187 x(im,k)=x(j,k) 188 x(j,k)=temp 189 enddo 190 nperm=-nperm 191 check=scl(j) 192 scl(j)=scl(im) 193 scl(im)=check 194 endif 195 iperm(j)=im 196 if(j.ne.n) then 197 if (x(j,j).ne.0) then 198 temp=(1.d0,0.d0)/x(j,j) 199 do i=j+1,n 200 x(i,j)=x(i,j)*temp 201 enddo 202 else 203 nperm=0 204 return 205 endif 206 endif 207 enddo 208 return 209 end 210 211 subroutine baksubluc(x,v,n,nm,iperm) 212! Solves a linear equation x'*w=v for w. X is the lower-upper 213! decomposition of x' from subroutine loupdc. Vector v is replaced by 214! the solutions w as output. 215 implicit none 216 integer n,nm,iperm(nm),i,j,ii,jj 217 double complex x(nm,nm),v(nm),temp 218 ii=0 219 do i=1,n 220 temp=v(iperm(i)) 221 v(iperm(i))=v(i) 222 if (ii.ne.0) then 223 do j=ii,i-1 224 temp=temp-x(i,j)*v(j) 225 enddo 226 elseif (temp.ne.0.d0) then 227 ii=i 228 endif 229 v(i)=temp 230 enddo 231 do i=n,1,-1 232 temp=v(i) 233 do j=i+1,n 234 temp=temp-x(i,j)*v(j) 235 enddo 236 v(i)=temp/x(i,i) 237 enddo 238 return 239 end 240 241 242 subroutine inversec(x,y,n,nm) 243! Returns matrix y as the inverse of matrix x, both 244! of dimension n x n. 245 implicit none 246 integer n,nm 247 double complex x(nm,nm),y(nm,nm),v(nm) 248 integer i,j,k,ii,jj,kk 249 integer iperm(nm),nperm 250! 251 double complex ll(nm,nm),uu(nm,nm),AA(nm,nm) 252 do i=1,n 253 do j=1,n 254 if (i.eq.j) then 255 y(i,j)=1.d0 256 else 257 y(i,j)=0.d0 258 endif 259 enddo 260 enddo 261 call loupdc(x,n,nm,iperm,nperm) 262 do j=1,n 263 call baksubluc(x,y(1,j),n,nm,iperm) 264 enddo 265 return 266 end 267 268 269