Project¶

Solution:

The documentation of the final project can be downloaded here `here`.

1. Python Implementation of Fortran LU Solver¶

Solution:

 ``` 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 119 120 121 122 123 124 125``` ```import os import glob import numpy as np import matplotlib.pyplot as plt import pylab import os os.chdir('../LinAlg') global threshold threshold=1e-14 def make_make(): if os.path.exists('LU.exe')==True: cmd1 = 'make clean' os.system(cmd1) cmd2='make' os.system(cmd2) def make_matrix(i): global b_py global mat_py global index global n index=i mat_master=([1,1,-2,1,2,1,-1,-2,1],[4,3,2,1,3,4,3,2,2,3,4,3,1,2,3,4],[1,-1,2,-3,-1,3,-4,7,1,-3,7,-10,-1,3,-7,14]) bmat_master=([1,0,1],[1,1,-1,-1],[0,2,-2,-8]) mat=mat_master[index] bmat=bmat_master[index] b_py=np.array(bmat) dim=len(mat) n=int(np.sqrt(dim)) x2=np.array(mat) mat_py=np.transpose(np.reshape(x2,(n,n))) if os.path.exists('matrix.init')==True: number=len(glob.glob('matrix.init*')) cmd1='cp matrix.init'+' matrix.init.'+repr(number) os.system(cmd1) cmd3='touch matrix.init' os.system(cmd3) f=open('matrix.init','w') for i in range(0,dim): f.write(repr(mat[i])+"\n") f.close() if os.path.exists('b.init')==True: number=len(glob.glob('b.init*')) cmd5='cp b.init'+' b.init.'+repr(number) os.system(cmd5) cmd4='touch b.init' os.system(cmd4) f=open('b.init','w') for i in range(0,n): f.write(repr(bmat[i])+"\n") f.close() global mtype print "please type 1 for LU without Pivoting and 2 for LU with Pivoting" mtype=int(raw_input('> ')) f=open('method.init','w') f.write(repr(mtype)) f.close() if os.path.exists('x.dat')==True: number=len(glob.glob('*.dat')) cmd='cp x.dat'+' x_'+repr(number)+'.dat' os.system(cmd) def run_LU(): kobe='./LU.exe' os.system(kobe) def check(): os.chdir('../LinAlg') check1=np.dot(np.linalg.inv(mat_py),b_py) kk=np.loadtxt("x.dat") error=abs(sum(check1-kk)) if error

2. Fortran Module 1: linear_solve.F90¶

Solution:

 ``` 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``` ```program LU use readin use sanity use LU_decomp use forward_solve use backward_solve use write_data implicit none real, allocatable, dimension(:) :: temp1,temp2,tt,b,x,d,b2,d2,x2,x3 real,allocatable,dimension(:,:) ::A,U,L,P,U2,L2 real,dimension(100)::A2 real::threshold,s integer::i,j,K,m,ind,n call readmatrix(A,U,L,P,b,x,d,n,ind) A=A U=U L=L P=P b=b x=x d=n ind=ind call sanity_print(A,b,n) call LU_compute(n,ind,A,P,U,L,b,d,L2,U2,b2,d2) L=L2 U=U2 b=b2 d=d2 call forward(n,ind,b,d,L,x2) x=x2 call backward(n,x,U,x3) x=x3 call printx(n,x) end program LU ```

Fortran Module 2: LU_decomp¶

Solution:

 ``` 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 119``` ```module LU_decomp implicit none contains subroutine LU_compute(n,ind,A,P,U,L,b,d,L2,U2,b2,d2) integer,intent(in)::ind,n real,allocatable,dimension(:,:)::A,U,L,P real,allocatable,dimension(:)::b,d real,allocatable,dimension(:,:),intent(out) ::U2,L2 real,allocatable,dimension(:),intent(out) ::b2,d2 real,allocatable,dimension(:) :: temp1,temp2,tt integer::i,j,K,m !initialize L matrix do K=1,n P(K,K)=1 enddo L=P if(ind==2) then !initialize L matrix do K=1,n P(K,K)=1 enddo L=P !LU Step do K=1,n-1 tt=maxloc(abs(U(K:,K))) m=tt(1)+(K-1) if(m>K) then temp1=U(K,K:n) temp2=U(m,K:n) U(K,K:n)=temp2 U(m,K:n)=temp1 temp1=L(K,1:(K-1)) temp2=L(m,1:(K-1)) L(K,1:(K-1))=temp2 L(m,1:(K-1))=temp1 temp1=P(K,:) temp2=P(m,:) P(K,:)=temp2 P(m,:)=temp1 endif if (A(K,K)==0) then exit endif do i=K+1,n L(i,K)=U(i,K)/U(K,K) enddo do j=K+1,n do i=K+1,n U(i,j)=U(i,j)-L(i,K)*U(K,j) enddo enddo enddo !Setting lower half U to zero since it's upper Triangular Matrix do j=1,n-1 do i=j+1,n U(i,j)=0 enddo enddo d=matmul(P,b) endif if(ind==1) then !LU Step do K=1,n-1 if (A(K,K)==0) then exit endif do i=K+1,n L(i,K)=U(i,K)/U(K,K) enddo do j=K+1,n do i=K+1,n U(i,j)=U(i,j)-L(i,K)*U(K,j) enddo enddo enddo !Setting lower half U to zero since it's upper Triangular Matrix do j=1,n-1 do i=j+1,n U(i,j)=0 enddo enddo endif allocate(L2(n,n)) allocate(U2(n,n)) allocate(b2(n)) allocate(d2(n)) L2=L U2=U b2=b d2=d end subroutine LU_compute end module LU_decomp ```

Solution:

 ``` 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``` ```module readin implicit none contains subroutine readmatrix(A,U,L,P,b,x,d,n,ind) implicit none real,dimension(100)::A2 real,allocatable,dimension(:,:),intent(out) ::A,U,L,P real, allocatable, dimension(:),intent(out) :: b,x,d integer,intent(out)::ind,n integer::i,j,K,m,stat open(15,file="method.init") read(15,*) ind open(11,file="matrix.init") do j=1,20 read(11,*,IOSTAT=stat) A2(j) IF(IS_IOSTAT_END(stat)) then exit endif enddo n=sqrt(real(j-1)) allocate(A(n,n)) allocate(L(n,n)) allocate(U(n,n)) allocate(P(n,n)) A=reshape(A2(1:(j-1)),(/n,n/)) U=A allocate(b(n)) allocate(x(n)) allocate(d(n)) open(12,file="b.init") do i=1,n read(12,*) b(i) enddo end subroutine readmatrix end module readin ```

Fortran Module 4: write_data.F90¶

Solution:

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23``` ```module write_data implicit none contains subroutine printx(n,x) integer,intent(in)::n real,allocatable,dimension(:),intent(in) ::x integer::i print *, "Answer X is" do i=1,n write (*,*) (x(i)) enddo open(21,file='x.dat') do i=1,n write(21,*) x(i) enddo end subroutine printx end module ```

Fortran Module 5: forward_solve.F90¶

Solution:

 ``` 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``` ```module forward_solve implicit none contains subroutine forward(n,ind,b,d,L,x2) integer,intent(in)::n,ind real,allocatable,dimension(:)::b,d real,allocatable,dimension(:,:),intent(in)::L real,allocatable,dimension(:),intent(out) ::x2 real::s integer::i,j allocate(x2(n)) if(ind==2) then x2(1)=d(1) do i=2,n s=0.0 do j=1,i-1 s=s+L(i,j)*x2(j) enddo x2(i)=d(i)-s enddo endif if(ind==1) then x2(1)=b(1) do i=2,n s=0 do j=1,i-1 s=s+L(i,j)*x2(j) enddo x2(i)=b(i)-s enddo endif end subroutine forward end module forward_solve ```

Fortran Module 6: backward_solve.F90¶

Solution:

 ``` 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``` ```module backward_solve implicit none contains subroutine backward(n,x,U,x3) integer,intent(in)::n real,allocatable,dimension(:,:),intent(in)::U real,allocatable,dimension(:),intent(in) ::x real,allocatable,dimension(:),intent(out) ::x3 real::s integer::i,j allocate(x3(n)) x3=x x3(n)=x3(n)/U(n,n) do i=n-1,1,-1 s=0.0 do j=i+1, n s=s+U(i,j)*x3(j) enddo x3(i)=(x3(i)-s)/U(i,i) enddo end subroutine backward end module backward_solve ```

Fortran Module 7: write_to_screen.F90¶

Solution:

 ``` 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``` ```module sanity implicit none contains subroutine sanity_print(A,b,n) real,allocatable,dimension(:,:),intent(in) ::A real, allocatable, dimension(:),intent(in) :: b integer,intent(in)::n integer::i,j print *, "Matrix A is" do i=1,n write (*,*) (A(i,j),j=1,n) enddo print *, "b is" do i=1,n write (*,*) (b(i)) enddo end subroutine sanity_print end module ```

Fortran Module 8: makefile¶

Solution:

 ``` 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``` ```FC= gfortran FFLAGS_DEBUG = -Wall -Wextra -Wimplicit-interface -fPIC -fmax-errors=1 -g -fcheck=all -fbacktrace FFLAGS_OPT = -ggdb -O3 -fdefault-real-8 -fdefault-double-8 -ffree-line-length-none -Wuninitialized EXE_FILE = LU.exe OBJS=linear_solve.o \ read_data.o\ write_to_screen.o\ LU_decomp.o\ forward_solve.o\ backward_solve.o\ write_data.o\ \$(EXE_FILE) : \$(OBJS) @\$(FC) \$(FFLAGS_OPT) \$(OBJS) -o \$(EXE_FILE) @echo "code is now linking..." %.o: %.F90 \$(FC) \$(FFLAGS_OPT) -c \$< clean: @rm -f *.o *.mod *~ LU.exe # debug: clean debug: FFLAGS_OPT = \$(FFLAGS_DEBUG) debug: \$(EXE_FILE) linear_solve.o: read_data.o write_to_screen.o LU_decomp.o forward_solve.o backward_solve.o write_data.o ```