From bd2af001b055c7c81d4a6adf5152bef4eee4715f Mon Sep 17 00:00:00 2001 From: william-dawson Date: Tue, 12 Sep 2023 09:14:56 +0900 Subject: [PATCH] Automatic Convergence Monitoring (#218) * Convergence monitor feature * lint * A more sensible convergence check for pseudo inverse * pure functions * Improve estimate of the gap by using power method * Fix power bounds guess * Improve power bounds implementation * Check consistent of power bounds estimates * Fix printing of estimate gap * Fix window bug --- Source/C/SolverParameters_c.h | 1 + Source/CPlusPlus/SolverParameters.cc | 5 + Source/CPlusPlus/SolverParameters.h | 3 + Source/Fortran/CMakeLists.txt | 1 + Source/Fortran/ConvergenceMonitorModule.F90 | 193 ++++++++++++++++++ Source/Fortran/DensityMatrixSolversModule.F90 | 68 +++--- Source/Fortran/EigenBoundsModule.F90 | 115 +++++++---- Source/Fortran/EigenSolversModule.F90 | 22 +- Source/Fortran/GeometryOptimizationModule.F90 | 20 +- Source/Fortran/InverseSolversModule.F90 | 36 ++-- Source/Fortran/LinearSolversModule.F90 | 15 +- Source/Fortran/RootSolversModule.F90 | 14 +- Source/Fortran/SignSolversModule.F90 | 18 +- Source/Fortran/SolverParametersModule.F90 | 40 +++- Source/Fortran/SquareRootSolversModule.F90 | 35 ++-- Source/Wrapper/SolverParametersModule_wrp.F90 | 24 ++- 16 files changed, 450 insertions(+), 160 deletions(-) create mode 100644 Source/Fortran/ConvergenceMonitorModule.F90 diff --git a/Source/C/SolverParameters_c.h b/Source/C/SolverParameters_c.h index 8f1bf02f..95f993a0 100644 --- a/Source/C/SolverParameters_c.h +++ b/Source/C/SolverParameters_c.h @@ -8,6 +8,7 @@ void SetParametersBeVerbose_wrp(int *ih_this, const bool *new_value); void SetParametersThreshold_wrp(int *ih_this, const double *new_value); void SetParametersLoadBalance_wrp(int *ih_this, const int *ih_permutation); void SetParametersStepThreshold_wrp(int *ih_this, const double *new_value); +void SetParametersMonitorConvergence_wrp(int *ih_this, const bool *new_value); void DestructSolverParameters_wrp(int *ih_this); #endif diff --git a/Source/CPlusPlus/SolverParameters.cc b/Source/CPlusPlus/SolverParameters.cc index 282fe14b..5ef1f6dc 100644 --- a/Source/CPlusPlus/SolverParameters.cc +++ b/Source/CPlusPlus/SolverParameters.cc @@ -43,4 +43,9 @@ void SolverParameters::SetLoadBalance(const Permutation &new_value) { void SolverParameters::SetStepThreshold(double new_value) { SetParametersStepThreshold_wrp(ih_this, &new_value); } + +//////////////////////////////////////////////////////////////////////////////// +void SolverParameters::SetMonitorConvergence(bool new_value) { + SetParametersMonitorConvergence_wrp(ih_this, &new_value); +} } // namespace NTPoly diff --git a/Source/CPlusPlus/SolverParameters.h b/Source/CPlusPlus/SolverParameters.h index fab640a5..d023c416 100644 --- a/Source/CPlusPlus/SolverParameters.h +++ b/Source/CPlusPlus/SolverParameters.h @@ -32,6 +32,9 @@ class SolverParameters { //! Thresholds for step size searches. //!\param new_value void SetStepThreshold(double new_value); + //! Whether to automatically monitor convergence + //!\param new_value + void SetMonitorConvergence(bool new_value); ~SolverParameters(); private: diff --git a/Source/Fortran/CMakeLists.txt b/Source/Fortran/CMakeLists.txt index 11e002b0..73248817 100644 --- a/Source/Fortran/CMakeLists.txt +++ b/Source/Fortran/CMakeLists.txt @@ -3,6 +3,7 @@ set(Fsrc AnalysisModule.F90 ChebyshevSolversModule.F90 CholeskySolversModule.F90 + ConvergenceMonitorModule.F90 DataTypesModule.F90 DensityMatrixSolversModule.F90 DMatrixModule.F90 diff --git a/Source/Fortran/ConvergenceMonitorModule.F90 b/Source/Fortran/ConvergenceMonitorModule.F90 new file mode 100644 index 00000000..e0daf25f --- /dev/null +++ b/Source/Fortran/ConvergenceMonitorModule.F90 @@ -0,0 +1,193 @@ +!> A module for monitoring convergence of an iterative algorithim. +!! In basic mode, we monitor that the last value isn't below the tight cutoff. +!! In automatic mode we monitor the following conditions: +!! o The last value can't be negative +!! o The moving average is within an order of magnitude +!! o The value isn't above the loose cutoff +MODULE ConvergenceMonitor + USE DataTypesModule, ONLY : NTREAL + USE LoggingModule, ONLY : EnterSubLog, ExitSubLog, WriteElement, & + & WriteListElement + IMPLICIT NONE + PRIVATE + !> Monitor convergence with a moving average + TYPE, PUBLIC :: Monitor_t + !> The small window + REAL(NTREAL), DIMENSION(:), ALLOCATABLE :: win_short + !> The large window + REAL(NTREAL), DIMENSION(:), ALLOCATABLE :: win_long + !> The number of values that have been added + INTEGER :: nval + !> We aren't converged if the average isn't below this. + REAL(NTREAL) :: loose_cutoff + !> We definitely are converged if the last value is below this. + REAL(NTREAL) :: tight_cutoff + !> Whether to do automatic exit + LOGICAL :: automatic + END TYPE Monitor_t +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + PUBLIC :: ConstructMonitor + PUBLIC :: DestructMonitor + PUBLIC :: AppendValue + PUBLIC :: CheckConverged +CONTAINS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Construct the monitor. + PURE SUBROUTINE ConstructMonitor(this, short_len_in, long_len_in, & + & loose_cutoff_in, tight_cutoff_in, automatic_in) + !> The monitor to construct. + TYPE(Monitor_t), INTENT(INOUT) :: this + !> The length of the short window (default: 3) + INTEGER, INTENT(IN), OPTIONAL :: short_len_in + !> The length of the long window (default: 6) + INTEGER, INTENT(IN), OPTIONAL :: long_len_in + !> If the average is greater than this than we aren't + !! converged (default: 0.01) + REAL(NTREAL), INTENT(IN), OPTIONAL :: loose_cutoff_in + !> If the last value is less than this, we definitely are converged + !! (default: 1e-8). + REAL(NTREAL), INTENT(IN), OPTIONAL :: tight_cutoff_in + !> Whether to use automatic mode (default: .true.) + LOGICAL, INTENT(IN), OPTIONAL :: automatic_in + + CALL DestructMonitor(this) + + !> Allocate Windows + IF (PRESENT(short_len_in)) THEN + ALLOCATE(this%win_short(short_len_in)) + ELSE + ALLOCATE(this%win_short(3)) + END IF + this%win_short = 0 + + IF (PRESENT(long_len_in)) THEN + ALLOCATE(this%win_long(long_len_in)) + ELSE + ALLOCATE(this%win_long(6)) + END IF + this%win_long = 0 + + IF (PRESENT(loose_cutoff_in)) THEN + this%loose_cutoff = loose_cutoff_in + ELSE + this%loose_cutoff = 1E-2_NTREAL + END IF + + IF (PRESENT(tight_cutoff_in)) THEN + this%tight_cutoff = tight_cutoff_in + ELSE + this%tight_cutoff = 1E-8_NTREAL + END IF + + IF (PRESENT(automatic_in)) THEN + this%automatic = automatic_in + ELSE + this%automatic = .TRUE. + END IF + + this%nval = 0 + END SUBROUTINE ConstructMonitor +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Destruct the monitor. + PURE SUBROUTINE DestructMonitor(this) + !> The monitor to destruct. + TYPE(Monitor_t), INTENT(INOUT) :: this + IF (ALLOCATED(this%win_short)) DEALLOCATE(this%win_short) + IF (ALLOCATED(this%win_long)) DEALLOCATE(this%win_long) + END SUBROUTINE DestructMonitor +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Add a value to the window + PURE SUBROUTINE AppendValue(this, val) + !> Add a value to the window + TYPE(Monitor_t), INTENT(INOUT) :: this + !> Value to add + REAL(NTREAL), INTENT(IN) :: val + !! Local Variables + INTEGER :: II + + !! Shift + DO II = 1, SIZE(this%win_short) -1 + this%win_short(II) = this%win_short(II + 1) + END DO + DO II = 1, SIZE(this%win_long) - 1 + this%win_long(II) = this%win_long(II + 1) + END DO + + !! Append + this%win_short(SIZE(this%win_short)) = val + this%win_long(SIZE(this%win_long)) = val + this%nval = this%nval + 1 + END SUBROUTINE AppendValue +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Check if Convergence Has Been Achieved + FUNCTION CheckConverged(this, be_verbose) RESULT(conv) + !> The window to test + TYPE(Monitor_t), INTENT(IN) :: this + !> True if converged. + LOGICAL :: conv + !> True if we should print convergence info + LOGICAL :: be_verbose + !! Local Variables + REAL(NTREAL) :: avg_short, avg_long + REAL(NTREAL) :: last, last2 + + !! Basic Check + last = this%win_short(SIZE(this%win_short)) + last2 = this%win_short(SIZE(this%win_short) - 1) + IF (be_verbose) THEN + CALL WriteListElement(key = "Convergence", VALUE = last) + END IF + IF (ABS(last) .GT. this%tight_cutoff) THEN + conv = .FALSE. + ELSE + conv = .TRUE. + CALL EnterSubLog + CALL WriteElement(key = "Trigger", VALUE = "Tight Criteria") + CALL ExitSubLog + END IF + IF (.NOT. this%automatic .OR. conv) RETURN + + !! Automatic disabled + conv = .TRUE. + + !! First check that we've seen enough values to make a judgement. + IF (this%nval .LT. SIZE(this%win_long)) conv = .FALSE. + + !! Compute Averages + avg_short = SUM(this%win_short) / SIZE(this%win_short) + avg_long = SUM(this%win_long) / SIZE(this%win_long) + IF (be_verbose) THEN + CALL EnterSubLog + CALL WriteElement(key = "Avg Short", VALUE = avg_short) + CALL WriteElement(key = "Avg Long", VALUE = avg_long) + CALL ExitSubLog + END IF + + !! Now check that the two windows are within an order of magnitude + IF (.NOT. (10 * avg_short .GT. avg_long .AND. & + & avg_short / 10 .LT. avg_long)) THEN + conv = .FALSE. + END IF + + !! The last convergence value should also be within an order of magnitude + !! of the long average. + IF (.NOT. (10 * last .GT. avg_long .AND. last / 10 .LT. avg_long)) THEN + conv = .FALSE. + END IF + + !! No convergence if the lastest value is negative + IF (last .LT. 0) conv = .FALSE. + + !! No convergence if the absolute value is decreasing in magnitude + IF (ABS(last) < ABS(last2)) conv = .FALSE. + + !! No convergence if the average value is greater than the cutoff + IF (avg_long .GT. this%loose_cutoff) conv = .FALSE. + + IF (conv) THEN + CALL EnterSubLog + CALL WriteElement(key = "Trigger", VALUE = "Automatic") + CALL ExitSubLog + END IF + END FUNCTION CheckConverged +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +END MODULE ConvergenceMonitor diff --git a/Source/Fortran/DensityMatrixSolversModule.F90 b/Source/Fortran/DensityMatrixSolversModule.F90 index ddea7e9f..27c9e7bf 100644 --- a/Source/Fortran/DensityMatrixSolversModule.F90 +++ b/Source/Fortran/DensityMatrixSolversModule.F90 @@ -1,6 +1,7 @@ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> A Module For Solving Quantum Chemistry Systems using Purification. MODULE DensityMatrixSolversModule + USE ConvergenceMonitor, ONLY : ConstructMonitor, CheckConverged, AppendValue USE DataTypesModule, ONLY : NTREAL, MPINTREAL USE EigenBoundsModule, ONLY : GershgorinBounds USE FermiOperatorModule, ONLY : ComputeDenseFOE @@ -64,7 +65,6 @@ SUBROUTINE PM(H, ISQ, trace, K, & REAL(NTREAL), DIMENSION(:), ALLOCATABLE :: sigma_array REAL(NTREAL) :: trace_value REAL(NTREAL) :: trace_value2 - REAL(NTREAL) :: norm_value REAL(NTREAL) :: energy_value, energy_value_old !! For computing the chemical potential REAL(NTREAL) :: zero_value, midpoint, interval_a, interval_b @@ -79,6 +79,9 @@ SUBROUTINE PM(H, ISQ, trace, K, & ELSE CALL ConstructSolverParameters(params) END IF + CALL ConstructMonitor(params%monitor, & + & automatic_in = params%monitor_convergence, & + & tight_cutoff_in=params%converge_diff) IF (params%be_verbose) THEN CALL WriteHeader("Density Matrix Solver") @@ -143,7 +146,6 @@ SUBROUTINE PM(H, ISQ, trace, K, & CALL EnterSubLog END IF II = 1 - norm_value = params%converge_diff + 1.0_NTREAL energy_value = 0.0_NTREAL DO II = 1, params%max_iterations !! Compute X_k2 @@ -190,18 +192,15 @@ SUBROUTINE PM(H, ISQ, trace, K, & !! Energy value based convergence energy_value_old = energy_value CALL DotMatrix(X_k, WH, energy_value) - norm_value = ABS(energy_value - energy_value_old) + CALL AppendValue(params%monitor, energy_value - energy_value_old) + IF (CheckConverged(params%monitor, params%be_verbose)) EXIT IF (params%be_verbose) THEN - CALL WriteListElement(key = "Convergence", VALUE = norm_value) CALL EnterSubLog CALL WriteElement("Energy Value", VALUE = energy_value) CALL ExitSubLog END IF - IF (norm_value .LE. params%converge_diff) THEN - EXIT - END IF END DO total_iterations = II - 1 IF (params%be_verbose) THEN @@ -310,7 +309,6 @@ SUBROUTINE TRS2(H, ISQ, trace, K, & REAL(NTREAL) :: e_min, e_max REAL(NTREAL), DIMENSION(:), ALLOCATABLE :: sigma_array REAL(NTREAL) :: trace_value - REAL(NTREAL) :: norm_value REAL(NTREAL) :: energy_value, energy_value_old !! For computing the chemical potential REAL(NTREAL) :: zero_value, midpoint, interval_a, interval_b @@ -325,6 +323,9 @@ SUBROUTINE TRS2(H, ISQ, trace, K, & ELSE CALL ConstructSolverParameters(params) END IF + CALL ConstructMonitor(params%monitor, & + & automatic_in = params%monitor_convergence, & + & tight_cutoff_in=params%converge_diff) IF (params%be_verbose) THEN CALL WriteHeader("Density Matrix Solver") @@ -375,7 +376,6 @@ SUBROUTINE TRS2(H, ISQ, trace, K, & CALL EnterSubLog END IF II = 1 - norm_value = params%converge_diff + 1.0_NTREAL energy_value = 0.0_NTREAL DO II = 1, params%max_iterations !! Compute Sigma @@ -402,18 +402,14 @@ SUBROUTINE TRS2(H, ISQ, trace, K, & !! Energy value based convergence energy_value_old = energy_value CALL DotMatrix(X_k, WH, energy_value) - norm_value = ABS(energy_value - energy_value_old) + CALL AppendValue(params%monitor, energy_value - energy_value_old) + IF (CheckConverged(params%monitor, params%be_verbose)) EXIT IF (params%be_verbose) THEN - CALL WriteListElement(key = "Convergence", VALUE = norm_value) CALL EnterSubLog CALL WriteElement("Energy Value", VALUE = energy_value) CALL ExitSubLog END IF - - IF (norm_value .LE. params%converge_diff) THEN - EXIT - END IF END DO total_iterations = II - 1 IF (params%be_verbose) THEN @@ -514,7 +510,6 @@ SUBROUTINE TRS4(H, ISQ, trace, K, & !! Local Variables REAL(NTREAL) :: e_min, e_max REAL(NTREAL), DIMENSION(:), ALLOCATABLE :: sigma_array - REAL(NTREAL) :: norm_value REAL(NTREAL) :: energy_value, energy_value_old !! For computing the chemical potential REAL(NTREAL) :: zero_value, midpoint, interval_a, interval_b @@ -531,6 +526,9 @@ SUBROUTINE TRS4(H, ISQ, trace, K, & ELSE CALL ConstructSolverParameters(params) END IF + CALL ConstructMonitor(params%monitor, & + & automatic_in = params%monitor_convergence, & + & tight_cutoff_in=params%converge_diff) IF (params%be_verbose) THEN CALL WriteHeader("Density Matrix Solver") @@ -584,7 +582,6 @@ SUBROUTINE TRS4(H, ISQ, trace, K, & CALL EnterSubLog END IF II = 1 - norm_value = params%converge_diff + 1.0_NTREAL energy_value = 0.0_NTREAL DO II = 1, params%max_iterations !! Compute X_k2 @@ -630,18 +627,14 @@ SUBROUTINE TRS4(H, ISQ, trace, K, & !! Energy value based convergence energy_value_old = energy_value CALL DotMatrix(X_k, WH, energy_value) - norm_value = ABS(energy_value - energy_value_old) + CALL AppendValue(params%monitor, energy_value - energy_value_old) + IF (CheckConverged(params%monitor, params%be_verbose)) EXIT IF (params%be_verbose) THEN - CALL WriteListElement(key = "Convergence", VALUE = norm_value) CALL EnterSubLog - CALL WriteElement(key = "Energy Value", VALUE = energy_value) + CALL WriteElement("Energy Value", VALUE = energy_value) CALL ExitSubLog END IF - - IF (norm_value .LE. params%converge_diff) THEN - EXIT - END IF END DO total_iterations = II - 1 IF (params%be_verbose) THEN @@ -756,7 +749,6 @@ SUBROUTINE HPCP(H, ISQ, trace, K, & REAL(NTREAL) :: mu REAL(NTREAL), DIMENSION(:), ALLOCATABLE :: sigma_array REAL(NTREAL) :: trace_value - REAL(NTREAL) :: norm_value REAL(NTREAL) :: energy_value, energy_value_old !! For computing the chemical potential REAL(NTREAL) :: zero_value, midpoint, interval_a, interval_b @@ -772,6 +764,9 @@ SUBROUTINE HPCP(H, ISQ, trace, K, & ELSE CALL ConstructSolverParameters(params) END IF + CALL ConstructMonitor(params%monitor, & + & automatic_in = params%monitor_convergence, & + & tight_cutoff_in=params%converge_diff) IF (params%be_verbose) THEN CALL WriteHeader("Density Matrix Solver") @@ -840,7 +835,6 @@ SUBROUTINE HPCP(H, ISQ, trace, K, & END IF II = 1 - norm_value = params%converge_diff + 1.0_NTREAL energy_value = 0.0_NTREAL DO II = 1, params%max_iterations !! Compute the hole matrix DH @@ -852,7 +846,6 @@ SUBROUTINE HPCP(H, ISQ, trace, K, & CALL MatrixMultiply(D1, DH, DDH, & & threshold_in = params%threshold, memory_pool_in = pool) CALL MatrixTrace(DDH, trace_value) - norm_value = ABS(trace_value) !! Compute D2DH CALL MatrixMultiply(D1, DDH, D2DH, & @@ -874,18 +867,14 @@ SUBROUTINE HPCP(H, ISQ, trace, K, & !! Energy value based convergence energy_value_old = energy_value CALL DotMatrix(D1, WH, energy_value) - norm_value = ABS(energy_value - energy_value_old) + CALL AppendValue(params%monitor, energy_value - energy_value_old) + IF (CheckConverged(params%monitor, params%be_verbose)) EXIT IF (params%be_verbose) THEN - CALL WriteListElement(key = "Convergence", VALUE = norm_value) CALL EnterSubLog CALL WriteElement("Energy Value", VALUE = energy_value) CALL ExitSubLog END IF - - IF (norm_value .LE. params%converge_diff) THEN - EXIT - END IF END DO total_iterations = II - 1 IF (params%be_verbose) THEN @@ -990,7 +979,6 @@ SUBROUTINE ScaleAndFold(H, ISQ, trace, K, & REAL(NTREAL) :: e_min, e_max REAL(NTREAL) :: Beta, BetaBar, alpha REAL(NTREAL) :: trace_value - REAL(NTREAL) :: norm_value REAL(NTREAL) :: energy_value, energy_value_old !! Temporary Variables TYPE(MatrixMemoryPool_p) :: pool @@ -1002,6 +990,9 @@ SUBROUTINE ScaleAndFold(H, ISQ, trace, K, & ELSE CALL ConstructSolverParameters(params) END IF + CALL ConstructMonitor(params%monitor, & + & automatic_in = params%monitor_convergence, & + & tight_cutoff_in=params%converge_diff) IF (params%be_verbose) THEN CALL WriteHeader("Density Matrix Solver") @@ -1053,7 +1044,6 @@ SUBROUTINE ScaleAndFold(H, ISQ, trace, K, & CALL EnterSubLog END IF II = 1 - norm_value = params%converge_diff + 1.0_NTREAL energy_value = 0.0_NTREAL DO II = 1, params%max_iterations !! Determine the path @@ -1081,18 +1071,14 @@ SUBROUTINE ScaleAndFold(H, ISQ, trace, K, & energy_value_old = energy_value CALL DotMatrix(X_k, WH, energy_value) energy_value = 2.0_NTREAL*energy_value - norm_value = ABS(energy_value - energy_value_old) + CALL AppendValue(params%monitor, energy_value - energy_value_old) + IF (CheckConverged(params%monitor, params%be_verbose)) EXIT IF (params%be_verbose) THEN - CALL WriteListElement(key = "Convergence", VALUE = norm_value) CALL EnterSubLog CALL WriteElement("Energy Value", VALUE = energy_value) CALL ExitSubLog END IF - - IF (norm_value .LE. params%converge_diff) THEN - EXIT - END IF END DO IF (params%be_verbose) THEN CALL ExitSubLog diff --git a/Source/Fortran/EigenBoundsModule.F90 b/Source/Fortran/EigenBoundsModule.F90 index 2fa3e1c8..7c14a4d7 100644 --- a/Source/Fortran/EigenBoundsModule.F90 +++ b/Source/Fortran/EigenBoundsModule.F90 @@ -1,6 +1,7 @@ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> A module for computing estimates of the bounds of the spectrum of a matrix. MODULE EigenBoundsModule + USE ConvergenceMonitor, ONLY : ConstructMonitor, CheckConverged, AppendValue USE DataTypesModule, ONLY : NTREAL, MPINTREAL USE LoggingModule, ONLY : EnterSubLog, ExitSubLog, WriteElement, & & WriteListElement, WriteHeader @@ -64,11 +65,13 @@ SUBROUTINE PowerBounds(this, max_value, solver_parameters_in) !> The parameters for this calculation. TYPE(SolverParameters_t), INTENT(IN), OPTIONAL :: solver_parameters_in !! Handling Optional Parameters - TYPE(SolverParameters_t) :: param + TYPE(SolverParameters_t) :: params !! Local Data + REAL(NTREAL), DIMENSION(3) :: ritz_values + REAL(NTREAL), DIMENSION(3) :: aitken_values + REAL(NTREAL) :: num, den TYPE(Matrix_ps) :: vector, vector2 REAL(NTREAL) :: scale_value - REAL(NTREAL) :: norm_value TYPE(TripletList_r) :: temp_list TYPE(Triplet_r) :: temp_triplet INTEGER :: II @@ -76,16 +79,19 @@ SUBROUTINE PowerBounds(this, max_value, solver_parameters_in) !! Optional Parameters IF (PRESENT(solver_parameters_in)) THEN - CALL CopySolverParameters(solver_parameters_in, param) + CALL CopySolverParameters(solver_parameters_in, params) ELSE - CALL ConstructSolverParameters(param) - param%max_iterations = 10 + CALL ConstructSolverParameters(params) + params%max_iterations = 10 END IF + CALL ConstructMonitor(params%monitor, & + & automatic_in = params%monitor_convergence, & + & tight_cutoff_in=params%converge_diff) - IF (param%be_verbose) THEN + IF (params%be_verbose) THEN CALL WriteHeader("Power Bounds Solver") CALL EnterSubLog - CALL PrintParameters(param) + CALL PrintParameters(params) END IF !! Diagonal matrices serve as vectors. @@ -94,57 +100,84 @@ SUBROUTINE PowerBounds(this, max_value, solver_parameters_in) !! Guess Vector CALL ConstructTripletList(temp_list) - IF (this%process_grid%global_rank .EQ. 0) THEN - temp_triplet%index_row = 1 - temp_triplet%index_column = 1 - temp_triplet%point_value = 1.0_NTREAL - CALL AppendToTripletList(temp_list, temp_triplet) + IF (this%start_row .EQ. 1) THEN + DO II = this%start_column, this%end_column - 1 + temp_triplet%index_row = 1 + temp_triplet%index_column = II + temp_triplet%point_value = & + & 1.0_NTREAL / vector%actual_matrix_dimension + CALL AppendToTripletList(temp_list, temp_triplet) + END DO END IF - CALL FillMatrixFromTripletList(vector, temp_list) + CALL FillMatrixFromTripletList(vector, temp_list, & + & preduplicated_in=.TRUE., prepartitioned_in=.TRUE.) !! Iterate - IF (param%be_verbose) THEN + ritz_values = 0 + aitken_values = 0 + IF (params%be_verbose) THEN CALL WriteHeader("Iterations") CALL EnterSubLog END IF II = 1 - norm_value = param%converge_diff + 1.0_NTREAL - DO II = 1, param%max_iterations - IF (param%be_verbose .AND. II .GT. 1) THEN - CALL WriteListElement(key = "Convergence", VALUE = norm_value) - END IF - + DO II = 1, params%max_iterations !! x = Ax CALL MatrixMultiply(this, vector, vector2, & - & threshold_in = param%threshold, memory_pool_in = pool) + & threshold_in = params%threshold, memory_pool_in = pool) + + !! Compute the ritz value + CALL DotMatrix(vector, vector, scale_value) + CALL DotMatrix(vector, vector2, max_value) + max_value = max_value / scale_value + !! x = x/||x|| scale_value = 1.0 / MatrixNorm(vector2) CALL ScaleMatrix(vector2, scale_value) - - !! Check if Converged - CALL IncrementMatrix(vector2, vector, alpha_in = -1.0_NTREAL) - norm_value = MatrixNorm(vector) - CALL CopyMatrix(vector2, vector) - IF (norm_value .LE. param%converge_diff) THEN - EXIT + !! Aitken's Extrapolation + ritz_values(1) = ritz_values(2) + ritz_values(2) = ritz_values(3) + ritz_values(3) = max_value + aitken_values(1) = aitken_values(2) + aitken_values(2) = aitken_values(3) + IF (II .GE. 3) THEN + num = ritz_values(3)*ritz_values(1) - ritz_values(2)**2 + den = ritz_values(3) - 2 * ritz_values(2) + ritz_values(1) + !! Avoid division by zero + IF (ABS(den) .GT. 1E-14_NTREAL) THEN + aitken_values(3) = num / den + ELSE + aitken_values(3) = ritz_values(3) + END IF + ELSE + aitken_values(3) = ritz_values(3) END IF + + !! Check if Converged - pass the negative value because we're looking + !! for the largest eigenvalue value. + CALL AppendValue(params%monitor, & + & - (aitken_values(3) - aitken_values(2))) + IF (CheckConverged(params%monitor, params%be_verbose)) THEN + !! Make sure the two estimates vaguely agree + IF(ABS(aitken_values(3) - ritz_values(3)) .LT. & + & params%monitor%loose_cutoff) THEN + EXIT + END IF + END IF + IF (params%be_verbose) THEN + CALL EnterSubLog + CALL WriteElement(key="Estimate", VALUE=ritz_values(3)) + CALL WriteElement(key="Aitken Estimate", VALUE=aitken_values(3)) + CALL ExitSubLog + END IF + END DO - IF (param%be_verbose) THEN + max_value = aitken_values(3) + IF (params%be_verbose) THEN CALL ExitSubLog CALL WriteElement(key = "Total Iterations", VALUE = II - 1) - END IF - - !! Compute The Largest Eigenvalue - CALL DotMatrix(vector, vector, scale_value) - CALL MatrixMultiply(this, vector, vector2, & - & threshold_in = param%threshold, memory_pool_in = pool) - CALL DotMatrix(vector, vector2, max_value) - max_value = max_value / scale_value - - IF (param%be_verbose) THEN - CALL WriteElement(key = "Max Eigen Value", VALUE = max_value) + CALL WriteElement(key = "Max Eigen Value", VALUE = aitken_values(3)) CALL ExitSubLog END IF @@ -152,7 +185,7 @@ SUBROUTINE PowerBounds(this, max_value, solver_parameters_in) CALL DestructMatrix(vector) CALL DestructMatrix(vector2) CALL DestructMatrixMemoryPool(pool) - CALL DestructSolverParameters(param) + CALL DestructSolverParameters(params) END SUBROUTINE PowerBounds !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! END MODULE EigenBoundsModule diff --git a/Source/Fortran/EigenSolversModule.F90 b/Source/Fortran/EigenSolversModule.F90 index 58ee794b..459b726d 100644 --- a/Source/Fortran/EigenSolversModule.F90 +++ b/Source/Fortran/EigenSolversModule.F90 @@ -8,7 +8,8 @@ MODULE EigenSolversModule #if EIGENEXA USE EigenExaModule, ONLY : EigenExa_s #endif - USE LoggingModule, ONLY : EnterSubLog, ExitSubLog, WriteHeader, WriteElement + USE LoggingModule, ONLY : EnterSubLog, ExitSubLog, WriteHeader, & + & WriteComment, WriteElement USE PMatrixMemoryPoolModule, ONLY : MatrixMemoryPool_p, & & DestructMatrixMemoryPool USE PSMatrixAlgebraModule, ONLY : IncrementMatrix, MatrixMultiply, & @@ -176,10 +177,22 @@ SUBROUTINE EstimateGap(H, K, chemical_potential, gap, solver_parameters_in) CALL PrintParameters(params) END IF + !! Use the power method to get the smallest eigenvalue. This works + !! because we project out the rydberg states, but if not we fall + !! back to the Gershgorin Bounds + CALL MatrixMultiply(K, H, KH, & + & threshold_in = params%threshold, memory_pool_in = pool) + CALL WriteElement("Estimate Minimum") + CALL EnterSubLog + CALL PowerBounds(KH, e_min, params) + CALL ExitSubLog + IF (e_min .GT. 0_NTREAL) THEN + CALL WriteComment("Switching to GershgorinBounds") + CALL GershgorinBounds(H, e_min, e_max) + END IF + CALL WriteElement("Estimated e_min", VALUE = e_min) + !! Shift the matrix so the HOMO is the largest magnitude eigenvalue - CALL GershgorinBounds(H, e_min, e_max) - CALL WriteElement("Gershgorin e_min", VALUE = e_min) - CALL WriteElement("Gershgorin e_max", VALUE = e_max) CALL ConstructEmptyMatrix(ShiftH, H) CALL FillMatrixIdentity(ShiftH) CALL ScaleMatrix(ShiftH, -e_min) @@ -194,6 +207,7 @@ SUBROUTINE EstimateGap(H, K, chemical_potential, gap, solver_parameters_in) e_max = e_max + e_min CALL WriteElement("HOMO Estimate", VALUE = e_max) gap = 2.0_NTREAL * (chemical_potential - e_max) + CALL WriteElement("Gap Estimate", VALUE = gap) !! Cleanup IF (params%be_verbose) THEN diff --git a/Source/Fortran/GeometryOptimizationModule.F90 b/Source/Fortran/GeometryOptimizationModule.F90 index 02cfdaea..6803e405 100644 --- a/Source/Fortran/GeometryOptimizationModule.F90 +++ b/Source/Fortran/GeometryOptimizationModule.F90 @@ -1,6 +1,7 @@ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> A Module For Geometry Optimization MODULE GeometryOptimizationModule + USE ConvergenceMonitor, ONLY : ConstructMonitor, CheckConverged, AppendValue USE DataTypesModule, ONLY : NTREAL USE LoadBalancerModule, ONLY : PermuteMatrix, UndoPermuteMatrix USE LoggingModule, ONLY : EnterSubLog, ExitSubLog, WriteHeader, & @@ -56,6 +57,9 @@ SUBROUTINE PurificationExtrapolate(PreviousDensity, Overlap, trace, & ELSE CALL ConstructSolverParameters(params) END IF + CALL ConstructMonitor(params%monitor, & + & automatic_in = params%monitor_convergence, & + & tight_cutoff_in=params%converge_diff) IF (params%be_verbose) THEN CALL WriteHeader("Density Matrix Extrapolator") @@ -108,23 +112,21 @@ SUBROUTINE PurificationExtrapolate(PreviousDensity, Overlap, trace, & CALL IncrementMatrix(WorkingDensity, NewDensity, 2.0_NTREAL) END IF - !! Check Convergence + !! Convergence Measure CALL IncrementMatrix(NewDensity, WorkingDensity, -1.0_NTREAL) norm_value = MatrixNorm(WorkingDensity) + !! Xn = Xn+1 + CALL CopyMatrix(NewDensity, WorkingDensity) + + !! Check Convergence + CALL AppendValue(params%monitor, norm_value) + IF (CheckConverged(params%monitor, params%be_verbose)) EXIT IF (params%be_verbose) THEN - CALL WriteListElement(key = "Convergence", VALUE = norm_value) CALL EnterSubLog CALL WriteElement(key = "Trace", VALUE = trace_value) CALL ExitSubLog END IF - - IF (norm_value .LE. params%converge_diff) THEN - EXIT - END IF - - !! Xn = Xn+1 - CALL CopyMatrix(NewDensity, WorkingDensity) END DO IF (params%be_verbose) THEN CALL ExitSubLog diff --git a/Source/Fortran/InverseSolversModule.F90 b/Source/Fortran/InverseSolversModule.F90 index 702b47e7..45c1e0a8 100644 --- a/Source/Fortran/InverseSolversModule.F90 +++ b/Source/Fortran/InverseSolversModule.F90 @@ -1,6 +1,7 @@ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> A Module For Computing The Inverse of a Matrix. MODULE InverseSolversModule + USE ConvergenceMonitor, ONLY : ConstructMonitor, CheckConverged, AppendValue USE DataTypesModule, ONLY : NTREAL USE EigenSolversModule, ONLY : DenseMatrixFunction USE LoadBalancerModule, ONLY : PermuteMatrix, UndoPermuteMatrix @@ -49,6 +50,9 @@ SUBROUTINE Invert(InputMat, OutputMat, solver_parameters_in) ELSE CALL ConstructSolverParameters(params) END IF + CALL ConstructMonitor(params%monitor, & + & automatic_in = params%monitor_convergence, & + & tight_cutoff_in=params%converge_diff) IF (params%be_verbose) THEN CALL WriteHeader("Inverse Solver") @@ -117,9 +121,8 @@ SUBROUTINE Invert(InputMat, OutputMat, solver_parameters_in) CALL IncrementMatrix(Temp2, OutputMat, & & threshold_in = params%threshold) - IF (norm_value .LE. params%converge_diff) THEN - EXIT - END IF + CALL AppendValue(params%monitor, norm_value) + IF (CheckConverged(params%monitor, params%be_verbose)) EXIT END DO IF (params%be_verbose) THEN CALL ExitSubLog @@ -205,6 +208,9 @@ SUBROUTINE PseudoInverse(InputMat, OutputMat, solver_parameters_in) ELSE CALL ConstructSolverParameters(params) END IF + CALL ConstructMonitor(params%monitor, & + & automatic_in = params%monitor_convergence, & + & tight_cutoff_in=params%converge_diff) IF (params%be_verbose) THEN CALL WriteHeader("Inverse Solver") @@ -249,31 +255,23 @@ SUBROUTINE PseudoInverse(InputMat, OutputMat, solver_parameters_in) II = 1 norm_value = params%converge_diff + 1.0_NTREAL DO II = 1,params%max_iterations - IF (params%be_verbose .AND. II .GT. 1) THEN - CALL WriteListElement(key = "Convergence", VALUE = norm_value) - END IF - CALL MatrixMultiply(OutputMat, BalancedMat, Temp1, & & threshold_in = params%threshold, memory_pool_in = pool) - CALL MatrixMultiply(Temp1, OutputMat, Temp2, alpha_in = -1.0_NTREAL, & - & threshold_in = params%threshold, memory_pool_in = pool) - !! Save a copy of the last inverse matrix - CALL CopyMatrix(OutputMat, Temp1) + !! Convergence measures + CALL CopyMatrix(Identity, Temp2) + CALL IncrementMatrix(Temp1, Temp2, -1.0_NTREAL) + norm_value = MatrixNorm(Temp2) + CALL MatrixMultiply(Temp1, OutputMat, Temp2, alpha_in = -1.0_NTREAL, & + & threshold_in = params%threshold, memory_pool_in = pool) CALL ScaleMatrix(OutputMat, 2.0_NTREAL) CALL IncrementMatrix(Temp2, OutputMat, & & threshold_in = params%threshold) !! Check if Converged - CALL IncrementMatrix(OutputMat, Temp1, -1.0_NTREAL) - norm_value = MatrixNorm(Temp1) - - !! Sometimes the first few values don't change so much, so that's why - !! I added the outer counter check - IF (norm_value .LE. params%converge_diff .AND. II .GT. 3) THEN - EXIT - END IF + CALL AppendValue(params%monitor, norm_value) + IF (CheckConverged(params%monitor, params%be_verbose)) EXIT END DO IF (params%be_verbose) THEN CALL ExitSubLog diff --git a/Source/Fortran/LinearSolversModule.F90 b/Source/Fortran/LinearSolversModule.F90 index 314fc4b1..cfbcf459 100644 --- a/Source/Fortran/LinearSolversModule.F90 +++ b/Source/Fortran/LinearSolversModule.F90 @@ -4,6 +4,7 @@ MODULE LinearSolversModule USE CholeskyModule, ONLY : ConstructRankLookup, AppendToVector, & & BroadcastVector, ConstructDiag, DotAllHelper, DotAllPivoted, & & GatherMatrixColumn, GetPivot, UnpackCholesky + USE ConvergenceMonitor, ONLY : ConstructMonitor, CheckConverged, AppendValue USE DataTypesModule, ONLY : NTREAL, MPINTREAL USE DMatrixModule, ONLY : Matrix_ldr, DestructMatrix, & & ConstructMatrixDFromS @@ -59,6 +60,9 @@ SUBROUTINE CGSolver(AMat, XMat, BMat, solver_parameters_in) ELSE CALL ConstructSolverParameters(params) END IF + CALL ConstructMonitor(params%monitor, & + & automatic_in = params%monitor_convergence, & + & tight_cutoff_in=params%converge_diff) !! Print out parameters IF (params%be_verbose) THEN @@ -107,13 +111,6 @@ SUBROUTINE CGSolver(AMat, XMat, BMat, solver_parameters_in) END IF norm_value = params%converge_diff + 1.0_NTREAL DO II = 1, params%max_iterations - IF (params%be_verbose .AND. II .GT. 1) THEN - CALL WriteListElement(key = "Convergence", VALUE = norm_value) - END IF - IF (norm_value .LE. params%converge_diff) THEN - EXIT - END IF - !! Compute the Step Size CALL MatrixMultiply(ABalanced, PMat, QMat, & & threshold_in = params%threshold, memory_pool_in = pool) @@ -151,6 +148,10 @@ SUBROUTINE CGSolver(AMat, XMat, BMat, solver_parameters_in) CALL ScaleMatrix(PMat, step_size) CALL IncrementMatrix(RMat, PMat) + !! Check Exit Condition + CALL AppendValue(params%monitor, norm_value) + IF (CheckConverged(params%monitor, params%be_verbose)) EXIT + END DO IF (params%be_verbose) THEN CALL ExitSubLog diff --git a/Source/Fortran/RootSolversModule.F90 b/Source/Fortran/RootSolversModule.F90 index 85731796..918bc935 100644 --- a/Source/Fortran/RootSolversModule.F90 +++ b/Source/Fortran/RootSolversModule.F90 @@ -1,6 +1,7 @@ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> A Module For Computing General Matrix Roots. MODULE RootSolversModule + USE ConvergenceMonitor, ONLY : ConstructMonitor, CheckConverged, AppendValue USE DataTypesModule, ONLY : NTREAL USE EigenBoundsModule, ONLY : GershgorinBounds USE InverseSolversModule, ONLY : Invert @@ -181,7 +182,7 @@ SUBROUTINE ComputeInverseRootImplemention(InputMat, OutputMat, root, params) !> Which inverse root to compute. INTEGER, INTENT(IN) :: root !> Parameters for the solver. - TYPE(SolverParameters_t), INTENT(IN) :: params + TYPE(SolverParameters_t), INTENT(INOUT) :: params !! Local Matrices TYPE(Matrix_ps) :: SqrtMat, FthrtMat TYPE(Matrix_ps) :: IdentityMat @@ -200,6 +201,11 @@ SUBROUTINE ComputeInverseRootImplemention(InputMat, OutputMat, root, params) INTEGER :: JJ TYPE(MatrixMemoryPool_p) :: pool + !! Setup the monitor + CALL ConstructMonitor(params%monitor, & + & automatic_in = params%monitor_convergence, & + & tight_cutoff_in=params%converge_diff) + IF (params%be_verbose) THEN CALL WriteHeader("Root Solver") CALL EnterSubLog @@ -291,9 +297,9 @@ SUBROUTINE ComputeInverseRootImplemention(InputMat, OutputMat, root, params) & alpha_in=-1.0_NTREAL) norm_value = MatrixNorm(Temp) - IF (norm_value .LE. params%converge_diff) THEN - EXIT - END IF + !! Check Exit Condition + CALL AppendValue(params%monitor, norm_value) + IF (CheckConverged(params%monitor, params%be_verbose)) EXIT END DO IF (params%be_verbose) THEN CALL ExitSubLog diff --git a/Source/Fortran/SignSolversModule.F90 b/Source/Fortran/SignSolversModule.F90 index 5eaf44f1..294f8919 100644 --- a/Source/Fortran/SignSolversModule.F90 +++ b/Source/Fortran/SignSolversModule.F90 @@ -1,6 +1,7 @@ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> A Module For Computing The Matrix Sign Function. MODULE SignSolversModule + USE ConvergenceMonitor, ONLY : ConstructMonitor, CheckConverged, AppendValue USE DataTypesModule, ONLY : NTREAL USE EigenBoundsModule, ONLY : GershgorinBounds USE EigenSolversModule, ONLY : DenseMatrixFunction @@ -152,7 +153,7 @@ SUBROUTINE CoreComputation(InMat, OutMat, params, needs_transpose) !> Output of the routine. TYPE(Matrix_ps), INTENT(INOUT) :: OutMat !> Parameters for the solver. - TYPE(SolverParameters_t), INTENT(IN) :: params + TYPE(SolverParameters_t), INTENT(INOUT) :: params !> Whether we need to perform transposes in this routine (for polar). LOGICAL, INTENT(IN) :: needs_transpose !! Local Matrices @@ -169,6 +170,11 @@ SUBROUTINE CoreComputation(InMat, OutMat, params, needs_transpose) REAL(NTREAL) :: norm_value INTEGER :: II + !! Setup the monitor + CALL ConstructMonitor(params%monitor, & + & automatic_in = params%monitor_convergence, & + & tight_cutoff_in=params%converge_diff) + !! Construct All The Necessary Matrices CALL ConstructEmptyMatrix(Identity, InMat) CALL ConstructEmptyMatrix(Temp1, InMat) @@ -199,9 +205,6 @@ SUBROUTINE CoreComputation(InMat, OutMat, params, needs_transpose) II = 1 norm_value = params%converge_diff + 1.0_NTREAL iterate: DO II = 1, params%max_iterations - IF (params%be_verbose .AND. II .GT. 1) THEN - CALL WriteListElement(key = "Convergence", VALUE = norm_value) - END IF !! Update Scaling Factors alpha_k = MIN(SQRT(3.0_NTREAL / (1.0_NTREAL + xk + xk**2)), alpha) @@ -230,9 +233,10 @@ SUBROUTINE CoreComputation(InMat, OutMat, params, needs_transpose) norm_value = MatrixNorm(OutMat) CALL CopyMatrix(Temp2, OutMat) - IF (norm_value .LE. params%converge_diff) THEN - EXIT - END IF + !! Check Exit Condition + CALL AppendValue(params%monitor, norm_value) + IF (CheckConverged(params%monitor, params%be_verbose)) EXIT + END DO iterate IF (params%be_verbose) THEN CALL ExitSubLog diff --git a/Source/Fortran/SolverParametersModule.F90 b/Source/Fortran/SolverParametersModule.F90 index 35901df0..ae089445 100644 --- a/Source/Fortran/SolverParametersModule.F90 +++ b/Source/Fortran/SolverParametersModule.F90 @@ -1,6 +1,7 @@ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> A Module For Storing The Parameters For Iterative Solvers. MODULE SolverParametersModule + USE ConvergenceMonitor, ONLY : Monitor_t, DestructMonitor USE DataTypesModule, ONLY : NTREAL USE LoggingModule, ONLY : EnterSubLog, ExitSubLog, WriteElement, & & WriteHeader @@ -25,6 +26,10 @@ MODULE SolverParametersModule TYPE(Permutation_t) :: BalancePermutation !> Thresholds for step size searches. REAL(NTREAL) :: step_thresh + !> Whether to do an automatic convergence detection. + LOGICAL :: monitor_convergence + !> The convergence monitor + TYPE(Monitor_t) :: monitor END TYPE SolverParameters_t !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! PUBLIC :: ConstructSolverParameters @@ -35,6 +40,7 @@ MODULE SolverParametersModule PUBLIC :: SetParametersBeVerbose PUBLIC :: SetParametersLoadBalance PUBLIC :: SetParametersStepThreshold + PUBLIC :: SetParametersMonitorConvergence PUBLIC :: PrintParameters PUBLIC :: DestructSolverParameters !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -46,7 +52,7 @@ MODULE SolverParametersModule !> Construct a data type which stores iterative solver parameters. SUBROUTINE ConstructSolverParameters(this, converge_diff_in, threshold_in, & & max_iterations_in, be_verbose_in, BalancePermutation_in, & - & step_thresh_in) + & step_thresh_in, monitor_convergence_in) !> The parameters to construct. TYPE(SolverParameters_t), INTENT(INOUT) :: this !> Converge_diff_in the difference between iterations to consider @@ -62,6 +68,8 @@ SUBROUTINE ConstructSolverParameters(this, converge_diff_in, threshold_in, & TYPE(Permutation_t), INTENT(IN), OPTIONAL :: BalancePermutation_in !> Step size for differential equation solvers. REAL(NTREAL), INTENT(IN), OPTIONAL :: step_thresh_in + !> Whether to do automatic convergence monitoring (default = True). + LOGICAL, INTENT(IN), OPTIONAL :: monitor_convergence_in CALL DestructSolverParameters(this) @@ -97,6 +105,11 @@ SUBROUTINE ConstructSolverParameters(this, converge_diff_in, threshold_in, & ELSE this%step_thresh = step_thresh_in END IF + IF (.NOT. PRESENT(monitor_convergence_in)) THEN + this%monitor_convergence = .TRUE. + ELSE + this%monitor_convergence = monitor_convergence_in + END IF END SUBROUTINE ConstructSolverParameters !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE CopySolverParameters(paramA, paramB) @@ -110,7 +123,7 @@ SUBROUTINE CopySolverParameters(paramA, paramB) END SUBROUTINE CopySolverParameters !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Set the value of the convergence difference. - PURE SUBROUTINE SetParametersConvergeDiff(this,new_value) + PURE SUBROUTINE SetParametersConvergeDiff(this, new_value) !> The parameter object. TYPE(SolverParameters_t), INTENT(INOUT) :: this !> Value to set it to. @@ -120,7 +133,7 @@ PURE SUBROUTINE SetParametersConvergeDiff(this,new_value) END SUBROUTINE SetParametersConvergeDiff !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Set the value of the max iterations. - PURE SUBROUTINE SetParametersMaxIterations(this,new_value) + PURE SUBROUTINE SetParametersMaxIterations(this, new_value) !> The parameter object. TYPE(SolverParameters_t), INTENT(INOUT) :: this !> Value to set it to. @@ -130,7 +143,7 @@ PURE SUBROUTINE SetParametersMaxIterations(this,new_value) END SUBROUTINE SetParametersMaxIterations !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Set the value of the threshold. - PURE SUBROUTINE SetParametersThreshold(this,new_value) + PURE SUBROUTINE SetParametersThreshold(this, new_value) !> The parameter object. TYPE(SolverParameters_t), INTENT(INOUT) :: this !> Value to set it to. @@ -140,7 +153,7 @@ PURE SUBROUTINE SetParametersThreshold(this,new_value) END SUBROUTINE SetParametersThreshold !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Set the value of the verbosity. - PURE SUBROUTINE SetParametersBeVerbose(this,new_value) + PURE SUBROUTINE SetParametersBeVerbose(this, new_value) !> The parameter object. TYPE(SolverParameters_t), INTENT(INOUT) :: this !> Value to set it to. @@ -150,7 +163,7 @@ PURE SUBROUTINE SetParametersBeVerbose(this,new_value) END SUBROUTINE SetParametersBeVerbose !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Set the value of the load balance. - PURE SUBROUTINE SetParametersLoadBalance(this,new_value) + PURE SUBROUTINE SetParametersLoadBalance(this, new_value) !> The parameter object. TYPE(SolverParameters_t), INTENT(INOUT) :: this !> Value to set it to. @@ -161,7 +174,7 @@ PURE SUBROUTINE SetParametersLoadBalance(this,new_value) END SUBROUTINE SetParametersLoadBalance !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Set the value of the step threshold. - PURE SUBROUTINE SetParametersStepThreshold(this,new_value) + PURE SUBROUTINE SetParametersStepThreshold(this, new_value) !> The parameter object. TYPE(SolverParameters_t), INTENT(INOUT) :: this !> Value to set it to. @@ -169,6 +182,16 @@ PURE SUBROUTINE SetParametersStepThreshold(this,new_value) this%step_thresh = new_value END SUBROUTINE SetParametersStepThreshold +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Set the value of the monitoring of convergence + PURE SUBROUTINE SetParametersMonitorConvergence(this, new_value) + !> The parameter object. + TYPE(SolverParameters_t), INTENT(INOUT) :: this + !> Value to set it to. + LOGICAL, INTENT(IN) :: new_value + + this%monitor_convergence = new_value + END SUBROUTINE SetParametersMonitorConvergence !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Print out the iterative solver parameter values. SUBROUTINE PrintParameters(this) @@ -189,6 +212,8 @@ SUBROUTINE PrintParameters(this) & VALUE = this%max_iterations) CALL WriteElement(key = "Step Threshold", & & VALUE = this%step_thresh) + CALL WriteElement(key = "Monitor Convergence", & + & VALUE = this%monitor_convergence) CALL ExitSubLog END SUBROUTINE PrintParameters !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -198,6 +223,7 @@ PURE SUBROUTINE DestructSolverParameters(this) TYPE(SolverParameters_t), INTENT(INOUT) :: this CALL DestructPermutation(this%BalancePermutation) + CALL DestructMonitor(this%monitor) END SUBROUTINE DestructSolverParameters !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! END MODULE SolverParametersModule diff --git a/Source/Fortran/SquareRootSolversModule.F90 b/Source/Fortran/SquareRootSolversModule.F90 index b4c68c45..82c85472 100644 --- a/Source/Fortran/SquareRootSolversModule.F90 +++ b/Source/Fortran/SquareRootSolversModule.F90 @@ -1,6 +1,7 @@ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> A Module For Computing The Square Root of a Matrix. MODULE SquareRootSolversModule + USE ConvergenceMonitor, ONLY : ConstructMonitor, CheckConverged, AppendValue USE DataTypesModule, ONLY : NTREAL USE EigenBoundsModule, ONLY : GershgorinBounds USE EigenSolversModule, ONLY : DenseMatrixFunction @@ -167,7 +168,7 @@ SUBROUTINE SquareRootSelector(InputMat, OutputMat, params, & !> The Matrix computed. TYPE(Matrix_ps), INTENT(INOUT) :: OutputMat !> Parameters about how to solve. - TYPE(SolverParameters_t),INTENT(IN) :: params + TYPE(SolverParameters_t),INTENT(INOUT) :: params !> True if we are computing the inverse square root. LOGICAL, INTENT(IN) :: compute_inverse !> The polynomial degree to use (optional, default = 5) @@ -200,7 +201,7 @@ SUBROUTINE NewtonSchultzISROrder2(InMat, OutMat, params, compute_inverse) !> Mat^-1/2 or Mat^1/2. TYPE(Matrix_ps), INTENT(INOUT) :: OutMat !> Parameters for the solver - TYPE(SolverParameters_t), INTENT(IN) :: params + TYPE(SolverParameters_t), INTENT(INOUT) :: params !> Whether to compute the inverse square root. LOGICAL, INTENT(IN) :: compute_inverse !! Local Variables @@ -215,6 +216,11 @@ SUBROUTINE NewtonSchultzISROrder2(InMat, OutMat, params, compute_inverse) REAL(NTREAL) :: norm_value TYPE(MatrixMemoryPool_p) :: mpool + !! Setup the monitor + CALL ConstructMonitor(params%monitor, & + & automatic_in = params%monitor_convergence, & + & tight_cutoff_in=params%converge_diff) + IF (params%be_verbose) THEN CALL WriteHeader("Newton Schultz Inverse Square Root") CALL EnterSubLog @@ -295,13 +301,10 @@ SUBROUTINE NewtonSchultzISROrder2(InMat, OutMat, params, compute_inverse) & threshold_in = params%threshold, memory_pool_in = mpool) CALL ScaleMatrix(SquareRootMat, SQRT(lambda)) - IF (params%be_verbose) THEN - CALL WriteElement(key = "Convergence", VALUE = norm_value) - END IF + !! Check Exit Condition + CALL AppendValue(params%monitor, norm_value) + IF (CheckConverged(params%monitor, params%be_verbose)) EXIT - IF (norm_value .LE. params%converge_diff) THEN - EXIT - END IF END DO IF (params%be_verbose) THEN CALL ExitSubLog @@ -343,7 +346,7 @@ SUBROUTINE NewtonSchultzISRTaylor(InMat, OutMat, params, & !> Mat^-1/2 or Mat^1/2. TYPE(Matrix_ps), INTENT(INOUT) :: OutMat !> Parameters for the solver. - TYPE(SolverParameters_t), INTENT(IN) :: params + TYPE(SolverParameters_t), INTENT(INOUT) :: params !> Order of polynomial to use. INTEGER, INTENT(IN) :: taylor_order !> Whether to compute the inverse square root or not. @@ -362,6 +365,11 @@ SUBROUTINE NewtonSchultzISRTaylor(InMat, OutMat, params, & REAL(NTREAL) :: norm_value TYPE(MatrixMemoryPool_p) :: mpool + !! Setup the monitor + CALL ConstructMonitor(params%monitor, & + & automatic_in = params%monitor_convergence, & + & tight_cutoff_in=params%converge_diff) + IF (params%be_verbose) THEN CALL WriteHeader("Newton Schultz Inverse Square Root") CALL EnterSubLog @@ -481,13 +489,10 @@ SUBROUTINE NewtonSchultzISRTaylor(InMat, OutMat, params, & CALL MatrixMultiply(Temp, X_k, SquareRootMat, & & threshold_in = params%threshold, memory_pool_in = mpool) - IF (params%be_verbose) THEN - CALL WriteListElement(key = "Convergence", VALUE = norm_value) - END IF + !! Check Exit Condition + CALL AppendValue(params%monitor, norm_value) + IF (CheckConverged(params%monitor, params%be_verbose)) EXIT - IF (norm_value .LE. params%converge_diff) THEN - EXIT - END IF END DO IF (params%be_verbose) THEN CALL ExitSubLog diff --git a/Source/Wrapper/SolverParametersModule_wrp.F90 b/Source/Wrapper/SolverParametersModule_wrp.F90 index fbd34dd0..dba41f45 100644 --- a/Source/Wrapper/SolverParametersModule_wrp.F90 +++ b/Source/Wrapper/SolverParametersModule_wrp.F90 @@ -23,6 +23,7 @@ MODULE SolverParametersModule_wrp PUBLIC :: SetParametersBeVerbose_wrp PUBLIC :: SetParametersLoadBalance_wrp PUBLIC :: SetParametersStepThreshold_wrp + PUBLIC :: SetParametersMonitorConvergence_wrp CONTAINS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Construct the iterat solver parameters. SUBROUTINE ConstructSolverParameters_wrp(ih_this) & @@ -47,7 +48,7 @@ SUBROUTINE DestructSolverParameters_wrp(ih_this) & END SUBROUTINE DestructSolverParameters_wrp !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Set the value of the convergence difference. - SUBROUTINE SetParametersConvergeDiff_wrp(ih_this,new_value) & + SUBROUTINE SetParametersConvergeDiff_wrp(ih_this, new_value) & & BIND(c,name="SetParametersConvergeDiff_wrp") INTEGER(kind=c_int), INTENT(INOUT) :: ih_this(SIZE_wrp) REAL(NTREAL), INTENT(IN) :: new_value @@ -58,7 +59,7 @@ SUBROUTINE SetParametersConvergeDiff_wrp(ih_this,new_value) & END SUBROUTINE SetParametersConvergeDiff_wrp !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Set the value of the max iterations. - SUBROUTINE SetParametersMaxIterations_wrp(ih_this,new_value) & + SUBROUTINE SetParametersMaxIterations_wrp(ih_this, new_value) & & BIND(c,name="SetParametersMaxIterations_wrp") INTEGER(kind=c_int), INTENT(INOUT) :: ih_this(SIZE_wrp) INTEGER(kind=c_int), INTENT(IN) :: new_value @@ -69,7 +70,7 @@ SUBROUTINE SetParametersMaxIterations_wrp(ih_this,new_value) & END SUBROUTINE SetParametersMaxIterations_wrp !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Set the value of the threshold. - SUBROUTINE SetParametersThreshold_wrp(ih_this,new_value) & + SUBROUTINE SetParametersThreshold_wrp(ih_this, new_value) & & BIND(c,name="SetParametersThreshold_wrp") INTEGER(kind=c_int), INTENT(INOUT) :: ih_this(SIZE_wrp) REAL(NTREAL), INTENT(IN) :: new_value @@ -80,7 +81,7 @@ SUBROUTINE SetParametersThreshold_wrp(ih_this,new_value) & END SUBROUTINE SetParametersThreshold_wrp !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Set the value of the verbosity. - SUBROUTINE SetParametersBeVerbose_wrp(ih_this,new_value) & + SUBROUTINE SetParametersBeVerbose_wrp(ih_this, new_value) & & BIND(c,name="SetParametersBeVerbose_wrp") INTEGER(kind=c_int), INTENT(INOUT) :: ih_this(SIZE_wrp) LOGICAL(kind=c_bool), INTENT(IN) :: new_value @@ -104,14 +105,25 @@ SUBROUTINE SetParametersLoadBalance_wrp(ih_this,ih_new_value) & END SUBROUTINE SetParametersLoadBalance_wrp !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Set the value of the step threshold. - SUBROUTINE SetParametersStepThreshold_wrp(ih_this,new_value) & + SUBROUTINE SetParametersStepThreshold_wrp(ih_this, new_value) & & BIND(c,name="SetParametersStepThreshold_wrp") INTEGER(kind=c_int), INTENT(INOUT) :: ih_this(SIZE_wrp) REAL(NTREAL), INTENT(IN) :: new_value TYPE(SolverParameters_wrp) :: h_this h_this = TRANSFER(ih_this,h_this) - CALL SetParametersStepThreshold(h_this%DATA,new_value) + CALL SetParametersStepThreshold(h_this%DATA, new_value) END SUBROUTINE SetParametersStepThreshold_wrp +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Set the value of the convergence monitoring. + SUBROUTINE SetParametersMonitorConvergence_wrp(ih_this, new_value) & + & BIND(c,name="SetParametersMonitorConvergence_wrp") + INTEGER(kind=c_int), INTENT(INOUT) :: ih_this(SIZE_wrp) + LOGICAL(kind=c_bool), INTENT(IN) :: new_value + TYPE(SolverParameters_wrp) :: h_this + + h_this = TRANSFER(ih_this, h_this) + CALL SetParametersMonitorConvergence(h_this%DATA, LOGICAL(new_value)) + END SUBROUTINE SetParametersMonitorConvergence_wrp !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! END MODULE SolverParametersModule_wrp