DOUG 0.2

pcg_robust.f90

Go to the documentation of this file.
00001 
00003 
00006 module pcgRobust_mod
00007   use RealKind
00008   use RobustCoarseMtx_mod
00009   use globals
00010   use Vect_mod
00011   use SpMtx_op_block
00012   use SpMtx_operation
00013 
00014   implicit none
00015 
00018   type RobustPreconditionMtx
00019      type(SpMtx), pointer :: I(:)
00020      type(SpMtx), pointer :: H(:)
00021      type(SpMtx), pointer :: AI(:)
00022 
00023      integer, pointer :: subsolve_ids(:)
00024   end type RobustPreconditionMtx
00025 
00026 contains
00027 
00028   function RobustPreconditionMtx_new() result (C)
00029     type(RobustPreconditionMtx) :: C
00030 
00031     C%I => NULL()
00032     C%H => NULL()
00033     C%AI => NULL()
00034     C%subsolve_ids => NULL()
00035   end function RobustPreconditionMtx_new
00036 
00038   subroutine pcg_forRCS (A, b, x)
00039     implicit none
00040 
00041     type(SumOfInversedSubMtx),intent(inout) :: A !< System matrix
00042     real(kind=rk),dimension(:),pointer :: b !< RHS
00043     real(kind=rk), pointer :: x(:) !< Solution
00044 
00045     real(kind=rk) :: tol   ! Tolerance
00046     integer       :: maxit ! Max number of iterations
00047 
00048     ! @todo: check which are unnecesary
00049     real(kind=rk),dimension(:),pointer :: r, p, q, m, z, b_k
00050     real(kind=rk) :: rho_curr, rho_prev, init_norm, res_norm
00051     real(kind=rk) :: alpha, beta, tmp, ratio_norm
00052     integer :: iter,nfreds
00053     integer :: ierr
00054 
00055     type(RobustPreconditionMtx) :: C
00056 
00057     C = RobustPreconditionMtx_new()
00058 
00059     write(stream,'(/a)') 'PCG for the sum of inversed submatrices:'
00060 
00061     tol = sctls%solve_tolerance
00062     maxit = sctls%solve_maxiters
00063 
00064     x = 0.0
00065 
00066     nfreds=size(b)
00067     allocate(r(nfreds))
00068     allocate(p(nfreds))
00069     allocate(q(nfreds))
00070     allocate(m(nfreds))
00071     allocate(z(nfreds))
00072     call SOISMtx_pmvm(r,A,x)
00073 
00074     r = b - r
00075     if (sctls%verbose>1) then
00076        write(stream,*) 'initial r = ', r
00077     end if
00078 
00079     init_norm = Vect_dot_product(r,r)
00080     if (init_norm == 0.0) init_norm = 1.0
00081 
00082     write(stream,'(a,i6,a,e8.2)') 'maxit = ',maxit,', tol = ',tol
00083     write(stream,'(a,e22.15)') 'init_norm = ', dsqrt(init_norm)
00084 
00085     iter = 0
00086     ratio_norm = 1.0_rk
00087     !x = 0.0
00088 
00089     do while((ratio_norm > tol*tol).and.(iter <= maxit))
00090        iter = iter + 1
00091 
00092        ! no preconditioner for now
00093        !z = r
00094 
00095        call precondition_forRCS(z, A, C, r)
00096        if (sctls%verbose > 10) then
00097           print *, "z="
00098           print "(10E10.2)", z
00099        end if
00100 
00101        ! compute current rho
00102        rho_curr = Vect_dot_product(r,z)
00103        if (iter == 1) then
00104           p = z
00105        else
00106           beta = rho_curr / rho_prev
00107           p = z + beta * p
00108        end if
00109        call SOISMtx_pmvm(q,A,p)
00110        ! compute alpha
00111        alpha = rho_curr / Vect_dot_product(p,q)
00112        x = x + alpha * p
00113        r = r - alpha * q
00114        rho_prev = rho_curr
00115        ! check
00116        res_norm = Vect_dot_product(r,r)
00117        ratio_norm = dsqrt(res_norm) / dsqrt(init_norm)
00118 !!$       ratio_norm = res_norm / init_norm
00119        if (ismaster()) write(stream, '(i5,a,e22.15)') iter,': res_norm=',dsqrt(res_norm)
00120     end do
00121     
00122   end subroutine pcg_forRCS
00123 
00125   subroutine precondition_forRCS(y, A, C, x)
00126     real(kind=rk), intent(out) :: y(:)
00127     real(kind=rk), pointer :: x(:)
00128     type(SumOfInversedSubMtx), intent(inout) :: A
00129     type(RobustPreconditionMtx), intent(inout) :: C
00130 
00131     integer :: nAggr, iAggr, m, n
00132     real(kind=rk), pointer :: yl(:), xl(:), xh(:), yh(:), yTemp(:)
00133 
00134     if (.not.associated(C%I)) call initialize(C, A)
00135 
00136     y = 0._rk
00137     allocate(yTemp(size(y)))
00138 
00139     nAggr = size(A%Ai)
00140     allocate(xl(size(x)), yl(size(y))) ! @todo: shrink
00141     m = maxval(C%H%nrows)
00142     n = maxval(C%H%ncols)
00143     allocate(xh(m), yh(n))
00144 
00145     do iAggr=1,nAggr
00146        ! Restrict
00147        call SpMtx_Ax(xl, A%R(iAggr), x, dozero=.TRUE.)
00148 
00149        if (sctls%verbose > 9) then
00150           print *, "xl="
00151           print "(10E10.2)", xl
00152        end if
00153        ! Apply B_j^(-1)
00154        call SpMtx_Ax(yl, A%Ai(iAggr), xl, dozero=.TRUE.)
00155 
00156        call SpMtx_Ax(xh, C%AI(iAggr), xl, transp=.TRUE., dozero=.TRUE.)
00157        if (sctls%verbose > 9) then
00158           print *, "xh="
00159           print "(10E10.2)", xh
00160        end if
00161 
00162        yh = 0.
00163        call sparse_singlesolve(C%subsolve_ids(iAggr), yh, xh, C%H(iAggr)%ncols, &
00164             C%H(iAggr)%nnz, C%H(iAggr)%indi, C%H(iAggr)%indj, C%H(iAggr)%val)
00165 
00166        if (sctls%verbose > 9) then
00167           print *, "yh="
00168           print "(10E10.2)", yh
00169        end if
00170 
00171        call SpMtx_Ax(xl, C%AI(iAggr), yh, dozero=.TRUE.)
00172 
00173        if (sctls%verbose > 9) then
00174           print *, "xl2="
00175           print "(10E10.2)", xl
00176        end if
00177 
00178        yl = yl - xl
00179 
00180        ! extend
00181        call SpMtx_Ax(yTemp, A%R(iAggr), yl, transp=.TRUE., dozero=.TRUE.)
00182 
00183        if (sctls%verbose > 9) then
00184           print *, "yTemp="
00185           print "(10E10.2)", yTemp
00186        end if
00187        y = y + yTemp
00188     end do
00189     
00190 
00191     deallocate(xl, yl, xh, yh, yTemp)
00192 
00193   contains
00194     subroutine initialize(C, A)
00195       type(RobustPreconditionMtx), intent(inout) :: C
00196       type(SumOfInversedSubMtx), intent(inout) :: A
00197       
00198       integer :: iAgr, nAgr, iOtherAgr
00199       type(SpMtx) :: Itemp, Htemp, Akl
00200 
00201       if (sctls%verbose>0) then
00202          write (stream, *) "Initializing C preconditioner"
00203       end if
00204 
00205       nAgr = size(A%Ai)
00206 
00207       allocate(C%subsolve_ids(nAgr))
00208       C%subsolve_ids = 0
00209 
00210       allocate(C%I(nAgr), C%H(nAgr), C%AI(nAgr))
00211       
00212       ! find Is (intersection of coarse supports)
00213       do iAgr=1,nAgr
00214          C%I(iAgr)=SpMtx_NewInit(0)
00215          do iOtherAgr=1,nAgr
00216             if (iAgr==iOtherAgr) cycle
00217 
00218             Itemp = SpMtx_AB2(A%R(iAgr), A%R(iOtherAgr), BT=.TRUE.)
00219             call SpMtx_addBlock(C%I(iAgr), Itemp, D_ADDBLOCK_OPERATION_COLS)
00220          end do
00221          if (sctls%verbose > 5) then
00222             print *, "I for", iAgr, ": nrows, ncols=", C%I(iAgr)%nrows, C%I(iAgr)%ncols
00223             call SpMtx_printRaw(C%I(iAgr))
00224          end if
00225       end do
00226       
00227       ! find Hs
00228       do iAgr=1,nAgr
00229          ! find A_{kl}
00230          Akl = SpMtx_newInit(0)
00231          do iOtherAgr=1,nAgr
00232             if (iAgr==iOtherAgr) cycle
00233 
00234             call SpMtx_addBlock(Akl, A%Ai(iOtherAgr), D_ADDBLOCK_OPERATION_DIAG)
00235          end do
00236          if (sctls%verbose > 9) then
00237             write (stream, *) "Akl for", iAgr, ": ", Akl%nrows, Akl%ncols
00238             call SpMtx_printRaw(Akl)
00239          end if
00240 
00241          C%AI(iAgr) = SpMtx_AB2(A%Ai(iAgr), C%I(iAgr))
00242          if (sctls%verbose > 7) then
00243             write (stream, *) "AI for", iAgr, ": ", C%AI(iAgr)%nrows, C%AI(iAgr)%ncols
00244             call SpMtx_printRaw(C%AI(iAgr))
00245          end if
00246 
00247          !print *, "mult2", iAgr
00248          Htemp = SpMtx_AB2(C%I(iAgr), C%AI(iAgr), AT=.TRUE.)
00249          if (sctls%verbose > 9) then
00250             write (stream, *) "Htemp for", iAgr, ": ", Htemp%nrows, Htemp%ncols
00251             call SpMtx_printRaw(Htemp)
00252          end if
00253 
00254          C%H(iAgr) = SpMtx_add(Akl, Htemp, 1.0_rk, 1.0_rk)
00255 
00256          !print *, "Akl for", iAgr, ": ", Akl%nrows, Akl%ncols
00257          !call SpMtx_printRaw(Akl)
00258          if (sctls%verbose > 5) then
00259             write (stream, *) "Hkl for", iAgr, ": ", C%H(iAgr)%nrows, C%H(iAgr)%ncols
00260             call SpMtx_printRaw(C%H(iAgr))
00261          end if
00262 
00263          call SpMtx_destroy(Akl)
00264          call SpMtx_destroy(Htemp)
00265 
00266          !call SpMtx_printRaw(C%H(iAgr))
00267       end do
00268 
00269       if (sctls%verbose>0) then
00270          write (stream, *) "Finished initializing C preconditioner"
00271       end if
00272     end subroutine initialize
00273   end subroutine precondition_forRCS
00274 
00275 end module pcgRobust_mod