DOUG 0.2

Fact.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 !-------------------------------------------------------
00026 module Fact_class
00027   use globals
00028   use DOUG_utils
00029 
00030 #include<doug_config.h>
00031 
00032 #ifdef D_WANT_MUMPS_YES
00033   include 'dmumps_struc.h'
00034 #endif
00035 
00036   !--------------------------------------------------------------------
00039   integer, parameter :: D_NULLSOLVER = -1 !< internal
00040   integer, parameter :: D_UMFPACK4   = 1  !< UMFPACK v4
00041   integer, parameter :: D_MUMPS      = 2  !< MUMPS package
00042 
00043   !--------------------------------------------------------------------
00048   type Fact
00049     integer :: solver_type !< solver type (D_MUMPS, D_UMFPACK4, etc)
00050     integer :: nfreds !< number of freedom
00051 #ifdef D_WANT_UMFPACK4_YES
00052     integer(kind=pointerk) :: umfpack4_id !< internal UMFPACK4 factorization object handle
00053 #endif
00054 #ifdef D_WANT_MUMPS_YES
00055     type (DMUMPS_STRUC), pointer :: mumps_id !< internal MUMPS factorization object handle
00056 #endif
00057   end type Fact
00058 
00059 contains
00060 
00061   !--------------------------------------------------------------------
00064   function Fact_New(solver_type, nfreds, nnz, indi, indj, val) result(fakt)
00065     implicit none
00066     integer, intent(in) :: solver_type !< what solver to use (D_UMFPACK4, etc)
00067     integer, intent(in) :: nfreds !< number of freedoms in system (matrix rows, columns)
00068     integer, intent(in) :: nnz !< number of non-zero elements in matrix
00069     integer, dimension(:), intent(inout), target :: indi !< matrix rows indices
00070     integer, dimension(:), intent(inout), target :: indj !< matrix column indices
00071     real(kind=rk), dimension(:), intent(in), target :: val !< matrix element values
00072     type(Fact) :: fakt !< factorization object
00073     !-----------------------------
00074     integer :: i, j, status, nz, n
00075     integer(kind=pointerk) :: symbolic
00076 #ifdef D_WANT_UMFPACK4_YES
00077     integer, dimension(:), allocatable :: Ap, Ai, Amap
00078     real(kind=rk), dimension(:), allocatable :: Av
00079     real(kind=rk), dimension(:) :: info90(90),control(20)
00080 #endif
00081 #ifdef D_WANT_MUMPS_YES
00082     type (DMUMPS_STRUC), pointer :: mumps_id => NULL()
00083 #endif
00084 
00085     nz = nnz
00086     n = nfreds
00087     ! If this is degenerate case, use nullsolver
00088     if (nz <= 0 .or. n <= 0) then
00089        fakt%solver_type = D_NULLSOLVER
00090        fakt%nfreds = nfreds
00091        return
00092     end if
00093     if (solver_type == D_UMFPACK4) then
00094 #ifdef D_WANT_UMFPACK4_YES
00095        ! convert input data
00096        indi = indi-1 ! ugly, we mess with input data
00097        indj = indj-1
00098        allocate(Av(nz))
00099        allocate(Ap(n+1))
00100        allocate(Ai(nz))
00101        allocate(Amap(nz))
00102        call umf4triplet2col(n,n,nz,indi,indj,val, &
00103           Ap,Ai,Av,Amap,status)
00104        if (status /= 0) then
00105           write(stream,*)'n,nz:',n,nz
00106           call DOUG_abort('[Fact_New] umf4triplet2col failed',-1)
00107        end if
00108 
00109        ! set default parameters
00110        call umf4def(control)
00111 
00112        ! pre-order and symbolic analysis
00113        symbolic = 0
00114        call umf4sym(n,n,Ap,Ai,Av,symbolic,control,info90)
00115        if (sctls%verbose > 3) then
00116           write(stream,70) info90(1), info90(16),            &
00117             (info90(21) * info90(4)) / 2**20,                &
00118             (info90(22) * info90(4)) / 2**20,                &
00119             info90(23), info90(24), info90(25)                
00120 70        format ('symbolic analysis:',/,                    &
00121             '   status:  ', f5.0, /,                         &
00122             '   time:    ', e10.2, ' (sec)'/,                &
00123             '   estimates (upper bound) for numeric LU:', /, &
00124             '   size of LU:    ', e10.2, ' (MB)', /,         &
00125             '   memory needed: ', e10.2, ' (MB)', /,         &
00126             '   flop count:    ', e10.2, /                   &
00127             '   nnz (L):       ', f10.0, /                   &
00128             '   nnz (U):       ', f10.0)
00129        end if
00130        if (info90(1)<0) then
00131           call DOUG_abort('[Fact_New] error occurred in umf4sym',-1)
00132        end if
00133 
00134        ! numeric factorization
00135        call umf4num(Ap,Ai,Av,symbolic,fakt%umfpack4_id,control,info90)
00136        if (sctls%verbose > 3) then
00137           write(stream,80) info90(1), info90(66),    &
00138             (info90(41) * info90(4)) / 2**20,        &
00139             (info90(42) * info90(4)) / 2**20,        &
00140             info90(43), info90(44), info90(45)        
00141 80        format ('numeric factorization:',/,        &
00142             '   status:  ', f5.0, /,                 &
00143             '   time:    ', e10.2, /,                &
00144             '   actual numeric LU statistics:', /,   &
00145             '   size of LU:    ', e10.2, ' (MB)', /, &
00146             '   memory needed: ', e10.2, ' (MB)', /, &
00147             '   flop count:    ', e10.2, /           &
00148             '   nnz (L):       ', f10.0, /           &
00149             '   nnz (U):       ', f10.0)
00150        end if
00151        ! check umf4num error condition
00152        if (info90(1) < 0) then
00153           write(stream,*)'umf4num info90(1) is:',info90(1)
00154           call umf4pinf(control,info90)
00155           call DOUG_abort('[Fact_New] error occurred in umf4num',-1)
00156        end if
00157 
00158        ! free the symbolic analysis
00159        call umf4fsym(symbolic)
00160        deallocate(Amap)
00161        deallocate(Ai)
00162        deallocate(Ap)
00163        deallocate(Av)
00164 #else
00165        call DOUG_abort('[Fact_New] UMFPACK4 support not compiled in!', -1)
00166 #endif
00167     else if (solver_type == D_MUMPS) then
00168 #ifdef D_WANT_MUMPS_YES
00169        ! Initialize
00170        allocate(fakt%mumps_id)
00171        mumps_id=>fakt%mumps_id
00172        mumps_id%comm=MPI_COMM_WORLD
00173        mumps_id%sym=0
00174        mumps_id%par=1
00175        mumps_id%job=-1
00176        call dmumps(mumps_id)
00177        if (sctls%verbose<3) then
00178           mumps_id%icntl(1)=0
00179           mumps_id%icntl(2)=0
00180           mumps_id%icntl(3)=0
00181        end if
00182        mumps_id%icntl(4)=sctls%verbose ! seems to be sensible 1-1 mapping
00183        if (sctls%verbose>3) mumps_id%icntl(11)=1
00184 
00185        ! pre-order and symbolic analysis and factorization
00186        mumps_id%n=n
00187        mumps_id%nz=nz
00188        mumps_id%irn=>indi
00189        mumps_id%jcn=>indj
00190        mumps_id%a=>val
00191        mumps_id%job=4
00192        call dmumps(mumps_id)
00193 
00194        if (sctls%verbose > 3) then
00195           write(stream, 90) mumps_id%infog(1),         &
00196             (mumps_id%infog(9) + mumps_id%infog(10)) / 2**20,        &
00197             mumps_id%infog(17) / 2**20,        &
00198             mumps_id%rinfog(3), mumps_id%infog(20)        
00199 90        format ('numeric factorization:',/,        &
00200             '   status:  ', f5.0, /,                 &
00201             '   actual numeric LU statistics:', /,   &
00202             '   size of LU:    ', e10.2, ' (MB)', /, &
00203             '   memory needed: ', e10.2, ' (MB)', /, &
00204             '   flop count:    ', e10.2, /           &
00205             '   nnz (total):   ', f10.0)
00206        end if
00207        ! check MUMPS error condition
00208        if (mumps_id%infog(1) < 0) then
00209           call DOUG_abort('[Fact_New] Error occurred in factorization', -1)
00210        end if
00211 #else
00212        call DOUG_abort('[Fact_New] MUMPS support not compiled in!', -1)
00213 #endif
00214     else
00215        call DOUG_abort('[Fact_New] illegal solver type!', -1)
00216     end if
00217     fakt%solver_type = solver_type
00218     fakt%nfreds = nfreds
00219   end function Fact_New
00220 
00221   !--------------------------------------------------------------------
00224   subroutine Fact_Destroy(fakt)
00225     type(Fact) :: fakt !< factorization object to destroy
00226 
00227     if (fakt%solver_type == D_UMFPACK4) then
00228 #ifdef D_WANT_UMFPACK4_YES
00229        call umf4fnum(fakt%umfpack4_id)
00230        fakt%umfpack4_id = 0
00231 #endif
00232     else if (fakt%solver_type == D_MUMPS) then
00233 #ifdef D_WANT_MUMPS_YES
00234        fakt%mumps_id%job = -2 ! free resources
00235        call dmumps(fakt%mumps_id)
00236        deallocate(fakt%mumps_id)
00237        fakt%mumps_id => NULL()
00238 #endif
00239     end if
00240     fakt%solver_type = D_NULLSOLVER
00241   end subroutine Fact_Destroy
00242 
00243   !--------------------------------------------------------------------
00246   subroutine Fact_solve(fakt, rhs, sol)
00247     type(Fact) :: fakt !< factorization object (factorized matrix)
00248     real(kind=rk), dimension(:), pointer :: sol !< solution vector
00249     real(kind=rk), dimension(:), pointer :: rhs !< right-hand side vector
00250     !---------------------
00251 #ifdef D_WANT_UMFPACK4_YES
00252     integer :: sys
00253     real(kind=rk), dimension(:) :: info90(90),control(20)
00254 #endif
00255 #ifdef D_WANT_MUMPS_YES
00256     type (DMUMPS_STRUC), pointer :: mumps_id => NULL()
00257 #endif
00258 
00259     if (fakt%solver_type == D_UMFPACK4) then
00260 #ifdef D_WANT_UMFPACK4_YES
00261       sys = 0
00262       call umf4sol(sys, sol, rhs, fakt%umfpack4_id, control, info90)
00263       if (info90(1) < 0) then
00264          call DOUG_abort('[Fact_solve] error occurred while solving',-1)
00265       end if
00266 #endif
00267     else if (fakt%solver_type == D_MUMPS) then
00268 #ifdef D_WANT_MUMPS_YES
00269       mumps_id => fakt%mumps_id
00270       allocate(mumps_id%rhs(fakt%nfreds))
00271       mumps_id%rhs = rhs
00272       mumps_id%job = 3
00273       call dmumps(mumps_id)
00274       sol = mumps_id%rhs
00275       deallocate(mumps_id%rhs)
00276       if (mumps_id%infog(1) < 0) then
00277          call DOUG_abort('[Fact_solve] error occurred while solving',-1)
00278       end if
00279 #endif
00280     end if
00281   end subroutine Fact_solve
00282 
00283 end module Fact_class