1

我正在使用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;

}
4

1 回答 1

1

ViennaCL 目前不支持复数。主要的技术原因是 OpenCL 本身并不支持复数。虽然通过实数算术模拟复数算术当然是可能的，但我们不愿意走这条路，并且（错误地？）希望复数的标准化很快就会到来。

于 2015-09-07T07:09:14.593 回答