DOUG 0.2

RobustCoarseMtx.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 
00024 
00027 module RobustCoarseMtx_mod
00028   use SpMtx_class
00029   use SpMtx_util
00030   use SpMtx_op_AB
00031   use subsolvers
00032   use SpMtx_op_Ax
00033 
00034   implicit none
00035 
00042   type SumOfInversedSubMtx
00043      type(SpMtx), pointer :: A !< Original matrix
00044      type(SpMtx), dimension(:), pointer :: Ai !< Submatrices \f$A_i\f$
00045      type(SpMtx), dimension(:), pointer :: R !< Restriction matrices for the submatrices
00046 
00047      integer, dimension(:), pointer :: subsolve_ids !< factorisation IDs (UMFPACK) for \e Ai submatrices
00048   end type SumOfInversedSubMtx
00049 
00050 contains
00053   function CoarseProjectionMtxsBuild(A,R,naggr) result (B)
00054     type(SpMtx), intent(inout) :: A
00055     type(SpMtx), intent(in) :: R !< restrict matrix from usual aggregate with smoothing
00056     integer,intent(in) :: naggr !< number of aggregates
00057     type(SumOfInversedSubMtx) :: B
00058 
00059     B = getConstrainedEnergyMatrix(A, R)
00060         
00061     ! ----- routines ---------
00062     contains
00063       function getConstrainedEnergyMatrix(A, R) result (B)
00064         type(SpMtx), intent(inout), target :: A
00065         type(SpMtx), intent(in) :: R
00066         type(SumOfInversedSubMtx) :: B
00067 
00069         type Array
00070            integer, pointer :: values(:)
00071         end type Array
00072 
00074         integer, dimension(:), pointer :: n_aggr_nodes
00075 
00076         type(Array), dimension(:), pointer :: aggr_nodes
00077 
00078         integer :: i, iNode, j
00079         integer :: coarse_node, fine_node, last_index
00080         integer, allocatable :: last_nodes(:) !< helper array: number of added nodes
00081         type(SpMtx) :: AiTemp
00082 
00083         ! For now hack: aggregate node numbers (incl. overlap nodes) is extracted from 
00084         ! the standard case restriction matrix.
00085 
00086         if(sctls%verbose>4) then
00087            print *, "Restrict : nrows, ncols=", R%nrows, R%ncols
00088            call SpMtx_printRaw(R)
00089         end if
00090 
00091         allocate(n_aggr_nodes(naggr))
00092         n_aggr_nodes = 0
00093         do i=1, R%nnz
00094            coarse_node = R%indi(i)
00095            n_aggr_nodes(coarse_node) = n_aggr_nodes(coarse_node) + 1
00096         end do
00097         allocate(aggr_nodes(naggr))
00098         do iNode=1,naggr
00099            allocate(aggr_nodes(iNode)%values(n_aggr_nodes(iNode)))
00100         end do
00101         allocate(last_nodes(naggr))
00102         last_nodes = 0
00103         do i=1, R%nnz
00104            coarse_node = R%indi(i)
00105            fine_node = R%indj(i)
00106            last_nodes(coarse_node) = last_nodes(coarse_node)+1
00107            last_index = last_nodes(coarse_node)
00108            aggr_nodes(coarse_node)%values(last_index) = fine_node
00109         end do
00110         
00111         ! Now set the result
00112         B%A => A
00113         allocate(B%R(naggr))
00114         do i=1,naggr
00115           B%R(i) = SpMtx_NewInit(n_aggr_nodes(i), nrows=n_aggr_nodes(i), ncols=A%nrows)
00116           forall(j=1:n_aggr_nodes(i)) B%R(i)%indi(j) = j
00117           B%R(i)%indj = aggr_nodes(i)%values
00118           B%R(i)%val = 1.0
00119           !B%R(i)%arrange_type = D_SpMtx_ARRNG_ROWS
00120           if(sctls%verbose>4) then
00121              print *, "R for", i, ": nrows, ncols=", B%R(i)%nrows, B%R(i)%ncols
00122              call SpMtx_printRaw(B%R(i))
00123           end if
00124         end do
00125 
00126         ! calculate submatrices
00127         allocate(B%Ai(naggr))
00128         do i=1,naggr
00129            AiTemp = SpMtx_AB2(B%R(i), B%A)
00130            B%Ai(i) = SpMtx_AB2(AiTemp, B%R(i), BT=.TRUE.) 
00131         end do
00132         
00133         ! finally allocate ID array
00134         allocate(B%subsolve_ids(naggr))
00135         B%subsolve_ids = 0
00136       end function getConstrainedEnergyMatrix
00137 
00138   end function CoarseProjectionMtxsBuild
00139 
00141   subroutine SOISMtx_pmvm(y,A,x)
00142     type(SumOfInversedSubMtx), intent(inout) :: A !< System matrix
00143     real(kind=rk),dimension(:), pointer :: x !< Vector
00144     !type(Mesh),   intent(in) :: M ! Mesh
00145     real(kind=rk),dimension(:),pointer :: y !< Result
00146 
00147     real(kind=rk),dimension(:),pointer :: xi, yi ! Temporary vectors
00148     real(kind=rk),dimension(:),pointer :: yTemp
00149     integer :: i
00150     logical :: factorised
00151 
00152     y = 0.0
00153     allocate(xi(size(x)))! may be smaller, i.e. max(A%Ai%nrows)
00154     allocate(yi(size(y))) 
00155     allocate(yTemp(size(y)))
00156 
00157     ! y <---(+)---  xTemp <---R^t--- yi <---Ai^(-1)--- xi <---R--- x
00158     
00159     do i=1,size(A%Ai)
00160        factorised = A%subsolve_ids(i)>0
00161        call SpMtx_Ax(xi, A%R(i), x, dozero=.TRUE.)
00162        call sparse_singlesolve(A%subsolve_ids(i), yi, xi, A%Ai(i)%nrows, &
00163             A%Ai(i)%nnz, A%Ai(i)%indi, A%Ai(i)%indj, A%Ai(i)%val)
00164        if(.not.factorised.and.A%subsolve_ids(i)>0) then
00165           A%Ai(i)%indi = A%Ai(i)%indi+1 ! hack, because singlesovle decr by 1
00166           A%Ai(i)%indj = A%Ai(i)%indj+1 ! hack, because singlesovle decr by 1
00167        end if
00168        call SpMtx_Ax(yTemp, A%R(i), yi, dozero=.TRUE., transp=.TRUE.)
00169        y = y+yTemp
00170     end do
00171 
00172     deallocate(xi, yi, yTemp)
00173   end subroutine SOISMtx_pmvm
00174 
00176   subroutine RobustRestrictMtxBuild(A,g,R)
00177     type(SumOfInversedSubMtx), intent(inout) :: A !< system matrix
00178     real(kind=rk), pointer :: g(:) !< solution of the equation \f$ B g = 1 \f$
00179     type(SpMtx), intent(out) :: R !< restriction matrix
00180 
00181     real(kind=rk), pointer :: gi(:), qi(:)
00182     integer :: i, nnz_cur, from, to
00183 
00184     R = SpMtx_newInit(sum(A%R%nnz))
00185 
00186     allocate(gi(size(g)), qi(size(g)))
00187     
00188     nnz_cur = 0
00189     R%nrows = size(A%Ai)
00190     R%ncols = A%R(1)%ncols
00191     do i=1,size(A%Ai)
00192        call SpMtx_Ax(gi, A%R(i), g, dozero=.TRUE.)
00193        call sparse_singlesolve(A%subsolve_ids(i), qi, gi, A%Ai(i)%nrows, &
00194             A%Ai(i)%nnz, A%Ai(i)%indi, A%Ai(i)%indj, A%Ai(i)%val)
00195        !print *, "qi=", qi(1:A%Ai(i)%ncols)
00196        from = nnz_cur + 1
00197        to = nnz_cur + A%R(i)%nnz
00198        R%indi(from:to) = i
00199        R%indj(from:to) = A%R(i)%indj
00200        R%val(from:to) = qi(1:A%R(i)%nnz)
00201        nnz_cur = nnz_cur + A%R(i)%nnz
00202     end do
00203 
00204     if(sctls%verbose>4) then
00205        print *, "R_robust", i, ": nrows, ncols=", R%nrows, R%ncols
00206        call SpMtx_printRaw(R)
00207     end if
00208     
00209     deallocate(gi, qi)
00210   end subroutine RobustRestrictMtxBuild
00211 
00212 end module RobustCoarseMtx_mod