|
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 module CoarseMtx_mod 00023 00024 use DOUG_utils 00025 use RealKind 00026 use SpMtx_class 00027 use Aggregate_mod 00028 use SpMtx_op_AB 00029 use SpMtx_op_Ax 00030 use globals 00031 use SpMtx_operation 00032 00033 use Mesh_class 00034 00035 implicit none 00036 00037 #include<doug_config.h> 00038 00039 #ifdef D_COMPLEX 00040 #define float complex 00041 #else 00042 #define float real 00043 #endif 00044 00045 ! Type(SpMtx), save :: Restrict !,Interp 00046 00047 contains 00048 00050 subroutine IntRestBuild(A,aggr,Restrict,A_ghost) 00051 implicit none 00052 Type(SpMtx), intent(in) :: A !< our fine level matrix 00053 Type(Aggrs), intent(in) :: aggr 00054 Type(SpMtx), intent(out) :: Restrict !< Our restriction matrix 00055 Type(SpMtx),intent(in),optional :: A_ghost !< additional part to the matrix 00056 integer :: nagr,nz,nagrnodes ! # aggregated nodes (there can be some isol.) 00057 integer, dimension(:), allocatable :: indi,indj 00058 integer :: i,j,k 00059 integer :: smoothers=0, Tnrows 00060 Type(SpMtx) :: S,S2,T,SS,SSS,SSSS 00061 float(kind=rk),dimension(:),allocatable :: diag,val 00062 !!!float(kind=rk) :: omega=0.66666666666667_rk 00063 !!!float(kind=rk) :: omega2=0.66666666666667_rk 00064 float(kind=rk) :: omega=0.33333333333333_rk 00065 !!float(kind=rk) :: omega2=0.66666666666667_rk 00066 !!float(kind=rk) :: omega=0.66666666666667_rk 00067 float(kind=rk) :: omega2=0.33333333333333_rk 00068 !float(kind=rk) :: omega=1.33333333333333_rk 00069 !float(kind=rk) :: omega2=1.33333333333333_rk 00070 real(kind=rk),dimension(:),pointer :: sol,rhs 00071 !logical,parameter :: smoothall=.true. 00072 logical,parameter :: smoothall=.false. 00073 00074 if (sctls%verbose>1) write(stream,*) 'Building restriction matrix' 00075 if (sctls%smoothers>sctls%radius1) then 00076 write(stream,*) '***NB! reducing smoothers to radius1=',sctls%radius1 00077 sctls%smoothers=sctls%radius1 00078 endif 00079 Tnrows = A%nrows 00080 if (present(A_ghost)) Tnrows = max(Tnrows, A_ghost%nrows) 00081 smoothers=sctls%smoothers 00082 if (smoothers==0) then 00083 nz=aggr%starts(aggr%nagr+1)-1 ! is actually A%nrows-nisolated 00084 nagr=aggr%nagr 00085 allocate(indi(nz)) 00086 do i=1,nagr 00087 j=aggr%starts(i) 00088 k=aggr%starts(i+1)-1 00089 indi(j:k)=i 00090 enddo 00091 Restrict = SpMtx_newInit( & 00092 nnz=nz, & ! non-overlapping simple case 00093 nblocks=1, & 00094 nrows=nagr, & 00095 ncols=Tnrows, & 00096 symmstruct=.false., & 00097 symmnumeric=.false., & 00098 indi=indi, & 00099 indj=aggr%nodes, & 00100 arrange_type=D_SpMtx_ARRNG_NO, & 00101 M_bound=aggr%starts & 00102 ) 00103 !Restrict%mtx_bbe(2,2)=A%aggr%starts(A%aggr%nagr_local+1)-1 00104 Restrict%val(1:nz)=1.0_rk 00105 !do i=1,nz ! trying averages instead... 00106 ! k=indi(i) ! aggrnumber 00107 ! j=A%aggr%starts(k+1)-A%aggr%starts(k) 00108 ! print *,i,' j=',j 00109 ! Restrict%val(i)=1.0_rk/j 00110 !enddo 00111 00112 !! The transpose 00113 !Interp = SpMtx_New() 00114 !Interp = SpMtx_Init( & 00115 ! nnz=nz, & ! non-overlapping simple case 00116 ! nblocks=1, & 00117 ! nrows=A%nrows, & ! should match for sparse Mult eg. 00118 ! ncols=nagr, & 00119 ! symmstruct=.false., & 00120 ! symmnumeric=.false., & 00121 ! indi=A%aggr%nodes, & 00122 ! indj=indi, & 00123 ! arrange_type=D_SpMtx_ARRNG_NO, & 00124 ! M_bound=A%aggr%starts & 00125 ! ) 00126 !Interp%val(1:nz)=1.0_rk 00127 deallocate(indi) 00128 elseif (smoothers>0) then ! smoothen: 00129 ! build Restrict: 00130 nz=aggr%starts(aggr%nagr+1)-1 ! is actually A%nrows-nisolated 00131 nagr=aggr%nagr 00132 !write(stream,*) "IntRestBuild aggr%nagr", aggr%nagr 00133 allocate(indi(nz)) 00134 do i=1,nagr 00135 j=aggr%starts(i) 00136 k=aggr%starts(i+1)-1 00137 indi(j:k)=i 00138 enddo 00139 T = SpMtx_newInit( & 00140 nnz=nz, & ! non-overlapping simple case 00141 nblocks=1, & 00142 nrows=nagr, & 00143 ncols=Tnrows, & ! should match for sparse Mult eg. 00144 symmstruct=.false., & 00145 symmnumeric=.false., & 00146 indi=indi, & 00147 indj=aggr%nodes, & ! what todo with indi? 00148 arrange_type=D_SpMtx_ARRNG_NO, & 00149 M_bound=aggr%starts & 00150 ) 00151 T%val=1.0_rk 00152 deallocate(indi) 00153 ! Build smoother 00154 allocate(diag(max(A%nrows,A%ncols))) 00155 diag=0.0; 00156 do i=1,A%nnz 00157 if (A%indi(i)==A%indj(i)) then 00158 diag(A%indi(i))=diag(A%indi(i))+A%val(i) 00159 elseif (.not.(A%strong(i).or.smoothall)) then 00160 ! A^epsilon (filtered matrix) has diagonal with added weak connections 00161 diag(A%indi(i))=diag(A%indi(i))+A%val(i) 00162 endif 00163 enddo 00164 if (present(A_ghost).and.associated(A_ghost%indi)) then 00165 do i=1,A_ghost%nnz 00166 if (A_ghost%indi(i)==A_ghost%indj(i)) then 00167 diag(A_ghost%indi(i))=A_ghost%val(i) 00168 endif 00169 enddo 00170 nz=A%nnz+A_ghost%nnz 00171 else 00172 nz=A%nnz 00173 endif 00174 allocate(indi(nz),indj(nz),val(nz)) 00175 nz=0 00176 do i=1,A%nnz 00177 if (A%indi(i)==A%indj(i)) then 00178 nz=nz+1 00179 indi(nz)=A%indi(i) 00180 indj(nz)=A%indj(i) 00181 val(nz)=1.0_rk-omega 00182 elseif (A%strong(i).or.smoothall) then 00183 nz=nz+1 00184 indi(nz)=A%indi(i) 00185 indj(nz)=A%indj(i) 00186 val(nz)=-omega/diag(A%indi(i))*A%val(i) 00187 else 00188 nz=nz+1 00189 indi(nz)=A%indi(i) 00190 indj(nz)=A%indj(i) 00191 val(nz)=0.0_rk 00192 if (sctls%verbose>4) then 00193 write(stream,*)'not strong:',A%indi(i),A%indj(i) 00194 endif 00195 endif 00196 enddo 00197 if (present(A_ghost).and.associated(A_ghost%indi)) then 00198 do i=1,A_ghost%nnz 00199 if (A_ghost%indi(i)==A_ghost%indj(i)) then 00200 nz=nz+1 00201 indi(nz)=A_ghost%indi(i) 00202 indj(nz)=A_ghost%indj(i) 00203 val(nz)=1.0_rk-omega/diag(A_ghost%indi(i))*A_ghost%val(i) 00204 ! todo: to think if these are still needed?: 00205 ! We might actually look at the strength of the opposite conn.in A 00206 elseif (smoothall.or.& 00207 (associated(A_ghost%strong).and.A_ghost%strong(i))) then 00208 nz=nz+1 00209 indi(nz)=A_ghost%indi(i) 00210 indj(nz)=A_ghost%indj(i) 00211 val(nz)=-omega/diag(A_ghost%indi(i))*A_ghost%val(i) 00212 elseif (associated(A_ghost%strong)) then 00213 nz=nz+1 00214 indi(nz)=A_ghost%indi(i) 00215 indj(nz)=A_ghost%indj(i) 00216 val(nz)=0.0_rk 00217 if (sctls%verbose>4) then 00218 write(stream,*)'not strong:',A_ghost%indi(i),A_ghost%indj(i) 00219 endif 00220 endif 00221 enddo 00222 endif 00223 S = SpMtx_newInit( & 00224 nnz=nz, & ! non-overlapping simple case 00225 nblocks=1, & 00226 nrows=Tnrows, & 00227 ncols=Tnrows, & ! should match for sparse Mult eg. 00228 symmstruct=.false., & 00229 symmnumeric=.false., & 00230 indi=indi, & 00231 indj=indj, & 00232 val=val, & 00233 arrange_type=D_SpMtx_ARRNG_NO ) 00234 deallocate(val,indj,indi) 00235 deallocate(diag) 00236 if (smoothers>=2) then 00237 ! Build smoother2 00238 allocate(diag(max(A%nrows,A%ncols))) 00239 do i=1,A%nnz 00240 if (A%indi(i)==A%indj(i)) then 00241 diag(A%indi(i))=A%val(i) 00242 elseif (.not.(A%strong(i).or.smoothall)) then 00243 ! A^epsilon (filtered matrix) has diagonal with added weak connections 00244 diag(A%indi(i))=diag(A%indi(i))+A%val(i) 00245 endif 00246 enddo 00247 if (present(A_ghost).and.associated(A_ghost%indi)) then 00248 do i=1,A_ghost%nnz 00249 if (A_ghost%indi(i)==A_ghost%indj(i)) then 00250 diag(A_ghost%indi(i))=A_ghost%val(i) 00251 endif 00252 enddo 00253 nz=A%nnz+A_ghost%nnz 00254 else 00255 nz=A%nnz 00256 endif 00257 allocate(indi(nz),indj(nz),val(nz)) 00258 nz=0 00259 do i=1,A%nnz 00260 if (A%indi(i)==A%indj(i)) then 00261 nz=nz+1 00262 indi(nz)=A%indi(i) 00263 indj(nz)=A%indj(i) 00264 val(nz)=1.0_rk-omega2 00265 elseif (A%strong(i).or.smoothall) then 00266 nz=nz+1 00267 indi(nz)=A%indi(i) 00268 indj(nz)=A%indj(i) 00269 val(nz)=-omega2/diag(A%indi(i))*A%val(i) 00270 else 00271 nz=nz+1 00272 indi(nz)=A%indi(i) 00273 indj(nz)=A%indj(i) 00274 val(nz)=0.0_rk 00275 if (sctls%verbose>4) then 00276 write(stream,*)'not strong:',A%indi(i),A%indj(i) 00277 endif 00278 endif 00279 enddo 00280 if (present(A_ghost).and.associated(A_ghost%indi)) then 00281 do i=1,A_ghost%nnz 00282 if (A_ghost%indi(i)==A_ghost%indj(i)) then 00283 nz=nz+1 00284 indi(nz)=A_ghost%indi(i) 00285 indj(nz)=A_ghost%indj(i) 00286 val(nz)=1.0_rk-omega2/diag(A_ghost%indi(i))*A_ghost%val(i) 00287 ! We might actually look at the strength of the opposite conn.in A 00288 elseif (smoothall.or.& 00289 (associated(A_ghost%strong).and.A_ghost%strong(i))) then 00290 nz=nz+1 00291 indi(nz)=A_ghost%indi(i) 00292 indj(nz)=A_ghost%indj(i) 00293 val(nz)=-omega2/diag(A_ghost%indi(i))*A_ghost%val(i) 00294 elseif (associated(A_ghost%strong)) then 00295 nz=nz+1 00296 indi(nz)=A_ghost%indi(i) 00297 indj(nz)=A_ghost%indj(i) 00298 val(nz)=0.0_rk 00299 if (sctls%verbose>4) then 00300 write(stream,*)'not strong:',A_ghost%indi(i),A_ghost%indj(i) 00301 endif 00302 endif 00303 enddo 00304 endif 00305 S2 = SpMtx_newInit( & 00306 nnz=nz, & ! non-overlapping simple case 00307 nblocks=1, & 00308 nrows=Tnrows, & 00309 ncols=Tnrows, & ! should match for sparse Mult eg. 00310 symmstruct=.false., & 00311 symmnumeric=.false., & 00312 indi=indi, & 00313 indj=indj, & 00314 val=val, & 00315 arrange_type=D_SpMtx_ARRNG_NO ) 00316 deallocate(val,indj,indi) 00317 deallocate(diag) 00318 SS=SpMtx_AB(A=S2, & 00319 B=S, & 00320 AT=.false., & 00321 BT=.false.) 00322 if (smoothers==3) then 00323 SSS=SpMtx_AB(A=SS, & 00324 B=S2, & 00325 AT=.false., & 00326 BT=.false.) 00327 Restrict = SpMtx_AB(A=T, & 00328 B=SSS, & 00329 AT=.false., & 00330 BT=.false.) 00331 call SpMtx_Destroy(SSS) 00332 elseif (smoothers==4) then 00333 SSSS=SpMtx_AB(A=SS, & 00334 B=SS, & 00335 AT=.false., & 00336 BT=.false.) 00337 Restrict = SpMtx_AB(A=T, & 00338 B=SSSS, & 00339 AT=.false., & 00340 BT=.false.) 00341 call SpMtx_Destroy(SSSS) 00342 elseif (smoothers==5) then 00343 SSSS=SpMtx_AB(A=SS, & 00344 B=SS, & 00345 AT=.false., & 00346 BT=.false.) 00347 SSS=SpMtx_AB(A=SSSS, & 00348 B=S2, & 00349 AT=.false., & 00350 BT=.false.) 00351 Restrict = SpMtx_AB(A=T, & 00352 B=SSS, & 00353 AT=.false., & 00354 BT=.false.) 00355 call SpMtx_Destroy(SSSS) 00356 call SpMtx_Destroy(SSS) 00357 elseif (smoothers==6) then 00358 SSSS=SpMtx_AB(A=SS, & 00359 B=SS, & 00360 AT=.false., & 00361 BT=.false.) 00362 SSS=SpMtx_AB(A=SSSS, & 00363 B=SS, & 00364 AT=.false., & 00365 BT=.false.) 00366 Restrict = SpMtx_AB(A=T, & 00367 B=SSS, & 00368 AT=.false., & 00369 BT=.false.) 00370 call SpMtx_Destroy(SSSS) 00371 call SpMtx_Destroy(SSS) 00372 elseif (smoothers==7) then 00373 SSSS=SpMtx_AB(A=SS, & !4 00374 B=SS, & 00375 AT=.false., & 00376 BT=.false.) 00377 SSS=SpMtx_AB(A=SSSS, & !6 00378 B=SS, & 00379 AT=.false., & 00380 BT=.false.) 00381 call SpMtx_Destroy(SSSS) !7 00382 SSSS=SpMtx_AB(A=SSS, & 00383 B=S2, & 00384 AT=.false., & 00385 BT=.false.) 00386 call SpMtx_Destroy(SSS) 00387 Restrict = SpMtx_AB(A=T, & 00388 B=SSSS, & 00389 AT=.false., & 00390 BT=.false.) 00391 call SpMtx_Destroy(SSSS) 00392 elseif (smoothers==8) then 00393 SSSS=SpMtx_AB(A=SS, & !4 00394 B=SS, & 00395 AT=.false., & 00396 BT=.false.) 00397 SSS=SpMtx_AB(A=SSSS, & !8 00398 B=SSSS, & 00399 AT=.false., & 00400 BT=.false.) 00401 Restrict = SpMtx_AB(A=T, & 00402 B=SSS, & 00403 AT=.false., & 00404 BT=.false.) 00405 call SpMtx_Destroy(SSSS) 00406 call SpMtx_Destroy(SSS) 00407 else ! smoothers==2 00408 if (smoothers>8) then 00409 write(stream,*) ' Warning: smoothers ',smoothers, & 00410 ' not supported, performing 2' 00411 endif 00412 Restrict = SpMtx_AB(A=T, & 00413 B=SS, & 00414 AT=.false., & 00415 BT=.true.) 00416 endif 00417 call SpMtx_Destroy(SS) 00418 else ! i.e. smoothers==1 : 00419 !write(stream,*)'bbb building Restrict...' 00420 !write(stream,*)'A%ncols===',A%ncols,maxval(A%indj) 00421 if (sctls%verbose>3) write(stream,*) "Restrict = T*S" 00422 Restrict = SpMtx_AB(A=T, & 00423 B=S, & 00424 AT=.false., & 00425 BT=.true.) 00426 !write(stream,*)'bbb done building Restrict...' 00427 endif 00428 call SpMtx_Destroy(S) 00429 if (smoothers>=2) then 00430 call SpMtx_Destroy(S2) 00431 endif 00432 call SpMtx_Destroy(T) 00433 elseif (smoothers==-1) then ! use exact smoother through solves on aggr, 00434 ! which are non-overlapping: 00435 !moved this to IntRestBuild2... 00436 ! 00437 !! build Restrict: 00438 !nz=A%aggr%starts(A%aggr%nagr+1)-1 ! is actually A%nrows-nisolated 00439 !nagr=A%aggr%nagr 00440 !allocate(indi(nz)) 00441 !do i=1,nagr 00442 ! j=A%aggr%starts(i) 00443 ! k=A%aggr%starts(i+1)-1 00444 ! indi(j:k)=i 00445 !enddo 00446 !Restrict = SpMtx_New() 00447 !Restrict = SpMtx_Init( & 00448 ! nnz=nz, & ! non-overlapping simple case 00449 ! nblocks/home/eero/share/AMG/Hetero1.txt=1, & 00450 ! nrows=nagr, & 00451 ! ncols=A%nrows, & ! should match for sparse Mult eg. 00452 ! symmstruct=.false., & 00453 ! symmnumeric=.false., & 00454 ! indi=indi, & 00455 ! indj=A%aggr%nodes, & ! what todo with indi? 00456 ! arrange_type=D_SpMtx_ARRNG_NO, & 00457 ! M_bound=A%aggr%starts & 00458 ! ) 00459 !deallocate(indi) 00460 !! Build smoother 00461 !allocate(rhs(A%nrows)) 00462 !allocate(sol(A%nrows)) 00463 !rhs=1.0_rk 00464 !sol=0.0_rk 00465 !call exact_sparse_multismoother(sol,A,rhs) 00466 !!print *,'sol=',sol 00467 !!stop 00468 !!Restrict%val=1.0_rk 00469 !Restrict%val=sol 00470 !!Restrict%val=1.0_rk/sol 00471 !deallocate(sol) 00472 !deallocate(rhs) 00473 endif 00474 end subroutine IntRestBuild 00475 00476 subroutine CoarseMtxBuild(A,AC,Restrict,ninner,A_ghost) 00477 Type(SpMtx),intent(inout) :: A ! the fine level matrix 00478 Type(SpMtx),intent(inout) :: AC ! coarse level matrix 00479 Type(SpMtx), intent(inout) :: Restrict ! the restriction matrix 00480 integer,intent(in) :: ninner !< number of inner nodes 00481 Type(SpMtx),intent(in),optional :: A_ghost !< additional part to the matrix 00482 Type(SpMtx) :: T,TT,RT !temporary matrix 00483 integer,dimension(:),pointer :: indi,indj 00484 real(kind=rk),dimension(:),pointer :: val 00485 integer :: i,nz 00486 00487 if (sctls%verbose>1) write(stream,*) 'Building coarse matrix' 00488 ! we need to work with a copy to preserve the structure and ordering of A: 00489 if (present(A_ghost).and.associated(A_ghost%indi)) then 00490 nz=A%nnz+A_ghost%nnz 00491 TT=SpMtx_newInit(nz) 00492 TT%indi(1:A%nnz)=A%indi 00493 TT%indj(1:A%nnz)=A%indj 00494 TT%val(1:A%nnz)=A%val 00495 TT%indi(A%nnz+1:nz)=A_ghost%indi 00496 TT%indj(A%nnz+1:nz)=A_ghost%indj 00497 TT%val(A%nnz+1:nz)=A_ghost%val 00498 TT%nrows=maxval(TT%indi) 00499 TT%ncols=maxval(TT%indj) 00500 else 00501 TT=SpMtx_Copy(A) 00502 endif 00503 T = SpMtx_AB(A=TT, & 00504 B=Restrict, & 00505 AT=.false., & 00506 BT=.true.) 00507 call SpMtx_Destroy(TT) 00508 RT = SpMtx_Copy(Restrict) 00509 if (sctls%input_type==DISTRIBUTION_TYPE_ASSEMBLED.or.& 00510 sctls%input_type==DISTRIBUTION_TYPE_STRUCTURED) then 00511 call KeepGivenColumnIndeces(RT,(/(i,i=1,ninner)/),.TRUE.) 00512 end if 00513 AC = SpMtx_AB(A=RT,B=T) 00514 call SpMtx_Destroy(RT) 00515 call SpMtx_Destroy(T) 00516 end subroutine CoarseMtxBuild 00517 00518 end module CoarseMtx_mod
1.7.3-20110217