Data Science & Machine Learning Basics
9 min
Friday, September 20, 2024
LU Factorization of a Full rank Matrix using Fortran
Data Science
Machine Learning
Computational Mathematics
Algorithmic Complexity
Programming
Computer Science
! This program factors a full rank matrix A into lower triangular (or trapezoidal) L and upper
! triangular matrix U
program LU_decompostion
implicit none
!############# #################List of main program variable##################################
integer::m,n
! m is the # of rows of the matrix that we are working with
! n is the # of columns that we are working with
doubleprecision, allocatable,dimension(:,:)::A,A1
! A is the working matrix, A1 is the original matrix preserved to check correctness of factoring
integer,allocatable,dimension(:):: P,Q
! P, Q are the row and column permutation vectors for partial and complete pivoting
character::method
!################################################################################################
! ############################ Open an Input Data File###########################################
open(unit=1,file="data.txt")
! ###############################################################################################
! ################################# Read m, n of the matrix A ###################################
write(*,*)"Input the number of rows of the matrix A, m"
read(*,*) m
write(*,*)"Input the number of columns of the matrix A, n"
read(*,*)n
! ##############################################################################################
! ########################### Allocate Space ###################################################
allocate(A(m,n),A1(m,n),P(m),Q(n))
! ##############################################################################################
! Create the matrix A
print*,
call matrixA(m,n,A,A1)
print*,
!##################################### Choose the method #######################################
print*,"What method you want to apply?"
print*,"For No Pivot input: N"
print*, "For Partial Pivot input: P"
print*, "For Complete Pivot input: C"
read(*,*) method
!###############################################################################################
! ############################### Execution of the methods #####################################
IF(method.eq."C".or.method.eq."c") THEN
print*, "Complete Pivoting method has been selected"
print*,
call completePivot(m,n,A,A1,P,Q)
ELSE IF(method.eq."P".or.method.eq."p") then
print*,"Partial Pivoting method has been selected"
print*,
call partialPivot(m,n,A,A1,P)
ELSE IF (method.eq."N".or.method.eq."n") then
print*, "No Pivoting method has been selected"
print*,
call noPivot(m,n,A,A1)
END IF
end program
subroutine matrixA(m,n,A,A1)
integer,intent(in)::m,n
doubleprecision,dimension(m,n),intent(inout)::A,A1
integer::i
print*,
print*, "This is the provided working matrix"
print*
do i=1,m
read(1,*)A(i,:)
A1(i,:)=A(i,:)
print*,A(i,:)
end do
do i=1,n
IF(A(i,i)==0) then
print*,"A 0 entry was found in the main diagonal."
print*, "Therefore, pivoting is a must required"
else if(A(i,i).lt.0.0001) then
print*, "WARNING!! Diagonal Element is too small."
END IF
end do
end subroutine
subroutine completePivot(m,n,A,A1,P,Q)
integer,intent(in)::m,n
doubleprecision,dimension(m,n),intent(inout)::A,A1
integer,dimension(m),intent(out)::P
integer,dimension(n),intent(out)::Q
integer::i,j,k,row,col
doubleprecision::temp
do k=1,n
call max_val(A,m,n,k,row,col)
do i=k,n
temp=A(i,col)
A(i,col)=A(i,k)
A(i,k)=temp
end do
Q(k)=col
do j=k,n
temp=A(k,j)
A(k,j)=A(row,j)
A(row,j)=temp
end do
P(k)=row
A(k+1:n,k)=A(k+1:n,k)/A(k,k)
do j=k+1,n
do i=k+1,n
A(i,j)=A(i,j)-A(i,k)*A(k,j)
end do
end do
end do
print*,
print*,"------------------------------------------------------"
print*," Complete Pivot A=LU factorized array"
print*,"-------------------------------------------------------"
print*,
do i=1,m
print*,A(i,:)
end do
print*,
print*, "Permutation vector P=(",(P(i),i=1,m-1),")"
print*,
print*, "Permutation vector Q=(",(Q(i),i=1,n-1),")"
print*,
print*,"*******************************************************"
print*," Checking Correctness of the factorization"
print*,"*******************************************************"
print*,
call CheckingCompletePivot(m,n,A,A1)
end subroutine
subroutine CheckingCompletePivot(m,n,A,A1)
integer,intent(in)::m,n
doubleprecision,dimension(m,n),intent(inout)::A,A1
doubleprecision,dimension(n,n):: U
doubleprecision,dimension(m,n)::L,A2
integer::i,j
do i=1,m
do j=1,n
if(i.le.j)then
L(i,j)=0
else if (i.gt.j) then
L(i,j)=A(i,j)
end if
L(i,i)=1
end do
end do
do i=1,n
do j=1,n
if(i.le.j)then
U(i,j)=A(i,j)
else
U(i,j)=0.0
end if
end do
end do
print*,
print*, '----------------------------------------------'
print*," Complete Pivot Upper triangular matrix U"
print*, '----------------------------------------------'
print*,
do i=1,n
print*,U(i,:)
end do
print*,
print*, '----------------------------------------------'
print*," Complete Pivot Lower triangular matrix L"
print*, '----------------------------------------------'
print*,
do i=1,m
print*,L(i,:)
end do
A2=matmul(L,U)
print*,
print*, '----------------------------------------------'
print*," Product of L U="
print*, '----------------------------------------------'
print*,
do i=1,m
print*,A2(i,:)
end do
print*,
print*, '----------------------------------------------'
print*, " Factoring Accuracy with the Frobenius-Norm"
print*, '----------------------------------------------'
print*,
print*,"Relative Error=",Frobenius(m,n,A1-A2)/Frobenius(m,n,A1)
print*
print*, "Growth Factor=", Frobenius(m,n,matmul(abs(L),abs(U)))/Frobenius(m,n,A1)
print*,
print*, '----------------------------------------------'
print*, " Factoring Accuracy with the 1-Norm"
print*, '----------------------------------------------'
print*,
print*,"Relative Error=",norm1(m,n,A1-A2)/norm1(m,n,A1)
print*
print*, "Growth Factor=", norm1(m,n,matmul(abs(L),abs(U)))/norm1(m,n,A1)
print*,
print*, '----------------------------------------------'
print*, " Factoring Accuracy with the Infinity-Norm"
print*, '----------------------------------------------'
print*,
print*,"Relative Error=",infinityNorm(m,n,A1-A2)/infinityNorm(m,n,A1)
print*
print*, "Growth Factor=", infinityNorm(m,n,matmul(abs(L),abs(U)))/infinityNorm(m,n,A1)
end subroutine
subroutine partialPivot(m,n,A,A1,P)
integer,intent(in)::m,n
doubleprecision,dimension(m,n),intent(inout)::A,A1
integer,dimension(m),intent(out)::P
integer::i,j,k,row
doubleprecision::temp
do k=1,n
call max_valP(A,m,n,k,row)
do j=k,n
temp=A(k,j)
A(k,j)=A(row,j)
A(row,j)=temp
end do
P(k)=row
A(k+1:n,k)=A(k+1:n,k)/A(k,k)
do j=k+1,n
do i=k+1,n
A(i,j)=A(i,j)-A(i,k)*A(k,j)
end do
end do
end do
print*,
print*,"------------------------------------------------------"
print*," Partial Pivot A=LU factorized array"
print*,"-------------------------------------------------------"
print*,
do i=1,m
print*,A(i,:)
end do
print*,
print*, "Permutation vector P=(",(P(i),i=1,m-1),")"
print*,
print*,"*******************************************************"
print*," Checking Correctness of the factorization"
print*,"*******************************************************"
print*,
call CheckingPartialPivot(m,n,A,A1)
end subroutine
subroutine CheckingPartialPivot(m,n,A,A1)
integer,intent(in)::m,n
doubleprecision,dimension(m,n),intent(inout)::A,A1
doubleprecision,dimension(n,n):: U
doubleprecision,dimension(m,n)::L,A2
integer::i,j
do i=1,m
do j=1,n
if(i.le.j)then
L(i,j)=0
else if (i.gt.j) then
L(i,j)=A(i,j)
end if
L(i,i)=1
end do
end do
do i=1,n
do j=1,n
if(i.le.j)then
U(i,j)=A(i,j)
else
U(i,j)=0.0
end if
end do
end do
print*,
print*, '----------------------------------------------'
print*," Partial Pivot Upper triangular matrix U"
print*, '----------------------------------------------'
print*,
do i=1,n
print*,U(i,:)
end do
print*,
print*, '----------------------------------------------'
print*," Partial Pivot Lower triangular matrix L"
print*, '----------------------------------------------'
print*,
do i=1,m
print*,L(i,:)
end do
A2=matmul(L,U)
print*,
print*, '----------------------------------------------'
print*," Partial Pivot Product of L U="
print*, '----------------------------------------------'
print*,
do i=1,m
print*,A2(i,:)
end do
print*,
print*, '----------------------------------------------'
print*, " Factoring Accuracy with the Frobenius-Norm"
print*, '----------------------------------------------'
print*,
print*,"Relative Error=",Frobenius(m,n,A1-A2)/Frobenius(m,n,A1)
print*
print*, "Growth Factor=", Frobenius(m,n,matmul(abs(L),abs(U)))/Frobenius(m,n,A1)
print*,
print*, '----------------------------------------------'
print*, " Factoring Accuracy with the 1-Norm"
print*, '----------------------------------------------'
print*,
print*,"Relative Error=",norm1(m,n,A1-A2)/norm1(m,n,A1)
print*
print*, "Growth Factor=", norm1(m,n,matmul(abs(L),abs(U)))/norm1(m,n,A1)
print*,
print*, '----------------------------------------------'
print*, " Factoring Accuracy with the Infinity-Norm"
print*, '----------------------------------------------'
print*,
print*,"Relative Error=",infinityNorm(m,n,A1-A2)/infinityNorm(m,n,A1)
print*
print*, "Growth Factor=", infinityNorm(m,n,matmul(abs(L),abs(U)))/infinityNorm(m,n,A1)
end subroutine
subroutine noPivot(m,n,A,A1)
integer,intent(in)::m,n
doubleprecision,dimension(m,n),intent(inout)::A,A1
integer::i,j,k
do k=1,n
A(k+1:n,k)=A(k+1:n,k)/A(k,k)
do j=k+1,n
do i=k+1,n
A(i,j)=A(i,j)-A(i,k)*A(k,j)
end do
end do
end do
print*,
print*,"------------------------------------------------------"
print*," No Pivot A=LU factorized array"
print*,"-------------------------------------------------------"
print*,
do i=1,m
print*,A(i,:)
end do
print*,
print*,"*******************************************************"
print*," Checking Correctness of the factorization"
print*,"*******************************************************"
print*,
call CheckingNoPivot(m,n,A,A1)
end subroutine
subroutine CheckingNoPivot(m,n,A,A1)
integer,intent(in)::m,n
doubleprecision,dimension(m,n),intent(inout)::A,A1
doubleprecision,dimension(n,n):: U
doubleprecision,dimension(m,n)::L,A2
integer::i,j
do i=1,m
do j=1,n
if(i.le.j)then
L(i,j)=0
else if (i.gt.j) then
L(i,j)=A(i,j)
end if
L(i,i)=1
end do
end do
do i=1,n
do j=1,n
if(i.le.j)then
U(i,j)=A(i,j)
else
U(i,j)=0.0
end if
end do
end do
print*,
print*, '----------------------------------------------'
print*," No Pivot Upper triangular matrix U"
print*, '----------------------------------------------'
print*,
do i=1,n
print*,U(i,:)
end do
print*,
print*, '----------------------------------------------'
print*," No Pivot Lower triangular matrix L"
print*, '----------------------------------------------'
print*,
do i=1,m
print*,L(i,:)
end do
A2=matmul(L,U)
print*,
print*, '----------------------------------------------'
print*," No Pivot Product of L U="
print*, '----------------------------------------------'
print*,
do i=1,m
print*,A2(i,:)
end do
print*,
print*, '----------------------------------------------'
print*, " Factoring Accuracy with the Frobenius Norm"
print*, '----------------------------------------------'
print*,
print*,"Relative Error=",Frobenius(m,n,A1-A2)/Frobenius(m,n,A1)
print*
print*, "Growth Factor=", Frobenius(m,n,matmul(abs(L),abs(U)))/Frobenius(m,n,A1)
print*,
print*, '----------------------------------------------'
print*, " Factoring Accuracy with the 1-Norm"
print*, '----------------------------------------------'
print*,
print*,"Relative Error=",norm1(m,n,A1-A2)/norm1(m,n,A1)
print*
print*, "Growth Factor=", norm1(m,n,matmul(abs(L),abs(U)))/norm1(m,n,A1)
print*,
print*, '----------------------------------------------'
print*, " Factoring Accuracy with the Infinity-Norm"
print*, '----------------------------------------------'
print*,
print*,"Relative Error=",infinityNorm(m,n,A1-A2)/infinityNorm(m,n,A1)
print*
print*, "Growth Factor=", infinityNorm(m,n,matmul(abs(L),abs(U)))/infinityNorm(m,n,A1)
end subroutine
subroutine max_val(A,m,n,k,row,col)
implicit none
integer,intent(in)::m,n,k
integer,intent(out)::row,col
doubleprecision,dimension(m,n),intent(inout)::A
doubleprecision::maximum
integer::i,j
maximum=maxval(A(k:m,k:n))
do i=k,m
do j=k,n
if(A(i,j)==maximum) then
row=i
col=j
goto 100
end if
end do
end do
100 end subroutine
subroutine max_valP(A,m,n,k,row)
implicit none
integer,intent(in)::m,n,k
integer,intent(out)::row
doubleprecision,dimension(m,n),intent(inout)::A
doubleprecision::maximum
integer::i,j
maximum=maxval(A(k:m,k))
do i=k,m
if(A(i,k)==maximum) then
row=i
goto 101
end if
end do
101 end subroutine
function norm1(m,n,A)
integer::m,n,i,j
doubleprecision,dimension(m,n)::A
doubleprecision,dimension(n):: normvector
doubleprecision::normval
normval=0
do j=1,n
do i=1,m
normval=normval+abs(A(i,j))
end do
normvector(j)=normval
end do
norm1=maxval(normvector(1:n))
return
end function
function infinityNorm(m,n,A)
integer::m,n,i,j
doubleprecision,dimension(m,n)::A
doubleprecision,dimension(n):: normvector
doubleprecision:: normval
normval=0
do j=1,n
do i=1,m
normval=normval+abs(A(i,j))
normvector(j)=normval
end do
end do
infinityNorm=maxval(normvector(1:m))
return
end function
function Frobenius(m,n,A)
integer::m,n,i,j
doubleprecision,dimension(m,n)::A
doubleprecision:: normval
normval=0
do i=1,m
do j=1,n
normval=normval+(abs(A(i,j)))**2
end do
end do
Frobenius=sqrt(normval)
return
end functionNo Pivot
Input Matrix (Stored in a data file in the same directory where the program file .f90 located)
8 2 9
4 9 4
6 7 9
Comand prompt:
Input the number of rows of the matrix A, m
3
Input the number of columns of the matrix A, n
3
This is the provided working matrix
8.0000000000000000 2.0000000000000000 9.0000000000000000
4.0000000000000000 9.0000000000000000 4.0000000000000000
6.0000000000000000 7.0000000000000000 9.0000000000000000
What method you want to apply?
For No Pivot input: N
For Partial Pivot input: P
For Complete Pivot input: C
n
No Pivoting method has been selected
------------------------------------------------------
No Pivot A=LU factorized array
-------------------------------------------------------
8.0000000000000000 2.0000000000000000 9.0000000000000000
0.50000000000000000 8.0000000000000000 -0.50000000000000000
0.75000000000000000 0.68750000000000000 2.5937500000000000
*******************************************************
Checking Correctness of the factorization
*******************************************************
----------------------------------------------
No Pivot Upper triangular matrix U
----------------------------------------------
8.0000000000000000 2.0000000000000000 9.0000000000000000
0.0000000000000000 8.0000000000000000 -0.50000000000000000
0.0000000000000000 0.0000000000000000 2.5937500000000000
----------------------------------------------
No Pivot Lower triangular matrix L
----------------------------------------------
1.0000000000000000 0.0000000000000000 0.0000000000000000
0.50000000000000000 1.0000000000000000 0.0000000000000000
0.75000000000000000 0.68750000000000000 1.0000000000000000
----------------------------------------------
No Pivot Product of L U=
----------------------------------------------
8.0000000000000000 2.0000000000000000 9.0000000000000000
4.0000000000000000 9.0000000000000000 4.0000000000000000
6.0000000000000000 7.0000000000000000 9.0000000000000000
----------------------------------------------
Factoring Accuracy with the Frobenius Norm
----------------------------------------------
Relative Error= 0.00000000
Growth Factor= 1.02520537
----------------------------------------------
Factoring Accuracy with the 1-Norm
----------------------------------------------
Relative Error= 0
Growth Factor= 1
----------------------------------------------
Factoring Accuracy with the Infinity-Norm
----------------------------------------------
Relative Error= 0
Growth Factor= 1Partial Pivot
Input Matrix (Stored in a data file in the same directory where the program file .f90 located)
1 2 4
2 1 3
3 2 4
Comand Prompt:
Input the number of rows of the matrix A, m
3
Input the number of columns of the matrix A, n
3
This is the provided working matrix
1.0000000000000000 2.0000000000000000 4.0000000000000000
2.0000000000000000 1.0000000000000000 3.0000000000000000
3.0000000000000000 2.0000000000000000 4.0000000000000000
What method you want to apply?
For No Pivot input: N
For Partial Pivot input: P
For Complete Pivot input: C
p
Partial Pivoting method has been selected
------------------------------------------------------
Partial Pivot A=LU factorized array
-------------------------------------------------------
3.0000000000000000 2.0000000000000000 4.0000000000000000
0.66666666666666663 1.3333333333333335 2.6666666666666670
0.33333333333333331 -0.24999999999999992 1.0000000000000000
Permutation vector P=( 3 3 )
*******************************************************
Checking Correctness of the factorization
*******************************************************
----------------------------------------------
Partial Pivot Upper triangular matrix U
----------------------------------------------
3.0000000000000000 2.0000000000000000 4.0000000000000000
0.0000000000000000 1.3333333333333335 2.6666666666666670
0.0000000000000000 0.0000000000000000 1.0000000000000000
----------------------------------------------
Partial Pivot Lower triangular matrix L
----------------------------------------------
1.0000000000000000 0.0000000000000000 0.0000000000000000
0.66666666666666663 1.0000000000000000 0.0000000000000000
0.33333333333333331 -0.24999999999999992 1.0000000000000000
----------------------------------------------
Partial Pivot Product of L U=
----------------------------------------------
3.0000000000000000 2.0000000000000000 4.0000000000000000
2.0000000000000000 2.6666666666666670 5.3333333333333339
1.0000000000000000 0.33333333333333337 1.6666666666666667
----------------------------------------------
Factoring Accuracy with the Frobenius-Norm
----------------------------------------------
Relative Error= 0.618016541
Growth Factor= 1.11492395
----------------------------------------------
Factoring Accuracy with the 1-Norm
----------------------------------------------
Relative Error= 0
Growth Factor= 1
----------------------------------------------
Factoring Accuracy with the Infinity-Norm
----------------------------------------------
Relative Error= 0
Growth Factor= 1Complete Pivot
Input Matrix (Stored in a data file in the same directory where the program file .f90 located)
2 3 4
4 7 5
4 9 5
Command Prompt:
Input the number of rows of the matrix A, m
3
Input the number of columns of the matrix A, n
3
This is the provided working matrix
2.0000000000000000 3.0000000000000000 4.0000000000000000
4.0000000000000000 7.0000000000000000 5.0000000000000000
4.0000000000000000 9.0000000000000000 5.0000000000000000
What method you want to apply?
For No Pivot input: N
For Partial Pivot input: P
For Complete Pivot input: C
c
Complete Pivoting method has been selected
------------------------------------------------------
Complete Pivot A=LU factorized array
-------------------------------------------------------
9.0000000000000000 4.0000000000000000 5.0000000000000000
0.77777777777777779 2.3333333333333335 0.66666666666666674
0.33333333333333331 0.47619047619047616 0.57142857142857140
Permutation vector P=( 3 3 )
Permutation vector Q=( 2 3 )
*******************************************************
Checking Correctness of the factorization
*******************************************************
----------------------------------------------
Complete Pivot Upper triangular matrix U
----------------------------------------------
9.0000000000000000 4.0000000000000000 5.0000000000000000
0.0000000000000000 2.3333333333333335 0.66666666666666674
0.0000000000000000 0.0000000000000000 0.57142857142857140
----------------------------------------------
Complete Pivot Lower triangular matrix L
----------------------------------------------
1.0000000000000000 0.0000000000000000 0.0000000000000000
0.77777777777777779 1.0000000000000000 0.0000000000000000
0.33333333333333331 0.47619047619047616 1.0000000000000000
----------------------------------------------
Product of L U=
----------------------------------------------
9.0000000000000000 4.0000000000000000 5.0000000000000000
7.0000000000000000 5.4444444444444446 4.5555555555555554
3.0000000000000000 2.4444444444444442 2.5555555555555554
----------------------------------------------
Factoring Accuracy with the Frobenius-Norm
----------------------------------------------
Relative Error= 0.683437467
Growth Factor= 1.00393677
----------------------------------------------
Factoring Accuracy with the 1-Norm
----------------------------------------------
Relative Error= 0
Growth Factor= 1
----------------------------------------------
Factoring Accuracy with the Infinity-Norm
----------------------------------------------
Relative Error= 0
Growth Factor= 1Share on
You may also like
No matching items
Citation
BibTeX citation:
@online{islam2021,
author = {Islam, Rafiq},
title = {LU {Factorization} of a {Full} Rank {Matrix} Using {Fortran}},
date = {2021-11-09},
url = {https://mrislambd.github.io/posts/lu/},
langid = {en}
}
For attribution, please cite this work as:
Islam, Rafiq. 2021. “LU Factorization of a Full Rank Matrix Using
Fortran.” November 9, 2021. https://mrislambd.github.io/posts/lu/.