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/SGMatrix.h> 00015 #include <shogun/mathematics/eigen3.h> 00016 #include <shogun/mathematics/linalg/linop/DenseMatrixOperator.h> 00017 #include <shogun/mathematics/linalg/linsolver/DirectLinearSolverComplex.h> 00018 00019 using namespace Eigen; 00020 00021 namespace shogun 00022 { 00023 00024 CDirectLinearSolverComplex::CDirectLinearSolverComplex() 00025 : CLinearSolver<complex128_t, float64_t>(), 00026 m_type(DS_QR_NOPERM) 00027 { 00028 SG_GCDEBUG("%s created (%p)\n", this->get_name(), this) 00029 } 00030 00031 CDirectLinearSolverComplex::CDirectLinearSolverComplex(EDirectSolverType type) 00032 : CLinearSolver<complex128_t, float64_t>(), 00033 m_type(type) 00034 { 00035 SG_GCDEBUG("%s created (%p)\n", this->get_name(), this) 00036 } 00037 00038 CDirectLinearSolverComplex::~CDirectLinearSolverComplex() 00039 { 00040 SG_GCDEBUG("%s destroyed (%p)\n", this->get_name(), this) 00041 } 00042 00043 SGVector<complex128_t> CDirectLinearSolverComplex::solve( 00044 CLinearOperator<complex128_t>* A, SGVector<float64_t> b) 00045 { 00046 SGVector<complex128_t> x(b.vlen); 00047 00048 REQUIRE(A, "Operator is NULL!\n"); 00049 REQUIRE(A->get_dimension()==b.vlen, "Dimension mismatch!\n"); 00050 00051 CDenseMatrixOperator<complex128_t> *op= 00052 dynamic_cast<CDenseMatrixOperator<complex128_t>*>(A); 00053 REQUIRE(op, "Operator is not CDenseMatrixOperator<complex128_t, float64_t> type!\n"); 00054 00055 SGMatrix<complex128_t> mat_A=op->get_matrix_operator(); 00056 Map<MatrixXcd> map_A(mat_A.matrix, mat_A.num_rows, mat_A.num_cols); 00057 Map<VectorXd> map_b(b.vector, b.vlen); 00058 Map<VectorXcd> map_x(x.vector, x.vlen); 00059 00060 switch (m_type) 00061 { 00062 case DS_LLT: 00063 { 00064 LLT<MatrixXcd> llt(map_A); 00065 map_x=llt.solve(map_b.cast<complex128_t>()); 00066 00067 // checking for success 00068 if (llt.info()==NumericalIssue) 00069 SG_WARNING("Matrix is not Hermitian positive definite!\n"); 00070 } 00071 break; 00072 case DS_QR_NOPERM: 00073 map_x=map_A.householderQr().solve(map_b.cast<complex128_t>()); 00074 break; 00075 case DS_QR_COLPERM: 00076 map_x=map_A.colPivHouseholderQr().solve(map_b.cast<complex128_t>()); 00077 break; 00078 case DS_QR_FULLPERM: 00079 map_x=map_A.fullPivHouseholderQr().solve(map_b.cast<complex128_t>()); 00080 break; 00081 case DS_SVD: 00082 map_x=map_A.jacobiSvd(ComputeThinU|ComputeThinV).solve(map_b.cast<complex128_t>()); 00083 break; 00084 } 00085 00086 return x; 00087 } 00088 00089 } 00090 #endif // HAVE_EIGEN3