DOUG 0.2

CoarseMtx.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 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