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<threshold:
        print "The answers computed from Fortran matches with Python Package"
    else:
        print "The error is too large and is equal to"
        print error         
    fig=plt.figure()
    ax=fig.add_subplot(1,1,1)
    ax.set_aspect('equal')
    plt.imshow(mat_py,interpolation='nearest',cmap=plt.cm.ocean)
    plt.colorbar()
    pylab.savefig('A_'+repr(index)+'.png')
    plt.show()
 
    fig=plt.figure()
    ax=fig.add_subplot(1,1,1)
    ax.set_aspect('equal')
    xx=np.array([kk,]*n).transpose()
    plt.imshow(xx,interpolation='nearest',cmap=plt.cm.ocean)
    plt.colorbar()
    pylab.savefig('x_'+repr(index)+'method_'+repr(mtype)+'.png')
    plt.show() 

    fig=plt.figure()
    ax=fig.add_subplot(1,1,1)
    ax.set_aspect('equal')
    bb=np.array([b_py,]*n).transpose()
    plt.imshow(bb,interpolation='nearest',cmap=plt.cm.ocean)
    plt.colorbar()
    pylab.savefig('b_'+repr(index)+'.png')
    plt.show()    
    

if __name__ == '__main__':
    for i in range(0,3):
        make_make()
        make_matrix(i)
        run_LU()
        check()
        
       
    
    
    
    
    
    
    
    
    
    

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

Fortran Module 3: read_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
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