c c c Sample_Code_1 c Solves A System of Linear Equations by LU Decomposition c A . X = B c LU Decomposition Through Subroutine ludcmp c Forward and Backward Substitution Through Subroutine lubksb c c IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (N = 4) ! N = Number of Equations c DOUBLE PRECISION A(N,N),B(N) ! Matrix A and RHS Vector B DOUBLE PRECISION dummy ! condition number of the matrix integer ipvt(N) ! Integer vector containing number of bivit columns OPEN(10,FILE='sample_code_1_input.txt') OPEN(11,FILE='sample_code_1_output.txt') c c initialize matrix A and vector B c do i = 1, N B(i) = 0.d0 do j = 1, N A(i,j)=0.d0 enddo enddo c c read input matrix A c do i = 1, N do j = 1, N read(10,*) A(i,j) enddo enddo c c read RHS vector B c do i = 1, N read(10,*) B(i) enddo c c LU-decompose matrix A c Note that matrix A will contains on return the L and U matrixes (original is destroyed) c call ludcmp(A,N,N,ipvt,dummy) c c Solve the equations by backsubstitution c Note that vector B will contains the unknowns X (original is destroyed) c call lubksb(A,N,N,ipvt,B) c....................................................................... do i= 1, N write(* ,'(e24.14)') B(i) write(11,'(e24.14)') B(i) enddo end c SUBROUTINE ludcmp(a,n,np,indx,d) IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER n,np,indx(n),NMAX DOUBLE PRECISION d,a(np,np),TINY PARAMETER (NMAX=500,TINY=1.0e-20) INTEGER i,imax,j,k DOUBLE PRECISION aamax,dum,sum,vv(NMAX) d=1. do 12 i=1,n aamax=0.d0 do 11 j=1,n if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j)) 11 continue if (aamax.eq.0.d0) pause 'singular matrix in ludcmp' vv(i)=1./aamax 12 continue do 19 j=1,n do 14 i=1,j-1 sum=a(i,j) do 13 k=1,i-1 sum=sum-a(i,k)*a(k,j) 13 continue a(i,j)=sum 14 continue aamax=0.d0 do 16 i=j,n sum=a(i,j) do 15 k=1,j-1 sum=sum-a(i,k)*a(k,j) 15 continue a(i,j)=sum dum=vv(i)*abs(sum) if (dum.ge.aamax) then imax=i aamax=dum endif 16 continue if (j.ne.imax)then do 17 k=1,n dum=a(imax,k) a(imax,k)=a(j,k) a(j,k)=dum 17 continue d=-d vv(imax)=vv(j) endif indx(j)=imax if(a(j,j).eq.0.d0)a(j,j)=TINY if(j.ne.n)then dum=1./a(j,j) do 18 i=j+1,n a(i,j)=a(i,j)*dum 18 continue endif 19 continue return END c SUBROUTINE lubksb(a,n,np,indx,b) IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER n,np,indx(n) DOUBLE PRECISION a(np,np),b(n) INTEGER i,ii,j,ll DOUBLE PRECISION sum ii=0 do 12 i=1,n ll=indx(i) sum=b(ll) b(ll)=b(i) if (ii.ne.0)then do 11 j=ii,i-1 sum=sum-a(i,j)*b(j) 11 continue else if (sum.ne.0.d0) then ii=i endif b(i)=sum 12 continue do 14 i=n,1,-1 sum=b(i) do 13 j=i+1,n sum=sum-a(i,j)*b(j) 13 continue b(i)=sum/a(i,i) 14 continue return END