SuperLU support added

CMakeLists.txt

@@ -7,6 +7,13 @@ set(EMSCHUR3D_VERSION_NOQUOTES "${EMSCHUR3D_VERSION_MAJOR}.${EMSCHUR3D_VERSION_M
 
 option ( LEMMA_MODULE_EMSCHUR3D TRUE )
 
+find_package( SuperLU )
+if (SUPERLU_FOUND)
+	message( STATUS "SuperLU was found" )
+	add_compile_options(-DHAVE_SUPERLU)
+	INCLUDE_DIRECTORIES(${SUPERLU_INCLUDES})
+endif()
+
 if ( LEMMA_VTK6_SUPPORT OR LEMMA_VTK7_SUPPORT OR LEMMA_VTK8_SUPPORT AND LEMMA_MODULE_EMSCHUR3D ) 
 
 	configure_file (
@@ -32,7 +39,7 @@ if ( LEMMA_VTK6_SUPPORT OR LEMMA_VTK7_SUPPORT OR LEMMA_VTK8_SUPPORT AND LEMMA_MO
 
 	# Linking
 	target_link_libraries(emschur3d ${VTK_LIBRARIES})
-
+	target_link_libraries(emschur3d ${SUPERLU_LIBRARIES})
 
 	# Testing
 	if (LEMMA_ENABLE_TESTING)

examples/EMSchur3D-vtk.cpp

@@ -114,8 +114,8 @@ int main( int argc, char** argv ) {
     // And solve
 
     // Use BiCGSTAB Diagonal preconditioner
-    auto EM3D = EMSchur3D< Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::IncompleteLUT<Complex> > >::NewSP();
-    //auto EM3D = EMSchur3D< Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, Eigen::ColMajor> > >::NewSP();
+    //auto EM3D = EMSchur3D< Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::IncompleteLUT<Complex> > >::NewSP();
+    auto EM3D = EMSchur3D< Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, Eigen::ColMajor> > >::NewSP();
 
     // LS CG
     //auto EM3D = EMSchur3D< Eigen::LeastSquaresConjugateGradient<Eigen::SparseMatrix<Complex, Eigen::ColMajor> > >::NewSP();

examples/EMSchur3D.cpp

@@ -29,7 +29,11 @@ using namespace Lemma;
 int main( int argc, char** argv ) {
 
     // BiCGSTAB Diagonal preconditioner
-    auto EM3D = EMSchur3D< Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::IncompleteLUT<Complex> > >::NewSP();
+    //auto EM3D = EMSchur3D< Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::IncompleteLUT<Complex> > >::NewSP();
+    //auto EM3D = EMSchur3D< Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, Eigen::ColMajor> > >::NewSP();
+
+    // SUPERLU
+    auto EM3D = EMSchur3D< Eigen::SuperLU<Eigen::SparseMatrix<Complex, Eigen::ColMajor> > >::NewSP();
 
     if (argc < 3) {
         std::cout << "EMSchur3D  <rgrid>  <1dmod> <3dmod>  <aemsurvey> " << std::endl;

include/EMSchur3D.h

@@ -24,6 +24,10 @@
 
 #include "EMSchur3DBase.h"
 #include "bicgstab.h"
+
+#ifdef HAVE_SUPERLU
+#include "Eigen/SuperLUSupport"
+#endif
 //#include "CSymSimplicialCholesky.h"
 
 namespace Lemma {
@@ -54,7 +58,6 @@ namespace Lemma {
          */
         static std::shared_ptr< EMSchur3D > NewSP() {
             return std::make_shared< EMSchur3D<Solver> >( ctor_key() );
-            //return std::make_shared< EMSchur3D< Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, Eigen::ColMajor> > > >( ctor_key() ) ;
         }
 
         /** Default protected constructor, use New */
@@ -205,7 +208,7 @@ namespace Lemma {
         logio << "solving RHS for source " << isource << std::endl;
 
         // TODO, this is stupid, try and get rid of this copy!
-        Eigen::SparseMatrix<Complex>  Cc  = Cvec[iw];
+        //Eigen::SparseMatrix<Complex>  Cc  = Cvec[iw];
 
         jsw_timer timer;
         jsw_timer timer2;
@@ -232,12 +235,13 @@ namespace Lemma {
 //         /* END EXPERIMENTAL */
 
         VectorXcr ADiv = D*A;  // ADiv == RHS == D C^I Se
-        VectorXcr Error = ((Cc.selfadjointView<Eigen::Lower>()*A).array() - Se.array());
+        //VectorXcr Error = ((Cc.selfadjointView<Eigen::Lower>()*A).array() - Se.array());
+        VectorXcr Error = ((Cvec[iw].selfadjointView<Eigen::Lower>()*A).array() - Se.array());
         logio << "|| Div(A) || = " << ADiv.norm()
               << "\tInital solution error="<<   Error.norm()  // Iteritive info
-              << "\tSolver reported error="<<   CSolver[iw].error()  // Iteritive info
+//              << "\tSolver reported error="<<   CSolver[iw].error()  // Iteritive info
               << "\ttime " << timer.end() / 60. << " [m]   "
-              <<  CSolver[iw].iterations() << "  iterations"
+//              <<  CSolver[iw].iterations() << "  iterations"
               << std::endl;
 
         //VectorXcr ADivMAC = ADiv.array() * MAC.array().cast<Complex>();
@@ -247,14 +251,13 @@ namespace Lemma {
         // Solve for Phi
         logio << "Solving for Phi " << std::flush;
         timer.begin();
-        tol = 1e-18;
+        tol = 1e-20;
         int success(2);
 
         success = implicitbicgstab(D, idx, ms, ADiv, Phi, CSolver[iw], max_it, tol, errorn, iter_done, logio);
         //Phi.array() *= MAC.array().cast<Complex>(); // remove phi from air regions
 
         /* Restart if necessary */
-/*
         int nrestart(1);
         // TODO send MAC to implicitbicgstab?
         while (success == 2 && nrestart < 18 && iter_done > 1) {
@@ -262,7 +265,6 @@ namespace Lemma {
             //Phi.array() *= MAC.array().cast<Complex>(); // remove phi from air regions
             nrestart += 1;
         }
-*/
 
         logio << "Implicit BiCGStab solution in " << iter_done << " iterations."
                 << " with error " << std::setprecision(8) << std::scientific << errorn << std::endl;
@@ -292,20 +294,23 @@ namespace Lemma {
               //<<  " with error " << errorn << "\t";
 
         // Report error of solutions
-        Error = ((Cc.selfadjointView<Eigen::Lower>()*A).array() + E.array() - Se.array());
-        //logio << "\tsolution error " << Error.norm()
-        //      << std::fixed << std::setprecision(2) << "\ttime " << timer.end()/60. << "\ttotal time " << timer2.end()/60. << std::endl;
+        //Error = ((Cc.selfadjointView<Eigen::Lower>()*A).array() + E.array() - Se.array());
+        Error = ((Cvec[iw].selfadjointView<Eigen::Lower>()*A).array() + E.array() - Se.array());
+
         //      << "\tSolver reported error="<<   CSolver[iw].error()  // Iteritive info
         //      << "\ttime " << timer.end() / 60. << " [m]   "
         //      <<  CSolver[iw].iterations() << "  iterations"
 
 
         logio << "|| Div(A) || = " << ADiv2.norm()
-              << "\tInital solution error="<<   Error.norm()  // Iteritive info
-              << "\tSolver reported error="<<   CSolver[iw].error()  // Iteritive info
+              << "\tSolution error="<<   Error.norm()  // Iteritive info
+//              << "\tSolver reported error="<<   CSolver[iw].error()  // Iteritive info
               << "\ttime " << timer.end() / 60. << " [m]   "
-              <<  CSolver[iw].iterations() << "  iterations"
-              << std::endl;
+//              <<  CSolver[iw].iterations() << "  iterations"
+              << std::endl << std::endl;
+
+        logio << std::fixed << std::setprecision(2) << "\ttime " << timer.end()/60. << "\ttotal time " << timer2.end()/60. << std::endl;
+
         logio.close();
 
         //////////////////////////////////////
@@ -356,7 +361,7 @@ namespace Lemma {
         }
     }
 
-    #ifdef HAVE_SUPERLUMT
+    #ifdef HAVE_SUPERLU
     template<>
     void EMSchur3D< Eigen::SuperLU<Eigen::SparseMatrix<Complex, Eigen::ColMajor> > >::BuildCDirectSolver() {
 
@@ -367,8 +372,8 @@ namespace Lemma {
             timer.begin();
 
             /* SuperLU */
-            //CSolver[iw].options().DiagPivotThresh = 0.01;
-            //CSolver[iw].options().SymmetricMode = YES;
+            CSolver[iw].options().DiagPivotThresh = 0.0;
+            CSolver[iw].options().SymmetricMode = YES;
             //CSolver[iw].options().ColPerm = MMD_AT_PLUS_A;
             //CSolver[iw].options().Trans = NOTRANS;
             //CSolver[iw].options().ConditionNumber = NO;
@@ -381,14 +386,14 @@ namespace Lemma {
             //std::cout << "\tCondition Number: " << CSolver[iw].options().ConditionNumber << std::endl;
 
             /*  Complex system */
-            std::cout << "SuperLU_MT pattern analyzing C_" << iw << ",";
+            std::cout << "SuperLU pattern analyzing C_" << iw << ",";
             std::cout.flush();
             CSolver[iw].analyzePattern( Cvec[iw].selfadjointView< Eigen::Lower>() );
             std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
 
             // factorize
             timer.begin();
-            std::cout << "SuperLU_MT factorising C_" << iw << ", ";
+            std::cout << "SuperLU factorising C_" << iw << ", ";
             std::cout.flush();
             CSolver[iw].factorize( Cvec[iw].selfadjointView< Eigen::Lower>() );
             std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
@@ -514,20 +519,22 @@ namespace Lemma {
         CSolver = new Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::IncompleteLUT<Complex> > [Omegas.size()];
         for (int iw=0; iw<Omegas.size(); ++iw) {
             Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
-            CSolver[iw].preconditioner().setDroptol(1e-5);      // 1e-12
-            //CSolver[iw].preconditioner().setFillfactor(5e1);     // 1e2
+            CSolver[iw].preconditioner().setDroptol(1e-6);       //1e-5);      // 1e-12
+            CSolver[iw].preconditioner().setFillfactor(5e1);     // 1e2
             jsw_timer timer;
             timer.begin();
             /*  Complex system */
             std::cout << "BiCGSTAB(ILU) pattern analyzing C_" << iw << ",";
             std::cout.flush();
             CSolver[iw].analyzePattern( Csym );
+            //CSolver[iw].analyzePattern( Cvec[iw]);
             std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
             /* factorize */
             timer.begin();
             std::cout << "BiCGSTAB(ILU) factorising C_" << iw << ", ";
             std::cout.flush();
             CSolver[iw].factorize( Csym );
+            //CSolver[iw].factorize( Cvec[iw] );
             std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
         }
     }

include/bicgstab.h

@@ -55,428 +55,25 @@
 using namespace Eigen;
 using namespace Lemma;
 
-//typedef Eigen::VectorXcd VectorXcr;
 typedef Eigen::SparseMatrix<Complex> SparseMat;
 
 
-// On Input
-// A = Matrix
-// B = Right hand side
-// X = initial guess, and solution
-// maxit = maximum Number of iterations
-// tol = error tolerance
-// On Output
-// X real solution vector
-// errorn = Real error norm
-int bicgstab(const SparseMat &A, const SparseMat &M, const VectorXcr &b, VectorXcr &x,
-                int &max_it, Real &tol, Real &errorn, int &iter_done,
-                const bool& banner = true) {
-
-    Complex omega, rho, rho_1, alpha, beta;
-    Real bnrm2, eps, errmin;
-    int n, iter; //, istat;
-
-    // Determine size of system and init vectors
-    n = x.size();
-    VectorXcr r(n);
-    VectorXcr r_tld(n);
-    VectorXcr p(n);
-    VectorXcr v(n);
-    VectorXcr p_hat(n);
-    VectorXcr s(n);
-    VectorXcr s_hat(n);
-    VectorXcr t(n);
-    VectorXcr xmin(n);
-
-    if (banner) {
-        std::cout << "Start BiCGStab, memory needed: "
-                  <<  (sizeof(Complex)*(9+2)*n/(1024.*1024*1024)) << " [Gb]\n";
-    }
-
-    // Initialise
-    iter_done = 0;
-    v.setConstant(0.); // not necessary I don't think
-    t.setConstant(0.);
-    eps = 1e-100;
-
-    bnrm2 = b.norm();
-    if (bnrm2 == 0) {
-        x.setConstant(0.0);
-        errorn = 0;
-        std::cerr << "Trivial case of Ax = b, where b is 0\n";
-        return (0);
-    }
-
-    // If there is an initial guess
-    if ( x.norm() ) {
-        r = b - A.selfadjointView<Eigen::Upper>()*x;
-        //r = b - A*x;
-    } else {
-        r = b;
-    }
-
-    errorn = r.norm() / bnrm2;
-    omega = 1.;
-    r_tld = r;
-    errmin = 1e30;
-
-    // Get down to business
-    for (iter=0; iter<max_it; ++iter) {
-
-        rho = r_tld.dot(r);
-        if ( abs(rho) < eps) return (0);
-
-        if (iter > 0) {
-            beta = (rho/rho_1) * (alpha/omega);
-            p = r.array() + beta*(p.array()-omega*v.array()).array();
-        } else {
-            p = r;
-        }
-
-        // Use pseudo inverse to get approximate answer
-        //#pragma omp sections
-        p_hat = M*p;
-        //v = A*p_hat; // TODO double check
-        v = A.selfadjointView<Eigen::Upper>()*p_hat; // TODO double check
-
-        alpha = rho / r_tld.dot(v);
-        s = r.array() - alpha*v.array();
-        errorn = s.norm()/bnrm2;
-
-        if (errorn < tol && iter > 1) {
-            x.array() += alpha*p_hat.array();
-            return (0);
-        }
-
-        s_hat = M*s;
-        t = A.selfadjointView<Eigen::Upper>()*s_hat;
-        //t = A*s_hat;
-
-        omega = t.dot(s)  / t.dot(t);
-        x.array() += alpha*p_hat.array() + omega*s_hat.array();
-        r = s.array() - omega*t.array();
-        errorn = r.norm() / bnrm2;
-        iter_done = iter;
-
-        if (errorn < errmin) {
-            // remember the model with the smallest norm
-            errmin = errorn;
-            xmin = x;
-        }
-
-        if ( errorn <= tol ) return (0);
-        if ( abs(omega) < eps ) return (0);
-        rho_1 = rho;
-
-    }
-    return (0);
-}
-
-template <typename Preconditioner>
-bool preconditionedBiCGStab(const SparseMat &A, const Preconditioner &M,
-        const Ref< VectorXcr const > b,
-        Ref <VectorXcr > x,
-        const int &max_it, const Real &tol,
-        Real &errorn, int &iter_done) {
-
-    Complex omega, rho, rho_1, alpha, beta;
-    Real bnrm2, eps;
-    int n, iter;
-    Real tol2 = tol*tol;
-
-    // Determine size of system and init vectors
-    n = x.size();
-
-    VectorXcd r(n);
-    VectorXcd r_tld(n);
-    VectorXcd p(n);
-    VectorXcd s(n);
-    VectorXcd s_hat(n);
-    VectorXcd p_hat(n);
-    VectorXcd v = VectorXcr::Zero(n);
-    VectorXcd t = VectorXcr::Zero(n);
-
-    //std::cout << "Start BiCGStab, memory needed: "
-    //          <<  (sizeof(Complex)*(8+2)*n/(1024.*1024)) / (1024.) << " [Gb]\n";
-
-    // Initialise
-    iter_done = 0;
-    eps = 1e-100;
-
-    bnrm2 = b.squaredNorm();
-    if (bnrm2 == 0) {
-        x.setConstant(0.0);
-        errorn = 0;
-        std::cerr << "Trivial case of Ax = b, where b is 0\n";
-        return (false);
-    }
 
-    // If there is an initial guess
-    if ( x.squaredNorm() ) {
-        r = b - A.selfadjointView<Eigen::Upper>()*x;
-    } else {
-        r = b;
-    }
-
-    errorn = r.squaredNorm() / bnrm2;
-    omega = 1.;
-    r_tld = r;
-
-    // Get down to business
-    for (iter=0; iter<max_it; ++iter) {
-
-        rho = r_tld.dot(r);
-        if (abs(rho) < eps) {
-            std::cerr << "arbitrary orthogonality issue in bicgstab\n";
-            std::cerr << "consider eigen restarting\n";
-            return (false);
-        }
-
-        if (iter > 0) {
-            beta = (rho/rho_1) * (alpha/omega);
-            p = r + beta*(p-omega*v);
-        } else {
-            p = r;
-        }
-
-        p_hat = M.solve(p);
-        v.noalias() = A.selfadjointView<Eigen::Upper>()*p_hat;
-
-        alpha = rho / r_tld.dot(v);
-        s = r - alpha*v;
-        errorn = s.squaredNorm()/bnrm2;
-
-        if (errorn < tol2 && iter > 1) {
-            x = x + alpha*p_hat;
-            errorn = std::sqrt(errorn);
-            return (true);
-        }
-
-        s_hat = M.solve(s);
-        t.noalias() = A.selfadjointView<Eigen::Upper>()*s_hat;
-
-        omega = t.dot(s)  / t.dot(t);
-        x += alpha*p_hat + omega*s_hat;
-        r = s - omega*t;
-        errorn = r.squaredNorm() / bnrm2;
-        iter_done = iter;
-
-        if ( errorn <= tol2 || abs(omega) < eps) {
-            errorn = std::sqrt(errorn);
-            return (true);
-        }
-
-        rho_1 = rho;
-    }
-    return (false);
-}
-
-template <typename Preconditioner>
-bool preconditionedSCBiCG(const SparseMat &A, const Preconditioner &M,
-        const Ref< VectorXcr const > b,
-        Ref <VectorXcr > x,
-        const int &max_iter, const Real &tol,
-        Real &errorn, int &iter_done) {
-
-    Real resid;
-    VectorXcr p, z, q;
-    Complex alpha, beta, rho, rho_1;
-
-    Real normb = b.norm( );
-    VectorXcr r = b - A*x;
-
-    if (normb == 0.0) normb = 1;
-
-    if ((resid = r.norm( ) / normb) <= tol) {
-        errorn = resid;
-        iter_done = 0;
-        return 0;
-    }
-
-    for (int i = 1; i <= max_iter; i++) {
-        z = M.solve(r);
-        rho = r.dot(z);
-
-        if (i == 1)  p = z;
-        else {
-            beta = rho / rho_1;
-            p = z + beta * p;
-        }
-
-        q = A*p;
-        alpha = rho / p.dot(q);
-
-        x += alpha * p;
-        r -= alpha * q;
-        std::cout << "resid\t" << resid << std::endl;
-        if ((resid = r.norm( ) / normb) <= tol) {
-            errorn = resid;
-            iter_done = i;
-            return 0;
-        }
-
-        rho_1 = rho;
-    }
-
-    errorn = resid;
-
-    return (false);
-}
-
-
-/** \internal Low-level conjugate gradient algorithm
-  * \param mat The matrix A
-  * \param rhs The right hand side vector b5
-  * \param x On input and initial solution, on output the computed solution.
-  * \param precond A preconditioner being able to efficiently solve for an
-  *                approximation of Ax=b (regardless of b)
-  * \param iters On input the max number of iteration, on output the number of performed iterations.
-  * \param tol_error On input the tolerance error, on output an estimation of the relative error.
-  */
-template<typename Rhs, typename Dest, typename Preconditioner>
-EIGEN_DONT_INLINE
-void conjugateGradient(const SparseMat& mat, const Rhs& rhs, Dest& x,
-                        const Preconditioner& precond, int& iters,
-                        typename Dest::RealScalar& tol_error)
-{
-  using std::sqrt;
-  using std::abs;
-  typedef typename Dest::RealScalar RealScalar;
-  typedef typename Dest::Scalar Scalar;
-  typedef Matrix<Scalar,Dynamic,1> VectorType;
-
-  RealScalar tol = tol_error;
-  int maxIters = iters;
-
-  int n = mat.cols();
-
-  VectorType residual = rhs - mat.selfadjointView<Eigen::Upper>() * x; //initial residual
-
-  RealScalar rhsNorm2 = rhs.squaredNorm();
-  if(rhsNorm2 == 0)
-  {
-    x.setZero();
-    iters = 0;
-    tol_error = 0;
-    return;
-  }
-  RealScalar threshold = tol*tol*rhsNorm2;
-  RealScalar residualNorm2 = residual.squaredNorm();
-  if (residualNorm2 < threshold)
-  {
-    iters = 0;
-    tol_error = sqrt(residualNorm2 / rhsNorm2);
-    return;
-  }
-
-  VectorType p(n);
-  p = precond.solve(residual);      //initial search direction
-
-  VectorType z(n), tmp(n);
-  RealScalar absNew = numext::real(residual.dot(p));  // the square of the absolute value of r scaled by invM
-  int i = 0;
-  while(i < maxIters)
-  {
-    tmp.noalias() = mat.selfadjointView<Eigen::Upper>() * p;              // the bottleneck of the algorithm
-
-    Scalar alpha = absNew / p.dot(tmp);   // the amount we travel on dir
-    x += alpha * p;                       // update solution
-    residual -= alpha * tmp;              // update residue
-
-    residualNorm2 = residual.squaredNorm();
-    if(residualNorm2 < threshold)
-      break;
-
-    z = precond.solve(residual);          // approximately solve for "A z = residual"
-
-    RealScalar absOld = absNew;
-    absNew = numext::real(residual.dot(z));     // update the absolute value of r
-    RealScalar beta = absNew / absOld;            // calculate the Gram-Schmidt value used to create the new search direction
-    p = z + beta * p;                             // update search direction
-    i++;
-  }
-  tol_error = sqrt(residualNorm2 / rhsNorm2);
-  iters = i;
-}
-
-// // Computes implicit
-// VectorXcr implicitDCInvBPhi (const SparseMat& D, const SparseMat& C,
-//                         const SparseMat& B, const SparseMat& MC,
-//                         const VectorXcr& Phi, Real& tol,
-//                         int& max_it) {
-//     int iter_done(0);
-//     Real errorn(0);
-//     VectorXcr b = B*Phi;
-//     VectorXcr y = VectorXcr::Zero(C.rows()) ; // = C^1*b;
-//     bicgstab(C, MC, b, y, max_it, tol, errorn, iter_done, false);
-//     //std::cout << "Temp " << errorn << std::endl;
-//     return  D*y;
-// }
-
-// Computes implicit
-VectorXcr implicitDCInvBPhi (const SparseMat& D, const SparseMat& C,
-                        const VectorXcr& ioms, const SparseMat& MC,
-                        const VectorXcr& Phi, Real& tol,
-                        int& max_it) {
-    int iter_done(0);
-    Real errorn(0);
-    VectorXcr b = (ioms).asDiagonal() * (D.transpose()*Phi);
-    VectorXcr y = VectorXcr::Zero(C.rows()) ; // = C^1*b;
-    bicgstab(C, MC, b, y, max_it, tol, errorn, iter_done, false);
-    //std::cout << "Temp " << errorn << std::endl;
-    max_it = iter_done;
-    return  D*y;
-}
-
-// Computes implicit
-template <typename Preconditioner>
-VectorXcr implicitDCInvBPhi2 (const SparseMat& D, const SparseMat& C,
-                        const Ref<VectorXcr const> ioms, const Preconditioner& solver,
-                        const Ref<VectorXcr const> Phi, Real& tol,
-                        int& max_it) {
-
-    VectorXcr b = (ioms).asDiagonal() * (D.transpose()*Phi);
-    VectorXcr y = VectorXcr::Zero(C.rows()) ; // = C^1*b;
-
-    // Home Made
-    //int iter_done(0);
-    //Real errorn(0);
-    //preconditionedBiCGStab(C, solver, b, y, max_it, tol, errorn, iter_done); //, false); // Jacobi M
-    //max_it = iter_done;
-
-    // Eigen BiCGStab
-    Eigen::BiCGSTAB<SparseMatrix<Complex> > BiCG;
-    BiCG.compute( C ); // TODO move this out of this loop!
-    y = BiCG.solve(b);
-    max_it = BiCG.iterations();
-    tol = BiCG.error();
-
-    // Direct
-/*
-    std::cout << "Computing LLT" << std::endl;
-    Eigen::SimplicialLLT<SparseMatrix<Complex>, Eigen::Upper, Eigen::AMDOrdering<int> >  LLT;
-    LLT.compute(C);
-    max_it = 1;
-    std::cout << "Computed LLT" << std::endl;
-    y = LLT.solve(b);
-*/
-
-    return  D*y;
-}
-
-// Computes implicit
-//template <typename Solver>
+// Computes implicit calculation
 template < typename Solver >
-inline VectorXcr implicitDCInvBPhi3 (const SparseMat& D, const Solver& solver,
+inline VectorXcr implicitDCInvBPhi3 (
+                        const SparseMat& D,
+                        const Solver& solver,
                         const Ref<VectorXcr const> ioms,
-                        const Ref<VectorXcr const> Phi, Real& tol,
-                        int& max_it) {
+                        const Ref<VectorXcr const> Phi,
+                        Real& tol,   // not used
+                        int& max_it  // not used
+                ) {
     VectorXcr b = (ioms).asDiagonal() * (D.transpose()*Phi);
     VectorXcr y = solver.solve(b);
-    //max_it = 0;
-    max_it = solver.iterations();
-    //errorn = solver.error();
+    //max_it = solver.iterations(); // actualy no need to pass this
     return  D*y;
+    //return  y;
 }
 
 
@@ -536,12 +133,19 @@ int implicitbicgstab(//const SparseMat& D,
 
     // Determine size of system and init vectors
     int n = idx.size();        // was phi.size();
+
+    std::cout << "BiCGStab SIZES  " << n << "\t" << phi.size() << "\t" << ioms.size() << std::endl;
+
     VectorXcr r(n);
     VectorXcr r_tld(n);
     VectorXcr p(n);
     VectorXcr s(n);
-    VectorXcr v = VectorXcr::Zero(n);
-    VectorXcr t = VectorXcr::Zero(n);
+
+    VectorXcr v = VectorXcr::Zero(ioms.size());
+    VectorXcr t = VectorXcr::Zero(ioms.size());
+
+//    VectorXcr vm1 = VectorXcr::Zero(ioms.size());
+//    VectorXcr tm1 = VectorXcr::Zero(ioms.size());
 
 //     TODO, refigure for implicit large system
 //     std::cout << "Start BiCGStab, memory needed: "
@@ -597,7 +201,6 @@ int implicitbicgstab(//const SparseMat& D,
         tol2 = tol;
 
         max_it2 = 500000;
-        //v = implicitDCInvBPhi2(D, C, ioms, solver, p, tol2, max_it2);
         ivmap(phi, p, idx);
         v = vmap(implicitDCInvBPhi3(D, solver, ioms, phi, tol2, max_it2), idx);
 
@@ -643,7 +246,6 @@ int implicitbicgstab(//const SparseMat& D,
               << max_it2+max_it2 << " iterations " << std::endl;
 
         // Check to see how progress is going
-
         if (errornold - errorn < 0) {
             logio << "Irregular non-monotonic (negative) convergence. Recommend restart. \n";
             ivmap( phi, phi2, idx );
@@ -668,275 +270,5 @@ int implicitbicgstab(//const SparseMat& D,
     return (0);
 }
 
-// On Input
-// A = Matrix
-// B = Right hand side
-// X = initial guess, and solution
-// maxit = maximum Number of iterations
-// tol = error tolerance
-// On Output
-// X real solution vector
-// errorn = Real error norm
-template < typename Solver >
-int implicitbicgstab_ei(const SparseMat&  D,
-                        const Ref< VectorXcr const > ioms,
-                        const Ref< VectorXcr const > rhs,
-                        Ref <VectorXcr> phi,
-                        Solver& solver,
-                        int &max_it, const Real &tol, Real &errorn, int &iter_done, ofstream& logio) {
-
-    logio << "using the preconditioned Eigen implicit solver" << std::endl;
-
-    Complex omega, rho, rho_1, alpha, beta;
-    Real tol2;
-    int  iter, max_it2,max_it1;
-
-    // Determine size of system and init vectors
-    int n = phi.size();
-    VectorXcr r(n);
-    VectorXcr r_tld(n);
-    VectorXcr p(n);
-    VectorXcr v(n);
-    VectorXcr s(n);
-    VectorXcr t(n);
-
-    // Initialise
-    iter_done = 0;
-    Real eps = 1e-100;
-
-    Real bnrm2 = rhs.norm();
-    if (bnrm2 == 0) {
-        phi.setConstant(0.0);
-        errorn = 0;
-        std::cerr << "Trivial case of Ax = b, where b is 0\n";
-        return (0);
-    }
-
-    // If there is an initial guess
-    if ( phi.norm() ) {
-        tol2 = tol;
-        max_it2 = 50000;
-        r = rhs - implicitDCInvBPhi3(D, solver, ioms, phi, tol2, max_it2);
-    } else {
-        r = rhs;
-    }
-
-    jsw_timer timer;
-
-    errorn = r.norm() / bnrm2;
-    omega = 1.;
-    r_tld = r;
-    Real errornold = 1e14;
-
-    // Get down to business
-    for (iter=0; iter<max_it; ++iter) {
-
-        timer.begin();
-
-        rho = r_tld.dot(r);
-        if (abs(rho) < eps) return (0);
-
-        if (iter > 0) {
-            beta = (rho/rho_1) * (alpha/omega);
-            p = r.array() + beta*(p.array()-omega*v.array()).array();
-        } else {
-            p = r;
-        }
-
-        tol2 = tol;
-        max_it2 = 500000;
-        v = implicitDCInvBPhi3(D, solver, ioms, p, tol2, max_it2);
-        max_it2 = solver.iterations();
-
-        alpha = rho / r_tld.dot(v);
-        s = r.array() - alpha*v.array();
-        errorn = s.norm()/bnrm2;
-
-        if (errorn < tol && iter > 1) {
-            phi.array() += alpha*p.array();
-            return (0);
-        }
-
-        tol2 = tol;
-        max_it1 = 500000;
-        t = implicitDCInvBPhi3(D, solver, ioms, s, tol2, max_it1);
-        max_it1 = solver.iterations();
-        omega = t.dot(s)  / t.dot(t);
-
-        r = s.array() - omega*t.array();
-        errorn = r.norm() / bnrm2;
-        iter_done = iter;
-
-        if (errorn <= tol ) return (0);
-        if (abs(omega) < eps) return (0);
-        rho_1 = rho;
-
-        logio << "iteration " << std::setw(4) << iter
-              << "\terrorn " << std::setw(6) << std::setprecision(4) << std::scientific << errorn
-              << "\timplicit iterations " << std::setw(5) << max_it1+max_it2
-              << "\ttime " << std::setw(10) << std::fixed << std::setprecision(2) << timer.end() << std::endl;
-
-        // Check to see how progress is going
-        if (errornold - errorn < 0) {
-            logio << "irregular (negative) convergence. Try again? \n";
-            return (2);
-        }
-
-        // only update phi if good things are happening
-        phi.array() += alpha*p.array() + omega*s.array();
-        errornold = errorn;
-
-    }
-    return (0);
-}
-
-
-// On Input
-// A = Matrix
-// B = Right hand side
-// X = initial guess, and solution
-// maxit = maximum Number of iterations
-// tol = error tolerance
-// On Output
-// X real solution vector
-// errorn = Real error norm
-int implicitbicgstabnt(const SparseMat& D,
-                     const SparseMat& C,
-                     const VectorXcr& ioms,
-                     const SparseMat& MC,
-                     Eigen::Ref< VectorXcr > rhs,
-                     VectorXcr& phi,
-                     int &max_it, Real &tol, Real &errorn, int &iter_done) {
-
-    Complex omega, rho, rho_1, alpha, beta;
-    Real errmin, tol2;
-    int  iter, max_it2;
-
-//     // Cholesky decomp
-//     SparseLLT<SparseMatrix<Complex>, Cholmod>
-//         CholC(SparseMatrix<Complex> (C.real()) );
-//     if(!CholC.succeeded()) {
-//         std::cerr << "decomposiiton failed\n";
-//         return EXIT_FAILURE;
-//     }
-
-    // Determine size of system and init vectors
-    int n = phi.size();
-    VectorXcr r(n);
-    VectorXcr r_tld(n);
-    VectorXcr p(n);
-    VectorXcr v(n);
-    //VectorXcr p_hat(n);
-    VectorXcr s(n);
-    //VectorXcr s_hat(n);
-    VectorXcr t(n);
-    VectorXcr xmin(n);
-
-//     TODO, refigure for implicit large system
-//     std::cout << "Start BiCGStab, memory needed: "
-//               <<  (sizeof(Complex)*(9+2)*n/(1024.*1024*1024)) << " [Gb]\n";
-
-    // Initialise
-    iter_done = 0;
-    v.setConstant(0.); // not necessary I don't think
-    t.setConstant(0.);
-    Real eps = 1e-100;
-
-    Real bnrm2 = rhs.norm();
-    if (bnrm2 == 0) {
-        phi.setConstant(0.0);
-        errorn = 0;
-        std::cerr << "Trivial case of Ax = b, where b is 0\n";
-        return (0);
-    }
-
-    // If there is an initial guess
-    if ( phi.norm() ) {
-        //r = rhs - A*phi;
-        tol2 = tol;
-        max_it2 = 50000;
-        std::cout << "Initial guess " << std::endl;
-        r = rhs - implicitDCInvBPhi(D, C, ioms, MC, phi, tol2, max_it2);
-        //r = rhs - implicitDCInvBPhi (D, C, B, CholC, phi, tol2, max_it2);
-    } else {
-        r = rhs;
-    }
-
-
-    errorn = r.norm() / bnrm2;
-    //std::cout << "Initial |r|  " << r.norm() << "\t" << errorn<< std::endl;
-    omega = 1.;
-    r_tld = r;
-    errmin = 1e30;
-    Real errornold = 1e6;
-    // Get down to business
-    for (iter=0; iter<max_it; ++iter) {
-
-        rho = r_tld.dot(r);
-        if (abs(rho) < eps) return (0);
-
-        if (iter > 0) {
-            beta = (rho/rho_1) * (alpha/omega);
-            p = r.array() + beta*(p.array()-omega*v.array()).array();
-        } else {
-            p = r;
-        }
-
-        // Use pseudo inverse to get approximate answer
-        //p_hat = p;
-        tol2  = std::max(1e-4*errorn, tol);
-        tol2 = tol;
-        max_it2 = 500000;
-        //v = A*p_hat;
-        v = implicitDCInvBPhi(D, C, ioms, MC, p, tol2, max_it2);
-        //v = implicitDCInvBPhi(D, C, B, CholC, p, tol2, max_it2);
-
-        alpha = rho / r_tld.dot(v);
-        s = r.array() - alpha*v.array();
-        errorn = s.norm()/bnrm2;
-
-        if (errorn < tol && iter > 1) {
-            phi.array() += alpha*p.array();
-            return (0);
-        }
-
-        // s_hat = M*s;
-        //tol2 = tol;
-        tol2  = std::max(1e-4*errorn, tol);
-        tol2 = tol;
-        max_it2 = 50000;
-        // t = A*s_hat;
-        t = implicitDCInvBPhi(D, C, ioms, MC, s, tol2, max_it2);
-        //t = implicitDCInvBPhi(D, C, B, CholC, s, tol2, max_it2);
-        omega = t.dot(s)  / t.dot(t);
-        phi.array() += alpha*p.array() + omega*s.array();
-        r = s.array() - omega*t.array();
-        errorn = r.norm() / bnrm2;
-        iter_done = iter;
-        if (errorn < errmin) {
-            // remember the model with the smallest norm
-            errmin = errorn;
-            xmin = phi;
-        }
-
-        if (errorn <= tol ) return (0);
-        if (abs(omega) < eps) return (0);
-        rho_1 = rho;
-
-        std::cout << "iteration " << std::setw(4) << iter << "\terrorn "  << std::setw(6) << std::scientific << errorn
-                  << "\timplicit iterations " << std::setw(5) << max_it2 << std::endl;
-        if (errornold - errorn < 1e-14) {
-            std::cout << "not making any progress. Giving up\n";
-            return (2);
-        }
-        errornold = errorn;
-
-    }
-    return (0);
-}
-
 #endif   // ----- #ifndef BICGSTAB_INC  -----
 
-//int bicgstab(const SparseMat &A, Eigen::SparseLU< Eigen::SparseMatrix<Complex, RowMajor> ,
-
-

src/EMSchur3DBase.cpp

@@ -26,8 +26,8 @@
 typedef Eigen::Triplet<Lemma::Complex> Tc;
 typedef Eigen::Triplet<Lemma::Real> Tr;
 
-#define UPPER 0
-#define LOWER 1 // 1=true, 0=false
+#define UPPER 0  // LOWER WAS 0
+#define LOWER 1  // 1=true, 0=false
 
 namespace Lemma {
 
@@ -107,8 +107,12 @@ namespace Lemma {
     void EMSchur3DBase::BuildC ( Real*** sigmax, Real*** sigmay, Real*** sigmaz, const int& iw) {
 
         Cvec[iw].resize( unx+uny+unz , unx+uny+unz );
+
+#if LOWER && UPPER
+        Cvec[iw].reserve(Eigen::VectorXi::Constant(unx+uny+unz, 7));   // Whole
+#else
         Cvec[iw].reserve(Eigen::VectorXi::Constant(unx+uny+unz, 4)); // Upper/Lower
-        //Cvec[iw].reserve(Eigen::VectorXi::Constant(unx+uny+unz, 7));   // Whole
+#endif
 
         //Cvec_s.resize( idx. )
         //CMMvec[iw].resize( unx+uny+unz , unx+uny+unz );