Changeset e083829


Ignore:
Timestamp:
03/05/11 13:48:20 (2 years ago)
Author:
Oleg Batrashev <ogbash@…>
Branches:
master, external
Children:
e8c8f58
Parents:
ea465d3
git-author:
Oleg Batrashev <ogbash@…> (03/05/11 13:48:20)
git-committer:
Oleg Batrashev <ogbash@…> (03/05/11 13:48:20)
Message:

Get rid of CoarseMtx_ and Restrict in pcg.

Location:
src
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • src/components/Preconditioner_base.F90

    rea465d3 re083829  
    4848 
    4949  ! -------- Coarse preconditioner 
     50  integer,parameter :: COARSE_PRECONDITIONER_TYPE_NONE=0, & 
     51       COARSE_PRECONDITIONER_TYPE_SMOOTH=1, & 
     52       COARSE_PRECONDITIONER_TYPE_GEOMETRIC=2, & 
     53       COARSE_PRECONDITIONER_TYPE_ROBUST=3 
    5054 
    5155  !> Base type for fine level preconditioner. 
    5256  type CoarsePreconditioner 
     57    integer :: type !< coarse preconditioner type 
    5358    type(CoarseData) :: cdat !<coarse data -- includes overlap 
    5459    type(CoarseData) :: cdat_vec !<coarse data -- w/o overlap, for vector collects 
     60    type(SpMtx) :: R !< Restriction matrix 
     61    type(SpMtx) :: AC !< Coarse matrix 
    5562 
    5663    ! implementations 
    … …  
    9198  end subroutine FinePreconditioner_InitAggrs 
    9299 
     100  function CoarsePreconditioner_New() result (CP) 
     101    type(CoarsePreconditioner) :: CP 
     102 
     103    CP%type = COARSE_PRECONDITIONER_TYPE_NONE 
     104 
     105  end function CoarsePreconditioner_New 
     106 
    93107end module Preconditioner_base_mod 

  • src/main/aggr.F90

    rea465d3 re083829  
    9292#endif 
    9393 
    94   type(SpMtx)    :: AC  !< coarse matrix 
    95   type(SpMtx)    :: Restrict !< Restriction matrix (for operation) 
    9694  type(SumOfInversedSubMtx) :: B_RCS !< B matrix for the Robust Coarse Spaces 
    9795  !type(SpMtx)    :: Rest_cmb !< Restriction matrix (for coarse matrix build) 
    … …  
    124122  integer,pointer :: aggrnum(:) 
    125123  integer :: nagr 
    126   ! Parallel coarse level 
    127   !type(CoarseData) :: CP%cdat --  
    128124 
    129125  ! Init DOUG 
    … …  
    157153    endif 
    158154     
     155    CP = CoarsePreconditioner_New() 
    159156    ! Testing coarse matrix and aggregation through it: 
    160157    if (numprocs>1) then 
     158      CP%type = COARSE_PRECONDITIONER_TYPE_SMOOTH 
    161159      allocate(aggrnum(D%mesh%nlf)) 
    162160      nagr = P%fAggr%inner%nagr 
    … …  
    170168      call SpMtx_unscale(D%A) 
    171169 
    172       call IntRestBuild(D%A,P%fAggr%inner,Restrict,D%A_ghost) 
    173       CS = CoarseSpace_Init(Restrict) 
     170      call IntRestBuild(D%A,P%fAggr%inner,CP%R,D%A_ghost) 
     171      CS = CoarseSpace_Init(CP%R) 
    174172      call CoarseData_Copy(CP%cdat,CP%cdat_vec) 
    175       call CoarseSpace_Expand(CS,Restrict,D%mesh,CP%cdat) 
    176       call CoarseMtxBuild(D%A,CP%cdat%LAC,Restrict,D%mesh%ninner,D%A_ghost) 
    177       call KeepGivenRowIndeces(Restrict, (/(i,i=1,P%fAggr%inner%nagr)/)) 
     173      call CoarseSpace_Expand(CS,CP%R,D%mesh,CP%cdat) 
     174      call CoarseMtxBuild(D%A,CP%cdat%LAC,CP%R,D%mesh%ninner,D%A_ghost) 
     175      call KeepGivenRowIndeces(CP%R, (/(i,i=1,P%fAggr%inner%nagr)/)) 
    178176 
    179177      if (sctls%verbose>3.and.CP%cdat%LAC%nnz<400) then 
    180         write(stream,*)'Restrict (local) is:==================' 
    181         call SpMtx_printRaw(A=Restrict) 
     178        write(stream,*)'CP%R (local) is:==================' 
     179        call SpMtx_printRaw(A=CP%R) 
    182180        write(stream,*)'A coarse (local) is:==================' 
    183181        call SpMtx_printRaw(A=CP%cdat%LAC) 
    … …  
    186184    else  
    187185      ! non-parallel 
    188       call IntRestBuild(D%A,P%fAggr%inner,Restrict) 
     186      call IntRestBuild(D%A,P%fAggr%inner,CP%R) 
    189187 
    190188      if (sctls%coarse_method<=1) then ! if not specified or ==1 
    191          call CoarseMtxBuild(D%A,AC,Restrict,D%mesh%ninner) 
     189        CP%type = COARSE_PRECONDITIONER_TYPE_SMOOTH 
     190        call CoarseMtxBuild(D%A,CP%AC,CP%R,D%mesh%ninner) 
    192191 
    193192      else if (sctls%coarse_method==2) then 
    194          ! use the Robust Coarse Spaces algorithm 
    195          B_RCS = CoarseProjectionMtxsBuild(D%A,Restrict,P%fAggr%inner%nagr) 
     193        ! use the Robust Coarse Spaces algorithm 
     194        CP%type = COARSE_PRECONDITIONER_TYPE_ROBUST 
     195 
     196         B_RCS = CoarseProjectionMtxsBuild(D%A,CP%R,P%fAggr%inner%nagr) 
    196197         allocate(rhs_1(D%A%nrows)) 
    197198         allocate(g(D%A%ncols)) 
    … …  
    204205         end if 
    205206 
    206          call RobustRestrictMtxBuild(B_RCS,g,Restrict) 
    207          call CoarseMtxBuild(D%A,AC,Restrict,D%mesh%ninner) 
     207         call RobustRestrictMtxBuild(B_RCS,g,CP%R) 
     208         call CoarseMtxBuild(D%A,CP%AC,CP%R,D%mesh%ninner) 
    208209 
    209210      else 
    … …  
    212213      endif 
    213214            
    214       if (sctls%verbose>3.and.AC%nnz<400) then 
     215      if (sctls%verbose>3.and.CP%AC%nnz<400) then 
    215216        write(stream,*)'A coarse is:==================' 
    216         call SpMtx_printRaw(A=AC) 
     217        call SpMtx_printRaw(A=CP%AC) 
    217218      endif 
    218219    endif 
    … …  
    220221    ! coarse aggregates 
    221222    if (numprocs==1) then 
    222       call Partitionings_CreateCoarse(P,D,AC) 
     223      call Partitionings_CreateCoarse(P,D,CP%AC) 
    223224      !call Aggrs_readFile_coarse(P%cAggr, "aggregates.txt") 
    224225 
    … …  
    283284     ! Preconditioned conjugate gradient 
    284285     t1 = MPI_WTIME() 
    285      if (sctls%levels==2) then 
    286        write(stream,*)'calling pcg_weigs with coarse matrix' 
    287        call pcg_weigs(A=D%A,b=D%rhs,x=xl,Msh=D%mesh,finePrec=FP,coarsePrec=CP,it=it,cond_num=cond_num, & 
    288             A_interf_=D%A_ghost, & 
    289             CoarseMtx_=AC,Restrict=Restrict, & 
     286     write(stream,*)'calling pcg_weigs' 
     287     call pcg_weigs(A=D%A,b=D%rhs,x=xl,Msh=D%mesh,& 
     288          finePrec=FP,coarsePrec=CP,& 
     289          it=it,cond_num=cond_num, A_interf_=D%A_ghost, & 
    290290            refactor_=.true.) 
    291      else 
    292        write(stream,*)'calling pcg_weigs' 
    293        call pcg_weigs(D%A, D%rhs, xl, D%mesh,FP,CP,it,cond_num, & 
    294             A_interf_=D%A_ghost, refactor_=.true.) 
    295      endif 
    296291     t=MPI_WTIME()-t1 
    297292     write(stream,*) 'time spent in pcg():',t 

  • src/main/geom.F90

    rea465d3 re083829  
    5656#endif 
    5757 
    58   type(SpMtx)    :: AC  ! Coarse matrix 
    59   type(SpMtx)    :: Restrict ! Restriction matrix 
    60  
    6158  float(kind=rk), dimension(:), pointer :: xl ! local solution vector 
    6259  float(kind=rk), dimension(:), pointer :: x  ! global solution on master 
    … …  
    132129    if (ismaster()) then 
    133130      if (sctls%verbose>0) write (stream,*) "Building coarse grid" 
    134  
     131      CP = CoarsePreconditioner_New() 
     132      CP%type = COARSE_PRECONDITIONER_TYPE_GEOMETRIC 
     133       
    135134      call CreateCoarse(D%mesh,C) 
    136135 
    … …  
    159158 
    160159      if (sctls%verbose>0) write (stream,*) "Creating Restriction matrix" 
    161       call CreateRestrict(LC,D%mesh,Restrict) 
     160      call CreateRestrict(LC,D%mesh,CP%R) 
    162161 
    163162      if (sctls%verbose>1) write (stream,*) "Cleaning Restriction matrix" 
    164       call CleanCoarse(LC,Restrict,D%mesh) 
     163      call CleanCoarse(LC,CP%R,D%mesh) 
    165164 
    166165      if (sctls%verbose>0)  write (stream,*) "Building coarse matrix" 
    167       call CoarseMtxBuild(D%A,CP%cdat%LAC,Restrict,D%mesh%ninner) 
     166      call CoarseMtxBuild(D%A,CP%cdat%LAC,CP%R,D%mesh%ninner) 
    168167 
    169168      if (sctls%verbose>1) write (stream, *) "Stripping the restriction matrix" 
    170       call StripRestrict(D%mesh,Restrict) 
     169      call StripRestrict(D%mesh,CP%R) 
    171170 
    172171      if (sctls%verbose>0) write (stream,*) "Transmitting local-to-global maps" 
    … …  
    205204     t1 = MPI_WTIME() 
    206205 
    207      if (sctls%input_type==DCTL_INPUT_TYPE_ELEMENTAL .and. & 
    208                         sctls%levels==2) then 
    209        call pcg_weigs(A=D%A,b=D%rhs,x=xl,Msh=D%mesh,finePrec=FP,coarsePrec=CP,it=it,cond_num=cond_num, & 
    210                       A_interf_=D%A_ghost,CoarseMtx_=AC,Restrict=Restrict, & 
    211                       refactor_=.true.) 
    212 !                     refactor_=.true., cdat_=CP%cdat) 
    213      else 
    214        call pcg_weigs(A=D%A,b=D%rhs,x=xl,Msh=D%mesh,finePrec=FP,coarsePrec=CP,it=it,cond_num=cond_num, & 
    215                         A_interf_=D%A_ghost,refactor_=.true.) 
    216      endif 
     206     call pcg_weigs(A=D%A,b=D%rhs,x=xl,Msh=D%mesh,finePrec=FP,coarsePrec=CP,it=it,cond_num=cond_num, & 
     207          A_interf_=D%A_ghost, & 
     208          refactor_=.true.) 
    217209 
    218210     write(stream,*) 'time spent in pcg():',MPI_WTIME()-t1 
    … …  
    274266 
    275267  if (sctls%input_type==DCTL_INPUT_TYPE_ELEMENTAL .and. sctls%levels==2) then 
    276       call SpMtx_Destroy(AC) 
    277       call SpMtx_Destroy(Restrict) 
     268      call SpMtx_Destroy(CP%AC) 
     269      call SpMtx_Destroy(CP%R) 
    278270!      call SpMtx_Destroy(Res_aux) 
    279271      call SendData_Destroy(CP%cdat%send) 

  • src/solvers/pcg.F90

    rea465d3 re083829  
    5151contains 
    5252 
    53   subroutine prec2Level(prepare,A,CP,sol,rhs,res,CoarseMtx_,Restrict,isFirstIter) 
     53  subroutine prec2Level(prepare,A,CP,sol,rhs,res,isFirstIter) 
    5454    use CoarseAllgathers 
    5555 
    … …  
    6060    real(kind=rk),dimension(:),pointer :: rhs !< right hand side 
    6161    real(kind=rk),dimension(:),pointer :: res !< residual vector 
    62     type(SpMtx) :: CoarseMtx_ !< Coarse matrix 
    63     type(SpMtx) :: Restrict   !< Restriction matrix 
    6462    logical :: isFirstIter 
    6563 
    … …  
    9189      endif 
    9290 
    93       call AllSendCoarseMtx(CP%cdat%LAC,CoarseMtx_,CP%cdat%lg_cfmap,& 
     91      call AllSendCoarseMtx(CP%cdat%LAC,CP%AC,CP%cdat%lg_cfmap,& 
    9492           CP%cdat%ngfc,CP%cdat%nprocs,CP%cdat%send) 
    95       call AllRecvCoarseMtx(CoarseMtx_,CP%cdat%send,add=add) ! Recieve it 
     93      call AllRecvCoarseMtx(CP%AC,CP%cdat%send,add=add) ! Recieve it 
    9694 
    9795    end subroutine prec2Level_exchangeMatrix 
    … …  
    106104        ! after coarse matrix is received allocate coarse vectors 
    107105        if (associated(crhs)) then 
    108           if (size(crhs)/=CoarseMtx_%ncols) then 
     106          if (size(crhs)/=CP%AC%ncols) then 
    109107            deallocate(csol) 
    110108            deallocate(crhs) 
    111             allocate(crhs(CoarseMtx_%ncols)) 
    112             allocate(csol(CoarseMtx_%nrows)) 
     109            allocate(crhs(CP%AC%ncols)) 
     110            allocate(csol(CP%AC%nrows)) 
    113111          endif 
    114112        else 
    115           allocate(crhs(CoarseMtx_%ncols)) 
    116           allocate(csol(CoarseMtx_%nrows)) 
     113          allocate(crhs(CP%AC%ncols)) 
     114          allocate(csol(CP%AC%nrows)) 
    117115        endif 
    118116        ! allocate memory for vector 
    … …  
    132130        if (sctls%verbose>6) write(stream,*) "Restricting into local coarse vector", size(clrhs) 
    133131        ! Send coarse vector 
    134         call SpMtx_Ax(clrhs,Restrict,rhs,dozero=.true.) ! restrict <RA> 
     132        call SpMtx_Ax(clrhs,CP%R,rhs,dozero=.true.) ! restrict <RA> 
    135133        if (CP%cdat_vec%active) then 
    136134          call AllSendCoarseVector(clrhs,CP%cdat_vec%nprocs,CP%cdat_vec%cdisps,& 
    … …  
    146144      if (.not.CP%cdat%active) then ! 1 processor case 
    147145        if (sctls%method>1.and.sctls%method/=5) then ! multiplicative Schwarz 
    148           call SpMtx_Ax(crhs,Restrict,res,dozero=.true.) ! restriction 
     146          call SpMtx_Ax(crhs,CP%R,res,dozero=.true.) ! restriction 
    149147        else 
    150           call SpMtx_Ax(crhs,Restrict,rhs,dozero=.true.) ! restriction 
     148          call SpMtx_Ax(crhs,CP%R,rhs,dozero=.true.) ! restriction 
    151149        endif 
    152150      end if 
    … …  
    167165      if (isFirstIter) then 
    168166        write (stream,*) & 
    169              'factorising coarse matrix of size',CoarseMtx_%nrows, & 
    170              ' and nnz:',CoarseMtx_%nnz 
    171         CoarseMtx_%subsolve_id=0 
     167             'factorising coarse matrix of size',CP%AC%nrows, & 
     168             ' and nnz:',CP%AC%nnz 
     169        CP%AC%subsolve_id=0 
    172170      end if 
    173171 
    174172      ! Coarse solve 
    175       call sparse_singlesolve(CoarseMtx_%subsolve_id,csol,crhs,& 
    176            nfreds=CoarseMtx_%nrows, & 
    177            nnz=CoarseMtx_%nnz,        & 
    178            indi=CoarseMtx_%indi,      & 
    179            indj=CoarseMtx_%indj,      & 
    180            val=CoarseMtx_%val) 
     173      call sparse_singlesolve(CP%AC%subsolve_id,csol,crhs,& 
     174           nfreds=CP%AC%nrows, & 
     175           nnz=CP%AC%nnz,        & 
     176           indi=CP%AC%indi,      & 
     177           indj=CP%AC%indj,      & 
     178           val=CP%AC%val) 
    181179      if (isFirstIter) then 
    182         CoarseMtx_%indi=CoarseMtx_%indi+1 
    183         CoarseMtx_%indj=CoarseMtx_%indj+1 
     180        CP%AC%indi=CP%AC%indi+1 
     181        CP%AC%indj=CP%AC%indj+1 
    184182      end if 
    185183 
    186184      if (CP%cdat_vec%active) then 
    187185        call Vect_remap(csol,clrhs,CP%cdat%gl_cfmap,dozero=.true.) 
    188         call SpMtx_Ax(tmpsol,Restrict,clrhs,dozero=.true.,transp=.true.) 
     186        call SpMtx_Ax(tmpsol,CP%R,clrhs,dozero=.true.,transp=.true.) 
    189187 
    190188      elseif (CP%cdat%active) then 
    191189        call Vect_remap(csol,clrhs,CP%cdat%gl_cfmap,dozero=.true.) 
    192         call SpMtx_Ax(tmpsol,Restrict,clrhs,dozero=.true.,transp=.true.) 
     190        call SpMtx_Ax(tmpsol,CP%R,clrhs,dozero=.true.,transp=.true.) 
    193191      else 
    194         call SpMtx_Ax(tmpsol,Restrict,csol,dozero=.true.,transp=.true.) ! interpolation 
     192        call SpMtx_Ax(tmpsol,CP%R,csol,dozero=.true.,transp=.true.) ! interpolation 
    195193      endif 
    196194 
    … …  
    221219  !------------------------------- 
    222220  subroutine preconditioner(sol,A,rhs,M,finePrec,coarsePrec,& 
    223                A_ghost,CoarseMtx_,Restrict,refactor,bugtrack_) 
     221               A_ghost,refactor,bugtrack_) 
    224222    use CoarseAllgathers 
    225223    use CoarseMtx_mod 
    … …  
    235233                                              !! here for multiplicative Schwarz 
    236234    type(SpMtx),optional               :: A_ghost  !< matr@interf. 
    237     type(SpMtx),optional               :: CoarseMtx_ !< Coarse matrix 
    238     type(SpMtx),optional               :: Restrict   !< Restriction matrix 
    239235    logical,intent(inout),optional :: refactor 
    240236    logical,optional                   :: bugtrack_ 
    … …  
    268264       
    269265    if (sctls%levels>1) then 
    270       call prec2Level(.true.,A,coarsePrec,sol,rhs,res,CoarseMtx_,Restrict,isFirstIter) 
     266      call prec2Level(.true.,A,coarsePrec,sol,rhs,res,isFirstIter) 
    271267    end if 
    272268 
    … …  
    275271 
    276272    if (sctls%levels>1) then 
    277       call prec2Level(.false.,A,coarsePrec,sol,rhs,res,CoarseMtx_,Restrict,isFirstIter) 
     273      call prec2Level(.false.,A,coarsePrec,sol,rhs,res,isFirstIter) 
    278274    end if 
    279275 
    … …  
    300296  !-------------------------- 
    301297  subroutine pcg_weigs (A,b,x,Msh,finePrec,coarsePrec,it,cond_num,A_interf_,tol_,maxit_, & 
    302        x0_,solinf,resvects_,CoarseMtx_,Restrict,refactor_) 
     298       x0_,solinf,resvects_,refactor_) 
    303299    use CoarseAllgathers 
    304300 
    … …  
    322318    type(ConvInf),            intent(in out), optional :: solinf !< Solution statistics 
    323319    logical,                      intent(in), optional :: resvects_ !< Fill in the 'resvect' or not 
    324     type(SpMtx),optional                               :: CoarseMtx_ !< Coarse matrix 
    325     type(SpMtx),optional                               :: Restrict !< Restriction mtx 
    326320    logical,intent(in),optional                        :: refactor_ 
    327321     
    … …  
    429423                            coarsePrec=coarsePrec, & 
    430424                    A_ghost=A_interf_,  & 
    431                    CoarseMtx_=CoarseMtx_, & 
    432                     Restrict=Restrict,    & 
    433425                    refactor=refactor,   & 
    434426                    bugtrack_=bugtrack) 
Note: See TracChangeset for help on using the changeset viewer.