MODULE CG_Data_Types USE Precision USE Error_Handling USE System_Monitors IMPLICIT NONE PUBLIC::CG_Errors,CG_Timers,CG_Solver,SPD_Matrix,CG_System PRIVATE INTEGER,PARAMETER,PUBLIC::CG_free_solver=0,CG_inactive_solver=-1,CG_active_solver=1 INTEGER,PARAMETER,PUBLIC::CG_unconverged=-1,CG_converged=-2,CG_working=0 TYPE CG_Errors REAL(KIND=r_wp)::relative_error=100.0*EPSILON(1.0_r_wp),& residual_norm=-1.0_r_wp,largest_residual=-1.0_r_wp ENDTYPE TYPE CG_Timers INTEGER::preconditioning_timer=-1,dot_timer=-1,& vector_timer=-1,multiplication_timer=-1 ENDTYPE TYPE CG_Solver INTEGER::n_iterations=0,max_iterations=0,min_iterations=-1 INTEGER::status=CG_free_solver INTEGER::convergence=CG_unconverged TYPE(CG_Errors)::errors TYPE(CG_Timers)::timers INTEGER::log_unit=0 INTEGER(KIND=i_wp)::lb=1,ub=0 REAL(KIND=r_wp),DIMENSION(:),POINTER::cg_Vx=>NULL(),cg_Vy=>NULL(),& cg_Vz=>NULL(),cg_residuals=>NULL() REAL(KIND=r_wp),DIMENSION(:),POINTER::cg_x=>NULL() ENDTYPE TYPE SPD_Matrix INTEGER::PIN=1 ENDTYPE TYPE CG_System INTEGER::PIN=1 TYPE(SPD_Matrix),POINTER::matrix=>NULL() TYPE(CG_Solver),POINTER::solver=>NULL() REAL(KIND=r_wp),DIMENSION(:),POINTER::cg_b=>NULL() ENDTYPE END MODULE CG_Data_Types MODULE Conjugate_Gradient USE Precision USE Error_Handling USE System_Monitors USE CG_Data_Types USE Vector_Operations IMPLICIT NONE PUBLIC::PreconditionedConjugateGradient PRIVATE PUBLIC::CG_Errors,CG_Timers,CG_Solver,CG_System PUBLIC::CG_free_solver,CG_inactive_solver,CG_active_solver,& CG_unconverged,CG_converged,CG_working CONTAINS SUBROUTINE PreconditionedConjugateGradient(system,& MatrixVectorMultiplication,Preconditioner) IMPLICIT NONE TYPE(CG_System),INTENT(INOUT)::system OPTIONAL::Preconditioner INTERFACE SUBROUTINE MatrixVectorMultiplication(system) USE CG_Data_Types TYPE(CG_System),INTENT(INOUT)::system END SUBROUTINE MatrixVectorMultiplication SUBROUTINE Preconditioner(system) USE CG_Data_Types TYPE(CG_System),INTENT(INOUT)::system END SUBROUTINE Preconditioner END INTERFACE TYPE(CG_Solver),POINTER::solver REAL(KIND=r_wp)::relative_error_tolerance,norm_2_tolerance,norm_inf_tolerance,& eps,alpha,beta,gamma,rho_new,rho_old,& residual,norm_rhs INTEGER::iteration,p_timer,d_timer,v_timer,m_timer INTEGER(KIND=i_wp)::variable,n_ub,n_lb LOGICAL::log_convergence,use_preconditioning,time_p,time_d,time_v,time_m,& converged,check_2_norm,check_inf_norm solver=>system%solver n_lb=solver%lb n_ub=solver%ub solver%status=CG_active_solver solver%convergence=CG_working p_timer=(solver%timers%preconditioning_timer) time_p=(p_timer>0) d_timer=(solver%timers%dot_timer) time_d=(d_timer>0) v_timer=(solver%timers%vector_timer) time_v=(v_timer>0) m_timer=(solver%timers%multiplication_timer) time_m=(m_timer>0) use_preconditioning=PRESENT(Preconditioner) log_convergence=(solver%log_unit/=0) eps=10.0_r_wp*EPSILON(1.0_r_wp) relative_error_tolerance=MAX(10.0_r_wp*eps,solver%errors%relative_error) check_2_norm=(solver%errors%residual_norm>eps) IF(check_2_norm)norm_2_tolerance=MAX(10.0_r_wp*eps,solver%errors%residual_norm) check_inf_norm=(solver%errors%largest_residual>eps) IF(check_inf_norm)norm_inf_tolerance=MAX(10.0_r_wp*eps,solver%errors%largest_residual) IF(time_v)CALL StartTimer(v_timer) CALL VectorCopy(target=solver%cg_Vx(n_lb:n_ub),source=solver%cg_x(n_lb:n_ub)) CALL StopTimer(v_timer) IF(time_m)CALL StartTimer(m_timer) CALL MatrixVectorMultiplication(system) IF(time_m)CALL StopTimer(m_timer) IF(time_v)CALL StartTimer(v_timer) CALL VectorSubtraction(from=system%cg_b(n_lb:n_ub),what=solver%cg_Vy(n_lb:n_ub),difference=solver%cg_residuals(n_lb:n_ub)) IF(time_v)CALL StopTimer(v_timer) iteration=0 CG:DO iteration=iteration+1 IF(use_preconditioning)THEN IF(time_p)CALL StartTimer(p_timer) CALL Preconditioner(system) IF(time_p)CALL StopTimer(p_timer) IF(time_d)CALL StartTimer(d_timer) rho_new=DOT_PRODUCT(solver%cg_residuals(n_lb:n_ub),solver%cg_Vz(n_lb:n_ub)) IF(time_d)CALL StopTimer(d_timer) ELSE IF(time_d)CALL StartTimer(d_timer) rho_new=DOT_PRODUCT(solver%cg_residuals(n_lb:n_ub),solver%cg_residuals(n_lb:n_ub)) IF(time_d)CALL StopTimer(d_timer) END IF IF(iteration==1)norm_rhs=SQRT(rho_new) IF(use_preconditioning)THEN solver%errors%relative_error=SQRT(rho_new)/ABS(norm_rhs+eps) ELSE solver%errors%relative_error=SQRT(rho_new)/ABS(norm_rhs+eps) END IF IF(log_convergence)THEN WRITE(UNIT=solver%log_unit,FMT="(I10,E10.3)",ADVANCE="NO")& iteration,solver%errors%relative_error END IF IF(check_2_norm)THEN IF(.NOT.use_preconditioning)THEN solver%errors%residual_norm=SQRT(rho_new) ELSE solver%errors%residual_norm=SQRT(DOT_PRODUCT(solver%cg_residuals(n_lb:n_ub),solver%cg_residuals(n_lb:n_ub))) END IF IF(log_convergence)THEN WRITE(UNIT=solver%log_unit,FMT="(E10.3)",ADVANCE="NO")solver%errors%residual_norm END IF END IF IF(check_inf_norm)THEN solver%errors%largest_residual=0.0_r_wp DO variable=INT(LBOUND(solver%cg_residuals(n_lb:n_ub),DIM=1),KIND=i_wp),INT(UBOUND(solver%cg_residuals(n_lb:n_ub),DIM=1),KIN& &D=i_wp) residual=solver%cg_residuals(variable) IF(ABS(residual)>=solver%errors%largest_residual)solver%errors%largest_residual=ABS(residual) END DO IF(log_convergence)THEN WRITE(UNIT=solver%log_unit,FMT="(E10.3)",ADVANCE="NO")solver%errors%largest_residual END IF END IF TestConvergence:IF(iteration>solver%max_iterations)THEN solver%convergence=CG_unconverged EXIT CG ELSE IF(solver%errors%relative_errornorm_2_tolerance)converged=.FALSE. END IF IF(check_inf_norm)THEN IF(solver%errors%largest_residual>norm_inf_tolerance)converged=.FALSE. END IF IF(converged.AND.(iteration>solver%min_iterations))THEN solver%convergence=CG_converged EXIT CG END IF END IF TestConvergence IF(log_convergence)THEN WRITE(solver%log_unit,*) END IF IF(use_preconditioning)THEN IF(time_v)CALL StartTimer(v_timer) IF(iteration==1)THEN CALL VectorCopy(target=solver%cg_Vx(n_lb:n_ub),source=solver%cg_Vz(n_lb:n_ub)) ELSE beta=rho_new/rho_old CALL VectorAXPBY(beta=beta,y=solver%cg_Vx(n_lb:n_ub),x=solver%cg_Vz(n_lb:n_ub)) END IF IF(time_v)CALL StopTimer(v_timer) ELSE IF(time_v)CALL StartTimer(v_timer) IF(iteration==1)THEN CALL VectorCopy(target=solver%cg_Vx(n_lb:n_ub),source=solver%cg_residuals(n_lb:n_ub)) ELSE beta=rho_new/rho_old CALL VectorAXPBY(beta=beta,y=solver%cg_Vx(n_lb:n_ub),x=solver%cg_residuals(n_lb:n_ub)) END IF IF(time_v)CALL StopTimer(v_timer) END IF rho_old=rho_new IF(time_m)CALL StartTimer(m_timer) CALL MatrixVectorMultiplication(system) IF(time_m)CALL StopTimer(m_timer) IF(time_d)CALL StartTimer(d_timer) gamma=DOT_PRODUCT(solver%cg_Vx(n_lb:n_ub),solver%cg_Vy(n_lb:n_ub)) alpha=rho_new/gamma IF(time_d)CALL StopTimer(d_timer) IF(time_v)CALL StartTimer(v_timer) CALL VectorAXPBY(alpha=alpha,x=solver%cg_Vx(n_lb:n_ub),y=solver%cg_x(n_lb:n_ub)) CALL VectorAXPBY(alpha=-alpha,x=solver%cg_Vy(n_lb:n_ub),y=solver%cg_residuals(n_lb:n_ub)) IF(time_v)CALL StopTimer(v_timer) END DO CG solver%n_iterations=iteration-1 solver%status=CG_inactive_solver END SUBROUTINE PreconditionedConjugateGradient END MODULE Conjugate_Gradient