Solving a complex sparse linear system

[marco-butz]

Hi together,

I'm trying to solve a complex sparse linear system using the following code using Ginkgo 1.8.0:

#include <iostream>

#include <ginkgo/ginkgo.hpp>

int main(int argc, char* argv[])
{
    using GlobalIndexType = gko::int64;
    using LocalIndexType = gko::int32;
    using ValueType = std::complex<double>;  // Using complex128 precision
    using RealValueType = gko::remove_complex<ValueType>;
    using vec = gko::matrix::Dense<ValueType>;
    using real_vec = gko::matrix::Dense<RealValueType>;
    using solver = gko::solver::Cg<ValueType>;

    auto gpu = gko::CudaExecutor::create(0, gko::OmpExecutor::create());
    std::cout << "reading A" << std::endl;
    auto input_A = std::ifstream(argv[1]);
    auto A = gko::share(gko::read<gko::matrix::Csr<ValueType, GlobalIndexType>>(input_A, gpu));
    std::cout << "reading b" << std::endl;
    auto input_b = std::ifstream(argv[2]);
    auto b = gko::read<gko::matrix::Dense<ValueType>>(input_b, gpu);
    std::cout << "reading x" << std::endl;
    auto input_x = std::ifstream(argv[3]);
    auto x = gko::read<gko::matrix::Dense<ValueType>>(input_x, gpu);
    std::cout << "done reading matrices" << std::endl;

    write(std::cout, A);
    write(std::cout, b);
    write(std::cout, x);

    std::shared_ptr<const gko::log::Convergence<ValueType>> logger =
        gko::log::Convergence<ValueType>::create();

    auto solver_obj =
        solver::build()
            .with_criteria(
                gko::stop::Iteration::build().with_max_iters(10000u),
                gko::stop::ResidualNorm<>::build().with_reduction_factor(1e-25))
            .on(gpu);

    solver_obj->add_logger(logger);

    solver_obj->generate(A)->apply(b, x);

    auto res_norm = gko::clone(gpu->get_master(),
                               gko::as<real_vec>(logger->get_residual_norm()));

    std::cout << "\nFinal Res norm: " << res_norm->at(0, 0)
             << "\nIteration count: " << logger->get_num_iterations()
             << std::endl;

    auto x_master = gko::clone(gpu->get_master(), x);

    auto output_file = std::ofstream("/tmp/result.mtx");
    gko::write(output_file, gko::as<vec>(x_master.get()), gko::layout_type::array);
    output_file.close();

    write(std::cout, gko::as<vec>(x_master.get()));
}

My problem is, that the solution diverges when I use complex matrices. Testing it with real-valued matrices worked just fine.
The output for complex valued matrices:

./fdfdGinkgoStandalone test_mtx/simple_complex/A.mtx test_mtx/simple_complex/b.mtx test_mtx/simple_complex/x0.mtx
reading A
reading b
reading x
done reading matrices
%%MatrixMarket matrix coordinate complex general
4 4 16
1 1 2 1
1 2 3 -1
1 3 4 0
1 4 0 -1
2 1 3 0
2 2 1 2
2 3 5 -1
2 4 2 0
3 1 1 -1
3 2 -4 0
3 3 2 1
3 4 3 0
4 1 2 0
4 2 3 2
4 3 -1 1
4 4 4 1
%%MatrixMarket matrix array complex general
4 1
1 2
2 -1
3 1
4 -3
%%MatrixMarket matrix array complex general
4 1
0 0
0 0
0 0
0 0

Final Res norm: 2.14549e+16
Iteration count: 10000
%%MatrixMarket matrix array complex general
4 1
6.61227e+14 1.24639e+15
-2.85568e+15 -4.80371e+14
-4.46031e+14 7.13731e+14
8.18186e+14 -8.93078e+14

The output for real valued matrices:

./fdfdGinkgoStandalone test_mtx/simple/A.mtx test_mtx/simple/b.mtx test_mtx/simple/x0.mtxreading A
reading b
reading x
done reading matrices
%%MatrixMarket matrix coordinate complex general
19 19 147
1 1 3 0
1 2 1 0
1 3 1 0
1 5 1 0
1 9 1 0
1 17 1 0
2 1 1 0
2 2 5 0
2 3 1 0
2 4 1 0
2 6 1 0
2 10 1 0
2 18 1 0
3 1 1 0
3 2 1 0
3 3 7 0
3 4 1 0
3 5 1 0
3 7 1 0
3 11 1 0
3 19 1 0
4 2 1 0
4 3 1 0
4 4 11 0
4 5 1 0
4 6 1 0
4 8 1 0
4 12 1 0
5 1 1 0
5 3 1 0
5 4 1 0
5 5 13 0
5 6 1 0
5 7 1 0
5 9 1 0
5 13 1 0
6 2 1 0
6 4 1 0
6 5 1 0
6 6 17 0
6 7 1 0
6 8 1 0
6 10 1 0
6 14 1 0
7 3 1 0
7 5 1 0
7 6 1 0
7 7 19 0
7 8 1 0
7 9 1 0
7 11 1 0
7 15 1 0
8 4 1 0
8 6 1 0
8 7 1 0
8 8 23 0
8 9 1 0
8 10 1 0
8 12 1 0
8 16 1 0
9 1 1 0
9 5 1 0
9 7 1 0
9 8 1 0
9 9 29 0
9 10 1 0
9 11 1 0
9 13 1 0
9 17 1 0
10 2 1 0
10 6 1 0
10 8 1 0
10 9 1 0
10 10 31 0
10 11 1 0
10 12 1 0
10 14 1 0
10 18 1 0
11 3 1 0
11 7 1 0
11 9 1 0
11 10 1 0
11 11 37 0
11 12 1 0
11 13 1 0
11 15 1 0
11 19 1 0
12 4 1 0
12 8 1 0
12 10 1 0
12 11 1 0
12 12 41 0
12 13 1 0
12 14 1 0
12 16 1 0
13 5 1 0
13 9 1 0
13 11 1 0
13 12 1 0
13 13 43 0
13 14 1 0
13 15 1 0
13 17 1 0
14 6 1 0
14 10 1 0
14 12 1 0
14 13 1 0
14 14 47 0
14 15 1 0
14 16 1 0
14 18 1 0
15 7 1 0
15 11 1 0
15 13 1 0
15 14 1 0
15 15 53 0
15 16 1 0
15 17 1 0
15 19 1 0
16 8 1 0
16 12 1 0
16 14 1 0
16 15 1 0
16 16 59 0
16 17 1 0
16 18 1 0
17 1 1 0
17 9 1 0
17 13 1 0
17 15 1 0
17 16 1 0
17 17 61 0
17 18 1 0
17 19 1 0
18 2 1 0
18 10 1 0
18 14 1 0
18 16 1 0
18 17 1 0
18 18 67 0
18 19 1 0
19 3 1 0
19 11 1 0
19 15 1 0
19 17 1 0
19 18 1 0
19 19 71 0
%%MatrixMarket matrix array complex general
19 1
1 0
1 0
1 0
1 0
1 0
1 0
1 0
1 0
1 0
1 0
1 0
1 0
1 0
1 0
1 0
1 0
1 0
1 0
1 0
%%MatrixMarket matrix array complex general
19 1
0 0
0 0
0 0
0 0
0 0
0 0
0 0
0 0
0 0
0 0
0 0
0 0
0 0
0 0
0 0
0 0
0 0
0 0
0 0

Final Res norm: 3.48709e-27
Iteration count: 38
%%MatrixMarket matrix array complex general
19 1
0.252218 0
0.108645 0
0.0662811 0
0.0630433 0
0.0384088 0
0.0396536 0
0.0402648 0
0.0338935 0
0.0193098 0
0.0234653 0
0.0211499 0
0.0196413 0
0.0199151 0
0.0181674 0
0.0162722 0
0.0150714 0
0.0107016 0
0.0121141 0
0.0123025 0

I attached the mtx files I used for testing. Does anyone know what is the problem here?

real.zip
complex.zip

[MarcelKoch]

Hi @marco-butz, I tried running your complex matrix with our solver benchmark and it also didn't converge with CG. It did converge however when I used GMRES. Did you check if your complex matrix is supposed to be solvable with CG, e.g. in matlab or similar. For CG with complex matrices it's not enough if the matrix is symmetric, it has to be hermitian.

[marco-butz]

Can you share the exact code that you use here? With the setup supplied in setup.zip above and reduction factor 1e-15 I still get divergence:

./fdfdGinkgoStandalone test_mtx/fdfd_2D/A.mtx test_mtx/fdfd_2D/b.mtx test_mtx/fdfd_2D/x0.mtx test_mtx/fdfd_2D/result.mtx
reading A
reading b
reading x
done reading matrices

Final Res norm: 8.12927e-05
Iteration count: 10000

My exact code is here:

#include <iostream>

#include <ginkgo/ginkgo.hpp>

int main(int argc, char* argv[])
{
    using GlobalIndexType = gko::int64;
    using LocalIndexType = gko::int32;
    using ValueType = std::complex<double>;  // Using complex128 precision
    using RealValueType = gko::remove_complex<ValueType>;
    using vec = gko::matrix::Dense<ValueType>;
    using real_vec = gko::matrix::Dense<RealValueType>;
    using solver = gko::solver::Gmres<ValueType>;

    auto gpu = gko::CudaExecutor::create(0, gko::OmpExecutor::create());
    std::cout << "reading A" << std::endl;
    auto input_A = std::ifstream(argv[1]);
    auto A = gko::share(gko::read<gko::matrix::Csr<ValueType, GlobalIndexType>>(input_A, gpu));
    std::cout << "reading b" << std::endl;
    auto input_b = std::ifstream(argv[2]);
    auto b = gko::read<gko::matrix::Dense<ValueType>>(input_b, gpu);
    std::cout << "reading x" << std::endl;
    auto input_x = std::ifstream(argv[3]);
    auto x = gko::read<gko::matrix::Dense<ValueType>>(input_x, gpu);
    std::cout << "done reading matrices" << std::endl;

    std::shared_ptr<const gko::log::Convergence<ValueType>> logger =
        gko::log::Convergence<ValueType>::create();

    auto solver_obj =
        solver::build()
            .with_criteria(
                gko::stop::Iteration::build().with_max_iters(10000u),
                gko::stop::ResidualNorm<>::build().with_reduction_factor(1e-15))
            .on(gpu);

    solver_obj->add_logger(logger);

    solver_obj->generate(A)->apply(b, x);

    auto res_norm = gko::clone(gpu->get_master(),
                               gko::as<real_vec>(logger->get_residual_norm()));

    std::cout << "\nFinal Res norm: " << res_norm->at(0, 0)
             << "\nIteration count: " << logger->get_num_iterations()
             << std::endl;

    auto x_master = gko::clone(gpu->get_master(), x);

    auto output_file = std::ofstream("/tmp/result.mtx");
    gko::write(output_file, gko::as<vec>(x_master.get()), gko::layout_type::array);
    output_file.close();
}

As a reference I used scipy, which solves the problem correctly:

import sys
import numpy as np
from scipy.io import mmread, mmwrite
from scipy.sparse.linalg import spsolve
from scipy.sparse import csr_matrix

def main():
    if len(sys.argv) != 5:
        print("Usage: python solve_system.py <matrix_A.mtx> <vector_b.mtx> <vector_x.mtx> <solution_x.mtx>")
        sys.exit(1)

    # Read the arguments
    matrix_file = sys.argv[1]
    vector_b_file = sys.argv[2]
    vector_x_file = sys.argv[3]
    solution_file = sys.argv[4]

    # Read matrix A (complex128) from file in Matrix Market format
    print("Reading matrix A from", matrix_file)
    A = mmread(matrix_file).tocsc()  # Assuming Matrix Market format for sparse matrices
    A = A.astype(np.complex128)  # Ensure it's complex128

    # Read vector b (complex128) from file in Matrix Market format
    print("Reading vector b from", vector_b_file)
    b = mmread(vector_b_file)
    b = b.astype(np.complex128)
    b = np.ravel(b)  # Convert from a sparse matrix to a 1D array, in case it's stored as a matrix

    # Read initial guess vector x (complex128) from file in Matrix Market format
    print("Reading vector x from", vector_x_file)
    x = mmread(vector_x_file)
    x = x.astype(np.complex128)
    x = np.ravel(x)  # Convert from a sparse matrix to a 1D array, if needed

    # Solve the system A*x = b using spsolve
    print("Solving the system A*x = b")
    x_solution = spsolve(A, b)

    # Print the solution
    print("Solution x:")
    print(x_solution)

    # Save the solution vector to file in Matrix Market format
    print("Writing solution to", solution_file)
    x_solution_2d = x_solution[:, np.newaxis]  # Convert to 2D array for Matrix Market format
    mmwrite(solution_file, x_solution_2d)  # Save as a sparse matrix in Matrix Market format

if __name__ == "__main__":
    main()

I again attached the setup and the solution vectors obtained by ginkgo (incorrect) and scipy(correct).

Uploading FDFD_setup_solutions_scipy_ginkgo.zip…

Any help would be very much appreciated.

[marco-butz]

Upload didn't work:

FDFD_setup_solutions_scipy_ginkgo.zip

[marco-butz]

I also have a MAGMA-based solver, which solves the problem correctly (same solution as scipy):

#define EIGEN_USE_MKL_ALL

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <iostream>
#include <math.h>
#include <chrono>
#include <Eigen/SparseCholesky>
#include <Eigen/SparseLU>
#include <Eigen/Sparse>
#include <Eigen/Dense>
// #include <Eigen/PardisoSupport>
// #include <Eigen/IterativeLinearSolvers>
#include <unsupported/Eigen/SparseExtra>

#include <magma_v2.h>
#include <magmasparse.h>

#include <complex>
#include <vector>

class SimBridge {
  int maxIterations;
  Eigen::VectorXcd b;
  //Eigen::PardisoLU<Eigen::SparseMatrix<std::complex<double>, Eigen::ColMajor>> directSolver;
  Eigen::BiCGSTAB<Eigen::SparseMatrix<std::complex<double>, Eigen::ColMajor>> iterativeCPUSolver;

  magma_device_t device;
  magma_zopts opts;
  magma_queue_t queue;

  magma_z_matrix b_mag={Magma_DENSE}, db_mag={Magma_DENSE};

public:
  bool twoD;
  SimBridge(bool twoD, int maxIters, float relativeTargetResidual) {
    this->twoD = twoD;
    this->maxIterations = maxIters;

    if(!this->twoD) {
      std::cout << "Initializing MAGMA now!" << std::endl;
      int err = magma_init();
      std::cout << magma_strerror(err) << std::endl;
      magma_getdevice(&device);
      std::cout << "magma device: " << magma_getdevice_arch() << std::endl;
      std::cout << "magma device multiprocessor count: " << magma_getdevice_multiprocessor_count() << std::endl;
      magma_queue_create( device, &queue );

      memset(&opts, 0, sizeof(magma_zopts));

      opts.solver_par.solver = Magma_BICGSTAB; // Use BiCGSTAB solver
      opts.solver_par.restart = 10;            // Restart after this number of iterations
      opts.solver_par.maxiter = maxIters;          // Maximum number of iterations
      opts.solver_par.rtol = relativeTargetResidual;             // Relative tolerance for convergence
      //opts.solver_par.atol = 1e-4;            // Absolute tolerance for convergence
      opts.solver_par.verbose = 100;             // Verbosity level (increase for more output)
      opts.solver_par.num_eigenvalues = 0;     // Number of eigenvalues for generalized problems (set to 0 for non-generalized)
      opts.solver_par.info = 0;                // Solver information (useful for debugging)

      opts.precond_par.solver = Magma_NONE;    // No preconditioner by default
      opts.precond_par.levels = 0;             // Number of levels for preconditioner (if applicable)
      opts.precond_par.trisolver = Magma_CG;   // Inner iterative solver for preconditioner (if applicable)

      magma_zsolverinfo_init(&opts.solver_par, &opts.precond_par, queue);
    }

    std::cout << "SimBridge has been initialized" << std::endl;
  }

  void deinit() {
    std::cout << "Deinitializing simulation bridge." << std::endl;

    if(!this->twoD) {
      magma_zmfree(&db_mag, queue);
      magma_zmfree(&b_mag, queue);  // won't do anything as MAGMA does not own the data.
      magma_queue_destroy(queue);
      magma_finalize();
    }

    std::cout << "SimBridge has been destroyed" << std::endl;
  }

  ~SimBridge() {

  }

  void setSource(Eigen::VectorXcd sourceVec) {
    if(!this->twoD) {
      std::cout << "setting source vector" << std::endl;
      magma_zvset(sourceVec.rows(), 1, reinterpret_cast <magmaDoubleComplex*> (sourceVec.data()), &b_mag, queue);
      magma_z_vtransfer(b_mag, &db_mag, Magma_CPU, Magma_DEV, queue);
    } else {
      this->b = sourceVec;
    }
    //std::cout << "Done setting sources." << std::endl;
  }

  Eigen::VectorXcd solveMagma(Eigen::SparseMatrix<std::complex<double>, Eigen::RowMajor> A, Eigen::VectorXcd initialGuess, bool verbose);
  //Eigen::VectorXcd solveEigenDirect(Eigen::SparseMatrix<std::complex<double>, Eigen::ColMajor> A);
};

Eigen::VectorXcd SimBridge::solveMagma(Eigen::SparseMatrix<std::complex<double>, Eigen::RowMajor> A, Eigen::VectorXcd initialGuess, bool verbose) {
  if(this->twoD) {
    std::cout << "Iterative solver not supported for 2D!" << std::endl;
  }

  int err = 0;

  magma_z_matrix A_mag={Magma_CSR}, dA_mag={Magma_CSR};
  magma_z_matrix x_mag={Magma_DENSE};//, dx_mag={Magma_DENSE};

  err = magma_zcsrset(A.rows(), A.cols(), A.outerIndexPtr(), A.innerIndexPtr(), reinterpret_cast <magmaDoubleComplex*> (A.valuePtr()), &A_mag, queue);

  // Copy the system to the device
  err = magma_z_mtransfer(A_mag, &dA_mag, Magma_CPU, Magma_DEV, queue);

  // initialize an initial guess for the iteration vector
  magma_z_matrix initial_mag={Magma_DENSE}, dx_mag={Magma_DENSE};
  magma_zvset(initialGuess.size(), 1, reinterpret_cast <magmaDoubleComplex*> (initialGuess.data()), &initial_mag, queue);
  magma_z_vtransfer(initial_mag, &dx_mag, Magma_CPU, Magma_DEV, queue);

  err = magma_z_solver(dA_mag, db_mag, &dx_mag, &opts, queue);
  magma_zmtransfer(dx_mag, &x_mag, Magma_DEV, Magma_CPU, queue);

  if (verbose) {
    std::cout << "Solver status: " << magma_strerror(err) << std::endl;
    std::cout << "Solver iterations: " << opts.solver_par.numiter << std::endl;
    std::cout << "Solver initial residual: " << opts.solver_par.init_res << std::endl;
    std::cout << "Solver final residual: " << opts.solver_par.final_res << std::endl;
    std::cout << "Solver runtime: " << opts.solver_par.runtime << std::endl;
  }

  magma_queue_sync(queue);

  // Free the allocated memory...
  magma_zmfree(&dx_mag, queue);
  magma_zmfree(&dA_mag, queue);
  magma_zmfree(&A_mag, queue);

  //construct eigen matrix again:
  Eigen::Map<Eigen::VectorXcd> vector(reinterpret_cast <std::complex<double>*> (x_mag.val), x_mag.num_rows);

  return vector.eval();
}

int main(int argc, char* argv[])
{
  SimBridge bridge(false, 10000, 0.00001);

  const char* fileNameA = argv[1];
  const char* fileNameB = argv[2];
  const char* fileNameX0 = argv[3];
  const char* fileNameResult = argv[4];

  Eigen::SparseMatrix<std::complex<double>, Eigen::RowMajor> A;
  Eigen::loadMarket(A, fileNameA);

  std::cout << "A loaded successfully. Size: "
              << A.rows() << " x " << A.cols() << std::endl;

  Eigen::VectorXcd b;
  Eigen::Matrix<std::complex<double>, Eigen::Dynamic, Eigen::Dynamic> matB;
  Eigen::loadMarketDense(matB, fileNameB);

  std::cout << "b loaded successfully. Size: "
              << matB.rows() << " x " << matB.cols() << std::endl;

  b = matB.col(0);

  Eigen::VectorXcd x0;
  Eigen::Matrix<std::complex<double>, Eigen::Dynamic, Eigen::Dynamic> matX;
  Eigen::loadMarketDense(matX, fileNameX0);

  std::cout << "x0 loaded successfully. Size: "
              << matX.rows() << " x " << matX.cols() << std::endl;

  x0 = matX.col(0);

  bridge.setSource(b);

  Eigen::VectorXcd solution;

  solution = bridge.solveMagma(A, x0, true);

  std::cout << "solution: "
              << solution.rows() << " x " << solution.cols() << std::endl;

  for (int i = 0; i < solution.size(); ++i) {
    std::cout << "solution[" << i << "]: " << solution[i] << std::endl;
  }

  Eigen::saveMarketDense(solution, fileNameResult);

  bridge.deinit();

  return 0;
}

[nbeams]

Hi @marco-butz, I noticed you are using Ginkgo 1.8.0. We recently fixed a bug for complex numbers in GMRES (#1678). If I solve your problem with a Ginkgo version older than that PR, I see output similar to what you posted:

1.52026e-05
Number of iterations 10000
Convergence status false

But if I add the commit from the PR:

1.22981e-19
Number of iterations 2974
Convergence status true

You can cherry-pick that one commit (aa09fbd), or make the change manually -- it's just one line of code. Hopefully that will fix your problem. You could also update to the latest develop branch of Ginkgo and get other improvements, too, but Ginkgo is planning a 1.9.0 release soon, if you want to wait for that.

[yhmtsai]

I use the latest code from ginkgo develop, which includes the fix as @nbeams mentioned.

From your magma code, you seems to use Bicgstab not Cg and the relative residual norm check is 0.00001. Do I get it right?
Do you mean solving it correctly by Cg in Magma with 1e-25, too?
If not, maybe also try Bicgstab with relaxed residual norm requirement in ginkgo rather than Cg because Cg only works for hermitian matrix as @MarcelKoch mentioned.

[marco-butz]

Hello Guys, actually the issue mentioned by @nbeams fixed the error. Thank you very much!

[yhmtsai]

I have checked with python scipy cg on 4x4 complex matrix.
I got the similar result with ginkgo from your report

[ 6.59958364e+14 +1.24504733e+15j 
 -2.85191667e+15 -4.80762466e+14j
 -4.45717893e+14 +7.12675465e+14j 
 8.17466744e+14 -8.91666699e+14j]

could you share the solution from your CG with MKL?