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;

}

Did you try using it with [`std::complex`](http://en.cppreference.com/w/cpp/numeric/complex)? Did you get any particular errors or problems, we might be able to help you with? As is your question seems to be overly broad. — πάντα ῥεῖ, Sep 06 '15 at 21:16
Yes I did. I include the example above. I don't know if someone here can help me, but I will try anyway. ;) — nunodsousa, Sep 06 '15 at 21:26
And what were the errors/problems you observed? Put everything verbatim in your question please. — πάντα ῥεῖ, Sep 06 '15 at 21:38

score 1 · Accepted Answer · answered Sep 07 '15 at 07:09

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.

viennacl with complex numbers

1 Answers1