1 /* 2 XLiFE++ is an extended library of finite elements written in C++ 3 Copyright (C) 2014 Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin 4 5 This program is free software: you can redistribute it and/or modify 6 it under the terms of the GNU General Public License as published by 7 the Free Software Foundation, either version 3 of the License, or 8 (at your option) any later version. 9 This program is distributed in the hope that it will be useful, 10 but WITHOUT ANY WARRANTY; without even the implied warranty of 11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12 GNU General Public License for more details. 13 You should have received a copy of the GNU General Public License 14 along with this program. If not, see <http://www.gnu.org/licenses/>. 15 */ 16 17 /*! 18 \file SymDenseStorage.hpp 19 \authors D. Martin, E. Lunéville, NGUYEN Manh Ha, N. Kielbasiewicz 20 \since 21 jun 2011 21 \date 21 feb 2013 22 23 \brief Definition of the xlifepp::SymDenseStorage class 24 25 The xlifepp::SymDenseStorage class deals with dense storage of a symmetric matrix M (square matrix) 26 storing main diagonal, then lower triangular part by rows 27 diagonal : M11, M22...Mnn, 28 lower triangular part : M21, M31, M32 ..., Mi1, ..., Mii-1, Mn1, ..., Mnn-1 29 it inherits from xlifepp::DenseStorage 30 */ 31 32 #ifndef SYM_DENSE_STORAGE_HPP 33 #define SYM_DENSE_STORAGE_HPP 34 35 #include "config.h" 36 #include "DenseStorage.hpp" 37 38 namespace xlifepp 39 { 40 41 /*! 42 \class SymDenseStorage 43 handles dense storage of "symetric" matrix storing 44 the lower "triangular" part row by row 45 the diagonal is stored first, then the strict lower part 46 */ 47 48 class SymDenseStorage : public DenseStorage { 49 public : 50 51 // Constructors, Destructor 52 SymDenseStorage(string_t id="SymmDenseStorage"); //!< default constructor 53 SymDenseStorage(number_t, string_t id="SymmDenseStorage"); //!< constructor by access type, number of columns and rows (square matrix) ~SymDenseStorage()54 virtual ~SymDenseStorage() {} //!< virtual destructor clone() const55 SymDenseStorage* clone() const //! create a clone (virtual copy constructor, covariant) 56 {return new SymDenseStorage(*this);} 57 SymDenseStorage* toScalar(dimen_t, dimen_t); //!< create a new scalar SymDense storage from current SymDense storage and submatrix sizes 58 59 // size utilities lowerPartSize() const60 number_t lowerPartSize() const //! returns number of matrix entries in lower triangular part of matrix 61 { return (nbRows_ * (nbRows_ - 1)) / 2; } upperPartSize() const62 number_t upperPartSize() const //! returns number of matrix entries in upper triangular part of matrix 63 { return (nbRows_ * (nbRows_ - 1)) / 2; } size() const64 number_t size() const //! returns number of stored entries of matrix 65 { return (nbRows_ * (nbRows_ + 1)) / 2; } 66 67 // access operator 68 number_t pos(number_t i, number_t j, SymType s = _noSymmetry) const; //!< overloaded pos returns adress of entry (i,j) 69 void positions(const std::vector<number_t>&, const std::vector<number_t>&, 70 std::vector<number_t>&, bool errorOn = true, SymType = _noSymmetry ) const; //!< access to submatrix positions 71 std::vector<std::pair<number_t, number_t> > getCol(SymType s, number_t c, number_t r1=1, number_t r2=0) const; //!< get (row indices, adress) of col c in set [r1,r2] 72 std::vector<std::pair<number_t, number_t> > getRow(SymType s, number_t r, number_t c1=1, number_t c2=0) const; //!< get (col indices, adress) of row r in set [c1,c2] 73 74 /* template specializations 75 This rather COMPLICATED implementation of multMatrixVector is made 76 to simulate "template virtual functions" in abstract base class MatrixStorage 77 and thus to avoid "switch...case:" programming and "downcasting" in Matrix*Vector 78 overloaded operator* of class LargeMatrix. 79 Note that the call to template functions above is forced with the "<>" template 80 descriptor in the following specialized virtual functions (this is not recursion) */ 81 82 //---------------------------------------------------------------------------------------------- 83 // set diagonal and Matrix + Matrix 84 //---------------------------------------------------------------------------------------------- 85 //! Set value of Diagonal 86 template<typename T> 87 void setDiagValueSymDense(std::vector<T>& m, const T k); 88 //@{ 89 //! sets values of diagonal matrix (specializations) setDiagValue(std::vector<real_t> & m,const real_t k)90 void setDiagValue(std::vector<real_t>& m, const real_t k) 91 { setDiagValueSymDense<>(m, k);} setDiagValue(std::vector<complex_t> & m,const complex_t k)92 void setDiagValue(std::vector<complex_t>& m, const complex_t k) 93 { setDiagValueSymDense<>(m, k);} 94 //@} 95 96 // template Matrix + Matrix addition and specializations (see below) 97 template<typename M1, typename M2, typename R> 98 void addMatrixMatrix(const std::vector<M1>&, const std::vector<M2>&, std::vector<R>&) const; //!< templated row dense Matrix + Matrix 99 //@{ 100 //! Matrix + Matrix addMatrixMatrix(const std::vector<real_t> & m,const std::vector<real_t> & v,std::vector<real_t> & rv) const101 void addMatrixMatrix(const std::vector<real_t>& m, const std::vector<real_t>& v, std::vector<real_t>& rv) const 102 { addMatrixMatrix<>(m, v, rv);} addMatrixMatrix(const std::vector<real_t> & m,const std::vector<complex_t> & v,std::vector<complex_t> & rv) const103 void addMatrixMatrix(const std::vector<real_t>& m, const std::vector<complex_t>& v, std::vector<complex_t>& rv) const 104 { addMatrixMatrix<>(m, v, rv); } addMatrixMatrix(const std::vector<complex_t> & m,const std::vector<real_t> & v,std::vector<complex_t> & rv) const105 void addMatrixMatrix(const std::vector<complex_t>& m, const std::vector<real_t>& v, std::vector<complex_t>& rv) const 106 { addMatrixMatrix<>(m, v, rv); } addMatrixMatrix(const std::vector<complex_t> & m,const std::vector<complex_t> & v,std::vector<complex_t> & rv) const107 void addMatrixMatrix(const std::vector<complex_t>& m, const std::vector<complex_t>& v, std::vector<complex_t>& rv) const 108 { addMatrixMatrix<>(m, v, rv); } 109 //@} 110 111 //----------------------------------------------------------------------------------------------- 112 // print stuff 113 //----------------------------------------------------------------------------------------------- 114 void printEntries(std::ostream&, const std::vector<real_t>&, 115 number_t vb = 0, const SymType sym = _noSymmetry) const; //!< output row dense matrix of real scalars 116 void printEntries(std::ostream&, const std::vector<complex_t>&, 117 number_t vb = 0, const SymType sym = _noSymmetry) const; //!< output row dense matrix of complex scalars 118 void printEntries(std::ostream&, const std::vector< Vector<real_t> >&, 119 number_t vb = 0, const SymType sym = _noSymmetry) const; //!< output row dense matrix of real vectors 120 void printEntries(std::ostream&, const std::vector< Vector<complex_t> >&, 121 number_t vb = 0, const SymType sym = _noSymmetry) const; //!< output row dense matrix of complex vectors 122 void printEntries(std::ostream&, const std::vector< Matrix<real_t> >&, 123 number_t vb = 0, const SymType sym = _noSymmetry) const; //!< output row dense matrix of real matrices 124 void printEntries(std::ostream&, const std::vector< Matrix<complex_t> >&, 125 number_t vb = 0, const SymType sym = _noSymmetry) const; //!< output row dense matrix of complex matrices 126 127 //---------------------------------------------------------------------------------------------- 128 // template Matrix x Vector & Vector x Matrix multiplications and specializations 129 //---------------------------------------------------------------------------------------------- 130 template<typename M, typename V, typename R> 131 void multMatrixVector(const std::vector<M>&, const std::vector<V>&, std::vector<R>&, SymType) const; //!< templated sym dense Matrix x Vector 132 template<typename M, typename V, typename R> 133 void multVectorMatrix(const std::vector<M>&, const std::vector<V>&, std::vector<R>&, SymType) const; //!< templated Vector x sym dense Matrix 134 template<typename M, typename V, typename R> 135 void multMatrixVector(const std::vector<M>&, V*, R*, SymType) const; //!< templated Vector x sym dense Matrix (pointer form) 136 template<typename M, typename V, typename R> 137 void multVectorMatrix(const std::vector<M>&, V*, R*, SymType) const; //!< templated Vector x sym dense Matrix (pointer form) 138 139 //@{ 140 //! Matrix (Scalar) * Vector (Scalar) multMatrixVector(const std::vector<real_t> & m,const std::vector<real_t> & v,std::vector<real_t> & rv,SymType sym) const141 void multMatrixVector(const std::vector<real_t>& m, const std::vector<real_t>& v, std::vector<real_t>& rv, SymType sym) const 142 { multMatrixVector<>(m, v, rv, sym); } multMatrixVector(const std::vector<real_t> & m,const std::vector<complex_t> & v,std::vector<complex_t> & rv,SymType sym) const143 void multMatrixVector(const std::vector<real_t>& m, const std::vector<complex_t>& v, std::vector<complex_t>& rv, SymType sym) const 144 { multMatrixVector<>(m, v, rv, sym); } multMatrixVector(const std::vector<complex_t> & m,const std::vector<real_t> & v,std::vector<complex_t> & rv,SymType sym) const145 void multMatrixVector(const std::vector<complex_t>& m, const std::vector<real_t>& v, std::vector<complex_t>& rv, SymType sym) const 146 { multMatrixVector<>(m, v, rv, sym); } multMatrixVector(const std::vector<complex_t> & m,const std::vector<complex_t> & v,std::vector<complex_t> & rv,SymType sym) const147 void multMatrixVector(const std::vector<complex_t>& m, const std::vector<complex_t>& v, std::vector<complex_t>& rv, SymType sym) const 148 { multMatrixVector<>(m, v, rv, sym); } 149 //@} 150 151 //@{ 152 //! Matrix (Matrix) * Vector (Vector) multMatrixVector(const std::vector<Matrix<real_t>> & m,const std::vector<Vector<real_t>> & v,std::vector<Vector<real_t>> & rv,SymType sym) const153 void multMatrixVector(const std::vector< Matrix<real_t> >& m, const std::vector< Vector<real_t> >& v, std::vector< Vector<real_t> >& rv, SymType sym) const 154 { multMatrixVector<>(m, v, rv, sym); } multMatrixVector(const std::vector<Matrix<real_t>> & m,const std::vector<Vector<complex_t>> & v,std::vector<Vector<complex_t>> & rv,SymType sym) const155 void multMatrixVector(const std::vector< Matrix<real_t> >& m, const std::vector< Vector<complex_t> >& v, std::vector< Vector<complex_t> >& rv, SymType sym) const 156 { multMatrixVector<>(m, v, rv, sym); } multMatrixVector(const std::vector<Matrix<complex_t>> & m,const std::vector<Vector<real_t>> & v,std::vector<Vector<complex_t>> & rv,SymType sym) const157 void multMatrixVector(const std::vector< Matrix<complex_t> >& m, const std::vector< Vector<real_t> >& v, std::vector< Vector<complex_t> >& rv, SymType sym) const 158 { multMatrixVector<>(m, v, rv, sym); } multMatrixVector(const std::vector<Matrix<complex_t>> & m,const std::vector<Vector<complex_t>> & v,std::vector<Vector<complex_t>> & rv,SymType sym) const159 void multMatrixVector(const std::vector< Matrix<complex_t> >& m, const std::vector< Vector<complex_t> >& v, std::vector< Vector<complex_t> >& rv, SymType sym) const 160 { multMatrixVector<>(m, v, rv, sym); } 161 //@} 162 163 //@{ 164 //! Vector (Scalar) * Matrix (Scalar) multVectorMatrix(const std::vector<real_t> & m,const std::vector<real_t> & v,std::vector<real_t> & rv,SymType sym) const165 void multVectorMatrix(const std::vector<real_t>& m, const std::vector<real_t>& v, std::vector<real_t>& rv, SymType sym) const 166 { multVectorMatrix<>(m, v, rv, sym); } multVectorMatrix(const std::vector<real_t> & m,const std::vector<complex_t> & v,std::vector<complex_t> & rv,SymType sym) const167 void multVectorMatrix(const std::vector<real_t>& m, const std::vector<complex_t>& v, std::vector<complex_t>& rv, SymType sym) const 168 { multVectorMatrix<>(m, v, rv, sym); } multVectorMatrix(const std::vector<complex_t> & m,const std::vector<real_t> & v,std::vector<complex_t> & rv,SymType sym) const169 void multVectorMatrix(const std::vector<complex_t>& m, const std::vector<real_t>& v, std::vector<complex_t>& rv, SymType sym) const 170 { multVectorMatrix<>(m, v, rv, sym); } multVectorMatrix(const std::vector<complex_t> & m,const std::vector<complex_t> & v,std::vector<complex_t> & rv,SymType sym) const171 void multVectorMatrix(const std::vector<complex_t>& m, const std::vector<complex_t>& v, std::vector<complex_t>& rv, SymType sym) const 172 { multVectorMatrix<>(m, v, rv, sym); } 173 //@} 174 175 //@{ 176 //! Matrix (Matrix) * Vector (Vector) multVectorMatrix(const std::vector<Matrix<real_t>> & m,const std::vector<Vector<real_t>> & v,std::vector<Vector<real_t>> & rv,SymType sym) const177 void multVectorMatrix(const std::vector< Matrix<real_t> >& m, const std::vector< Vector<real_t> >& v, std::vector< Vector<real_t> >& rv, SymType sym) const 178 { multVectorMatrix<>(m, v, rv, sym); } multVectorMatrix(const std::vector<Matrix<real_t>> & m,const std::vector<Vector<complex_t>> & v,std::vector<Vector<complex_t>> & rv,SymType sym) const179 void multVectorMatrix(const std::vector< Matrix<real_t> >& m, const std::vector< Vector<complex_t> >& v, std::vector< Vector<complex_t> >& rv, SymType sym) const 180 { multVectorMatrix<>(m, v, rv, sym); } multVectorMatrix(const std::vector<Matrix<complex_t>> & m,const std::vector<Vector<real_t>> & v,std::vector<Vector<complex_t>> & rv,SymType sym) const181 void multVectorMatrix(const std::vector< Matrix<complex_t> >& m, const std::vector< Vector<real_t> >& v, std::vector< Vector<complex_t> >& rv, SymType sym) const 182 { multVectorMatrix<>(m, v, rv, sym); } multVectorMatrix(const std::vector<Matrix<complex_t>> & m,const std::vector<Vector<complex_t>> & v,std::vector<Vector<complex_t>> & rv,SymType sym) const183 void multVectorMatrix(const std::vector< Matrix<complex_t> >& m, const std::vector< Vector<complex_t> >& v, std::vector< Vector<complex_t> >& rv, SymType sym) const 184 { multVectorMatrix<>(m, v, rv, sym); } 185 //@} 186 187 //@{ 188 //! Matrix * Vector (pointer form) multMatrixVector(const std::vector<real_t> & m,real_t * vp,real_t * rp,SymType sym) const189 void multMatrixVector(const std::vector<real_t>& m, real_t* vp, real_t* rp, SymType sym) const 190 { multMatrixVector<>(m, vp, rp, sym); } multMatrixVector(const std::vector<real_t> & m,complex_t * vp,complex_t * rp,SymType sym) const191 void multMatrixVector(const std::vector<real_t>& m, complex_t* vp, complex_t* rp, SymType sym) const 192 { multMatrixVector<>(m, vp, rp, sym); } multMatrixVector(const std::vector<complex_t> & m,real_t * vp,complex_t * rp,SymType sym) const193 void multMatrixVector(const std::vector<complex_t>& m, real_t* vp, complex_t* rp, SymType sym) const 194 { multMatrixVector<>(m, vp, rp, sym); } multMatrixVector(const std::vector<complex_t> & m,complex_t * vp,complex_t * rp,SymType sym) const195 void multMatrixVector(const std::vector<complex_t>& m, complex_t* vp, complex_t* rp, SymType sym) const 196 { multMatrixVector<>(m, vp, rp, sym); } 197 //@} 198 199 //@{ 200 //! Vector * Matrix (pointer form) multVectorMatrix(const std::vector<real_t> & m,real_t * vp,real_t * rp,SymType sym) const201 void multVectorMatrix(const std::vector<real_t>& m, real_t* vp, real_t* rp, SymType sym) const 202 { multVectorMatrix<>(m, vp, rp, sym); } multVectorMatrix(const std::vector<real_t> & m,complex_t * vp,complex_t * rp,SymType sym) const203 void multVectorMatrix(const std::vector<real_t>& m, complex_t* vp, complex_t* rp, SymType sym) const 204 { multVectorMatrix<>(m, vp, rp, sym); } multVectorMatrix(const std::vector<complex_t> & m,real_t * vp,complex_t * rp,SymType sym) const205 void multVectorMatrix(const std::vector<complex_t>& m, real_t* vp, complex_t* rp, SymType sym) const 206 { multVectorMatrix<>(m, vp, rp, sym); } multVectorMatrix(const std::vector<complex_t> & m,complex_t * vp,complex_t * rp,SymType sym) const207 void multVectorMatrix(const std::vector<complex_t>& m, complex_t* vp, complex_t* rp, SymType sym) const 208 { multVectorMatrix<>(m, vp, rp, sym); } 209 //@} 210 211 /*--------------------------------------------------------------------------------------------------------------------- 212 triangular part matrix * vector 213 ------------------------------------------------------------------------------------------------------------------*/ 214 template<typename M, typename V, typename R> 215 void lowerMatrixVector(const std::vector<M>&, const std::vector<V>&, std::vector<R>&) const; 216 //@{ 217 //! lower Matrix (Scalar) * Vector (Scalar) lowerMatrixVector(const std::vector<real_t> & m,const std::vector<real_t> & v,std::vector<real_t> & rv,SymType s) const218 virtual void lowerMatrixVector(const std::vector<real_t>&m, const std::vector<real_t>&v, std::vector<real_t>&rv, SymType s) const 219 {lowerMatrixVector<>(m,v,rv);} lowerMatrixVector(const std::vector<real_t> & m,const std::vector<complex_t> & v,std::vector<complex_t> & rv,SymType s) const220 virtual void lowerMatrixVector(const std::vector<real_t>&m, const std::vector<complex_t>&v, std::vector<complex_t>&rv, SymType s) const 221 {lowerMatrixVector<>(m,v,rv);} lowerMatrixVector(const std::vector<complex_t> & m,const std::vector<real_t> & v,std::vector<complex_t> & rv,SymType s) const222 virtual void lowerMatrixVector(const std::vector<complex_t>&m, const std::vector<real_t>&v, std::vector<complex_t>&rv, SymType s) const 223 {lowerMatrixVector<>(m,v,rv);} lowerMatrixVector(const std::vector<complex_t> & m,const std::vector<complex_t> & v,std::vector<complex_t> & rv,SymType s) const224 virtual void lowerMatrixVector(const std::vector<complex_t>&m, const std::vector<complex_t>&v, std::vector<complex_t>&rv, SymType s) const 225 {lowerMatrixVector<>(m,v,rv);} 226 //@} 227 228 template<typename M, typename V, typename R> 229 void lowerD1MatrixVector(const std::vector<M>&, const std::vector<V>&, std::vector<R>&) const; 230 //@{ 231 //! lower Matrix diag 1 (Scalar) * Vector (Scalar) lowerD1MatrixVector(const std::vector<real_t> & m,const std::vector<real_t> & v,std::vector<real_t> & rv,SymType s) const232 virtual void lowerD1MatrixVector(const std::vector<real_t>& m, const std::vector<real_t>& v, std::vector<real_t>& rv, SymType s) const 233 {lowerD1MatrixVector<>(m,v,rv);} lowerD1MatrixVector(const std::vector<real_t> & m,const std::vector<complex_t> & v,std::vector<complex_t> & rv,SymType s) const234 virtual void lowerD1MatrixVector(const std::vector<real_t>& m, const std::vector<complex_t>& v, std::vector<complex_t>& rv, SymType s) const 235 {lowerD1MatrixVector<>(m,v,rv);} lowerD1MatrixVector(const std::vector<complex_t> & m,const std::vector<real_t> & v,std::vector<complex_t> & rv,SymType s) const236 virtual void lowerD1MatrixVector(const std::vector<complex_t>& m, const std::vector<real_t>& v, std::vector<complex_t>& rv, SymType s) const 237 {lowerD1MatrixVector<>(m,v,rv);} lowerD1MatrixVector(const std::vector<complex_t> & m,const std::vector<complex_t> & v,std::vector<complex_t> & rv,SymType s) const238 virtual void lowerD1MatrixVector(const std::vector<complex_t>& m, const std::vector<complex_t>& v, std::vector<complex_t>& rv, SymType s) const 239 {lowerD1MatrixVector<>(m,v,rv);} 240 //@} 241 242 template<typename M, typename V, typename R> 243 void upperMatrixVector(const std::vector<M>&, const std::vector<V>&, std::vector<R>&, SymType) const; 244 //@{ 245 //! upper Matrix (Scalar) * Vector (Scalar) upperMatrixVector(const std::vector<real_t> & m,const std::vector<real_t> & v,std::vector<real_t> & rv,SymType s) const246 virtual void upperMatrixVector(const std::vector<real_t>&m, const std::vector<real_t>&v, std::vector<real_t>&rv, SymType s) const 247 {upperMatrixVector<>(m,v,rv,s);} upperMatrixVector(const std::vector<real_t> & m,const std::vector<complex_t> & v,std::vector<complex_t> & rv,SymType s) const248 virtual void upperMatrixVector(const std::vector<real_t>&m, const std::vector<complex_t>&v, std::vector<complex_t>&rv, SymType s) const 249 {upperMatrixVector<>(m,v,rv,s);} upperMatrixVector(const std::vector<complex_t> & m,const std::vector<real_t> & v,std::vector<complex_t> & rv,SymType s) const250 virtual void upperMatrixVector(const std::vector<complex_t>&m, const std::vector<real_t>&v, std::vector<complex_t>&rv, SymType s) const 251 {upperMatrixVector<>(m,v,rv,s);} upperMatrixVector(const std::vector<complex_t> & m,const std::vector<complex_t> & v,std::vector<complex_t> & rv,SymType s) const252 virtual void upperMatrixVector(const std::vector<complex_t>&m, const std::vector<complex_t>&v, std::vector<complex_t>&rv, SymType s) const 253 {upperMatrixVector<>(m,v,rv,s);} 254 //@} 255 256 template<typename M, typename V, typename R> 257 void upperD1MatrixVector(const std::vector<M>&, const std::vector<V>&, std::vector<R>&, SymType s) const; 258 //@{ 259 //! upper Matrix diag 1 (Scalar) * Vector (Scalar) upperD1MatrixVector(const std::vector<real_t> & m,const std::vector<real_t> & v,std::vector<real_t> & rv,SymType s) const260 virtual void upperD1MatrixVector(const std::vector<real_t>&m, const std::vector<real_t>&v, std::vector<real_t>&rv, SymType s) const 261 {upperD1MatrixVector<>(m,v,rv,s);} upperD1MatrixVector(const std::vector<real_t> & m,const std::vector<complex_t> & v,std::vector<complex_t> & rv,SymType s) const262 virtual void upperD1MatrixVector(const std::vector<real_t>&m, const std::vector<complex_t>&v, std::vector<complex_t>&rv, SymType s) const 263 {upperD1MatrixVector<>(m,v,rv,s);} upperD1MatrixVector(const std::vector<complex_t> & m,const std::vector<real_t> & v,std::vector<complex_t> & rv,SymType s) const264 virtual void upperD1MatrixVector(const std::vector<complex_t>&m, const std::vector<real_t>&v, std::vector<complex_t>&rv, SymType s) const 265 {upperD1MatrixVector<>(m,v,rv,s);} upperD1MatrixVector(const std::vector<complex_t> & m,const std::vector<complex_t> & v,std::vector<complex_t> & rv,SymType s) const266 virtual void upperD1MatrixVector(const std::vector<complex_t>&m, const std::vector<complex_t>&v, std::vector<complex_t>&rv, SymType s) const 267 {upperD1MatrixVector<>(m,v,rv,s);} 268 //@} 269 270 template<typename M, typename V, typename R> 271 void diagonalMatrixVector(const std::vector<M>&, const std::vector<V>&, std::vector<R>&) const; 272 //@{ 273 //! diag Matrix * Vector (Scalar) diagonalMatrixVector(const std::vector<real_t> & m,const std::vector<real_t> & v,std::vector<real_t> & rv,SymType s) const274 virtual void diagonalMatrixVector(const std::vector<real_t>&m, const std::vector<real_t>&v, std::vector<real_t>&rv, SymType s) const 275 {diagonalMatrixVector<>(m,v,rv);} diagonalMatrixVector(const std::vector<real_t> & m,const std::vector<complex_t> & v,std::vector<complex_t> & rv,SymType s) const276 virtual void diagonalMatrixVector(const std::vector<real_t>&m, const std::vector<complex_t>&v, std::vector<complex_t>&rv, SymType s) const 277 {diagonalMatrixVector<>(m,v,rv);} diagonalMatrixVector(const std::vector<complex_t> & m,const std::vector<real_t> & v,std::vector<complex_t> & rv,SymType s) const278 virtual void diagonalMatrixVector(const std::vector<complex_t>&m, const std::vector<real_t>&v, std::vector<complex_t>&rv, SymType s) const 279 {diagonalMatrixVector<>(m,v,rv);} diagonalMatrixVector(const std::vector<complex_t> & m,const std::vector<complex_t> & v,std::vector<complex_t> & rv,SymType s) const280 virtual void diagonalMatrixVector(const std::vector<complex_t>&m, const std::vector<complex_t>&v, std::vector<complex_t>&rv, SymType s) const 281 {diagonalMatrixVector<>(m,v,rv);} 282 //@} 283 284 //---------------------------------------------------------------------------------------------- 285 // ldlt factorization stuff 286 //---------------------------------------------------------------------------------------------- 287 // template ldlt 288 template<typename T> 289 void ldlt(std::vector<T>& A, std::vector<T>& LD) const; //!< LDLt factorization with no pivoting (omp) 290 //@{ 291 //! LDLt specializations ldlt(std::vector<real_t> & A,std::vector<real_t> & LD,const SymType sym=_noSymmetry) const292 void ldlt(std::vector<real_t>& A, std::vector<real_t>& LD, const SymType sym = _noSymmetry) const 293 {ldlt<>(A, LD);} ldlt(std::vector<complex_t> & A,std::vector<complex_t> & LD,const SymType sym=_noSymmetry) const294 void ldlt(std::vector<complex_t>& A, std::vector<complex_t>& LD, const SymType sym = _noSymmetry) const 295 {ldlt<>(A, LD);} 296 //@} 297 298 //! Diagonal linear system solver : D x = b 299 template<typename M, typename V, typename X> 300 void diagonalSolver(const std::vector<M>& m, std::vector<V>& v, std::vector<X>& x) const; 301 //@{ 302 //! Specializations of diagonal linear solvers D x = v diagonalSolver(const std::vector<real_t> & m,std::vector<real_t> & v,std::vector<real_t> & x) const303 void diagonalSolver(const std::vector<real_t>& m, std::vector<real_t>& v, std::vector<real_t>& x) const 304 { diagonalSolver<>(m, v, x); } diagonalSolver(const std::vector<real_t> & m,std::vector<complex_t> & v,std::vector<complex_t> & x) const305 void diagonalSolver(const std::vector<real_t>& m, std::vector<complex_t>& v, std::vector<complex_t>& x) const 306 { diagonalSolver<>(m, v, x); } diagonalSolver(const std::vector<complex_t> & m,std::vector<real_t> & v,std::vector<complex_t> & x) const307 void diagonalSolver(const std::vector<complex_t>& m, std::vector<real_t>& v, std::vector<complex_t>& x) const 308 { diagonalSolver<>(m, v, x); } diagonalSolver(const std::vector<complex_t> & m,std::vector<complex_t> & v,std::vector<complex_t> & x) const309 void diagonalSolver(const std::vector<complex_t>& m, std::vector<complex_t>& v, std::vector<complex_t>& x) const 310 { diagonalSolver<>(m, v, x); } 311 //@} 312 313 template<typename M, typename V, typename X> 314 void lowerD1Solver(const std::vector<M>& m, std::vector<V>& v, std::vector<X>& x) const; 315 //@{ 316 //! Specializations of lower triangular part with unit diagonal linear solvers (I + L) x = v */ lowerD1Solver(const std::vector<real_t> & m,std::vector<real_t> & v,std::vector<real_t> & x) const317 void lowerD1Solver(const std::vector<real_t>& m, std::vector<real_t>& v, std::vector<real_t>& x) const 318 { lowerD1Solver<>(m, v, x); } lowerD1Solver(const std::vector<real_t> & m,std::vector<complex_t> & v,std::vector<complex_t> & x) const319 void lowerD1Solver(const std::vector<real_t>& m, std::vector<complex_t>& v, std::vector<complex_t>& x) const 320 { lowerD1Solver<>(m, v, x); } lowerD1Solver(const std::vector<complex_t> & m,std::vector<real_t> & v,std::vector<complex_t> & x) const321 void lowerD1Solver(const std::vector<complex_t>& m, std::vector<real_t>& v, std::vector<complex_t>& x) const 322 { lowerD1Solver<>(m, v, x); } lowerD1Solver(const std::vector<complex_t> & m,std::vector<complex_t> & v,std::vector<complex_t> & x) const323 void lowerD1Solver(const std::vector<complex_t>& m, std::vector<complex_t>& v, std::vector<complex_t>& x) const 324 { lowerD1Solver<>(m, v, x); } 325 //@} 326 327 //! Upper triangular with unit diagonal linear system solver : (I+U) x = b 328 template<typename M, typename V, typename X> 329 void upperD1Solver(const std::vector<M>& m, std::vector<V>& v, std::vector<X>& x, const SymType sym = _noSymmetry) const; 330 //@{ 331 //! Specializations of upper triangular part with unit diagonal linear solvers (I + U) x = v upperD1Solver(const std::vector<real_t> & m,std::vector<real_t> & v,std::vector<real_t> & x,const SymType sym=_noSymmetry) const332 void upperD1Solver(const std::vector<real_t>& m, std::vector<real_t>& v, std::vector<real_t>& x, const SymType sym = _noSymmetry) const 333 { upperD1Solver<>(m, v, x, sym); } upperD1Solver(const std::vector<real_t> & m,std::vector<complex_t> & v,std::vector<complex_t> & x,const SymType sym=_noSymmetry) const334 void upperD1Solver(const std::vector<real_t>& m, std::vector<complex_t>& v, std::vector<complex_t>& x, const SymType sym = _noSymmetry) const 335 { upperD1Solver<>(m, v, x, sym); } upperD1Solver(const std::vector<complex_t> & m,std::vector<real_t> & v,std::vector<complex_t> & x,const SymType sym=_noSymmetry) const336 void upperD1Solver(const std::vector<complex_t>& m, std::vector<real_t>& v, std::vector<complex_t>& x, const SymType sym = _noSymmetry) const 337 { upperD1Solver<>(m, v, x, sym); } upperD1Solver(const std::vector<complex_t> & m,std::vector<complex_t> & v,std::vector<complex_t> & x,const SymType sym=_noSymmetry) const338 void upperD1Solver(const std::vector<complex_t>& m, std::vector<complex_t>& v, std::vector<complex_t>& x, const SymType sym = _noSymmetry) const 339 { upperD1Solver<>(m, v, x, sym); } 340 //@} 341 342 //! upper triangular linear system solver : (D+U) x = b 343 template<typename M, typename V, typename X> 344 void upperSolver(const std::vector<M>& m, std::vector<V>& v, std::vector<X>& x, const SymType sym = _noSymmetry) const; 345 //@{ 346 //! specializations of upper triangular part linear solvers (D + U) x = v upperSolver(const std::vector<real_t> & m,std::vector<real_t> & v,std::vector<real_t> & x,const SymType sym=_noSymmetry) const347 void upperSolver(const std::vector<real_t>& m, std::vector<real_t>& v, std::vector<real_t>& x, const SymType sym = _noSymmetry) const 348 { upperSolver<>(m, v, x, sym); } upperSolver(const std::vector<real_t> & m,std::vector<complex_t> & v,std::vector<complex_t> & x,const SymType sym=_noSymmetry) const349 void upperSolver(const std::vector<real_t>& m, std::vector<complex_t>& v, std::vector<complex_t>& x, const SymType sym = _noSymmetry) const 350 { upperSolver<>(m, v, x, sym); } upperSolver(const std::vector<complex_t> & m,std::vector<real_t> & v,std::vector<complex_t> & x,const SymType sym=_noSymmetry) const351 void upperSolver(const std::vector<complex_t>& m, std::vector<real_t>& v, std::vector<complex_t>& x, const SymType sym = _noSymmetry) const 352 { upperSolver<>(m, v, x, sym); } upperSolver(const std::vector<complex_t> & m,std::vector<complex_t> & v,std::vector<complex_t> & x,const SymType sym=_noSymmetry) const353 void upperSolver(const std::vector<complex_t>& m, std::vector<complex_t>& v, std::vector<complex_t>& x, const SymType sym = _noSymmetry) const 354 { upperSolver<>(m, v, x, sym); } 355 //@} 356 357 //---------------------------------------------------------------------------------------------- 358 // UMFpack stuff 359 //---------------------------------------------------------------------------------------------- 360 //! conversion to umfpack format 361 template<typename M, typename OrdinalType> 362 void toUmfPack(const std::vector<M>& values, std::vector<OrdinalType>& colPointer, std::vector<OrdinalType>& rowIndex, std::vector<M>& mat, const SymType sym = _noSymmetry) const; 363 //@{ 364 //! specializations od umfpack conversion toUmfPack(const std::vector<real_t> & values,std::vector<int_t> & colPointer,std::vector<int_t> & rowIndex,std::vector<real_t> & mat,const SymType sym) const365 void toUmfPack(const std::vector<real_t>& values, std::vector<int_t>& colPointer, std::vector<int_t>& rowIndex, std::vector<real_t>& mat, const SymType sym) const 366 { toUmfPack<>(values,colPointer,rowIndex,mat, sym);} toUmfPack(const std::vector<complex_t> & values,std::vector<int_t> & colPointer,std::vector<int_t> & rowIndex,std::vector<complex_t> & mat,const SymType sym) const367 void toUmfPack(const std::vector<complex_t>& values, std::vector<int_t>& colPointer, std::vector<int_t>& rowIndex, std::vector<complex_t>& mat, const SymType sym) const 368 { toUmfPack<>(values,colPointer,rowIndex,mat, sym);} 369 //@} 370 }; 371 372 } // end of namespace xlifepp 373 374 #endif // SYM_DENSE_STORAGE_HPP 375 376