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

File:

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

src/solvers/pcg.F90

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