|
DOUG 0.2
|
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
1.7.3-20110217