viennacl complex numbers

[Nuno de Sousa]

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 veinal?

Above is an example of a working code with real numbers and a nonworking example with complex 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 number example

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;

}

[Karl Rupp]

Unfortunately, ViennaCL currently (version 1.7.1) does not support complex numbers.