Changes in src/solvers/pcg.F90 [b378423:08dfb5f]

File:

• src/solvers/pcg.F90 (modified) (2 diffs)

src/solvers/pcg.F90

@@ -41 +41 @@
 #endif
 
-  integer,save :: coarseid=0
   private :: &
        preconditioner, &
@@ -50 +49 @@
 
 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)