viennacl with complex numbers

Question

I'm using viennacl to solve a linear system of equations (AX = B) with the graphic card. Also, the code uses armadillo.

My system of equations have complex numbers. So the question is: Can I solve a system of equations (with complex numbers) using viennacl?

Above is an example of a working code with real numbers.

// System headers
#include <iostream>

// Armadillo headers (disable BLAS and LAPACK to avoid linking issues)
#define ARMA_DONT_USE_BLAS
#define ARMA_DONT_USE_LAPACK
#include <armadillo>
#include <complex>

#define VIENNACL_WITH_ARMADILLO 1

// ViennaCL headers
#include "viennacl/linalg/cg.hpp"
#include "viennacl/linalg/bicgstab.hpp"
#include "viennacl/linalg/gmres.hpp"
#include "viennacl/io/matrix_market.hpp"
#include "vector-io.hpp"

//using namespace arma;
using namespace viennacl::linalg;
using namespace std;

typedef arma::mat armat;
typedef arma::vec arvec;

typedef complex<double> dcmplx;

int main(void)
{

    int N = 500;
    armat A(N,N);
    A.randu();
    arvec B(N);
    B.randu();
    arvec X(N);
    arvec residual(N);

    viennacl::matrix<double> vcl_A(N, N);
    viennacl::vector<double> vcl_B(N);
    viennacl::vector<double> vcl_X(N);
    viennacl::vector<double> vcl_result(N);

    viennacl::copy(A, vcl_A);
    viennacl::copy(B, vcl_B);
    viennacl::copy(X, vcl_X);

    std::cout << "----- Running GMRES -----" << std::endl;
    vcl_X = viennacl::linalg::solve(vcl_A, vcl_B,               viennacl::linalg::gmres_tag());

    viennacl::copy(vcl_A, A);
    viennacl::copy(vcl_B, B);
    viennacl::copy(vcl_X, X);

    residual = A * X - B;
    cout << "Relative residual: " << norm(residual) / norm(B) << endl;
}

Complex version of the code:

#include <iostream>

// Armadillo headers (disable BLAS and LAPACK to avoid linking issues)
#define ARMA_DONT_USE_BLAS
#define ARMA_DONT_USE_LAPACK
#include <armadillo>
#include <complex>

#define VIENNACL_WITH_ARMADILLO 1

// ViennaCL headers
#include "viennacl/linalg/cg.hpp"
#include "viennacl/linalg/bicgstab.hpp"
#include "viennacl/linalg/gmres.hpp"
#include "viennacl/io/matrix_market.hpp"
#include "vector-io.hpp"


//using namespace arma;
using namespace viennacl::linalg;
using namespace std;

typedef arma::cx_mat armat;
typedef arma::cx_vec arvec;

typedef complex<double> dcmplx;

int main(void)
{

int N = 500;
armat A(N,N);
A.randu();
arvec B(N);
B.randu();
arvec X(N);
arvec residual(N);

viennacl::matrix<dcmplx> vcl_A(N, N);
viennacl::vector<dcmplx> vcl_B(N);
viennacl::vector<dcmplx> vcl_X(N);
viennacl::vector<dcmplx> vcl_result(N);

viennacl::copy(A, vcl_A);
viennacl::copy(B, vcl_B);
viennacl::copy(X, vcl_X);

std::cout << "----- Running GMRES -----" << std::endl;
vcl_X = viennacl::linalg::solve(vcl_A, vcl_B, viennacl::linalg::gmres_tag());

viennacl::copy(vcl_A, A);
viennacl::copy(vcl_B, B);
viennacl::copy(vcl_X, X);

residual = A * X - B;
cout << "Relative residual: " << norm(residual) / norm(B) << endl;



std::cout << "----- Running BiCGStab -----" << std::endl;
vcl_X = viennacl::linalg::solve(vcl_A, vcl_B, viennacl::linalg::bicgstab_tag());

viennacl::copy(vcl_A, A);
viennacl::copy(vcl_B, B);
viennacl::copy(vcl_X, X);

residual = A * X - B;
cout << "Relative residual: " << norm(residual) / norm(B) << endl;

std::cout << "----- Running CG -----" << std::endl;
vcl_X = viennacl::linalg::solve(vcl_A, vcl_B, viennacl::linalg::cg_tag());

viennacl::copy(vcl_A, A);
viennacl::copy(vcl_B, B);
viennacl::copy(vcl_X, X);

residual = A * X - B;
cout << "Relative residual: " << norm(residual) / norm(B) << endl;

}

Answer 1

ViennaCL currently does not support complex numbers. The primary technical reason is that OpenCL does not natively provide support for complex numbers. While emulating complex arithmetic via real arithmetic is certainly possible, we were reluctant to go down this path and (wrongly?) hoped that a standardization for complex will come soon.