sor_sprs.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Subroutine :  s o r _ s p r s
!
!> @file
!!
!! SOR solver for a system of linear equations lgs_a*lgs_x=lgs_b
!! [matrix in compressed row storage CRS,
!! represented by arrays lgs_a_value(values), lgs_a_index (indices)
!! and lgs_a_ptr (pointers)].
!!
!! @section Copyright
!!
!! Copyright 2009-2013 Ralf Greve, Tatsuru Sato
!!
!! @section License
!!
!! This file is part of SICOPOLIS.
!!
!! SICOPOLIS is free software: you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! SICOPOLIS is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with SICOPOLIS.  If not, see <http://www.gnu.org/licenses/>.
!<
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!-------------------------------------------------------------------------------
!> SOR solver for a system of linear equations lgs_a*lgs_x=lgs_b
!! [matrix in compressed row storage CRS,
!! represented by arrays lgs_a_value(values), lgs_a_index (indices)
!! and lgs_a_ptr (pointers)].
!<------------------------------------------------------------------------------
subroutine sor_sprs(lgs_a_value, lgs_a_index, lgs_a_diag_index, lgs_a_ptr, &
                    lgs_b_value, &
                    nnz, nmax, omega, eps_sor, lgs_x_value, ierr)

use sico_types

implicit none

integer(i4b),                     intent(in) :: nnz, nmax
real(dp),                         intent(in) :: omega, eps_sor
integer(i4b), dimension(nmax+1),  intent(in) :: lgs_a_ptr
integer(i4b), dimension(nnz),     intent(in) :: lgs_a_index
integer(i4b), dimension(nmax),    intent(in) :: lgs_a_diag_index
real(dp),     dimension(nnz),     intent(in) :: lgs_a_value
real(dp),     dimension(nmax),    intent(in) :: lgs_b_value

integer(i4b),                    intent(out) :: ierr
real(dp),     dimension(nmax), intent(inout) :: lgs_x_value

integer(i4b) :: iter
integer(i4b) :: nr, k
real(dp), allocatable, dimension(:) :: lgs_x_value_prev
real(dp) :: b_nr
logical :: flag_convergence

allocate(lgs_x_value_prev(nmax))

iter_loop : do iter=1, 1000

   lgs_x_value_prev = lgs_x_value

   do nr=1, nmax

      b_nr = 0.0_dp 

      do k=lgs_a_ptr(nr), lgs_a_ptr(nr+1)-1
         b_nr = b_nr + lgs_a_value(k)*lgs_x_value(lgs_a_index(k))
      end do

      lgs_x_value(nr) = lgs_x_value(nr) &
                        -omega*(b_nr-lgs_b_value(nr)) &
                              /lgs_a_value(lgs_a_diag_index(nr))

   end do

   flag_convergence = .true.
   do nr=1, nmax
      if (abs(lgs_x_value(nr)-lgs_x_value_prev(nr)) > eps_sor) then
         flag_convergence = .false.
         exit
      end if
   end do

   if (flag_convergence) then
      write(6,'(11x,a,i5)') 'sor_sprs: iter = ', iter
      ierr = 0   ! convergence criterion fulfilled
      deallocate(lgs_x_value_prev)
      return
   end if

end do iter_loop

write(6,'(a,i5)') 'sor_sprs: iter = ', iter
ierr = -1   ! convergence criterion not fulfilled
deallocate(lgs_x_value_prev)

end subroutine sor_sprs
!