SparseLU

examples/EMSchur3D-vtk.cpp

 
     //////////////////////////////////////////////////
     // And solve
+
     // Use BiCGSTAB Diagonal preconditioner
     //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();
+
+    // LS CG
+    //auto EM3D = EMSchur3D< Eigen::LeastSquaresConjugateGradient<Eigen::SparseMatrix<Complex, Eigen::RowMajor> > >::NewSP();
+
+    // CG...dangerous
     //auto EM3D = EMSchur3D< Eigen::ConjugateGradient<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower > >::NewSP();
+    //auto EM3D = EMSchur3D< Eigen::ConjugateGradient<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower, Eigen::IncompleteCholesky<Complex> > >::NewSP();
 
+    // LU direct
+    auto EM3D = EMSchur3D< Eigen::SparseLU<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::COLAMDOrdering<int> > >::NewSP();
 
     EM3D->SetRectilinearGrid( gimp->GetGrid() );
 

include/EMSchur3D.h

                 // << " in " << iter_done << " iterations"
               //<<  " with error " << errorn << "\t"
               << "\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;
+//              << "\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;
             timer.begin();
 
             CSolver[iw].isSymmetric(true);
-            CSolver[iw].setPivotThreshold(0.0);
+            //CSolver[iw].setPivotThreshold(0.0); // OK for symmetric complex systems with real and imaginary positive definite parts.
+            //                                  // but our imaginary part is negative definite
+            //http://www.ams.org/journals/mcom/1998-67-224/S0025-5718-98-00978-8/S0025-5718-98-00978-8.pdf
 
             /*  Complex system */
             std::cout << "SparseLU pattern analyzing C_" << iw << ",";
         }
     }
 
+    template<>
+    void EMSchur3D<   Eigen::ConjugateGradient<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower, Eigen::IncompleteCholesky<Complex> > > ::BuildCDirectSolver() {
+        CSolver = new Eigen::ConjugateGradient<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower, Eigen::IncompleteCholesky<Complex> > [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<IncompleteCholesky> 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<IncompleteCholesky> 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() {
 //         CSolver = new Eigen::PastixLLT<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower > [Omegas.size()];

include/bicgstab.h

     VectorXcr b = (ioms).asDiagonal() * (D.transpose()*Phi);
     VectorXcr y = solver.solve(b);
     //max_it = 0;
-    max_it = solver.iterations();
+    //max_it = solver.iterations();
     //errorn = solver.error();
     return  D*y;
 }