1 subroutine ludcmp(aa,n,np,indx,d) 2c--------------------------------------------------------------------- 3c 4c find the LU decomposition of a matrix. 5c 6c Numerical Recipes: p. 38 7c 8c--------------------------------------------------------------------- 9 include "common.h" 10 11 dimension aa(np,np), indx(n), vv(500) 12 13 d=1.0 14 do i=1,n 15 aamax=0.0 16 do j=1,n 17 if (abs(aa(i,j)).gt.aamax) aamax=abs(aa(i,j)) 18 enddo 19 vv(i)=1.0/aamax 20 enddo 21 22 do j=1,n 23 do i=1,j-1 24 sum=aa(i,j) 25 do k=1,i-1 26 sum=sum-aa(i,k)*aa(k,j) 27 enddo 28 aa(i,j)=sum 29 enddo 30 aamax = 0.0 31 do i=j,n 32 sum=aa(i,j) 33 do k=1,j-1 34 sum=sum-aa(i,k)*aa(k,j) 35 enddo 36 aa(i,j)=sum 37 dum=vv(i)*abs(sum) 38 if (dum.ge.aamax) then 39 imax=i 40 aamax=dum 41 endif 42 enddo 43 if (j.ne.imax) then 44 do k=1,n 45 dum=aa(imax,k) 46 aa(imax,k)=aa(j,k) 47 aa(j,k)=dum 48 enddo 49 d=-d 50 vv(imax)=vv(j) 51 endif 52 indx(j)=imax 53 if(j.ne.n)then 54 dum=1.0/aa(j,j) 55 do i=j+1,n 56 aa(i,j)=aa(i,j)*dum 57 enddo 58 endif 59 enddo 60 61 return 62 end 63 64 65 66 67 68