gauss3.f90

I suggest you either typed  this as it appears into the compiler, or as an easier way, just copy and paste it!

Note: I do not claim the program will 100% work, but I give it my best. If it doesn't, try tweaking it yourself, or email me about it, and I'll modify and post the improved version later. Good luck!

module linear_equations
implicit none
private
public :: gaussian_solve

contains
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subroutine gaussian_solve(a,b,error)
!solves Ax=b
!A is stored in a
!solution is put in b
!error indicates if error is found
!dummy arguments
real, dimension(:,:), intent(inout) :: a
real, dimension(:), intent(inout) ::b
integer, intent(out) :: error
!reduce the equations by gaussian elimination
call gaussian_elimination(a,b,error)

!if reduction is succesful, calculate solution by back subtitution
if (error==0) then 
call back_substitution(a,b,error)
end if
return
end subroutine gaussian_solve

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subroutine gaussian_elimination(a,b,error)
!performs gaussian eliminatiuon on a system of linear equations
!dummy arguments
!a contains the coefficients
!b contains the right hand side
real, dimension(:,:), intent(inout)::a
real, dimension(:), intent(inout)::b
integer, intent(out) ::error

!local variable
real, dimension(size(a,1))::temp_array !automatic array
integer, dimension(1) :: ksave
integer :: i,j,k,n
real ::temp, m

!validity check
n=size(a,1)
if (n==0) then
error=-1 !there is no prob to solve
return
end if
if (n/=size(a,2)) then
error=-2 !a is not square
return
end if
if (n/=size(b)) then
error=-3 !size of b does not match a
return
end if

!dimension of array r ok !&&&&&&&&&&&&& got prob!!! &&&&&&&&&&&&
error=0
do i=1, n-1
!find row with largest value of |a(j,i)|, j=i,...,n
ksave=maxloc(abs(a(i:n,i)))
!check whether largest |a(j,i)| is zero
k=ksave(1) +(i-1)
if (abs(a(k,i))<=1e-5) then
error =-4 !no solution poissible !************************************ error 4 **********
return
end if

!interchange row i and row k if necessary
if (k/=i) then
a(i,:) =a(k,:)
a(k,:)=temp_array
!interchanges corresponding elements of b
temp=b(i)
b(i)=b(k)
b(k)=temp
end if

!subtract multiples of row i from subsequent rows to
!zero all subsequent coefficients of x sub i
do j=i+1,n
m=a(j,i)/a(i,i)
a(j,:)=a(j,:)-m*a(i,:)
b(j)=b(j)-m*b(i)
end do
end do
return
end subroutine gaussian_elimination

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subroutine back_substitution(a,b,error)
!this conducts back substitution when equation has been reduced by gaussian wlimination
! dummy arguments
!array a contain coefficients
!b contains ride hand side coefficients and is solution
!error is non zero if an error occurs
real, dimension(:,:), intent(inout) ::a
real, dimension(:), intent(inout) ::b
integer, intent(out)::error

!local variables
real::sum
integer::i,j,n

error=0
n=size(b)

!solve for each variable in turn

do i=n,1,-1
!check for zero coefficients
!if (abs(a(i,i))<=1e-5) then
!error=-4 !************************************ error 4 **********
!return
!end if

sum=b(i)
do j=i+1,n
sum=sum-a(i,j)*b(j)
end do
b(i)=sum/a(i,i)
end do
end subroutine back_substitution

end module linear_equations

!***********************************************************************


program gaussian_elimination_proper
use linear_equations
implicit none
!defines coefficient and solve them


!allocatable array for coefficients
real, allocatable, dimension(:,:) ::a
real, allocatable, dimension(:) ::b
!size of array
integer::n
!loop variables and error flag
integer::i,j,error

!get size of prob
print*, "how many equations are there?"
read*, n

!allocate arrays
allocate(a(n,n),b(n))

!get corfficients
print*, "type coefficients for each equation in turn"
do i=1,n
read*, (a(i,j), j=1,n)
print*, "again"
read*, b(i)
end do

!attempt to solve system of equations
call gaussian_solve(a,b,error)

!check to see if there were any errors
if (error<=-1 .and. error >=-3) then
print*, "error in call to gaussian_solve"
else if (error==-4) then
print*, "system is degenerate"
else
print*, ""
print*, "solution is"
print'(1X, "x("i2,")=",f6.2)', (i,b(i),i=1,n)
end if

end program gaussian_elimination_proper