Hi, everyone! I use BiCGSTAB to solve a sparse matrix. This is my code:

Code: Select all

#include "wimplicitFD3D.h"

void implicitFinite3D(sucomplex *data, int xsize, int ysize, float dx, float dy, float *vel
   , float dz, float v0, float omega, int verbose)
{
   int size = xsize*ysize;
   sucomplex im(0, 1);
   
   sucomplex ccc(0,0);
   sucomplex dd, aim1j, aijm1, aij, aija1, aia1j;

   int n = size;

   VectorXcf x(n), d(n); SparseMatrix< std::complex<float>, ColMajor > A(n, n);

// fill A and d A * x = d

   for (int i = 0; i<xsize; i++)
   {
      for (int j = 0; j<ysize; j++)
      {
         float aaa = (vel[j*xsize + i] * vel[j*xsize + i] + vel[j*xsize + i] * v0 + v0*v0) / (4 * omega*omega);
         float bbb = (vel[j*xsize + i] - v0) / (2 * omega);
         ccc = aaa + im*bbb*dz;

         // fill A and d
         if (0 == i)
         {
            if (0 == j)
            {
               dd = (dx*dx*dy*dy - 2 * ccc*(dx*dx+dy*dy))*data[j*xsize + i] +
                  ccc*dx*dx*data[(j + 1)*xsize + i] + ccc*dy*dy*data[j*xsize + i + 1];

               aij = (dx*dx*dy*dy - 2 * aaa*(dx*dx + dy*dy));
               aija1 = aaa * dx*dx;
               aia1j = aaa * dy*dy;
               A.coeffRef(i, j).imag ( aij.i); A.coeffRef(i, j).real ( aij.r);
               A.coeffRef(i, j + 1).imag ( aija1.i); A.coeffRef(i, j + 1).real ( aija1.r);
               A.coeffRef(i + 1, j).imag ( aia1j.i); A.coeffRef(i + 1, j).real ( aia1j.r);
            }
            else if ((ysize-1) == j)
            {
               dd = ccc*dx*dx*data[(j - 1)*xsize + i] +
                  (dx*dx*dy*dy - 2 * ccc*(dx*dx + dy*dy))*data[j*xsize + i] +
                   ccc*dy*dy*data[j*xsize + i + 1];

               aij = (dx*dx*dy*dy - 2 * aaa*(dx*dx + dy*dy));
               aijm1 = aaa * dx*dx;
               aia1j = aaa * dy*dy;
               A.coeffRef(i, j).imag ( aij.i); A.coeffRef(i, j).real ( aij.r);
               A.coeffRef(i, j - 1).imag ( aijm1.i); A.coeffRef(i, j - 1).real ( aijm1.r);
               A.coeffRef(i + 1, j).imag ( aia1j.i); A.coeffRef(i + 1, j).real ( aia1j.r);
            }   
         }
         else if ((xsize-1) == i)
         {
            if (0 == j)
            {
               dd = (dx*dx*dy*dy - 2 * ccc*(dx*dx + dy*dy))*data[j*xsize + i] +
                  ccc*dx*dx*data[(j + 1)*xsize + i] + ccc*dy*dy*data[j*xsize + i - 1];

               aij = (dx*dx*dy*dy - 2 * aaa*(dx*dx + dy*dy));
               aija1 = aaa * dx*dx;
               aim1j = aaa * dy*dy;
               A.coeffRef(i, j).imag ( aij.i); A.coeffRef(i, j).real ( aij.r);
               A.coeffRef(i, j + 1).imag ( aija1.i); A.coeffRef(i, j + 1).real ( aija1.r);
               A.coeffRef(i - 1, j).imag ( aim1j.i); A.coeffRef(i - 1, j).real ( aim1j.r);
            }
            else if ((ysize - 1) == j)
            {
               dd = ccc*dx*dx*data[(j - 1)*xsize + i] +
                  (dx*dx*dy*dy - 2 * ccc*(dx*dx + dy*dy))*data[j*xsize + i] +
                  ccc*dy*dy*data[j*xsize + i - 1];

               aij = (dx*dx*dy*dy - 2 * aaa*(dx*dx + dy*dy));
               aijm1 = aaa * dx*dx;
               aim1j = aaa * dy*dy;
               A.coeffRef(i, j).imag ( aij.i); A.coeffRef(i, j).real ( aij.r);
               A.coeffRef(i, j - 1).imag ( aijm1.i); A.coeffRef(i, j - 1).real ( aijm1.r);
               A.coeffRef(i - 1, j).imag ( aim1j.i); A.coeffRef(i - 1, j).real ( aim1j.r);
            }
         }
         else
         {
            dd = ccc*dx*dx*data[(j - 1)*xsize + i] + ccc*dy*dy*data[j*xsize + i - 1] +
               (dx*dx*dy*dy - 2 * ccc*(dx*dx + dy*dy))*data[j*xsize + i] +
               ccc*dx*dx*data[(j + 1)*xsize + i] + ccc*dx*dx*data[j*xsize + i + 1];

            aij = (dx*dx*dy*dy - 2 * aaa*(dx*dx + dy*dy));
            aim1j = aaa * dy*dy;
            aijm1 = aaa * dx*dx;
            aia1j = aaa * dy*dy;
            aija1 = aaa * dx*dx;

            A.coeffRef(i, j).imag(aij.i); A.coeffRef(i, j).real(aij.r);
            A.coeffRef(i - 1, j).imag( aim1j.i); A.coeffRef(i - 1, j).real( aim1j.r);
            A.coeffRef(i, j + 1).imag ( aija1.i); A.coeffRef(i, j + 1).real ( aija1.r);
            A.coeffRef(i + 1, j).imag ( aia1j.i); A.coeffRef(i + 1, j).real ( aia1j.r);
         }
         d[j*xsize + i].imag (dd.i); d[j*xsize + i].real ( dd.r);
      }
   }

   setNbThreads(7);

   // BiCGSTAB solver
   BiCGSTAB<SparseMatrix<std::complex<float> >, IncompleteLUT<std::complex<float>  >  > solver;
   solver.preconditioner().setDroptol(0.001);
   solver.compute(A);
   x = solver.solve(d);


//   if (3 <= verbose )
//   {
//      float err = solver.error();
//      sf_warning("#iterations: %d", solver.iterations());
//      sf_warning("estimated error: %f ", err );
//   }

   for (int i = 0; i < xsize; i++)
   {
      for (int j = 0; j < ysize; j++)
      {
         data[j*xsize + i].i = x[j*xsize + i].imag(); data[j*xsize + i].r = x[j*xsize + i].real();
      }
   }

   return;
}

The header file is

Code: Select all

#ifndef WIMPLICITFD3D_H
#define WIMPLICITFD3D_H
#include <rsf.hh>

#include "sucomplex.h"
#include "wconstant.h"

#include "Eigen/IterativeLinearSolvers"
using namespace Eigen;

void implicitFinite3D(sucomplex *data, int xsize, int ysize, float dx, float dy, float *vel
   , float dz, float v0, float omega, int verbose);

#endif

And I get this running error :
eigen3.2.7/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h:171: const Eigen::internal::solve_retval<Eigen::IncompleteLUT<_Scalar>, Rhs> Eigen::IncompleteLUT<_Scalar>::solve(const Eigen::MatrixBase<OtherDerived>&) const [with Rhs = Eigen::Matrix<std::complex<float>, -1, 1>; _Scalar = std::complex<float>]: Assertion `m_isInitialized && "IncompleteLUT is not initialized."' failed.

I don't why. The BiCGSTAB solver can't solve the complex matrix? Please help me.
Thanks a lot!!!

You probably need to initialize the preconditioner before using it.
Something like:

Code: Select all

solver.preconditioner().setDroptol(0.001);
solver.preconditioner().compute(A); // should also check solver.preconditioner().info()
solver.compute(A);


I guess if you call solver.compute(A) several times you might not have to call it again as long as the matrix A is similar to the previous calls.

P.S.- I did not try to compile the code.

solver.compute(A) essentially calls "preconditioner.compute(A)", so the propblem is more likely that the ILUT preconditioner failed. This can check with:

if(solver.info()==Eigen::Success) ...

after solver.compute(A).

Could it be that you matrix contains empty rows or columns?