[问题]:有人可以帮助我在 Eigen-3.3.5 上使用 DiagonalPreconditioner 包装器来实现无矩阵求解器吗？
我正在尝试修改下面链接上的示例，以便我可以将 DiagonalPreconditioner 与 CG 稀疏交互式求解器结合使用。
[特征3.3.5] https://eigen.tuxfamily.org/dox/group__MatrixfreeSolverExample.html
在 Eigen-3.2.10 上，他们为预处理器提供了包装器，但我在修改此示例的求解器部分时陷入困境，并由于两个版本之间的重大变化而放弃使用它。 [特征3.2.10] http://eigen.tuxfamily.org/dox-3.2/group__MatrixfreeSolverExample.html
[背景]:我正在 Mex 上编译 cg 求解器，以便可以在 Matlab 上使用它。原因是，对于我现在正在工作的线性系统(LS)的大小，特征求解器比Matlab的pcg求解器更有效，250k x 250k(它可能大到1e6 x 1e6，但幸运的是它非常稀疏) 。此外，LS 的右侧(RHS)和左侧(LHS)的值每次迭代都会改变，但稀疏模式是固定的。出于这个原因，我认为 Eigen 将是一个不错的选择，因为可以分解“计算”步骤，如下所述:https://eigen.tuxfamily.org/dox/group__TopicSparseSystems.html
[操作系统]:CentOS 7
[代码]:如果有人注释与无矩阵求解器相关的部分并取消注释所指示的部分(我相信它已被很好地指示)，则下面的代码将起作用。
#include "mex.h"
#include "math.h"
#include "matrix.h"
#include <Eigen/Sparse>
#include <Eigen/Dense>
#include <Eigen/Core>
#include <vector>
#include <algorithm>
#include <iostream>
#include <omp.h>
#include <fstream>
//using namespace Eigen;
using namespace std;
typedef Eigen::Triplet<double> T;
class MatrixReplacement;
using Eigen::SparseMatrix;
namespace Eigen {
namespace internal {
// MatrixReplacement looks-like a SparseMatrix, so let's inherits its traits:
template<>
struct traits<MatrixReplacement> : public Eigen::internal::traits<Eigen::SparseMatrix<double> >
{};
}
}
// Example of a matrix-free wrapper from a user type to Eigen's compatible type
// For the sake of simplicity, this example simply wrap a Eigen::SparseMatrix.
class MatrixReplacement : public Eigen::EigenBase<MatrixReplacement> {
public:
// Required typedefs, constants, and method:
typedef double Scalar;
typedef double RealScalar;
typedef int StorageIndex;
enum {
ColsAtCompileTime = Eigen::Dynamic,
MaxColsAtCompileTime = Eigen::Dynamic,
IsRowMajor = false
};
Index rows() const { return mp_mat->rows(); }
Index cols() const { return mp_mat->cols(); }
//Index outerSize() const { return mp_mat->outerSize(); }
template<typename Rhs>
Eigen::Product<MatrixReplacement,Rhs,Eigen::AliasFreeProduct> operator*(const Eigen::MatrixBase<Rhs>& x) const {
return Eigen::Product<MatrixReplacement,Rhs,Eigen::AliasFreeProduct>(*this, x.derived());
}
// Custom API:
MatrixReplacement() : mp_mat(0) {}
void attachMyMatrix(const SparseMatrix<double> &mat) {
mp_mat = &mat;
}
const SparseMatrix<double> my_matrix() const { return *mp_mat; }
private:
const SparseMatrix<double> *mp_mat;
};
// Implementation of MatrixReplacement * Eigen::DenseVector though a specialization of internal::generic_product_impl:
namespace Eigen {
namespace internal {
template<typename Rhs>
struct generic_product_impl<MatrixReplacement, Rhs, SparseShape, DenseShape, GemvProduct> // GEMV stands for matrix-vector
: generic_product_impl_base<MatrixReplacement,Rhs,generic_product_impl<MatrixReplacement,Rhs> >
{
typedef typename Product<MatrixReplacement,Rhs>::Scalar Scalar;
template<typename Dest>
static void scaleAndAddTo(Dest& dst, const MatrixReplacement& lhs, const Rhs& rhs, const Scalar& alpha)
{
// This method should implement "dst += alpha * lhs * rhs" inplace,
// however, for iterative solvers, alpha is always equal to 1, so let's not bother about it.
assert(alpha==Scalar(1) && "scaling is not implemented");
EIGEN_ONLY_USED_FOR_DEBUG(alpha);
// Here we could simply call dst.noalias() += lhs.my_matrix() * rhs,
// but let's do something fancier (and less efficient):
for(Index i=0; i<lhs.cols(); ++i)
dst += rhs(i) * lhs.my_matrix().col(i);
}
};
}
}
// the gateway function
void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
//
Eigen::initParallel();
//input vars
double *A_ir;
double *A_jc;
double *A_val;
double *b;
double *x0;
int maxiter;
double tol;
double icholShift;
int nWorkers;
//output vars
double *x;
// long int numite;
// double resrel;
//temp vars
long int nrows;
long int nnz;
std::vector<T> tripletList;
//-----------------
// GetData
//-----------------
if(nrhs != 9){
mexErrMsgIdAndTxt("MEX:c_sol_pcg_eigen:rhs",
"This function takes too much input arguments.");
}
A_ir = mxGetPr(prhs[0]);
A_jc = mxGetPr(prhs[1]);
A_val = mxGetPr(prhs[2]);
b = mxGetPr(prhs[3]);
x0 = mxGetPr(prhs[4]);
maxiter = mxGetScalar(prhs[5]);
tol = mxGetScalar(prhs[6]);
icholShift = mxGetScalar(prhs[7]);
nrows = mxGetM(prhs[3]);
nnz = mxGetM(prhs[0]);
nWorkers = mxGetScalar(prhs[8]);
plhs[0] = mxCreateDoubleMatrix(nrows,1,mxREAL);
x = mxGetPr(plhs[0]);
Eigen::setNbThreads(nWorkers);
//-----------------
//calculations
//-----------------
//covert A_ir, A_jc, A_val to Eigen Sparse matrix
tripletList.reserve(nnz);
for(long int i=0; i<nnz; i++){
tripletList.push_back(T(A_ir[i]-1, A_jc[i]-1, A_val[i]));
}
Eigen::SparseMatrix<double> A_eigen(nrows, nrows);
A_eigen.setFromTriplets(tripletList.begin(), tripletList.end());
// call matrix-free API
MatrixReplacement AF_eigen; // <-comment here for matrix context
AF_eigen.attachMyMatrix(A_eigen); // <-comment here for matrix context
Eigen::Map<Eigen::VectorXd> b_eigen(b, nrows);
//CG solver of Eigen
Eigen::VectorXd x_eigen;
Eigen::Map<Eigen::VectorXd> x0_eigen(x0, nrows);
// solve problem
//Eigen::ConjugateGradient<SparseMatrix<double>, Eigen::Lower|Eigen::Upper,Eigen::DiagonalPreconditioner<double>> cg; // <-uncomment here for matrix context
// matrix free solver
Eigen::ConjugateGradient<MatrixReplacement, Eigen::Lower|Eigen::Upper,Eigen::DiagonalPreconditioner<double>> cg; // <-comment here for matrix context
cg.setTolerance(tol);
cg.setMaxIterations(maxiter);
//cg.compute(A_eigen); // <-uncomment here for matrix context
cg.compute(AF_eigen); // <-comment here for matrix context
x_eigen = cg.solveWithGuess(b_eigen, x0_eigen);
for(long int i=0; i<nrows; i++){
x[i] = x_eigen(i);
}
}
如果按原样编译，则会返回以下错误(我已删除本地计算机路径并替换为错误日志中的 path_to_eigen 或 path_to_mex_file):
Error using mex
In file included from
/path_to_eigen/include/eigen3/Eigen/IterativeLinearSolvers:39:0,
from
/path_to_eigen/include/eigen3/Eigen/Sparse:33,
from
/path_to_mex_file/mex/EigenPCGIcholOMPFREE.cpp:23:
/path_to_eigen/include/eigen3/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h:
In instantiation of ‘Eigen::DiagonalPreconditioner<_Scalar>&
Eigen::DiagonalPreconditioner<_Scalar>::factorize(const MatType&) [with MatType
= MatrixReplacement; _Scalar = double]’:
/path_to_eigen/include/eigen3/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h:84:27:
required from ‘Eigen::DiagonalPreconditioner<_Scalar>&
Eigen::DiagonalPreconditioner<_Scalar>::compute(const MatType&) [with MatType =
MatrixReplacement; _Scalar = double]’
/path_to_eigen/include/eigen3/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h:241:5:
required from ‘Derived& Eigen::IterativeSolverBase<Derived>::compute(const
Eigen::EigenBase<OtherDerived>&) [with MatrixDerived = MatrixReplacement;
Derived = Eigen::ConjugateGradient<MatrixReplacement, 3,
Eigen::DiagonalPreconditioner<double> >]’
/path_to_mex_file/EigenPCGIcholOMPFREE.cpp:175:24:
required from here
/path_to_eigen/include/eigen3/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h:68:26:
error: ‘const class MatrixReplacement’ has no member named ‘outerSize’
for(int j=0; j<mat.outerSize(); ++j)
~~~~^~~~~~~~~
/path_to_eigen/include/eigen3/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h:70:41:
error: no type named ‘InnerIterator’ in ‘class MatrixReplacement’
typename MatType::InnerIterator it(mat,j);
^~
/path_to_eigen/include/eigen3/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h:70:41:
error: no type named ‘InnerIterator’ in ‘class MatrixReplacement’
我能够通过添加来修复“outerSize”
索引outerSize() const { return mp_mat->outerSize(); } 下面 Index cols() const { return mp_mat->cols(); }。但是，我不知道如何解决 InnerIterator 问题。
我相信 Mex 不是问题，但这里是在 mex 上编译的 Matlab 文件:
MEXOPTS={'-v','-largeArrayDims','-DMEX','-DNDEBUG'};
MSSE42='CXXFLAGS=$CXXFLAGS -msse4.2';
STDCPP11='CXXFLAGS=$CXXFLAGS -fopenmp';
LDFLAGS = 'LDFLAGS=$LDFLAGS -fopenmp';
EIGEN_INC{1}='-I/path_to_eigen/include/eigen3'; // version 3.3.5
%% not sure if boost is need, I have being carrying this from other functions
BOOST_INC='-I/path_to_boost_include/include';
BOOST_LIB{1}='-L/path_to_local_boost_lib/boost/lib'; // version 1.67
BOOST_LIB{2}='-lboost_thread';
BOOST_LIB{3}='-lboost_system';
GMP_INC='-I/apps/cent7/gcc/6.3.0/gmp-6.1.0/include';
GMP_LIB{1}='-L/apps/cent7/gcc/6.3.0/gmp-6.1.0/lib';
GMP_LIB{2}='-lgmpxx';
GMP_LIB{3}='-lgmp';
MPFR_INC='-I/apps/cent7/gcc/6.3.0/mpfr-3.1.5/include';
MPFR_LIB{1}='-L/apps/cent7/gcc/6.3.0/mpfr-3.1.5/lib';
MPFR_LIB{2}='-lmpfr';
MPC_INC='-I/apps/cent7/gcc/6.3.0/mpc-1.0.3/include';
MPC_LIB{1}='-L/apps/cent7/gcc/6.3.0/mpc-1.0.3/lib';
MPC_LIB{2}='-lmpc';
mex( MEXOPTS{:}, MSSE42,STDCPP11,LDFLAGS,...
BOOST_INC, BOOST_LIB{:}, EIGEN_INC{:}, GMP_INC, GMP_LIB{:}, ...
MPFR_INC,MPFR_LIB{:}, MPC_INC, MPC_LIB{:}, ...
'EigenPCGIcholOMPFREE.cpp');
修改这个包装器对于我的“程序思维”来说太过分了。我希望你们能帮助我。
感谢您提前提供的帮助!
最佳答案
我知道这篇文章已经过时了，您可能不再需要答案，但我想在这里为像我这样通过 Google 找到这篇文章的人留下答案。我自己对 Eigen 库还是个新手，但希望以下解决方案比没有解决方案要好。
要实现对角预处理器，您显然需要某种方法从线性运算符中提取对角系数。如何完成此操作取决于无矩阵运算符类的实现细节，但更糟糕的是，您可以将基本 vector 传递给运算符(除了感兴趣的列中的 1 之外，全部为零)。这将返回线性运算符的矩阵表示形式的列 vector ，您可以从中提取对角线元素。这基本上相当于计算运算符的密集矩阵表示，只是不存储大部分值，因此除非没有其他选择，否则您真的不想这样做。根据这种方法的效率，最好只使用恒等预处理器。
查看Eigen::DiagonalPreconditioner的源代码，您会看到它使用MatrixReplacement::InnerIterator来提取对角系数，如下所示:
for(int j=0; j<mat.outerSize(); ++j)
{
typename MatType::InnerIterator it(mat,j);
while(it && it.index()!=j) ++it;
if(it && it.index()==j && it.value()!=Scalar(0))
m_invdiag(j) = Scalar(1)/it.value();
else
m_invdiag(j) = Scalar(1);
}
您需要在 MatrixReplacement 中创建一个迭代器类，它将返回运算符矩阵表示的指定行/列的对角线系数。我的实现如下:
class InnerIterator {
public:
InnerIterator(const MatrixReplacement& mat, Index row)
: mat(mat), has_val(false), row(row), col(row){ }
operator bool() { return row==col; }
InnerIterator& operator++(){
col++;
has_val = false;
return *this;
}
Index index(){ return col; }
Scalar value(){
if(!has_val){ // cache value since this function is called twice
stored_val = mat.diagonalCoefficient(row); // your implementation here
has_val = true;
}
return stored_val;
}
private:
const MatrixReplacement& mat;
bool has_val;
Index row, col;
Scalar stored_val;
};
关于c++ - Eigen-3.3.5 上的 DiagonalPreconditioner 包装器，用于无矩阵稀疏求解器 CG，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/51846432/