具有复数的维也纳

Question

我正在使用viennacl用图形卡求解线性方程组(AX=B(。此外，代码使用了armadillo。

我的方程组有复数。所以问题是：我能用viennacl解一个方程组(复数(吗？

上面是一个使用实数的工作代码示例。

// 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;
}

复杂版本的代码：

#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目前不支持复数。主要的技术原因是OpenCL本身并没有提供对复数的支持。虽然通过真实的算术来模拟复杂的算术当然是可能的，但我们不愿意走这条路，(错误地？(希望复杂的标准化很快就会到来。