AMS 209: Project: Fortran Code

 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
!!------------------------------------------------------------------
!! Final Project Type A – Numerical Linear Algebra:
!! Gaussian Elimination with and without Pivoting for Solving Linear
!! Systems
!!
!! This code is written by Geetanjali Rakshit for AMS 209.
!!
!! This module has one subroutine which reads parameters (dimension and equation
!! set from the init file 'setup.init'.
!!------------------------------------------------------------------

module setup_module

	contains

	subroutine setup_init(eq_set, n)
		implicit none
		integer, intent(out) :: eq_set, n

		open(unit=10, file="setup.init")
		read(10,*) eq_set
        	read(10,*) n
		close(10)

	end subroutine setup_init

end module setup_module

 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
!!------------------------------------------------------------------
!! Final Project Type A – Numerical Linear Algebra:
!! Gaussian Elimination with and without Pivoting for Solving Linear
!! Systems
!!
!! This code is written by Geetanjali Rakshit for AMS 209.
!!
!! This module has two subroutines which reads in matrix A
!! and vector b from 'A_i.dat' and 'b_i.dat' files.
!!------------------------------------------------------------------


module read_data
	integer :: i
    
	contains

	subroutine read_matrix(eq_set, A, n)
        	implicit none
        	integer, intent(in) :: eq_set, n
        	real (kind=8), dimension(n, n) :: A

		if (eq_set == 1) then
			open(unit=10, file="A_1.dat")
            		do i = 1, n
				read(10, *) A(i, :)
			end do
			close(10)

		elseif (eq_set == 2) then
			open(unit=10, file="A_2.dat")
            		do i = 1, n
				read(10, *) A(i, :)
			end do
			close(10)

		elseif (eq_set == 3) then
			open(unit=10, file="A_3.dat")
            		do i = 1, n
				read(10, *) A(i, :)
			end do
			close(10)

		else
			print *, "Incorrect value entered!"        
		end if

	end subroutine read_matrix

	
	subroutine read_vector(eq_set, b, n)
		implicit none
		integer, intent(in) :: eq_set, n
		real (kind=8), dimension(n) :: b

		if (eq_set == 1) then
			open(unit=10, file="b_1.dat")
			do i = 1, n
				read(10, *) b(i)
			end do
			close(10)

		elseif (eq_set == 2) then
			open(unit=10, file="b_2.dat")
            		do i = 1, n
				read(10, *) b(i)
			end do
			close(10)

		elseif (eq_set == 3) then
			open(unit=10, file="b_3.dat")
			do i = 1, n
				read(10, *) b(i)
			end do
			close(10)

		else
			print *, "Incorrect value entered!"   
		end if

	end subroutine read_vector
            
end module read_data

 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
!!------------------------------------------------------------------
!! Final Project Type A – Numerical Linear Algebra:
!! Gaussian Elimination with and without Pivoting for Solving Linear
!! Systems
!!
!! This code is written by Geetanjali Rakshit for AMS 209.
!!
!! This module has two subroutines: one prints the solution vector x 
!! on the screen and another which writes x to a file 'x_i.dat'.
!!------------------------------------------------------------------


module write_data
	integer :: i
    	
	contains

	subroutine show_x(x, n)
		implicit none
		integer, intent(in) :: n
		real (kind=8), dimension(n) :: x

		print *, "Solution x: "
		do i=1, n
			print *, x(i)
		end do

	end subroutine show_x


	subroutine write_x_to_file(eq_set, x, n)
		implicit none
		integer, intent(in) :: eq_set, n
		real (kind=8), dimension(n) :: x

		if (eq_set == 1) then
			open(unit=10, file="x_1.dat")
			do i = 1, n
				write(10, *) x(i)
			end do
			close(10)

		elseif (eq_set == 2) then
			open(unit=10, file="x_2.dat")
			do i = 1, n
				write(10, *) x(i)
			end do
			close(10)

		elseif (eq_set == 3) then
			open(unit=10, file="x_3.dat")
			do i = 1, n
				write(10, *) x(i)
			end do
			close(10)

		else
			print *, "Incorrect value entered!"    
		end if

	end subroutine write_x_to_file
            
end module write_data

  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
!!------------------------------------------------------------------
!! Fortran code for solving Ax=b using gaussian elimination
!! 
!! This code is written by Geetanjali Rakshit for AMS 209.
!!
!!
!! * Two methods of iteration:
!!   1. With partial pivoting
!!   2. Without partial pivoting
!!
!! * This routine is a driver routine which calls subroutines:
!!   
!!   linear_solve:
!!       |- read_init (from setup_module)
!!       |
!!       |- read_matrix (from read_data)
!!       |- read_vector (from read_data)
!!       |
!!       |- write_to_console (from write_to_screen)
!!       |
!!       |   / GE_without_pivoting (from LU_decomp)
!!       |- |
!!       |   \ GE_with_pivoting (from LU_decomp)
!!       |
!!       |- solve_for_y (from forward_solve)
!!       |
!!       |- solve_for_x (from backward_solve)
!!       |
!!       |- show_x (from write_data)
!!       |- write_x_to_file (from write_data)
!!
!!------------------------------------------------------------------

   
program linear_solve
	!! Begining of the real implementation of the driver
	!! Define usages of module variables and subroutines 

	use setup_module,          only : setup_init
	use read_data,	           only : read_matrix, read_vector
	use write_to_screen,       only : write_to_console
	use LU_decomp,             only : GE_without_pivoting, GE_with_pivoting
	use forward_solve,         only : solve_for_y
	use backward_solve,        only : solve_for_x
	use write_data,            only : show_x, write_x_to_file

	implicit none
	integer :: eq_set, n
	logical :: errorFlag
	real (kind=8), dimension(:,:), allocatable :: A
	real (kind=8), dimension(:,:), allocatable :: L
	real (kind=8), dimension(:,:), allocatable :: U
	real (kind=8), dimension(:), allocatable :: b
	real (kind=8), dimension(:), allocatable :: x
	real (kind=8), dimension(:), allocatable :: y
	character(1) :: pivot_flag

	!! Initialize parameters
  	call setup_init(eq_set, n)

	allocate(A(n,n))
	allocate(L(n,n))
	allocate(U(n,n))
	allocate(b(n))
	allocate(y(n))
	allocate(x(n))

	!! Read values for A and b from .dat files
	call read_matrix(eq_set, A, n)
	call read_vector(eq_set, b, n)

	!! Print A and b to screen
	call write_to_console(A, b, n)

	!! Choose method
	print *, "Use pivoting? (y/n)"
	read *, pivot_flag

	if (pivot_flag == "y") then
		!! Gaussian Elimination without pivoting
		call GE_without_pivoting(A, L, U, n)
	elseif (pivot_flag == "n") then
		!! Gaussian Elimination with pivoting 
		call GE_with_pivoting(A, L, U, b, n)
	else
		print *, "Error!!!"
		call exit(0)
	end if

	!! Forward substitution
	call solve_for_y(L, y, b, n)

	!! Backward substitution
	call solve_for_x(U, x, y, n)

	!! Print solution to screen
	call show_x(x, n)

	!! Write x to file
	call write_x_to_file(eq_set, x, n)

	deallocate(A)
	deallocate(L)
	deallocate(U)
	deallocate(b)
	deallocate(x)
	deallocate(y)

end program linear_solve

  

  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
!!------------------------------------------------------------------
!! Final Project Type A – Numerical Linear Algebra:
!! Gaussian Elimination with and without Pivoting for Solving Linear
!! Systems
!!
!! This code is written by Geetanjali Rakshit for AMS 209.
!!
!! This module has two subroutines that implement Gaussian
!! elimination with and without pivoting.
!!------------------------------------------------------------------


module LU_decomp

	contains
	
	subroutine GE_without_pivoting(A, L, U, n)
		implicit none
		integer, intent(in) :: n
		real (kind=8), dimension(n, n) :: A
		real (kind=8), dimension(n, n) :: L
		real (kind=8), dimension(n, n) :: U
		integer :: i, j, k

		!! Initialize L and U matrices 
		do i=1, n
			do j=1, n
				L(i,j) = 0
				U(i,j) = A(i,j)
				if (i == j) then
					L(i,j) = 1
				end if
			end do
		end do

		do k=1, (n-1)
			if (U(k,k) == 0) exit
			do i=(k+1), n
				L(i,k) = U(i,k)/U(k,k)
			end do

			do j=(k+1), n
				do i=(k+1), n
					U(i,j) = U(i,j) - (L(i,k)*U(k,j))
				end do
			end do
		end do

		do i=2, n
			do j=1, (i-1)
				U(i,j) = 0
			end do
		end do
	end subroutine GE_without_pivoting

	
	subroutine GE_with_pivoting(A, L, U, b, n)
		implicit none
		integer, intent(in) :: n
		real (kind=8) :: pivot, b_temp
		real (kind=8), dimension(n, n) :: A
		real (kind=8), dimension(n,n) :: L
		real (kind=8), dimension(n,n) :: U
		real (kind=8), dimension(n) :: b
		real (kind=8), dimension(n) :: LU_temp
		integer :: i, j, k, ind, p

		!! Initialize L and U matrices 
		do i=1, n
			do j=1, n
				L(i,j) = 0
				U(i,j) = A(i,j)
				if (i == j) then
					L(i,j) = 1
				end if
			end do
		end do

		do k=1, (n-1)
			ind = k
			pivot = abs(U(k,k))
            		!! Find pivot index
			do p=(k+1), n
				if (abs(U(p,k)) > pivot) then
					ind = p
					pivot = abs(U(p,k))
				end if
			end do

		!! Interchange rows
		if (ind /= k) then
			b_temp = b(ind)
			b(ind) = b(k)
			b(k) = b_temp
			do i=k, n
				LU_temp(i) = U(ind,i)
				U(ind,i) = U(k,i)
				U(k,i) = LU_temp(i)
			end do

			do i=1, (k-1)
				LU_temp(i) = L(ind,i)
				L(ind,i) = L(k,i)
				L(k,i) = LU_temp(i)
			end do
		end if

		if (U(k,k) /= 0) then
			do j=(k+1), n
				L(j,k) = U(j,k)/U(k,k)
				do i=k, n
					U(j,i) = U(j,i)-(L(j,k)*U(k,i))
				end do
			end do
		end if
	end do
	
	end subroutine GE_with_pivoting

end module LU_decomp

 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
!!------------------------------------------------------------------
!! Final Project Type A – Numerical Linear Algebra:
!! Gaussian Elimination with and without Pivoting for Solving Linear
!! Systems
!!
!! This code is written by Geetanjali Rakshit for AMS 209.
!!
!! This module has one subroutine that implements backward
!! substitution and solves Ux = y.
!!------------------------------------------------------------------


module backward_solve

	contains

	subroutine solve_for_x(U, x, y, n)
		implicit none
		integer, intent(in) :: n
		real (kind=8), dimension(n) :: x
		real (kind=8), dimension(n) :: y
		real (kind=8), dimension(n,n) :: U
		real (kind=8) :: temp
		integer :: j, k

		x(n) = y(n)/U(n,n)

		do k=(n-1), 1, -1
			temp = 0
			do j=(k+1), n
				temp = temp + (x(j)*U(k,j))
			end do
			x(k) = (y(k)-temp)/U(k,k)
		end do
            
	end subroutine solve_for_x

end module backward_solve

 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
!!------------------------------------------------------------------
!! Final Project Type A – Numerical Linear Algebra:
!! Gaussian Elimination with and without Pivoting for Solving Linear
!! Systems
!!
!! This code is written by Geetanjali Rakshit for AMS 209.
!!
!! This module has one subroutine that implements forward
!! substitution and solves Ly = b.
!!------------------------------------------------------------------


module forward_solve

	contains

	subroutine solve_for_y(L, y, b, n)

		implicit none
		integer, intent(in) :: n
		real (kind=8), dimension(n) :: y
		real (kind=8), dimension(n,n) :: L
		real (kind=8), dimension(n) :: b
		real (kind=8) :: temp
		integer :: j, k

        	y(1) = b(1)

		do k=2, n
			temp = 0
			do j=1, (k-1)
				temp = temp + (y(j)*L(k,j))
			end do
		y(k) = b(k) - temp
		end do

	end subroutine solve_for_y

end module forward_solve

 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
!!------------------------------------------------------------------
!! Final Project Type A – Numerical Linear Algebra:
!! Gaussian Elimination with and without Pivoting for Solving Linear
!! Systems
!!
!! This code is written by Geetanjali Rakshit for AMS 209.
!!
!! This module has one subroutine which writes matrix A
!! and vector b to screen.
!!------------------------------------------------------------------


module write_to_screen
	contains

	subroutine write_to_console(A, b, n)
		implicit none
		integer, intent(in) :: n
		real (kind=8), dimension(n, n) :: A
		real (kind=8), dimension(n) :: b
		integer :: i

		print *, "A: "
		do i=1, n
			print *, A(i,:)
		end do

		print *, "b: "
		do i=1, n
			print *, b(i)
		end do

	end subroutine write_to_console

end module write_to_screen