Changeset b77a46f


Ignore:
Timestamp:
12/11/10 12:29:18 (2 years ago)
Author:
Oleg Batrashev <ogbash@…>
Branches:
master, external, fix-prolong, refactor, refactor-subsolvers
Children:
026b2c5
Parents:
981f775
git-author:
Oleg Batrashev <ogbash@…> (12/11/10 12:29:18)
git-committer:
Oleg Batrashev <ogbash@…> (12/11/10 12:29:18)
Message:

Unsuccessful refactoring of coarse preconditioner.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/solvers/pcg.F90

    ra10a43b rb77a46f  
    359359  end subroutine preconditioner_1level 
    360360 
     361  subroutine prec2Level_exchangeMatrix(CoarseMtx_) 
     362    use CoarseAllgathers 
     363    type(SpMtx) :: CoarseMtx_ !< Coarse matrix 
     364    integer :: ol 
     365    logical :: add 
     366 
     367    ol=max(sctls%overlap,sctls%smoothers) 
     368    if (ol==0) then 
     369      add=.false. 
     370    else 
     371      add=.true. 
     372    endif 
     373 
     374    call AllSendCoarseMtx(cdat%LAC,CoarseMtx_,cdat%lg_cfmap,& 
     375         cdat%ngfc,cdat%nprocs,cdat%send) 
     376    call AllRecvCoarseMtx(CoarseMtx_,cdat%send,add=add) ! Recieve it 
     377 
     378  end subroutine prec2Level_exchangeMatrix 
     379 
    361380  !------------------------------- 
    362381  !> Make preconditioner 
    … …  
    385404    logical :: add,bugtrack 
    386405    real(kind=rk) :: t1 
     406    logical :: isFirstIter 
    387407 
    388408    t1 = MPI_WTime() 
    … …  
    394414    endif 
    395415    ! ---------------------------- 
     416    isFirstIter = .not.(present(refactor_).and.refactor_) 
    396417    if (sctls%method==0) then 
    397418      sol=rhs 
    … …  
    405426      if (.not.present(Restrict)) call DOUG_abort("Restriction matrix needs to be passed along with the coarse matrix!") 
    406427      if (cdat%active) then 
    407         ! Not the first iteration, so send the coarse vector here 
    408         if (.not.(present(refactor_).and.refactor_)) then 
    409           call SpMtx_Ax(clrhs,Restrict,rhs,dozero=.true.) ! restrict <RA> 
    410          ! And send 
    411 if (bugtrack)write(stream,*) "local coarse RHS is:",clrhs 
    412           if (cdat_vec%active) then 
    413             call AllSendCoarseVector(clrhs,cdat_vec%nprocs,cdat_vec%cdisps,& 
    414                       cdat_vec%send,useprev=.true.) 
    415           else 
    416             call AllSendCoarseVector(clrhs,cdat%nprocs,cdat%cdisps,& 
    417                       cdat%send,useprev=.true.) 
     428        if (isFirstIter) then 
     429          ! First iteration - send matrix 
     430          call prec2Level_exchangeMatrix(CoarseMtx_) 
     431 
     432          ! factorise coarse matrix 
     433          write (stream,*) & 
     434            'factorising coarse matrix of size',CoarseMtx_%nrows, & 
     435            ' and nnz:',CoarseMtx_%nnz 
     436          call free_spmtx_subsolves(CoarseMtx_) 
     437          allocate(CoarseMtx_%subsolve_ids(1)) 
     438          CoarseMtx_%subsolve_ids=0 
     439          CoarseMtx_%nsubsolves=1 
     440 
     441          ! Factorise the matrix 
     442          if (sctls%verbose>2) then 
     443            write(stream,*)'Global coarse matrix is:---------' 
     444            do i=1,CoarseMtx_%nnz 
     445               write(stream,*) CoarseMtx_%indi(i),& 
     446                     CoarseMtx_%indj(i),CoarseMtx_%val(i) 
     447            enddo 
     448            write(stream,*)'---------------------------------' 
    418449          endif 
    419         else ! First iteration - send matrix 
    420           call AllSendCoarseMtx(cdat%LAC,CoarseMtx_,cdat%lg_cfmap,& 
    421                                 cdat%ngfc,cdat%nprocs,cdat%send) 
     450          call factorise(CoarseMtx_%subsolve_ids(1), & 
     451               nfreds=CoarseMtx_%nrows,   & 
     452               nnz=CoarseMtx_%nnz,        & 
     453               indi=CoarseMtx_%indi,      & 
     454               indj=CoarseMtx_%indj,      & 
     455               val=CoarseMtx_%val) 
     456          CoarseMtx_%indi=CoarseMtx_%indi+1 
     457          CoarseMtx_%indj=CoarseMtx_%indj+1 
     458        end if 
     459 
     460        ! Send coarse vector 
     461        call SpMtx_Ax(clrhs,Restrict,rhs,dozero=.true.) ! restrict <RA> 
     462        if (cdat_vec%active) then 
     463          call AllSendCoarseVector(clrhs,cdat_vec%nprocs,cdat_vec%cdisps,& 
     464               cdat_vec%send,useprev=isFirstIter) 
     465        else 
     466          call AllSendCoarseVector(clrhs,cdat%nprocs,cdat%cdisps,& 
     467               cdat%send,useprev=isFirstIter) 
    422468        endif 
    423       endif 
    424  
    425       ol=max(sctls%overlap,sctls%smoothers) 
    426       if (sctls%input_type==DCTL_INPUT_TYPE_ELEMENTAL.or.numprocs>1) then 
    427           call sparse_multisolve(sol=sol,A=A,M=M,rhs=rhs,res=res, & 
    428                         A_interf_=A_interf_, & 
    429                         refactor=refactor_) !fine solves  
    430       else 
    431           call sparse_multisolve(sol=sol,A=A,M=M,rhs=rhs,res=res, & 
    432                         A_interf_=A_interf_,AC=CoarseMtx_, & 
    433                         refactor=refactor_,Restrict=Restrict) !fine solves  
    434 if (bugtrack)call Print_Glob_Vect(sol,M,'global sol===') 
    435       endif  
    436  
    437       if (sctls%levels>1) then 
    438         if (ol==0) then 
    439           add=.false. 
    440         else 
    441           add=.true. 
    442         endif 
    443         if (present(refactor_).and.refactor_) then !{ setup coarse solve: 
    444           ! We have a matrix send pending 
    445           if (cdat%active) then 
    446             !write(stream,*) "Waiting the Matrix!" 
    447             call AllRecvCoarseMtx(CoarseMtx_,cdat%send,add=add) ! Recieve it 
    448             !write(stream,*) "Got Matrix!" 
    449             allocate(clrhs(Restrict%nrows)) ! allocate memory for vector       
    450             call SpMtx_Ax(clrhs,Restrict,rhs,dozero=.true.) ! restrict <RA> 
    451             ! And send 
    452 if (bugtrack)write(stream,*) "(f) local Coarse RHS is:",clrhs 
    453             if (cdat_vec%active) then 
    454               call AllSendCoarseVector(clrhs,cdat_vec%nprocs,& 
    455                      cdat_vec%cdisps,cdat_vec%send) 
    456             else 
    457               call AllSendCoarseVector(clrhs,cdat%nprocs,& 
    458                      cdat%cdisps,cdat%send) 
    459             endif 
    460           endif 
     469 
     470        ! allocate coarse vectors 
     471        if (isFirstIter) then 
    461472          if (associated(crhs)) then 
    462473            if (size(crhs)/=CoarseMtx_%ncols) then 
    … …  
    474485            allocate(tmpsol(size(rhs))) 
    475486          endif 
    476  
    477           if (.not.cdat%active) then 
     487        end if 
     488 
     489        ! first level prec 
     490        if (sctls%input_type==DCTL_INPUT_TYPE_ELEMENTAL.or.numprocs>1) then 
     491          call sparse_multisolve(sol=sol,A=A,M=M,rhs=rhs,res=res, & 
     492                        A_interf_=A_interf_, & 
     493                        refactor=refactor_) !fine solves  
     494        else 
     495          call sparse_multisolve(sol=sol,A=A,M=M,rhs=rhs,res=res, & 
     496                        A_interf_=A_interf_,AC=CoarseMtx_, & 
     497                        refactor=refactor_,Restrict=Restrict) !fine solves  
     498        endif 
     499      end if ! cdat%active 
     500 
     501      if (sctls%levels>1) then 
     502        ol=max(sctls%overlap,sctls%smoothers) 
     503        if (isFirstIter) then 
     504          if (.not.cdat%active) then ! 1 processor case 
    478505            if (sctls%smoothers==-1) then 
    479506              allocate(tmpsol2(A%nrows)) 
    … …  
    489516            endif 
    490517          endif 
    491           !csol=0.0_rk 
    492           write (stream,*) & 
    493             'factorising coarse matrix of size',CoarseMtx_%nrows, & 
    494             ' and nnz:',CoarseMtx_%nnz 
    495           call free_spmtx_subsolves(CoarseMtx_) 
    496           allocate(CoarseMtx_%subsolve_ids(1)) 
    497           CoarseMtx_%subsolve_ids=0 
    498           CoarseMtx_%nsubsolves=1 
    499  
    500           !! TODO: force it to do it AFTER factorization, but before solve 
    501           !if (present(cdat)) then 
    502           if (cdat%active) then 
    503             ! Factorise the matrix 
    504             if (sctls%verbose>2) then 
    505               write(stream,*)'Global coarse matrix is:---------' 
    506               do i=1,CoarseMtx_%nnz 
    507                 write(stream,*) CoarseMtx_%indi(i),& 
    508                 CoarseMtx_%indj(i),CoarseMtx_%val(i) 
    509               enddo 
    510               write(stream,*)'---------------------------------' 
    511             endif 
    512             call factorise(CoarseMtx_%subsolve_ids(1), & 
    513                  nfreds=CoarseMtx_%nrows,   & 
    514                  nnz=CoarseMtx_%nnz,        & 
    515                  indi=CoarseMtx_%indi,      & 
    516                  indj=CoarseMtx_%indj,      & 
    517                  val=CoarseMtx_%val) 
    518             ! Recieve the vector for solve 
    519             if (cdat_vec%active) then 
    520               call AllRecvCoarseVector(crhs,cdat_vec%nprocs,& 
    521                      cdat_vec%cdisps,cdat_vec%glg_cfmap,cdat_vec%send) 
    522             else 
    523               call AllRecvCoarseVector(crhs,cdat%nprocs,& 
    524                      cdat%cdisps,cdat%glg_cfmap,cdat%send) 
    525             endif 
    526             !write(stream,*) "Got coarse vector!" 
    527             !call MPI_BARRIER(MPI_COMM_WORLD,i) 
     518        end if 
     519 
     520        !csol=0.0_rk 
     521        ! Recieve the vector for solve 
     522        if (cdat_vec%active) then 
     523          call AllRecvCoarseVector(crhs,cdat_vec%nprocs,& 
     524               cdat_vec%cdisps,cdat_vec%glg_cfmap,cdat_vec%send) 
     525        else 
     526          call AllRecvCoarseVector(crhs,cdat%nprocs,& 
     527               cdat%cdisps,cdat%glg_cfmap,cdat%send) 
     528        endif 
     529        !write(stream,*) "Got coarse vector!" 
     530        !call MPI_BARRIER(MPI_COMM_WORLD,i) 
     531 
     532        call sparse_singlesolve(CoarseMtx_%subsolve_ids(1),csol,crhs,nfreds=CoarseMtx_%nrows) 
     533        write (stream,*) 'coarse factorisation done!',CoarseMtx_%subsolve_ids(1) 
     534        if (bugtrack)write(stream,*) "(f) Coarse SOL is:",csol 
     535 
     536        if (cdat_vec%active) then 
     537          call Vect_remap(csol,clrhs,cdat_vec%gl_cfmap,dozero=.true.) 
     538          call SpMtx_Ax(tmpsol,Restrict,clrhs,dozero=.true.,transp=.true.) 
     539        elseif (cdat%active) then 
     540          call Vect_remap(csol,clrhs,cdat%gl_cfmap,dozero=.true.) 
     541          call SpMtx_Ax(tmpsol,Restrict,clrhs,dozero=.true.,transp=.true.) 
     542        else 
     543          if (sctls%smoothers==-1) then 
     544            call SpMtx_Ax(tmpsol2,Restrict,csol,dozero=.true.,transp=.true.) ! interpolation 
     545            tmpsol=0.0_rk 
     546            call exact_sparse_multismoother(tmpsol,A,tmpsol2) 
     547          else 
     548            call SpMtx_Ax(tmpsol,Restrict,csol,dozero=.true.,transp=.true.) ! interpolation 
    528549          endif 
    529 if (bugtrack) then 
    530  if (isnan(sum(crhs))) then 
    531   write(stream,*) "(f) global Coarse RHS is:",crhs 
    532   call DOUG_abort('Got NaN receive fromAllRecvCoarseVector...',838465) 
    533  endif 
    534 !else 
    535 ! write(stream,*) "(f) global Coarse RHS is:",crhs 
    536 endif 
    537  
    538  
    539           call sparse_singlesolve(CoarseMtx_%subsolve_ids(1),csol,crhs, & 
    540                  nfreds=CoarseMtx_%nrows,   & 
    541                  nnz=CoarseMtx_%nnz,        & 
    542                  indi=CoarseMtx_%indi,      & 
    543                  indj=CoarseMtx_%indj,      & 
    544                  val=CoarseMtx_%val) 
    545           write (stream,*) 'coarse factorisation done!',CoarseMtx_%subsolve_ids(1) 
    546 if (bugtrack)write(stream,*) "(f) Coarse SOL is:",csol 
    547  
    548           if (cdat_vec%active) then 
    549             CoarseMtx_%indi=CoarseMtx_%indi+1 
    550             CoarseMtx_%indj=CoarseMtx_%indj+1 
    551             call Vect_remap(csol,clrhs,cdat_vec%gl_cfmap,dozero=.true.) 
    552             call SpMtx_Ax(tmpsol,Restrict,clrhs,dozero=.true.,transp=.true.) 
    553           elseif (cdat%active) then 
    554             CoarseMtx_%indi=CoarseMtx_%indi+1 
    555             CoarseMtx_%indj=CoarseMtx_%indj+1 
    556             call Vect_remap(csol,clrhs,cdat%gl_cfmap,dozero=.true.) 
    557             call SpMtx_Ax(tmpsol,Restrict,clrhs,dozero=.true.,transp=.true.) 
     550        endif 
     551 
     552        if (.not.(cdat_vec%active)) then 
     553          if (sctls%smoothers==-1) then 
     554            tmpsol2=0.0_rk 
     555            call exact_sparse_multismoother(tmpsol2,A,rhs) 
     556            call SpMtx_Ax(crhs,Restrict,tmpsol2,dozero=.true.) ! restriction 
    558557          else 
    559             if (sctls%smoothers==-1) then 
    560               call SpMtx_Ax(tmpsol2,Restrict,csol,dozero=.true.,transp=.true.) ! interpolation 
    561               tmpsol=0.0_rk 
    562               call exact_sparse_multismoother(tmpsol,A,tmpsol2) 
    563             else 
    564               call SpMtx_Ax(tmpsol,Restrict,csol,dozero=.true.,transp=.true.) ! interpolation 
    565             endif 
     558            call SpMtx_Ax(crhs,Restrict,rhs,dozero=.true.) ! restriction 
    566559          endif 
    567 if (bugtrack)call Print_Glob_Vect(tmpsol,M,'tmpsol===',chk_endind=M%ninner) 
    568         else !}{ apply coarse solve: 
    569           if (cdat_vec%active) then 
    570             call AllRecvCoarseVector(crhs,cdat_vec%nprocs,cdat_vec%cdisps,& 
    571                                           cdat_vec%glg_cfmap,cdat_vec%send) 
    572           elseif (cdat%active) then 
    573             call AllRecvCoarseVector(crhs,cdat%nprocs,cdat%cdisps,& 
    574                                           cdat%glg_cfmap,cdat%send) 
    575           else 
    576             if (sctls%smoothers==-1) then 
    577               tmpsol2=0.0_rk 
    578               call exact_sparse_multismoother(tmpsol2,A,rhs) 
    579               call SpMtx_Ax(crhs,Restrict,tmpsol2,dozero=.true.) ! restriction 
    580             else 
    581               call SpMtx_Ax(crhs,Restrict,rhs,dozero=.true.) ! restriction 
    582             endif 
    583           endif 
    584           call sparse_singlesolve(CoarseMtx_%subsolve_ids(1),csol,crhs, & 
    585                  nfreds=CoarseMtx_%nrows) 
    586           if (cdat_vec%active) then 
    587             call Vect_remap(csol,clrhs,cdat_vec%gl_cfmap,dozero=.true.) 
    588             call SpMtx_Ax(tmpsol,Restrict,clrhs,dozero=.true.,transp=.true.) 
    589           elseif (cdat%active) then 
    590             call Vect_remap(csol,clrhs,cdat%gl_cfmap,dozero=.true.) 
    591             call SpMtx_Ax(tmpsol,Restrict,clrhs,dozero=.true.,transp=.true.) 
    592           else 
    593             if (sctls%smoothers==-1) then 
    594               call SpMtx_Ax(tmpsol2,Restrict,csol,dozero=.true.,transp=.true.) ! interpolation 
    595               tmpsol=0.0_rk 
    596               call exact_sparse_multismoother(tmpsol,A,tmpsol2) 
    597             else 
    598               call SpMtx_Ax(tmpsol,Restrict,csol,dozero=.true.,transp=.true.) ! interpolation 
    599 !              call Vect_print(tmpsol,"csol") 
    600             endif 
    601           endif 
    602         endif !} 
    603         if (sctls%method==1) then 
     560        endif 
     561      endif !} 
     562      if (sctls%method==1) then 
    604563!call Print_Glob_Vect(sol,M,'sol===',chk_endind=M%ninner) 
    605           sol(1:A%nrows)=sol(1:A%nrows)+tmpsol(1:A%nrows) 
    606         elseif (sctls%method==3) then ! fully multiplicative Schwarz 
    607           sol(1:A%nrows)=sol(1:A%nrows)+tmpsol(1:A%nrows) 
    608           ! calculate the residual: 
    609           call SpMtx_Ax(res,A,sol,dozero=.true.) !  
    610           res=rhs-res 
    611         elseif (sctls%method==2.and.ol>0) then ! fully multiplicative Schwarz 
    612           sol(1:A%nrows)=tmpsol(1:A%nrows) 
    613           ! calculate the residual: 
    614           call SpMtx_Ax(res,A,sol,dozero=.true.) !  
    615           res=rhs-res 
    616         endif 
     564        sol(1:A%nrows)=sol(1:A%nrows)+tmpsol(1:A%nrows) 
     565      elseif (sctls%method==3) then ! fully multiplicative Schwarz 
     566        sol(1:A%nrows)=sol(1:A%nrows)+tmpsol(1:A%nrows) 
     567        ! calculate the residual: 
     568        call SpMtx_Ax(res,A,sol,dozero=.true.) !  
     569        res=rhs-res 
     570      elseif (sctls%method==2.and.ol>0) then ! fully multiplicative Schwarz 
     571        sol(1:A%nrows)=tmpsol(1:A%nrows) 
     572        ! calculate the residual: 
     573        call SpMtx_Ax(res,A,sol,dozero=.true.) !  
     574        res=rhs-res 
    617575      endif 
    618576      if (sctls%method>1) then ! For multiplicative Schwarz method...: 
Note: See TracChangeset for help on using the changeset viewer.