Optimizing solvers

examples/EMSchur3D-vtk.cpp

@@ -116,6 +116,7 @@ int main( int argc, char** argv ) {
     //auto EM3D = EMSchur3D< Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::IncompleteLUT<Complex> > >::NewSP();
     //auto EM3D = EMSchur3D< Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, Eigen::RowMajor> > >::NewSP();
     auto EM3D = EMSchur3D< Eigen::LeastSquaresConjugateGradient<Eigen::SparseMatrix<Complex, Eigen::RowMajor> > >::NewSP();
+    //auto EM3D = EMSchur3D< Eigen::ConjugateGradient<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower > >::NewSP();
 
 
     EM3D->SetRectilinearGrid( gimp->GetGrid() );

include/EMSchur3D.h

@@ -215,7 +215,7 @@ namespace Lemma {
 
         /////////////////////////////////////////
         // Solve for RHS
-        CSolver[iw].setMaxIterations(10000);
+        //CSolver[iw].setMaxIterations(750);
         VectorXcr A = CSolver[iw].solve(Se);
 
 //         // Solve Real system instead
@@ -236,8 +236,9 @@ namespace Lemma {
         logio << "|| Div(A) || = " << ADiv.norm()
                 // << " in " << iter_done << " iterations"
               //<<  " with error " << errorn << "\t"
-              << "\tInital solution error "<<   Error.norm()  // Iteritive info
-              << "\ttime " << timer.end() / 60. << " [m]" << std::endl;
+              << "\tInital solution error="<<   Error.norm()  // Iteritive info
+              << "\tSolver reported error="<<   CSolver[iw].error()  // Iteritive info
+              << "\ttime " << timer.end() / 60. << " [m]   " <<  CSolver[iw].iterations() << "  iterations" << std::endl;
 
         //VectorXcr ADivMAC = ADiv.array() * MAC.array().cast<Complex>();
         //logio << "|| Div(A) || on MAC grid " << ADivMAC.norm() << std::endl;
@@ -499,8 +500,8 @@ namespace Lemma {
         CSolver = new Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::IncompleteLUT<Complex> > [Omegas.size()];
         for (int iw=0; iw<Omegas.size(); ++iw) {
             Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
-            CSolver[iw].preconditioner().setDroptol(1e-12);
-            CSolver[iw].preconditioner().setFillfactor(1e2);
+            CSolver[iw].preconditioner().setDroptol(1e-12);      // 1e-12
+            CSolver[iw].preconditioner().setFillfactor(1e1);     // 1e2
             jsw_timer timer;
             timer.begin();
             /*  Complex system */
@@ -559,26 +560,26 @@ namespace Lemma {
         }
     }
 
-//     template<>
-//     void EMSchur3D<   Eigen::ConjugateGradient<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower > > ::BuildCDirectSolver() {
-//         CSolver = new Eigen::ConjugateGradient<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower > [Omegas.size()];
-//         for (int iw=0; iw<Omegas.size(); ++iw) {
-//             //Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
-//             jsw_timer timer;
-//             timer.begin();
-//             /*  Complex system */
-//             std::cout << "ConjugateGradient pattern analyzing C_" << iw << ",";
-//             std::cout.flush();
-//             CSolver[iw].analyzePattern( Cvec[iw] );
-//             std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
-//             // factorize
-//             timer.begin();
-//             std::cout << "ConjugateGradient factorising C_" << iw << ", ";
-//             std::cout.flush();
-//             CSolver[iw].factorize( Cvec[iw] );
-//             std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
-//         }
-//     }
+    template<>
+    void EMSchur3D<   Eigen::ConjugateGradient<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower > > ::BuildCDirectSolver() {
+        CSolver = new Eigen::ConjugateGradient<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower > [Omegas.size()];
+        for (int iw=0; iw<Omegas.size(); ++iw) {
+            //Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
+            jsw_timer timer;
+            timer.begin();
+            /*  Complex system */
+            std::cout << "ConjugateGradient pattern analyzing C_" << iw << ",";
+            std::cout.flush();
+            CSolver[iw].analyzePattern( Cvec[iw] );
+            std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
+            // factorize
+            timer.begin();
+            std::cout << "ConjugateGradient factorising C_" << iw << ", ";
+            std::cout.flush();
+            CSolver[iw].factorize( Cvec[iw] );
+            std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
+        }
+    }
 
 //     template<>
 //     void EMSchur3D<   Eigen::PastixLLT<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower > > ::BuildCDirectSolver() {

include/bicgstab.h

@@ -473,8 +473,8 @@ inline VectorXcr implicitDCInvBPhi3 (const SparseMat& D, const Solver& solver,
                         int& max_it) {
     VectorXcr b = (ioms).asDiagonal() * (D.transpose()*Phi);
     VectorXcr y = solver.solve(b);
-    max_it = 0;
-    //max_it = solver.iterations();
+    //max_it = 0;
+    max_it = solver.iterations();
     //errorn = solver.error();
     return  D*y;
 }
@@ -639,7 +639,8 @@ int implicitbicgstab(//const SparseMat& D,
         logio << "iteration " << std::setw(3) << iter
               << "  errorn " << std::setw(6) << std::setprecision(4) << std::scientific << errorn
               //<< "\timplicit iterations " << std::setw(5) << max_it1+max_it2
-              << "  time " << std::setw(6) << std::fixed << std::setprecision(2) << timer.end() << std::endl;
+              << "  time " << std::setw(6) << std::fixed << std::setprecision(2) << timer.end()/60. << " [m]  "
+              << max_it2+max_it2 << " iterations " << std::endl;
 
         // Check to see how progress is going