Hey there, let's work through this problem together—it’s super frustrating when a fast solver spits out garbage results while the accurate one crawls along, especially when MATLAB is showing you what’s possible. Let’s break down what might be going wrong with Pardiso and how to get that fast, correct solution you need.

First, let’s diagnose the Pardiso issue #

Your matrix II is built from terms like J^T * J, which means it’s symmetric positive definite (SPD) (since J is invertible). Pardiso’s default settings might not be tuned for SPD matrices, which is likely why you’re getting near-zero garbage solutions. Here’s how to fix that:

1. Use the right Pardiso solver variant
   Stop using the generic Pardiso interface—opt for the SPD-specific PardisoLDLT or PardisoLLT instead. These variants know to leverage the matrix’s symmetry, which not only speeds things up but also avoids numerical mistakes from treating it as a general matrix. Try this code snippet:

   Eigen::initParallel();
   II.makeCompressed();
   
   // Use SPD-optimized Pardiso solver
   Eigen::PardisoLDLT<Eigen::SparseMatrix<double>> pardiso_solver;
   pardiso_solver.setVerbose(true); // Turn on debug output to catch errors
   
   if (pardiso_solver.compute(II) != Eigen::Success) {
       std::cerr << "Pardiso decomposition failed—check matrix for issues!" << std::endl;
       return;
   }
   
   Eigen::VectorXd Delta = pardiso_solver.solve(EE);
   if (pardiso_solver.info() != Eigen::Success) {
       std::cerr << "Pardiso solve step failed!" << std::endl;
       return;
   }
   
   // Verify the solution with residual norm
   double residual_norm = (II * Delta - EE).norm();
   std::cout << "Residual norm: " << residual_norm << std::endl;
   

   The verbose output will tell you if Pardiso hit issues like singular matrices or invalid entries, which could explain the bad solution.

2. Double-check matrix symmetry
   Your II should be symmetric (since U = JPSlice^T*JPSlice, V = JDRow^T*JDRow + WeightHH, and W^T pairs with W). If there’s any numerical drift breaking symmetry (e.g., from floating-point errors), explicitly enforce it to help Pardiso:

   // Symmetrize II to fix any small numerical discrepancies
   II = 0.5 * (II + II.transpose());
   II.makeCompressed();
   

Get MATLAB-like speed with SuiteSparse/CHOLMOD #

You mentioned trying SuiteSparse with no luck—but MATLAB’s \ operator uses CHOLMOD under the hood for SPD sparse matrices, so this should be your best bet for matching MATLAB’s performance. The issue is likely incorrect setup:

• Ensure you’re linking SuiteSparse properly during compilation. You’ll need to link against cholmod, colamd, amd, and other SuiteSparse libraries.
• Use Eigen’s CholmodLDLT interface directly:

  Eigen::CholmodLLT<Eigen::SparseMatrix<double>> cholmod_solver;
  cholmod_solver.compute(II);
  Eigen::VectorXd Delta = cholmod_solver.solve(EE);
  
  // Check for success
  if (cholmod_solver.info() != Eigen::Success) {
      std::cerr << "CHOLMOD failed—check matrix condition or library links!" << std::endl;
  }
  

  CHOLMOD automatically handles reordering to reduce fill-in (just like MATLAB), which is a big part of its speed.

Other quick fixes to try #

• Check for numerical instability: If your matrix II has a high condition number (common with J^T*J where J is ill-conditioned), try adding a small diagonal shift to stabilize it:

  II += 1e-8 * Eigen::SparseMatrix<double>::Identity(II.rows(), II.cols());
  

• Parallelize properly: Make sure Eigen, MKL, and SuiteSparse are compiled with OpenMP/MKL parallelism enabled. Your Eigen::initParallel(); is a good start, but verify that the solver libraries are using multiple cores.

Final note #

If you still get bad results, compute the residual norm ((II*Delta - EE).norm()) for both the SimplicialLDLT solution and the Pardiso/CHOLMOD solution. A small residual means the solver is working correctly, and your "correct" SimplicialLDLT solution might actually be wrong (unlikely, but worth checking). A large residual points to a solver configuration or matrix issue.

备注：内容来源于stack exchange，提问作者Yingyu.Wang