Changeset 5aaa704


Ignore:
Timestamp:
08/01/06 17:52:40 (7 years ago)
Author:
eero <eero@…>
Branches:
master, external, fix-prolong, refactor, refactor-ext, refactor-subsolvers
Children:
821f89d
Parents:
98817b2
git-author:
eero <eero@…> (08/01/06 17:52:40)
git-committer:
eero <eero@…> (08/01/06 17:52:40)
Message:

Multiplicative Schwarz implementation (in the case of 1 processor) -- test of concept, still need to optimise residual calculations for test of better performance than the additive Schwarz. Best so far to choose method=2 in case of smoothers==0 and method=4 in case of smoothers>0.

git-svn-id:  svn://dougdevel.org/doug/trunk@50 146eac83-2b0f-0410-9363-9ebcc3d868b9

File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/solvers/pcg.f90

    r792ba1b r5aaa704  
    217217  end subroutine pcg 
    218218 
    219  
    220   !------------------------------- 
    221   ! Setup preconditioner 
    222   !------------------------------- 
    223   subroutine setup_preconditioner(sol,A,rhs,CoarseMtx_, Restrict) 
    224     use subsolvers 
    225     use CoarseMtx_mod 
    226     implicit none 
    227     real(kind=rk),dimension(:),pointer :: sol 
    228     type(SpMtx),           intent(in)  :: A 
    229     real(kind=rk),dimension(:),pointer :: rhs 
    230     type(SpMtx),optional               :: CoarseMtx_ ! Coarse matrix 
    231     type(SpMtx),optional               :: Restrict ! Restriction matrix 
    232  
    233     ! ----- local: ------ 
    234     real(kind=rk),dimension(:),pointer :: csol,crhs,tmpsol 
    235     ! ---------------------------- 
    236     sol=0.0_rk 
    237     if (present(CoarseMtx_)) then 
    238       call sparse_multisolve(sol,A,rhs,CoarseMtx_) !fine solves 
    239       !call sparse_multisolve(sol,A,rhs) !fine solves 
    240       if (sctls%levels>1) then 
    241         if (coarseid==0) then ! setup coarse solve: 
    242           allocate(crhs(CoarseMtx_%ncols)) 
    243           allocate(csol(CoarseMtx_%nrows)) 
    244           allocate(tmpsol(A%nrows)) 
    245           call SpMtx_Ax(crhs,Restrict,rhs,dozero=.true.) ! restriction 
    246           !csol=0.0_rk 
    247           write (stream,*) & 
    248             'factorising coarse matrix of size',CoarseMtx_%nrows, & 
    249             ' and nnz:',CoarseMtx_%nnz 
    250           call sparse_singlesolve(coarseid,csol,crhs, & 
    251                  nfreds=CoarseMtx_%nrows,   & 
    252                  nnz=CoarseMtx_%nnz,        & 
    253                  indi=CoarseMtx_%indi,      & 
    254                  indj=CoarseMtx_%indj,      & 
    255                  val=CoarseMtx_%val) 
    256           write (stream,*) 'coarse factorisation done!' 
    257           !tmpsol=0.0_rk 
    258           call SpMtx_Ax(tmpsol,Restrict,csol,dozero=.true.,transp=.true.) ! interpolation 
    259           !call SpMtx_Ax(tmpsol,Interp,csol,dozero=.true.,transp=.false.) ! interpolation 
    260         else ! apply coarse solve: 
    261           call SpMtx_Ax(crhs,Restrict,rhs,dozero=.true.) ! restriction 
    262           !csol=0.0_rk 
    263           call sparse_singlesolve(coarseid,csol,crhs, & 
    264                  nfreds=CoarseMtx_%nrows) 
    265           !tmpsol=0.0_rk 
    266           call SpMtx_Ax(tmpsol,Restrict,csol,dozero=.true.,transp=.true.) ! interpolation 
    267           !call SpMtx_Ax(tmpsol,Interp,csol,dozero=.true.,transp=.false.) ! interpolation 
    268         endif 
    269         !sol=sol+tmpsol 
    270         sol(1:A%nrows)=sol(1:A%nrows)+tmpsol(1:A%nrows) 
    271       endif 
    272     else 
    273       call sparse_multisolve(sol,A,rhs) 
    274     endif 
    275   end subroutine setup_preconditioner 
    276  
    277219  !------------------------------- 
    278220  ! Make preconditioner 
    279221  !------------------------------- 
    280   subroutine preconditioner(sol,A,rhs,A_interf_,CoarseMtx_,Restrict,refactor_,cdat) 
     222  subroutine preconditioner(sol,A,rhs,M,A_interf_,CoarseMtx_,Restrict,refactor_,cdat) 
    281223    use subsolvers 
    282224    use CoarseMtx_mod 
    … …  
    287229    type(SpMtx)                        :: A 
    288230    real(kind=rk),dimension(:),pointer :: rhs 
     231    type(Mesh),intent(in)              :: M ! Mesh 
     232    real(kind=rk),dimension(:),pointer :: res ! -- residual vector, allocated 
     233                                              ! here for multiplicative Schwarz 
    289234    type(SpMtx),optional               :: A_interf_  ! matr@interf. 
    290235    type(SpMtx),optional               :: CoarseMtx_ ! Coarse matrix 
    291236    type(SpMtx),optional               :: Restrict ! Restriction matrix 
    292     logical,intent(in),optional :: refactor_ 
     237    logical,intent(inout),optional :: refactor_ 
    293238    type(CoarseData),intent(inout),optional :: cdat 
    294239    ! ----- local: ------ 
    … …  
    307252    if (present(CoarseMtx_)) then !{ 
    308253 
    309         if (.not.present(Restrict)) call DOUG_abort("Restriction matrix needs to be passed along with the coarse matrix!") 
     254      if (.not.present(Restrict)) call DOUG_abort("Restriction matrix needs to be passed along with the coarse matrix!") 
    310255         
    311     if (present(cdat)) then 
    312  
    313 !do i=1,MG%nlf 
    314 !   if (MG%inner_interf_fmask(i)/=0) &  
    315 !   write(stream,*) MG%lg_fmap(i), rhs(i) 
    316 !enddo 
    317  
    318  
     256      if (present(cdat)) then 
    319257        ! Not the first iteration, so send the coarse vector here 
    320258        if (.not.(present(refactor_).and.refactor_)) then 
    321          call SpMtx_Ax(clrhs,Restrict,rhs,dozero=.true.) ! restrict <RA> 
     259          call SpMtx_Ax(clrhs,Restrict,rhs,dozero=.true.) ! restrict <RA> 
    322260!         call Vect_print(rhs,"RHS") 
    323  
     261       
    324262         ! And send 
    325          call AllSendCoarseVector(clrhs,cdat%nprocs,cdat%cdisps,cdat%send, & 
     263          call AllSendCoarseVector(clrhs,cdat%nprocs,cdat%cdisps,cdat%send, & 
    326264                                  useprev=.true.) 
    327265        else ! First iteration - send matrix 
    328          call AllSendCoarseMtx(cdat%LAC,CoarseMtx_,cdat%lg_cfmap,& 
     266          call AllSendCoarseMtx(cdat%LAC,CoarseMtx_,cdat%lg_cfmap,& 
    329267                                cdat%ngfc,cdat%nprocs,cdat%send) 
    330          allocate(ctmp(cdat%ngfc)) 
    331  
     268          allocate(ctmp(cdat%ngfc)) 
     269       
    332270        endif 
    333271      endif 
    334272 
     273      if (sctls%method>1) then ! For multiplicative Schwarz method...: 
     274        allocate(res(size(rhs))) 
     275      endif 
    335276      if (sctls%input_type==DCTL_INPUT_TYPE_ELEMENTAL) then 
    336           call sparse_multisolve(sol=sol,A=a,rhs=rhs, & 
     277          call sparse_multisolve(sol=sol,A=A,M=M,rhs=rhs,res=res, & 
    337278                        A_interf_=A_interf_, & 
    338279                        refactor=refactor_) !fine solves  
    339280      else 
    340           call sparse_multisolve(sol=sol,A=a,rhs=rhs, & 
     281          call sparse_multisolve(sol=sol,A=A,M=M,rhs=rhs,res=res, & 
    341282                        A_interf_=A_interf_,AC=CoarseMtx_, & 
    342283                        refactor=refactor_,Restrict=Restrict) !fine solves  
     284!call Print_Glob_Vect(sol,M,'global sol===') 
    343285      endif  
    344286 
    … …  
    377319 
    378320          if (.not.present(cdat)) then 
    379           if (sctls%smoothers==-1) then 
    380             allocate(tmpsol2(A%nrows)) 
    381             tmpsol2=0.0_rk 
    382             call exact_sparse_multismoother(tmpsol2,A,rhs) 
    383             call SpMtx_Ax(crhs,Restrict,tmpsol2,dozero=.true.) ! restriction 
    384           else 
    385             call SpMtx_Ax(crhs,Restrict,rhs,dozero=.true.) ! restriction 
    386           endif 
     321            if (sctls%smoothers==-1) then 
     322              allocate(tmpsol2(A%nrows)) 
     323              tmpsol2=0.0_rk 
     324              call exact_sparse_multismoother(tmpsol2,A,rhs) 
     325              call SpMtx_Ax(crhs,Restrict,tmpsol2,dozero=.true.) ! restriction 
     326            else 
     327              if (sctls%method>1.and.sctls%method/=5) then ! multiplicative Schwarz 
     328                call SpMtx_Ax(crhs,Restrict,res,dozero=.true.) ! restriction 
     329              else 
     330                call SpMtx_Ax(crhs,Restrict,rhs,dozero=.true.) ! restriction 
     331              endif 
     332            endif 
    387333          endif 
    388334          !csol=0.0_rk 
    … …  
    412358          endif 
    413359 
    414           !call sparse_singlesolve(coarseid,csol,crhs, & 
    415360          call sparse_singlesolve(CoarseMtx_%subsolve_ids(1),csol,crhs, & 
    416361                 nfreds=CoarseMtx_%nrows,   & 
    … …  
    433378            call SpMtx_Ax(tmpsol,Restrict,clrhs,dozero=.true.,transp=.true.) ! RB 
    434379          else 
    435           if (sctls%smoothers==-1) then 
    436             call SpMtx_Ax(tmpsol2,Restrict,csol,dozero=.true.,transp=.true.) ! interpolation 
    437             tmpsol=0.0_rk 
    438             call exact_sparse_multismoother(tmpsol,A,tmpsol2) 
    439           else 
    440             call SpMtx_Ax(tmpsol,Restrict,csol,dozero=.true.,transp=.true.) ! interpolation 
    441           endif 
     380            if (sctls%smoothers==-1) then 
     381              call SpMtx_Ax(tmpsol2,Restrict,csol,dozero=.true.,transp=.true.) ! interpolation 
     382              tmpsol=0.0_rk 
     383              call exact_sparse_multismoother(tmpsol,A,tmpsol2) 
     384            else 
     385              call SpMtx_Ax(tmpsol,Restrict,csol,dozero=.true.,transp=.true.) ! interpolation 
     386            endif 
    442387          endif 
    443388        else ! apply coarse solve: 
    … …  
    450395 
    451396          else 
    452           if (sctls%smoothers==-1) then 
    453             tmpsol2=0.0_rk 
    454             call exact_sparse_multismoother(tmpsol2,A,rhs) 
    455             call SpMtx_Ax(crhs,Restrict,tmpsol2,dozero=.true.) ! restriction 
    456           else 
    457             call SpMtx_Ax(crhs,Restrict,rhs,dozero=.true.) ! restriction 
    458           endif 
     397            if (sctls%smoothers==-1) then 
     398              tmpsol2=0.0_rk 
     399              call exact_sparse_multismoother(tmpsol2,A,rhs) 
     400              call SpMtx_Ax(crhs,Restrict,tmpsol2,dozero=.true.) ! restriction 
     401            else 
     402              call SpMtx_Ax(crhs,Restrict,rhs,dozero=.true.) ! restriction 
     403            endif 
    459404          endif 
    460405          call sparse_singlesolve(CoarseMtx_%subsolve_ids(1),csol,crhs, & 
    … …  
    470415            call SpMtx_Ax(tmpsol,Restrict,clrhs,dozero=.true.,transp=.true.) ! RB 
    471416          else 
    472           if (sctls%smoothers==-1) then 
    473             call SpMtx_Ax(tmpsol2,Restrict,csol,dozero=.true.,transp=.true.) ! interpolation 
    474             tmpsol=0.0_rk 
    475             call exact_sparse_multismoother(tmpsol,A,tmpsol2) 
    476           else 
    477             call SpMtx_Ax(tmpsol,Restrict,csol,dozero=.true.,transp=.true.) ! interpolation 
    478 !            call Vect_print(tmpsol,"csol") 
    479           endif 
     417            if (sctls%smoothers==-1) then 
     418              call SpMtx_Ax(tmpsol2,Restrict,csol,dozero=.true.,transp=.true.) ! interpolation 
     419              tmpsol=0.0_rk 
     420              call exact_sparse_multismoother(tmpsol,A,tmpsol2) 
     421            else 
     422              call SpMtx_Ax(tmpsol,Restrict,csol,dozero=.true.,transp=.true.) ! interpolation 
     423!              call Vect_print(tmpsol,"csol") 
     424            endif 
    480425          endif 
    481426        endif 
     427        if (sctls%method==1) then 
     428          sol(1:A%nrows)=sol(1:A%nrows)+tmpsol(1:A%nrows) 
     429        elseif (sctls%method==3) then ! fully multiplicative Schwarz 
     430          sol(1:A%nrows)=sol(1:A%nrows)+tmpsol(1:A%nrows) 
     431          ! calculate the residual: 
     432          call SpMtx_Ax(res,A,sol,dozero=.true.) !  
     433          res=rhs-res 
     434        elseif (sctls%method==4) then ! fully multiplicative Schwarz 
     435          sol(1:A%nrows)=tmpsol(1:A%nrows) 
     436          ! calculate the residual: 
     437          call SpMtx_Ax(res,A,sol,dozero=.true.) !  
     438          res=rhs-res 
     439        endif 
     440      endif 
     441      if (sctls%method>1) then ! For multiplicative Schwarz method...: 
     442        refactor_=.false. 
     443        if (sctls%input_type==DCTL_INPUT_TYPE_ELEMENTAL) then 
     444            call sparse_multisolve(sol=sol,A=A,M=M,rhs=rhs,res=res, & 
     445                          A_interf_=A_interf_, & 
     446                          refactor=refactor_) !fine solves  
     447        else 
     448            call sparse_multisolve(sol=sol,A=A,M=M,rhs=rhs,res=res, & 
     449                          A_interf_=A_interf_,AC=CoarseMtx_, & 
     450                          refactor=refactor_,Restrict=Restrict) !fine solves  
     451        endif  
     452      endif 
     453      if ((sctls%method==2.or.sctls%method==5).and.sctls%levels>1) then  
     454        ! multiplicative on fine level, additive with coarse level:  
    482455        sol(1:A%nrows)=sol(1:A%nrows)+tmpsol(1:A%nrows) 
    483456      endif 
    484457    else !}{ 
    485       if (refactor_.and.present(A_interf_)) then 
     458      if (refactor_.and.present(A_interf_)) then!{ 
    486459        if (sctls%verbose>9) then 
    487           call SpMtx_printMat(A) 
     460          !call SpMtx_printMat(A) 
    488461          call SpMtx_printRaw(A) 
    489           call SpMtx_printMat(A_interf_) 
     462          !call SpMtx_printMat(A_interf_) 
    490463          call SpMtx_printRaw(A_interf_) 
    491464          !stop 
    … …  
    505478          !endif 
    506479          call SpMtx_arrange(A,D_SpMtx_ARRNG_ROWS,sort=.true.) 
    507    write(stream,*)'AAAAAAAAAAAAAAAA is:' 
    508    call SpMtx_printRaw(A) 
    509    call SpMtx_printRaw(A_interf_) 
    510    write(stream,*)'AAAAAAAAAAAAAAAA :' 
    511    call flush(stream) 
    512           call sparse_multisolve(sol=sol,A=A,rhs=rhs, & 
     480          call sparse_multisolve(sol=sol,A=A,M=M,rhs=rhs,res=res, & 
    513481                                 A_interf_=A_interf_, & 
    514482                                  refactor=refactor_) !fine solves  
     483          if (sctls%method>1) then ! For multiplicative Schwarz method...: 
     484            refactor_=.false. 
     485            call sparse_multisolve(sol=sol,A=A,M=M,rhs=rhs,res=res, & 
     486                                   A_interf_=A_interf_, & 
     487                                    refactor=refactor_) !fine solves  
     488          endif 
    515489          ! put the original structure and orders back: 
    516490          A%indi=A_tmp%indi 
    … …  
    523497          call SpMtx_Destroy(A_tmp) 
    524498        else 
    525           call sparse_multisolve(sol=sol,A=a,rhs=rhs, & 
     499          call sparse_multisolve(sol=sol,A=A,M=M,rhs=rhs,res=res, & 
     500                               A_interf_=A_interf_, & 
     501                                refactor=refactor_) !fine solves  
     502          if (sctls%method>1) then ! For multiplicative Schwarz method...: 
     503            refactor_=.false. 
     504            call sparse_multisolve(sol=sol,A=A,M=M,rhs=rhs,res=res, & 
     505                                 A_interf_=A_interf_, & 
     506                                  refactor=refactor_) !fine solves  
     507          endif 
     508        endif 
     509      elseif (refactor_.and.sctls%input_type==DCTL_INPUT_TYPE_ASSEMBLED) then!}{ 
     510        A_tmp=SpMtx_newInit(nnz=A%nnz,nblocks=A%nblocks,nrows=A%nrows,& 
     511                           ncols=A%ncols,& 
     512                           indi=A%indi,indj=A%indj,val=A%val,& 
     513                           arrange_type=A%arrange_type,M_bound=A%M_bound) 
     514        A%arrange_type=D_SpMtx_ARRNG_NO 
     515        call SpMtx_arrange(A,D_SpMtx_ARRNG_ROWS,sort=.true.) 
     516        call sparse_multisolve(sol=sol,A=A,M=M,rhs=rhs,res=res, & 
     517                                refactor=refactor_) !fine solves  
     518        if (sctls%method>1) then ! For multiplicative Schwarz method...: 
     519          refactor_=.false. 
     520          call sparse_multisolve(sol=sol,A=A,M=M,rhs=rhs,res=res, & 
     521                                  refactor=refactor_) !fine solves  
     522        endif 
     523        ! put the original structure and orders back: 
     524        A%indi=A_tmp%indi 
     525        A%indj=A_tmp%indj 
     526        A%val=A_tmp%val 
     527        A%M_bound=A_tmp%M_bound 
     528        A%nrows=A_tmp%nrows 
     529        A%ncols=A_tmp%ncols 
     530        A%arrange_type=A_tmp%arrange_type 
     531        call SpMtx_Destroy(A_tmp) 
     532      else!}{ 
     533        call sparse_multisolve(sol=sol,A=A,M=M,rhs=rhs,res=res, & 
     534                             A_interf_=A_interf_, & 
     535                              refactor=refactor_) !fine solves  
     536        if (sctls%method>1) then ! For multiplicative Schwarz method...: 
     537          refactor_=.false. 
     538          call sparse_multisolve(sol=sol,A=A,M=M,rhs=rhs,res=res, & 
    526539                               A_interf_=A_interf_, & 
    527540                                refactor=refactor_) !fine solves  
    528541        endif 
    529       elseif (refactor_.and.sctls%input_type==DCTL_INPUT_TYPE_ASSEMBLED) then 
    530           A_tmp=SpMtx_newInit(nnz=A%nnz,nblocks=A%nblocks,nrows=A%nrows,& 
    531                              ncols=A%ncols,& 
    532                              indi=A%indi,indj=A%indj,val=A%val,& 
    533                              arrange_type=A%arrange_type,M_bound=A%M_bound) 
    534           A%arrange_type=D_SpMtx_ARRNG_NO 
    535           call SpMtx_arrange(A,D_SpMtx_ARRNG_ROWS,sort=.true.) 
    536  write(stream,*)'BBBBBBBBBBBBAAAA is:' 
    537  write(stream,*)'A%nrows=',A%nrows 
    538  call SpMtx_printRaw(A,startnz=1,endnz=A%mtx_bbe(2,2)) 
    539  write(stream,*)'BBBBBBBBBBBBAAAA :' 
    540  call flush(stream) 
    541           call sparse_multisolve(sol=sol,A=A,rhs=rhs, & 
    542                                   refactor=refactor_) !fine solves  
    543           ! put the original structure and orders back: 
    544           A%indi=A_tmp%indi 
    545           A%indj=A_tmp%indj 
    546           A%val=A_tmp%val 
    547           A%M_bound=A_tmp%M_bound 
    548           A%nrows=A_tmp%nrows 
    549           A%ncols=A_tmp%ncols 
    550           A%arrange_type=A_tmp%arrange_type 
    551           call SpMtx_Destroy(A_tmp) 
    552       else 
    553         call sparse_multisolve(sol=sol,A=a,rhs=rhs, & 
    554                              A_interf_=A_interf_, & 
    555                               refactor=refactor_) !fine solves  
    556       endif 
     542      endif!} 
    557543    endif !} 
     544    if (sctls%method>1) then ! For multiplicative Schwarz method...: 
     545      if (associated(res)) deallocate(res) 
     546    endif 
    558547  end subroutine preconditioner 
    559548 
    … …  
    690679  
    691680!call Print_Glob_Vect(r,Msh,'global r===') 
     681!write(stream,*)'present(CoarseMtx_):',present(CoarseMtx_) 
    692682      call preconditioner(sol=z,          & 
    693683                            A=A,          & 
    694684                          rhs=r,          & 
     685                            M=Msh,        & 
    695686                    A_interf_=A_interf_,  & 
    696687                   CoarseMtx_=CoarseMtx_, & 
    … …  
    702693        call Add_common_interf(z,A,Msh) 
    703694      endif 
    704 !write(stream,*)'localz==:',z 
     695!call Print_Glob_Vect(z,Msh,'global z===') 
    705696!call Print_Glob_Vect(z,Msh,'global z===') 
    706697      ! compute current rho 
    … …  
    724715      ! check 
    725716      res_norm = Vect_dot_product(r,r) 
    726       write (stream,*) "Norm is", res_norm 
     717      !write (stream,*) "Norm is", res_norm 
    727718      ratio_norm = res_norm / init_norm 
    728719      if (ismaster()) & 
    … …  
    764755      if (alpha(i)==0.or.beta(i)<0) then 
    765756        print *,'Unable to compute eigenvalue number',i,' and ' 
    766         do j=1,it 
     757        do j=1,i 
    767758          print *,j,' alpha:',alpha(j),' beta:',beta 
    768759        enddo 
Note: See TracChangeset for help on using the changeset viewer.