hello!

When I use the BiCGSTAB to slove the 10000-dimensions sparse matrix, it worked out successfully. But when I use the BiCGSTAB to solve the 250000-dimensions spare matrix which is very similar to the 10000-dimensions spare matrix I solved before in structure, it turned out to be all "nan" value. Can you help me please? Thank you !

By the way, the Simplicilcholesky can solve both the 10000-dimensions and the 250000-dimensions, but I suppose the matrix is not SPD, so I want to examine the results with BiCGSTAB.

here is the code:

Code: Select all

/main.cpp/
#include <Eigen/Eigen>
#include <Eigen/Sparse>
#include <vector>
#include <iostream>
#include <cmath>

using namespace std;
using namespace Eigen;


typedef Eigen::SparseMatrix< double > SpMat;
typedef Eigen::Triplet<  double  > T;

void build(vector<T>& coefficients,Eigen::VectorXd& b, int n,double k,double h,int nx, int ny , Eigen::VectorXd& e);
 int n=500;
    int m=n*n;
    VectorXd x(m),b(m),e(m);
    vector<T> coefficients;

int main()
{
    for(int j=0;j<n;j++)
    {
        for(int i=0;i<n;i++)
            e(i+j*n)=1;
    }


    for(int j=0;j<n;j++)
    {
     for(int i=0;i<n;i++)
     {
         if((i<=25)||(i>=475)||(j<=25)||(j>=475))
         e(i+j*n)=6;
         if((i>=300)&&(i<=325)&&(j<=150))
         e(i+j*n)=6;
         if((j>=125)&&(j<=150)&&(i>=375))
         e(i+j*n)=6;
         if((j>=150)&&(j<=175)&&(i<=225))
         e(i+j*n)=6;
         if((j>=375)&&(j<=400)&&(i<=225))
         e(i+j*n)=6;
         if((i>=200)&&(i<=225)&&(j>=150)&&(j<=210))
         e(i+j*n)=6;
         if((i>=200)&&(i<=225)&&(j>=300)&&(j<=400))
         e(i+j*n)=6;
     }
    }


    double lenth=5;
    double k=2*3.1415*2.4*1000000000/300000000;
    double h=1.0*lenth/n;
    build(coefficients, b, n,k,h,100,300,e );
    SpMat A(m,m);
    A.setFromTriplets(coefficients.begin(), coefficients.end());
    BiCGSTAB<SpMat> chol;
    chol.compute(A);
    x = chol.solve(b);
    for(int i=0;i<=m-1;i++)
    cout<<x(i)<<" ";
    return 0;
}


/fuction.cpp/
#include <iostream>
#include <Eigen/Sparse>
#include <vector>

using namespace std;


typedef Eigen::SparseMatrix< double > SpMat;
typedef Eigen::Triplet<  double  > T;

void inser(int id, int i, int j, double w, vector<T>& coeffs,int n,
                       Eigen::VectorXd& b)
{

  int id1 = i+j*n;
        if(i==-1 || i==n) b(id) -=0;
  else  if(j==-1 || j==n) b(id) -=0;
  else  coeffs.push_back(T(id,id1,w));

}


void build(vector<T>& coefficients, Eigen::VectorXd& b, int n,double k,double h,int nx, int ny ,Eigen::VectorXd& e )
{
  b.setZero();
  for(int j=0; j<n; ++j)
  {
    for(int i=0; i<n; ++i)
    {
        int id = i+j*n;

      inser(id, i-1,j, -1, coefficients, n,b);
      inser(id, i+1,j, -1, coefficients, n,b);
      inser(id, i,j-1, -1, coefficients, n,b);
      inser(id, i,j+1, -1, coefficients, n,b);
      inser(id, i,j, 4-k*k*h*h*e(i+j*n), coefficients,n, b);
      if((i==nx)&&(j==ny))
      b(id)-=h*h;
    }
  }
}



The BiCGSTAB algorithm is not guaranteed to always converge, especially if the problem is not SPD, so that's not a good option to check the solution of a LDLT factorization. The best is to check the residual error:

Code: Select all

cout << "error: " << (A*x-b).norm()/b.norm() << endl;


Thank you for your reply!

Then，how can I solve a matrix which is not SPD ?

Give a try at SparseLU and SparseQR