Index: src/solvers/pcg.F90 =================================================================== --- src/solvers/pcg.F90 (revision b378423872c8daefa84453552e956c20d5b92073) +++ src/solvers/pcg.F90 (revision 08dfb5f36c9e00fbfadeb13ac323151abe3d4871) @@ -41,5 +41,4 @@ #endif - integer,save :: coarseid=0 private :: & preconditioner, & @@ -50,197 +49,4 @@ contains - - !-------------------------------------------- - ! Preconditioned conjugate gradient method - !-------------------------------------------- - subroutine pcg (A, b, x, Msh, tol_, maxit_, & - x0_, solinf, resvects_, CoarseMtx_) - implicit none - - type(SpMtx),intent(in out) :: A ! System matrix (sparse) - float(kind=rk),dimension(:),pointer :: b ! RHS - float(kind=rk),dimension(:),pointer :: x ! Solution - ! Mesh - aux data for Ax operation - type(Mesh), intent(in) :: Msh - - ! Optional arguments - ! Tolerance of the method - real(kind=rk), intent(in), optional :: tol_ - ! Max number of iterations - integer, intent(in), optional :: maxit_ - ! Initial guess - float(kind=rk), dimension(:), intent(in), optional :: x0_ - ! Solution statistics - type(ConvInf), intent(in out), optional :: solinf - ! Fill in the 'resvect' or not - logical, intent(in), optional :: resvects_ - type(SpMtx),optional :: CoarseMtx_ ! Coarse matrix - - real(kind=rk) :: tol ! Tolerance - integer :: maxit ! Max number of iterations - logical :: resvects=.false. - - !real(kind=rk), dimension(size(b)) :: r, p, q, m, z, b_k - real(kind=rk),dimension(:),pointer :: r, p, q, m, z, b_k - real(kind=rk) :: rho_curr, rho_prev, init_norm, res_norm - real(kind=rk) :: alpha, beta, tmp, ratio_norm - integer :: iter,nfreds - integer :: ierr - ! + test - float(kind=rk), dimension(:), pointer :: arr_copy, gv_tmp - -!!$ real(kind=rk) :: r_max - - write(stream,'(/a)') 'Preconditioned conjugate gradient:' - - if (size(b) /= size(x)) & - call DOUG_abort('[pcg] : SEVERE : size(b) /= size(x)',-1) - -!!$ if (present(x0_)) -!!$ ! use method with initial guess - - - nfreds=size(b) - allocate(r(nfreds)) - allocate(p(nfreds)) - allocate(q(nfreds)) - allocate(m(nfreds)) - allocate(z(nfreds)) - allocate(b_k(nfreds)) - - - tol = sctls%solve_tolerance - if (present(tol_)) & - tol = tol_ - - maxit = sctls%solve_maxiters - if (maxit==-1) maxit=100000 - if (present(maxit_)) & - maxit = maxit_ - - if (present(resvects_).and.(.not.present(solinf))) & - call DOUG_abort('[pcg] : SEVERE : "resvects_" must be given along'//& - ' with "solinf".',-1) - - ! Allocate convergence info structure - if (present(solinf).and.(.not.present(resvects_))) then - call ConvInf_Init(solinf) - else if (present(solinf)) then - if (present(resvects_).and.(resvects_.eqv.(.true.))) then - call ConvInf_Init(solinf, maxit) - resvects = .true. - end if - end if - - ! Initialise auxiliary data structures - ! to assist with pmvm - call pmvmCommStructs_init(A, Msh) - - ! + test - if (ismaster()) & - allocate(gv_tmp(Msh%ngf)) - allocate(arr_copy(Msh%nlf)) - - ! Build preconditioner - ! call preconditioner(A, m, CoarseMtx_) -!!$ ! + test -!!$ arr_copy = m -!!$ call Vect_newToOldPerm(arr_copy) -!!$ gv_tmp = 0.0_rk -!!$ call Vect_Gather(arr_copy, gv_tmp, Msh) -!!$ if (ismaster()) & -!!$ call Vect_Print(gv_tmp, 'm ::') - -!!$ ! + test -!!$ !b = 1.0_rk -!!$ arr_copy = b -!!$ call Vect_newToOldPerm(arr_copy) -!!$ call Vect_Print(arr_copy, 'b local (original ordering) ::') -!!$ gv_tmp = 0.0_rk -!!$ call Vect_Gather(arr_copy, gv_tmp, Msh) -!!$ if (ismaster()) & -!!$ call Vect_Print(gv_tmp, 'b global (original ordering) ::') - -!!$ ! + test -!!$ arr_copy = x -!!$ call Vect_newToOldPerm(arr_copy) -!!$ call Vect_Print(arr_copy, 'x local (original ordering) ::') -!!$ gv_tmp = 0.0_rk -!!$ call Vect_Gather(arr_copy, gv_tmp, Msh) -!!$ if (ismaster()) & -!!$ call Vect_Print(gv_tmp, 'x global (original ordering) ::') - - !r = b - SpMtx_pmvm(A, x, Msh) - call SpMtx_pmvm(r,A,x,Msh) - r = b - r - !call Vect_Print(r,'initial residual') - - ! + test -!!$ arr_copy = r -!!$ call Vect_newToOldPerm(arr_copy) -!!$ gv_tmp = 0.0_rk -!!$ call Vect_Gather(arr_copy, gv_tmp, Msh) -!!$ if (ismaster()) & -!!$ call Vect_Print(gv_tmp, 'r ::') - - init_norm = Vect_dot_product(r,r) - if (init_norm == 0.0) & - init_norm = 1.0 - - write(stream,'(a,i6,a,e8.2)') 'maxit = ',maxit,', tol = ',tol - write(stream,'(a,e22.15)') 'init_norm = ', dsqrt(init_norm) - - iter = 1 - ratio_norm = 1.0_rk - do while((ratio_norm > tol*tol).and.(iter <= maxit)) - call msolve(m,r,z) - !if (present(CoarseMtx_)) then - ! call preconditioner(z,A,r,CoarseMtx_) - ! !z=r - !else - ! call preconditioner(z,A,r) - !endif - - !call Vect_Print(z,'preconditioner') - - ! compute current rho - rho_curr = Vect_dot_product(r,z) - if (iter == 1) then - p = z - else - beta = rho_curr / rho_prev - p = z + beta * p - end if - !q = SpMtx_pmvm(A, p, Msh) - call SpMtx_pmvm(q,A,p,Msh) - ! compute alpha - alpha = rho_curr / Vect_dot_product(p,q) - x = x + alpha * p - r = r - alpha * q - rho_prev = rho_curr - ! check - res_norm = Vect_dot_product(r,r) - ratio_norm = dsqrt(res_norm) / dsqrt(init_norm) -!!$ ratio_norm = res_norm / init_norm - if (ismaster()) & - write(stream, '(i5,a,e22.15)') iter,': res_norm=',dsqrt(res_norm) - !write(stream, '(i5,a,e22.15,e22.15)') iter,': dsqrt(ratio_norm)=',& - !dsqrt(ratio_norm), dsqrt(res_norm) - iter = iter + 1 - end do - - ! Deallocate auxiliary data structures - ! helped to assist with pmvm - call pmvmCommStructs_destroy() - - deallocate(b_k) - deallocate(z) - deallocate(m) - deallocate(q) - deallocate(p) - deallocate(r) - if (associated(gv_tmp)) deallocate(gv_tmp) - if (associated(arr_copy)) deallocate(arr_copy) - end subroutine pcg subroutine prec1Level(DD,sol,A,rhs,A_ghost,refactor)