program main !==================================================================== ! A system of linear equations A*x=b ! Method: iterative SOR ! AG: November 2009 !==================================================================== implicit none integer, parameter :: n=3 double precision a(n,n), b(n), x(n), eps, omega integer i,j, iter ! matrix A data (a(1,i), i=1,3) / 8.0, 2.0, 4.0 / data (a(2,i), i=1,3) / 2.0, 6.0, 1.0 / data (a(3,i), i=1,3) / 1.0, 1.0, 8.0 / ! vector b data (b(i), i=1,3) / 2.0, 6.0, 4.0 / ! vector x (initial guess) data (x(i), i=1,3) / 0.0, 0.0, 0.0 / eps = 0.00001 omega =1.2 ! print a header and the original equations write (*,200) do i=1,n write (*,201) (a(i,j),j=1,n), b(i) end do ! print the trial solution write (*,205) write (*,201) (x(i),i=1,n) call gs_sor(a,b,x,omega,eps,n,iter) ! print solutions write (*,203) write (*,201) (x(i),i=1,n) write (*,204) iter 200 format (' The successive-over-relaxation (SOR) ' & ,/,/,' Matrix A and vector b') 201 format (6f12.5) 203 format (/,' Solutions x(n)') 204 format (/,' iterations = ',I5) 205 format (/,' Trial solutions x(n)') end subroutine gs_sor(a,b,x,omega,eps,n,iter) !========================================================== ! Solutions to a system of linear equations A*x=b ! Method: The successive-over-relaxation (SOR) ! Alex G. (November 2009) !---------------------------------------------------------- ! input ... ! a(n,n) - array of coefficients for matrix A ! b(n) - array of the right hand coefficients b ! x(n) - solutions (initial guess) ! n - number of equations (size of matrix A) ! omega - the over-ralaxation factor ! eps - convergence tolerance ! output ... ! x(n) - solutions ! iter - number of iterations to achieve the tolerance ! coments ... ! kmax - max number of allowed iterations !========================================================== implicit none integer, parameter::kmax=100 integer n double precision a(n,n), b(n), x(n) double precision c, omega, eps, delta, conv, sum integer i, j, k, iter, flag ! check if the system is diagonally dominant flag = 0 do i=1,n sum = 0.0 do j=1,n if(i == j) cycle sum = sum+abs(a(i,j)) end do if(abs(a(i,i)) < sum) flag = flag+1 end do if(flag >0) write(*,*) 'The system is NOT diagonally dominant' do k=1,kmax conv = 0.0 do i=1,n delta = b(i) do j=1,n delta = delta - a(i,j)*x(j) end do x(i) = x(i)+omega*delta/a(i,i) if(abs(delta) > conv) conv=abs(delta) end do if(conv < eps) exit end do iter = k if(k == kmax) write (*,*)'The system failed to converge' end subroutine gs_sor