Extrapolation of density matrix in ROM MVP solver

examples/PinnedH2O/job.ref

@@ -1,14 +1,14 @@
 #!/bin/tcsh
 #SBATCH -N 1
-#SBATCH -t 1:00:00
-#SBATCH -p pdebug
+#SBATCH -t 3:00:00
+#SBATCH -p pbatch
 
 date
 
 setenv OMP_NUM_THREADS 1
 #setenv KMP_DETERMINISTIC_REDUCTION 1
 
-set ncpus = 8
+set ncpus = 1
 set case = 2
 
 set maindir = /p/lustre2/cheung26/mgmol

examples/PinnedH2O/job.rom_3DOF

@@ -1,14 +1,14 @@
 #!/bin/tcsh
 #SBATCH -N 1
-#SBATCH -t 1:00:00
-#SBATCH -p pdebug
+#SBATCH -t 2:00:00
+#SBATCH -p pbatch
 
 date
 
 setenv OMP_NUM_THREADS 1
 #setenv KMP_DETERMINISTIC_REDUCTION 1
 
-set ncpus = 8
+set ncpus = 1
 set case = 2
 
 set maindir = /p/lustre2/cheung26/mgmol

examples/PinnedH2O/mgmol_rom_3DOF_test2.cfg

@@ -38,7 +38,7 @@ output_level=4
 [ROM]
 stage=online_pinned_H2O_3dof
 [ROM.offline]
-basis_file=/usr/workspace/nlrom/MGmol/PinnedH2O_3DOF/data_8/PinnedH2O_3DOF_orbitals_basis_2_2
+basis_file=/usr/workspace/nlrom/MGmol/PinnedH2O_3DOF/data_1/PinnedH2O_3DOF_orbitals_basis_2_2
 [ROM.basis]
 compare_md=false
 number_of_orbital_basis=34

src/Control.cc

@@ -112,7 +112,6 @@ Control::Control()
     out_restart_file_type             = -1;
     spread_radius                     = -1.;
     iprint_residual                   = -1;
-    override_restart                  = -1;
     dot_product_type                  = -1;
     spread_penalty_damping_           = -1;
     spread_penalty_alpha_             = -1;
@@ -386,7 +385,7 @@ void Control::sync(void)
         short_buffer[47] = out_restart_file_naming_strategy;
         short_buffer[48] = enforceVmass0;
         short_buffer[49] = dm_inner_steps;
-        short_buffer[50] = override_restart;
+        short_buffer[50] = -1;
         short_buffer[51] = fgmres_kim;
         short_buffer[52] = fgmres_maxits;
         short_buffer[53] = ilu_type;
@@ -599,41 +598,41 @@ void Control::sync(void)
     out_restart_file_naming_strategy = short_buffer[47];
     enforceVmass0                    = short_buffer[48];
     dm_inner_steps                   = short_buffer[49];
-    override_restart                 = short_buffer[50];
-    fgmres_kim                       = short_buffer[51];
-    fgmres_maxits                    = short_buffer[52];
-    ilu_type                         = short_buffer[53];
-    ilu_lof                          = short_buffer[54];
-    ilu_maxfil                       = short_buffer[55];
-    coloring_algo_                   = short_buffer[56];
-    diel_flag_                       = short_buffer[57];
-    poisson_pc_nu1                   = short_buffer[58];
-    poisson_pc_nu2                   = short_buffer[59];
-    poisson_pc_nlev                  = short_buffer[60];
-    system_charge_                   = short_buffer[61];
-    md_print_freq                    = short_buffer[62];
-    use_kernel_functions             = short_buffer[63];
-    ngpts_[0]                        = short_buffer[64];
-    ngpts_[1]                        = short_buffer[65];
-    ngpts_[2]                        = short_buffer[66];
-    computeCondGram_                 = short_buffer[67];
-    lrs_extrapolation                = short_buffer[68];
-    parallel_transport               = (bool)short_buffer[69];
-    with_spin_                       = (bool)short_buffer[70];
-    conv_criterion_                  = short_buffer[71];
-    load_balancing_max_iterations    = short_buffer[72];
-    load_balancing_modulo            = short_buffer[73];
-    write_clusters                   = short_buffer[74];
-    DM_solver_                       = short_buffer[75];
-    dm_algo_                         = short_buffer[80];
-    dm_approx_order                  = short_buffer[81];
-    dm_approx_ndigits                = short_buffer[82];
-    dm_approx_power_maxits           = short_buffer[83];
-    spread_penalty_type_             = short_buffer[84];
-    dm_use_old_                      = short_buffer[85];
-    max_electronic_steps_tight_      = short_buffer[86];
-    hartree_reset_                   = short_buffer[88];
-    MD_last_step_                    = short_buffer[89];
+    //...                 = short_buffer[50];
+    fgmres_kim                    = short_buffer[51];
+    fgmres_maxits                 = short_buffer[52];
+    ilu_type                      = short_buffer[53];
+    ilu_lof                       = short_buffer[54];
+    ilu_maxfil                    = short_buffer[55];
+    coloring_algo_                = short_buffer[56];
+    diel_flag_                    = short_buffer[57];
+    poisson_pc_nu1                = short_buffer[58];
+    poisson_pc_nu2                = short_buffer[59];
+    poisson_pc_nlev               = short_buffer[60];
+    system_charge_                = short_buffer[61];
+    md_print_freq                 = short_buffer[62];
+    use_kernel_functions          = short_buffer[63];
+    ngpts_[0]                     = short_buffer[64];
+    ngpts_[1]                     = short_buffer[65];
+    ngpts_[2]                     = short_buffer[66];
+    computeCondGram_              = short_buffer[67];
+    lrs_extrapolation             = short_buffer[68];
+    parallel_transport            = (bool)short_buffer[69];
+    with_spin_                    = (bool)short_buffer[70];
+    conv_criterion_               = short_buffer[71];
+    load_balancing_max_iterations = short_buffer[72];
+    load_balancing_modulo         = short_buffer[73];
+    write_clusters                = short_buffer[74];
+    DM_solver_                    = short_buffer[75];
+    dm_algo_                      = short_buffer[80];
+    dm_approx_order               = short_buffer[81];
+    dm_approx_ndigits             = short_buffer[82];
+    dm_approx_power_maxits        = short_buffer[83];
+    spread_penalty_type_          = short_buffer[84];
+    dm_use_old_                   = short_buffer[85];
+    max_electronic_steps_tight_   = short_buffer[86];
+    hartree_reset_                = short_buffer[88];
+    MD_last_step_                 = short_buffer[89];
     poisson_lap_type_ = static_cast<PoissonFDtype>(short_buffer[90]);
 
     numst    = int_buffer[0];
@@ -1787,7 +1786,6 @@ void Control::setOptions(const boost::program_options::variables_map& vm)
         // options not available in configure file
         lr_updates_type         = 0;
         precond_factor_computed = false;
-        override_restart        = 0;
         mix_pot                 = 1.;
         project_out_psd         = 0;
         multipole_order         = 1;

src/Control.h

@@ -562,7 +562,6 @@ class Control
     short out_restart_file_naming_strategy;
     short restart_file_type;
     short out_restart_file_type;
-    short override_restart;
 
     short verbose;
 

src/DensityMatrix.cc

@@ -22,12 +22,6 @@
 
 const double factor_kernel4dot = 10.;
 
-#define MGMOL_DENSITYMATRIX_FAIL(X)                                            \
-    {                                                                          \
-        std::cerr << "DensityMatrix failure:" << std::endl;                    \
-        std::cerr << "Error Message: " << X << std::endl;                      \
-    }
-
 #define MGMOL_DENSITYMATRIX_FAIL(X)                                            \
     {                                                                          \
         std::cerr << "DensityMatrix failure:" << std::endl;                    \

src/MGmol.cc

@@ -45,6 +45,7 @@
 #include "MGmol.h"
 #include "MLWFTransform.h"
 #include "MPIdata.h"
+#include "MVPSolver.h"
 #include "MasksSet.h"
 #include "Mesh.h"
 #include "OrbitalsPreconditioning.h"
@@ -901,6 +902,9 @@ void MGmol<OrbitalsType>::printTimers()
     AndersonMix<OrbitalsType>::update_tm().print(os_);
     proj_matrices_->printTimers(os_);
     ShortSightedInverse::printTimers(os_);
+    if (std::is_same<OrbitalsType, ExtendedGridOrbitals>::value)
+        MVPSolver<ExtendedGridOrbitals,
+            dist_matrix::DistMatrix<DISTMATDTYPE>>::printTimers(os_);
     VariableSizeMatrixInterface::printTimers(os_);
     DataDistribution::printTimers(os_);
     PackedCommunicationBuffer::printTimers(os_);

src/MGmol.h

@@ -229,7 +229,6 @@ class MGmol : public MGmolInterface
 
     void printEigAndOcc();
 
-    int readCoordinates(std::ifstream* tfile, const bool cell_relative);
     int readCoordinates(const std::string& filename, const bool cell_relative);
     double computeConstraintResidual(OrbitalsType& orbitals,
         const OrbitalsType& hphi, OrbitalsType& res, const bool print_residual,

src/MVPSolver.cc

@@ -398,12 +398,8 @@ int MVPSolver<OrbitalsType, MatrixType>::solve(OrbitalsType& orbitals)
 template <class OrbitalsType, class MatrixType>
 void MVPSolver<OrbitalsType, MatrixType>::printTimers(std::ostream& os)
 {
-    if (onpe0)
-    {
-        os << std::setprecision(2) << std::fixed << std::endl;
-        solve_tm_.print(os);
-        target_tm_.print(os);
-    }
+    solve_tm_.print(os);
+    target_tm_.print(os);
 }
 
 template class MVPSolver<LocGridOrbitals,

src/MVPSolver.h

@@ -16,7 +16,6 @@
 
 class Ions;
 class Electrostatic;
-class ProjectedMatrices2N;
 template <class OrbitalsType>
 class ProjectedMatrices;
 
@@ -67,7 +66,7 @@ class MVPSolver
     ~MVPSolver();
 
     int solve(OrbitalsType& orbitals);
-    void printTimers(std::ostream& os);
+    static void printTimers(std::ostream& os);
 };
 
 #endif

src/main.cc

@@ -58,7 +58,7 @@ int main(int argc, char** argv)
     /*
      * read runtime parameters
      */
-    std::string input_filename("");
+    std::string coords_filename("");
     std::string lrs_filename;
     std::string constraints_filename("");
 
@@ -70,7 +70,7 @@ int main(int argc, char** argv)
     // read from PE0 only
     if (MPIdata::onpe0)
     {
-        read_config(argc, argv, vm, input_filename, lrs_filename,
+        read_config(argc, argv, vm, coords_filename, lrs_filename,
             constraints_filename, total_spin, with_spin);
     }
 
@@ -89,18 +89,18 @@ int main(int argc, char** argv)
     int ret = ct.checkOptions();
     if (ret < 0) return ret;
 
-    mmpi.bcastGlobal(input_filename);
+    mmpi.bcastGlobal(coords_filename);
     mmpi.bcastGlobal(lrs_filename);
 
     // Enter main scope
     {
         MGmolInterface* mgmol;
         if (ct.isLocMode())
             mgmol = new MGmol<LocGridOrbitals>(global_comm, *MPIdata::sout,
-                input_filename, lrs_filename, constraints_filename);
+                coords_filename, lrs_filename, constraints_filename);
         else
             mgmol = new MGmol<ExtendedGridOrbitals>(global_comm, *MPIdata::sout,
-                input_filename, lrs_filename, constraints_filename);
+                coords_filename, lrs_filename, constraints_filename);
 
         mgmol->setup();
 

src/md.cc

@@ -340,7 +340,7 @@ void MGmol<OrbitalsType>::md(OrbitalsType** orbitals, Ions& ions)
 
     constraints_->printConstraints(os_);
 
-    if (ct.restart_info > 0 && !ct.override_restart)
+    if (ct.restart_info > 0)
     {
         if (onpe0) os_ << "Use restart file to initialize MD..." << std::endl;
         stepper->init(*h5f_file_);
@@ -406,9 +406,13 @@ void MGmol<OrbitalsType>::md(OrbitalsType** orbitals, Ions& ions)
         h5f_file_.reset();
     }
 
-    bool ROM_MVP = (ct.getROMOptions().rom_stage == ROMStage::ONLINE_PINNED_H2O_3DOF);
 #ifdef MGMOL_HAS_LIBROM
+    bool ROM_MVP = (ct.getROMOptions().rom_stage == ROMStage::ONLINE_PINNED_H2O_3DOF);
+    DensityMatrix<dist_matrix::DistMatrix<DISTMATDTYPE>>* dm = nullptr;
+    dist_matrix::DistMatrix<DISTMATDTYPE>* previous_dm = nullptr;
     // ROM - initialize orbitals and density matrix
+    // assume ct.Mehrstellen() nor ct.short_sighted nor use_replicated_matrix
+    // so that Line 267 of MGmol.cc is activated 
     if (ROM_MVP)
     {
         if (onpe0) os_ << "Setup ROM MVP solver..." << std::endl;
@@ -422,6 +426,11 @@ void MGmol<OrbitalsType>::md(OrbitalsType** orbitals, Ions& ions)
             = getProjectedMatrices();
         projmatrices->setDMuniform(ct.getNelSpin());
         projmatrices->printDM(os_);
+
+        std::shared_ptr<ProjectedMatrices<dist_matrix::DistMatrix<DISTMATDTYPE>>> projmatrices_downcast =
+            std::dynamic_pointer_cast<ProjectedMatrices<dist_matrix::DistMatrix<DISTMATDTYPE>>>(projmatrices);
+        dm = &(projmatrices_downcast->getDM());
+        previous_dm = new dist_matrix::DistMatrix<DISTMATDTYPE>(dm->getMatrix());
     }
 #endif
 
@@ -487,6 +496,9 @@ void MGmol<OrbitalsType>::md(OrbitalsType** orbitals, Ions& ions)
         if (ROM_MVP)
         {
             md_updateDMandEnergy_tm.start();
+            dist_matrix::DistMatrix<DISTMATDTYPE> tmp(dm->getMatrix());
+            dm->linearExtrapolate(*previous_dm);
+            *previous_dm = tmp;
             updateDMandEnergy(**orbitals, *ROM_ions, eks);
             md_updateDMandEnergy_tm.stop();
         }