DOUG 0.2

subsolvers.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 !!----------------------------------
00023 !! Preconditioned Conjugate Gradient
00024 !!----------------------------------
00025 module subsolvers
00026   use globals
00027   use DOUG_utils
00028   use Fact_class
00029   use SpMtx_class
00030   use SpMtx_arrangement
00031   use Decomposition_mod
00032   implicit none
00033 
00034 #include<doug_config.h>
00035 
00036 #ifdef D_COMPLEX
00037 #define float complex
00038 #else
00039 #define float real
00040 #endif
00041 
00042   integer :: nfacts=0 
00043   integer :: maxnfacts=0
00044 #ifdef D_WANT_UMFPACK4_YES
00045   integer :: subsolver=D_UMFPACK4
00046 #else
00047   integer :: subsolver=D_MUMPS
00048 #endif
00049   real(kind=rk) :: setuptime=0.0_rk, factorisation_time=0.0_rk, backsolve_time=0.0_rk
00050   type(Fact), dimension(:), pointer :: fakts => NULL()
00051 
00052 contains
00053 
00054   subroutine Solve_subdomains(sol,DD,subsolve_ids,rhs)
00055     type(Decomposition),intent(in) :: DD !< subdomains
00056     integer, dimension(:), pointer   :: subsolve_ids !< numeric object handles of (UMFPACK,...) factorizations
00057     real(kind=rk),dimension(:),pointer :: sol,rhs
00058 
00059     ! ----- local: ------
00060     integer :: sd,n,nsubsolves
00061 
00062     real(kind=rk),pointer :: subrhs(:),subsol(:)
00063 
00064     ! solve
00065     nsubsolves = size(subsolve_ids)
00066     allocate(subrhs(maxval(DD%subd(1:nsubsolves)%ninds)))
00067     allocate(subsol(maxval(DD%subd(1:nsubsolves)%ninds)))
00068     do sd=1,nsubsolves
00069       n=DD%subd(sd)%ninds
00070       subrhs(1:n)=rhs(DD%subd(sd)%inds(1:n))
00071       call Fact_solve(fakts(subsolve_ids(sd)),subrhs,subsol)
00072       sol(DD%subd(sd)%inds(1:n))=sol(DD%subd(sd)%inds(1:n))+subsol(1:n)
00073     enddo
00074     deallocate(subrhs,subsol)
00075 
00076   end subroutine Solve_subdomains
00077 
00078   subroutine Factorise_subdomains(DD,A,A_ghost,subsolve_ids)
00079     type(Decomposition),intent(in) :: DD
00080     type(SpMtx),intent(in) :: A
00081     type(SpMtx),optional :: A_ghost
00082     integer, dimension(:), pointer   :: subsolve_ids !< numeric object handles of (UMFPACK,...) factorizations
00083 
00084     integer :: cAggr, ol, i, nagr, nsubsolves
00085     double precision :: t1
00086 
00087     setuptime=0.0_rk
00088     factorisation_time=0.0_rk
00089     t1 = MPI_WTime()
00090     
00091     nsubsolves = size(DD%subd)
00092     allocate(subsolve_ids(nsubsolves))
00093     subsolve_ids = 0
00094 
00095     do i=1,nsubsolves
00096       setuptime=setuptime+(MPI_WTIME()-t1) ! switchoff clock
00097       call Factorise_subdomain(A,DD%subd(i)%inds,subsolve_ids(i),A_ghost)
00098       t1 = MPI_WTIME() ! switchon clock
00099     end do
00100 
00101     !deallocate(subrhs,subsol)
00102   end subroutine factorise_subdomains
00103 
00104   subroutine Factorise_subdomain(A,nodes,id,A_ghost)
00105     type(SpMtx),intent(in) :: A !< matrix
00106     integer,intent(in) :: nodes(:) !< subdomain nodes
00107     integer,intent(out) :: id
00108     type(SpMtx),intent(in),optional :: A_ghost !< ghost matrix values
00109     
00110     integer :: snnz,i,node,nnz,nnodes,nallnodes
00111     integer,dimension(:),allocatable :: floc,sindi,sindj
00112     real(kind=rk),dimension(:),allocatable :: sval
00113 
00114     if(present(A_ghost)) then
00115       nallnodes = max(A%nrows,A_ghost%nrows)
00116       allocate(floc(nallnodes))
00117       nnz = A%nnz+A_ghost%nnz
00118     else 
00119       nallnodes = A%nrows
00120       allocate(floc(nallnodes))
00121       nnz = A%nnz
00122     endif
00123     nnodes = size(nodes)
00124 
00125     ! create global to local array for the subdomain
00126     floc=0
00127     do i=1,nnodes
00128       node=nodes(i) ! node number
00129       floc(node)=i
00130     enddo
00131 
00132     ! now we have gathered information for subdomain
00133     allocate(sindi(nnz),sindj(nnz),sval(nnz))
00134     snnz=0
00135     do i=1,A%nnz
00136       ! todo:
00137       !   need to be avoided going it all through again and again?
00138       ! idea: to use linked-list arrays as in old DOUG
00139       if (A%indi(i)<=nallnodes.and.A%indj(i)<=nallnodes) then
00140         if (floc(A%indi(i))>0.and.floc(A%indj(i))>0) then ! wheather in the subdomain?
00141           snnz=snnz+1
00142           sindi(snnz)=floc(A%indi(i))
00143           sindj(snnz)=floc(A%indj(i))
00144           sval(snnz)=A%val(i)
00145         endif
00146       endif
00147     enddo
00148 
00149     ! do the same for the ghost matrix
00150     if(present(A_ghost)) then
00151       do i=1,A_ghost%nnz
00152         if (floc(A_ghost%indi(i))>0.and.floc(A_ghost%indj(i))>0) then ! wheather in the subdomain?
00153           snnz=snnz+1
00154           sindi(snnz)=floc(A_ghost%indi(i))
00155           sindj(snnz)=floc(A_ghost%indj(i))
00156           sval(snnz)=A_ghost%val(i)
00157         endif
00158       enddo
00159     end if
00160 
00161     ! factorise
00162     call factorise(id,nnodes,snnz,sindi,sindj,sval)
00163   end subroutine factorise_subdomain
00164 
00165   subroutine sparse_singlesolve(id,sol,rhs,nfreds,nnz,indi,indj,val,tot_nfreds,nnz_est)
00166     ! For adding a new factorised matrix id must be 0.
00167     !   id > 0 will be returned as a handle for the factors
00168     integer,intent(inout) :: id
00169     real(kind=rk),dimension(:),pointer :: sol,rhs
00170     integer,intent(in) :: nfreds
00171     integer,intent(in),optional :: nnz
00172     integer,dimension(:),intent(inout),optional :: indi,indj
00173     real(kind=rk),dimension(:),optional :: val
00174     integer,intent(in),optional :: tot_nfreds
00175     integer,intent(in),optional :: nnz_est ! estimate for aver. nnz
00176 
00177     call factorise_and_solve(id,sol,rhs,nfreds,nnz,indi,indj,val)
00178   end subroutine sparse_singlesolve
00179 
00180   subroutine factorise(id,nfreds,nnz,indi,indj,val)
00181     ! For adding a new factorised matrix id must be 0.
00182     !   id > 0 will be returned as a handle for the factors
00183     integer,intent(out) :: id
00184     integer,intent(in) :: nfreds
00185     integer,intent(in),optional :: nnz
00186     integer,dimension(:),intent(inout),optional :: indi,indj
00187     real(kind=rk),dimension(:),intent(in),optional :: val
00188 
00189     ! ---- local -----
00190     integer :: i,nz,n
00191     integer :: fakts_size
00192     type(Fact),dimension(:),pointer :: fakts_temp
00193     real(kind=rk) :: t1
00194 
00195     id=-1
00196 
00197       if (present(nnz)) then
00198         nz=nnz
00199       elseif (present(val)) then
00200         nz=size(val)
00201       else
00202         call DOUG_abort('[factorise_and_solve] unable to get nnz', -1)
00203       endif
00204       n=nfreds
00205       if (.not.present(indi)) call DOUG_abort('[factorise_and_solve] indi must be given',-1)
00206       if (.not.present(indj)) call DOUG_abort('[factorise_and_solve] indj must be given',-1)
00207       if (.not.present(val))  call DOUG_abort('[factorise_and_solve] val must be given',-1)
00208 
00209       fakts_size = 0
00210       if (associated(fakts)) fakts_size = size(fakts)
00211       do i=1,fakts_size
00212         if (fakts(i)%solver_type == 0) then
00213           id=i
00214           exit
00215         endif
00216       enddo
00217       if (id<=0) then
00218         allocate(fakts_temp(fakts_size*3/2+5))
00219         fakts_temp(1:size(fakts_temp))%solver_type=0
00220         if (associated(fakts)) fakts_temp(1:fakts_size)=fakts
00221         id=fakts_size+1
00222         if (associated(fakts)) deallocate(fakts)
00223         fakts=>fakts_temp
00224       else
00225         if (id>size(fakts)) call DOUG_abort('[factorise_and_solve] id is out of range', -1)
00226       endif
00227 
00228       t1=MPI_WTIME()      
00229       fakts(id)=Fact_New(subsolver, n, nz, indi, indj, val)
00230       factorisation_time=factorisation_time + MPI_WTIME()-t1
00231   end subroutine factorise
00232 
00233 
00234   subroutine factorise_and_solve(id,sol,rhs,nfreds,nnz,indi,indj,val)
00235     ! For adding a new factorised matrix id must be 0.
00236     !   id > 0 will be returned as a handle for the factors
00237 !use SpMtx_util
00238     integer,intent(inout) :: id
00239     real(kind=rk),dimension(:),pointer :: sol,rhs
00240     integer,intent(in) :: nfreds
00241     integer,intent(in),optional :: nnz
00242     integer,dimension(:),intent(inout),optional :: indi,indj
00243     real(kind=rk),dimension(:),intent(in),optional :: val
00244 
00245     ! ---- local -----
00246     integer :: i,nz,n
00247     integer :: fakts_size
00248     type(Fact),dimension(:),pointer :: fakts_temp
00249     real(kind=rk) :: t1
00250 
00251     ! ---- for check ----
00252     !logical,parameter :: check=.true.
00253     logical,parameter :: check=.false.
00254     logical :: docheck=.false.
00255     real(kind=rk),dimension(:),pointer :: chk
00256     real(kind=rk) :: rtmp
00257 
00258     if (check.and.id<=0) then
00259       docheck=.true.
00260     else
00261       docheck=.false.
00262     endif
00263     if (id<=0) then
00264       call factorise(id,nfreds,nnz,indi,indj,val)
00265       indi=indi+1
00266       indj=indj+1
00267     end if
00268 
00269     t1=MPI_WTIME()      
00270     call Fact_Solve(fakts(id), rhs, sol)
00271     backsolve_time=backsolve_time + MPI_WTIME()-t1
00272 
00273     if (docheck) then
00274       allocate(chk(size(rhs)))
00275       chk=0.0_rk
00276       do i=1,nnz
00277         chk(indi(i)+1)=chk(indi(i)+1)+val(i)*sol(indj(i)+1)
00278       enddo
00279       rtmp=0.0_rk
00280       do i=1,size(sol)
00281         rtmp=rtmp+abs(rhs(i)-chk(i))
00282       enddo
00283       deallocate(chk)
00284       write(stream,*)'******* subdomain solution error^2:',rtmp
00285     endif
00286   end subroutine factorise_and_solve
00287 
00288   function total_setup_time() result(t)
00289     real(kind=rk) :: t
00290     t=setuptime+factorisation_time
00291   end function total_setup_time
00292 
00293   function total_factorisation_time() result(t)
00294     real(kind=rk) :: t
00295     t=factorisation_time
00296   end function total_factorisation_time
00297 
00298   function total_backsolve_time() result(t)
00299     real(kind=rk) :: t
00300     t=backsolve_time
00301   end function total_backsolve_time
00302 
00303 end module subsolvers
00304