DOUG 0.2

stationary.F90

Go to the documentation of this file.
00001 ! DOUG - Domain decomposition On Unstructured Grids
00002 ! Copyright (C) 1998-2006 Faculty of Computer Science, University of Tartu and
00003 ! Department of Mathematics, University of Bath
00004 !
00005 ! This library is free software; you can redistribute it and/or
00006 ! modify it under the terms of the GNU Lesser General Public
00007 ! License as published by the Free Software Foundation; either
00008 ! version 2.1 of the License, or (at your option) any later version.
00009 !
00010 ! This library is distributed in the hope that it will be useful,
00011 ! but WITHOUT ANY WARRANTY; without even the implied warranty of
00012 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00013 ! Lesser General Public License for more details.
00014 !
00015 ! You should have received a copy of the GNU Lesser General Public
00016 ! License along with this library; if not, write to the Free Software
00017 ! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
00018 ! or contact the authors (University of Tartu, Faculty of Computer Science, Chair
00019 ! of Distributed Systems, Liivi 2, 50409 Tartu, Estonia, http://dougdevel.org,
00020 ! mailto:info(at)dougdevel.org)
00021 
00022 !!----------------------------------
00024 
00025 module stationary_mod
00026   use SpMtx_class
00027   use SpMtx_arrangement
00028 
00029   implicit none
00030 
00031 contains
00033   subroutine SymGaussSeidel(A,x,rhs,iter)
00034     type(SpMtx),intent(inout) :: A
00035     real(kind=rk),intent(in) :: rhs(:) !< right-hand side
00036     real(kind=rk),intent(inout) :: x(:) !< approximation
00037     integer,intent(in) :: iter !< number of iterations
00038 
00039     integer :: it,k,i,j
00040     real(kind=rk),allocatable :: diag(:)
00041 
00042     if (A%arrange_type/=D_SpMtx_ARRNG_ROWS) &
00043          call SpMtx_arrange(A, D_SpMtx_ARRNG_ROWS)
00044 
00045     ! get diagonals
00046     allocate(diag(size(x)))
00047     diag = 0
00048     do k=1,A%nnz
00049       i = A%indi(k)
00050       j = A%indj(k)
00051       if (i==j) then
00052         diag(i) = A%val(k)
00053       end if
00054     end do
00055 
00056     ! do the iterations
00057     do it=1,iter      
00058       ! forward
00059       do i=1,A%nrows
00060         ! trick: ignore row if diag==0
00061         if (diag(i)==0) cycle
00062 
00063         x(i) = rhs(i)
00064         do k=A%m_bound(i),A%m_bound(i+1)-1
00065           j = A%indj(k)
00066           if (i/=j) then
00067             x(i) = x(i) - A%val(k)*x(j)
00068           end if
00069         end do
00070         x(i) = x(i)/diag(i)
00071       end do
00072 
00073       ! backward
00074       do i=A%nrows,1,-1
00075         ! trick: ignore row if diag==0
00076         if (diag(i)==0) cycle
00077 
00078         x(i) = rhs(i)
00079         do k=A%m_bound(i),A%m_bound(i+1)-1
00080           j = A%indj(k)
00081           if (i/=j) then
00082             x(i) = x(i) - A%val(k)*x(j)
00083           end if
00084         end do
00085         x(i) = x(i)/diag(i)
00086       end do
00087     end do
00088   end subroutine SymGaussSeidel
00089 end module stationary_mod