SHOGUN
v3.2.0
|
00001 /* 00002 * This program is free software; you can redistribute it and/or modify 00003 * it under the terms of the GNU General Public License as published by 00004 * the Free Software Foundation; either version 3 of the License, or 00005 * (at your option) any later version. 00006 * 00007 * Written (W) 2013 Soumyajit De 00008 */ 00009 00010 #include <shogun/lib/config.h> 00011 00012 #ifdef HAVE_EIGEN3 00013 #include <shogun/lib/SGVector.h> 00014 #include <shogun/lib/SGSparseMatrix.h> 00015 #include <shogun/mathematics/eigen3.h> 00016 #include <shogun/mathematics/linalg/linop/SparseMatrixOperator.h> 00017 #include <shogun/mathematics/linalg/linsolver/DirectSparseLinearSolver.h> 00018 00019 using namespace Eigen; 00020 00021 namespace shogun 00022 { 00023 00024 CDirectSparseLinearSolver::CDirectSparseLinearSolver() 00025 : CLinearSolver<float64_t, float64_t>() 00026 { 00027 } 00028 00029 CDirectSparseLinearSolver::~CDirectSparseLinearSolver() 00030 { 00031 } 00032 00033 SGVector<float64_t> CDirectSparseLinearSolver::solve( 00034 CLinearOperator<float64_t>* A, SGVector<float64_t> b) 00035 { 00036 REQUIRE(A, "Operator is NULL!\n"); 00037 REQUIRE(A->get_dimension()==b.vlen, "Dimension mismatch!\n"); 00038 CSparseMatrixOperator<float64_t>* op 00039 =dynamic_cast<CSparseMatrixOperator<float64_t>*>(A); 00040 REQUIRE(op, "Operator is not SparseMatrixOperator type!\n"); 00041 00042 // creating eigen3 Sparse Matrix 00043 SGSparseMatrix<float64_t> sm=op->get_matrix_operator(); 00044 typedef SparseMatrix<float64_t> MatrixType; 00045 const MatrixType& m=EigenSparseUtil<float64_t>::toEigenSparse(sm); 00046 00047 // creating eigen3 maps for vectors 00048 SGVector<float64_t> x(m.cols()); 00049 x.set_const(0.0); 00050 Map<VectorXd> map_x(x.vector, x.vlen); 00051 Map<VectorXd> map_b(b.vector, b.vlen); 00052 00053 // using LLT to solve the system Ax=b 00054 SimplicialLLT<MatrixType> llt; 00055 llt.compute(m); 00056 map_x=llt.solve(map_b); 00057 00058 // checking for success 00059 if (llt.info()==NumericalIssue) 00060 SG_WARNING("Matrix is not Hermitian positive definite!\n"); 00061 00062 return x; 00063 } 00064 00065 } 00066 #endif // HAVE_EIGEN3