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 DualDenseStorage.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::DualDenseStorage class 24 25 The xlifepp::DualDenseStorage class deals with dense storage of a matrix M 26 storing main diagonal, then lower triangular part by rows and upper triangular part by columns 27 diagonal : M11, M22...Mpp, with p=min(m,n) 28 lower triangular part : M21, M31, M32 ..., Mi1, ..., Mii-1, Mm1, ..., Mmn-1 29 upper triangular part : M12, M13, M23 ..., M1j, ..., Mj-1j, M1n, ..., Mm-1n 30 it inherits from xlifepp::DenseStorage 31 */ 32 33 #ifndef DUAL_DENSE_STORAGE_HPP 34 #define DUAL_DENSE_STORAGE_HPP 35 36 #include "config.h" 37 #include "DenseStorage.hpp" 38 39 namespace xlifepp 40 { 41 42 /*! 43 \class DualDenseStorage 44 handles dense storage of matrix stored 45 row by row for the lower "triangular" part and column by column for the upper "triangular" 46 */ 47 48 class DualDenseStorage : public DenseStorage 49 { 50 public: 51 // Constructors, Destructor 52 DualDenseStorage(string_t id="DualDenseStorage"); //!< default constructor 53 DualDenseStorage(number_t, string_t id="DualDenseStorage"); //!< constructor by access type, number of columns and rows (square matrix) 54 DualDenseStorage(number_t, number_t, string_t id="DualDenseStorage"); //!< constructor by access type, number of columns and rows ~DualDenseStorage()55 virtual ~DualDenseStorage() {} //!< virtual destructor clone() const56 DualDenseStorage* clone() const //! create a clone (virtual copy constructor, covariant) 57 {return new DualDenseStorage(*this);} 58 DualDenseStorage* toScalar(dimen_t, dimen_t); //!< create a new scalar DualDense storage from current DualDense storage and submatrix sizes 59 60 61 // size utilities 62 number_t lowerPartSize() const; //!< returns number of matrix entries in lower triangular part of matrix 63 number_t upperPartSize() const; //!< returns number of matrix entries in upper triangular part of matrix 64 65 // access operator 66 number_t pos(number_t i, number_t j, SymType s = _noSymmetry) const; //!< overloaded pos returns adress of entry (i,j) 67 void positions(const std::vector<number_t>&, const std::vector<number_t>&, 68 std::vector<number_t>&, bool errorOn = true, SymType = _noSymmetry) const; //!< access to submatrix positions 69 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] 70 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] 71 72 73 /* template specializations 74 This rather COMPLICATED implementation of multMatrixVector is made 75 to simulate "template virtual functions" in abstract base class MatrixStorage 76 and thus to avoid "switch...case:" programming and "downcasting" in Matrix*Vector 77 overloaded operator* of class LargeMatrix. 78 Note that the call to template functions above is forced with the "<>" template 79 descriptor in the following specialized virtual functions (this is not recursion) */ 80 81 /*------------------------------------------------------------------------------------------------------------------ 82 input and output stuff 83 ------------------------------------------------------------------------------------------------------------------*/ 84 void printEntries(std::ostream&, const std::vector<real_t>&, 85 number_t vb = 0, const SymType sym = _noSymmetry) const; //!< output dual dense matrix of real scalars 86 void printEntries(std::ostream&, const std::vector<complex_t>&, 87 number_t vb = 0, const SymType sym = _noSymmetry) const; //!< output dual dense matrix of complex scalars 88 void printEntries(std::ostream&, const std::vector< Vector<real_t> >&, 89 number_t vb = 0, const SymType sym = _noSymmetry) const; //!< output dual dense matrix of real vectors 90 void printEntries(std::ostream&, const std::vector< Vector<complex_t> >&, 91 number_t vb = 0, const SymType sym = _noSymmetry) const; //!< output dual dense matrix of complex vectors 92 void printEntries(std::ostream&, const std::vector< Matrix<real_t> >&, 93 number_t vb = 0, const SymType sym = _noSymmetry) const; //!< output dual dense matrix of real matrices 94 void printEntries(std::ostream&, const std::vector< Matrix<complex_t> >&, 95 number_t vb = 0, const SymType sym = _noSymmetry) const; //!< output dual dense matrix of complex matrices 96 97 /*------------------------------------------------------------------------------------------------------------------ 98 basic algebra stuff 99 ------------------------------------------------------------------------------------------------------------------*/ 100 //! Set value of Diagonal 101 template<typename T> 102 void setDiagValueDualDense(std::vector<T>& m, const T k); 103 104 //@{ 105 //! set value of diagonal setDiagValue(std::vector<real_t> & m,const real_t k)106 void setDiagValue(std::vector<real_t>& m, const real_t k) 107 { setDiagValueDualDense<>(m, k);} setDiagValue(std::vector<complex_t> & m,const complex_t k)108 void setDiagValue(std::vector<complex_t>& m, const complex_t k) 109 { setDiagValueDualDense<>(m, k);} 110 //@} 111 112 // template Matrix + Matrix addition and specializations (see below) 113 template<typename M1, typename M2, typename R> 114 void addMatrixMatrix(const std::vector<M1>&, const std::vector<M2>&, std::vector<R>&) const; //!< templated row dense Matrix + Matrix 115 116 //@{ 117 //! Matrix+Matrix addMatrixMatrix(const std::vector<real_t> & m,const std::vector<real_t> & v,std::vector<real_t> & rv) const118 void addMatrixMatrix(const std::vector<real_t>& m, const std::vector<real_t>& v, std::vector<real_t>& rv) const 119 { addMatrixMatrix<>(m, v, rv);} addMatrixMatrix(const std::vector<real_t> & m,const std::vector<complex_t> & v,std::vector<complex_t> & rv) const120 void addMatrixMatrix(const std::vector<real_t>& m, const std::vector<complex_t>& v, std::vector<complex_t>& rv) const 121 { addMatrixMatrix<>(m, v, rv); } addMatrixMatrix(const std::vector<complex_t> & m,const std::vector<real_t> & v,std::vector<complex_t> & rv) const122 void addMatrixMatrix(const std::vector<complex_t>& m, const std::vector<real_t>& v, std::vector<complex_t>& rv) const 123 { addMatrixMatrix<>(m, v, rv); } addMatrixMatrix(const std::vector<complex_t> & m,const std::vector<complex_t> & v,std::vector<complex_t> & rv) const124 void addMatrixMatrix(const std::vector<complex_t>& m, const std::vector<complex_t>& v, std::vector<complex_t>& rv) const 125 { addMatrixMatrix<>(m, v, rv); } 126 //@} 127 128 129 /*------------------------------------------------------------------------------------------------------------------ 130 template Matrix x Vector & Vector x Matrix multiplications and specializations 131 ------------------------------------------------------------------------------------------------------------------*/ 132 template<typename M, typename V, typename R> 133 void multMatrixVector(const std::vector<M>&, const std::vector<V>&, std::vector<R>&) const; //!< templated row dense Matrix x Vector 134 template<typename M, typename V, typename R> 135 void multVectorMatrix(const std::vector<M>&, const std::vector<V>&, std::vector<R>&) const; //!< templated Vector x row dense Matrix 136 template<typename M, typename V, typename R> 137 void multMatrixVector(const std::vector<M>& m, V* vp, R* rp) const; //!< templated dual dense Matrix x Vector (pointer) 138 template<typename M, typename V, typename R> 139 void multVectorMatrix(const std::vector<M>& m, V* vp, R* rp) const; //!< templated Vector x dual dense Matrix (pointer) 140 141 //@{ 142 //! Matrix (Scalar) * Vector (Scalar) multMatrixVector(const std::vector<real_t> & m,const std::vector<real_t> & v,std::vector<real_t> & rv,SymType) const143 void multMatrixVector(const std::vector<real_t>& m, const std::vector<real_t>& v, std::vector<real_t>& rv, SymType) const 144 { multMatrixVector<>(m, v, rv); } multMatrixVector(const std::vector<real_t> & m,const std::vector<complex_t> & v,std::vector<complex_t> & rv,SymType) const145 void multMatrixVector(const std::vector<real_t>& m, const std::vector<complex_t>& v, std::vector<complex_t>& rv, SymType) const 146 { multMatrixVector<>(m, v, rv); } multMatrixVector(const std::vector<complex_t> & m,const std::vector<real_t> & v,std::vector<complex_t> & rv,SymType) const147 void multMatrixVector(const std::vector<complex_t>& m, const std::vector<real_t>& v, std::vector<complex_t>& rv, SymType) const 148 { multMatrixVector<>(m, v, rv); } multMatrixVector(const std::vector<complex_t> & m,const std::vector<complex_t> & v,std::vector<complex_t> & rv,SymType) const149 void multMatrixVector(const std::vector<complex_t>& m, const std::vector<complex_t>& v, std::vector<complex_t>& rv, SymType) const 150 { multMatrixVector<>(m, v, rv); } 151 //@} 152 153 //@{ 154 //! 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) const155 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 156 { multMatrixVector<>(m, v, rv); } multMatrixVector(const std::vector<Matrix<real_t>> & m,const std::vector<Vector<complex_t>> & v,std::vector<Vector<complex_t>> & rv,SymType sym) const157 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 158 { multMatrixVector<>(m, v, rv); } multMatrixVector(const std::vector<Matrix<complex_t>> & m,const std::vector<Vector<real_t>> & v,std::vector<Vector<complex_t>> & rv,SymType sym) const159 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 160 { multMatrixVector<>(m, v, rv); } multMatrixVector(const std::vector<Matrix<complex_t>> & m,const std::vector<Vector<complex_t>> & v,std::vector<Vector<complex_t>> & rv,SymType sym) const161 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 162 { multMatrixVector<>(m, v, rv); } 163 //@} 164 165 //@{ 166 //! Vector (Scalar) * Matrix (Scalar) multVectorMatrix(const std::vector<real_t> & m,const std::vector<real_t> & v,std::vector<real_t> & rv,SymType) const167 void multVectorMatrix(const std::vector<real_t>& m, const std::vector<real_t>& v, std::vector<real_t>& rv, SymType) const 168 { multVectorMatrix<>(m, v, rv); } multVectorMatrix(const std::vector<real_t> & m,const std::vector<complex_t> & v,std::vector<complex_t> & rv,SymType) const169 void multVectorMatrix(const std::vector<real_t>& m, const std::vector<complex_t>& v, std::vector<complex_t>& rv, SymType) const 170 { multVectorMatrix<>(m, v, rv); } multVectorMatrix(const std::vector<complex_t> & m,const std::vector<real_t> & v,std::vector<complex_t> & rv,SymType) const171 void multVectorMatrix(const std::vector<complex_t>& m, const std::vector<real_t>& v, std::vector<complex_t>& rv, SymType) const 172 { multVectorMatrix<>(m, v, rv); } multVectorMatrix(const std::vector<complex_t> & m,const std::vector<complex_t> & v,std::vector<complex_t> & rv,SymType) const173 void multVectorMatrix(const std::vector<complex_t>& m, const std::vector<complex_t>& v, std::vector<complex_t>& rv, SymType) const 174 { multVectorMatrix<>(m, v, rv); } 175 //@} 176 177 //@{ 178 //! 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) const179 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 180 { multVectorMatrix<>(m, v, rv); } multVectorMatrix(const std::vector<Matrix<real_t>> & m,const std::vector<Vector<complex_t>> & v,std::vector<Vector<complex_t>> & rv,SymType sym) const181 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 182 { multVectorMatrix<>(m, v, rv); } multVectorMatrix(const std::vector<Matrix<complex_t>> & m,const std::vector<Vector<real_t>> & v,std::vector<Vector<complex_t>> & rv,SymType sym) const183 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 184 { multVectorMatrix<>(m, v, rv); } multVectorMatrix(const std::vector<Matrix<complex_t>> & m,const std::vector<Vector<complex_t>> & v,std::vector<Vector<complex_t>> & rv,SymType sym) const185 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 186 { multVectorMatrix<>(m, v, rv); } 187 //@} 188 189 //@{ 190 //! Matrix * Vector (pointer form) multMatrixVector(const std::vector<real_t> & m,real_t * vp,real_t * rp,SymType sym) const191 void multMatrixVector(const std::vector<real_t>& m, real_t* vp, real_t* rp, SymType sym) const 192 { multMatrixVector<>(m, vp, rp); } multMatrixVector(const std::vector<real_t> & m,complex_t * vp,complex_t * rp,SymType sym) const193 void multMatrixVector(const std::vector<real_t>& m, complex_t* vp, complex_t* rp, SymType sym) const 194 { multMatrixVector<>(m, vp, rp); } multMatrixVector(const std::vector<complex_t> & m,real_t * vp,complex_t * rp,SymType sym) const195 void multMatrixVector(const std::vector<complex_t>& m, real_t* vp, complex_t* rp, SymType sym) const 196 { multMatrixVector<>(m, vp, rp); } multMatrixVector(const std::vector<complex_t> & m,complex_t * vp,complex_t * rp,SymType sym) const197 void multMatrixVector(const std::vector<complex_t>& m, complex_t* vp, complex_t* rp, SymType sym) const 198 { multMatrixVector<>(m, vp, rp); } 199 //@} 200 201 //@{ 202 //! Vector * Matrix (pointer form) multVectorMatrix(const std::vector<real_t> & m,real_t * vp,real_t * rp,SymType sym) const203 void multVectorMatrix(const std::vector<real_t>& m, real_t* vp, real_t* rp, SymType sym) const 204 { multVectorMatrix<>(m, vp, rp); } multVectorMatrix(const std::vector<real_t> & m,complex_t * vp,complex_t * rp,SymType sym) const205 void multVectorMatrix(const std::vector<real_t>& m, complex_t* vp, complex_t* rp, SymType sym) const 206 { multVectorMatrix<>(m, vp, rp); } multVectorMatrix(const std::vector<complex_t> & m,real_t * vp,complex_t * rp,SymType sym) const207 void multVectorMatrix(const std::vector<complex_t>& m, real_t* vp, complex_t* rp, SymType sym) const 208 { multVectorMatrix<>(m, vp, rp); } multVectorMatrix(const std::vector<complex_t> & m,complex_t * vp,complex_t * rp,SymType sym) const209 void multVectorMatrix(const std::vector<complex_t>& m, complex_t* vp, complex_t* rp, SymType sym) const 210 { multVectorMatrix<>(m, vp, rp); } 211 //@} 212 213 /*------------------------------------------------------------------------------------------------------------------ 214 triangular part matrix * vector 215 ------------------------------------------------------------------------------------------------------------------*/ 216 template<typename M, typename V, typename R> 217 void lowerMatrixVector(const std::vector<M>&, const std::vector<V>&, std::vector<R>&) const; 218 //@{ 219 //! 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) const220 virtual void lowerMatrixVector(const std::vector<real_t>&m, const std::vector<real_t>&v, std::vector<real_t>&rv, SymType s) const 221 {lowerMatrixVector<>(m,v,rv);} lowerMatrixVector(const std::vector<real_t> & m,const std::vector<complex_t> & v,std::vector<complex_t> & rv,SymType s) const222 virtual void lowerMatrixVector(const std::vector<real_t>&m, const std::vector<complex_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<real_t> & v,std::vector<complex_t> & rv,SymType s) const224 virtual void lowerMatrixVector(const std::vector<complex_t>&m, const std::vector<real_t>&v, std::vector<complex_t>&rv, SymType s) const 225 {lowerMatrixVector<>(m,v,rv);} lowerMatrixVector(const std::vector<complex_t> & m,const std::vector<complex_t> & v,std::vector<complex_t> & rv,SymType s) const226 virtual void lowerMatrixVector(const std::vector<complex_t>&m, const std::vector<complex_t>&v, std::vector<complex_t>&rv, SymType s) const 227 {lowerMatrixVector<>(m,v,rv);} 228 //@} 229 template<typename M, typename V, typename R> 230 void lowerD1MatrixVector(const std::vector<M>&, const std::vector<V>&, std::vector<R>&) const; 231 //@{ 232 //! 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) const233 virtual void lowerD1MatrixVector(const std::vector<real_t>& m, const std::vector<real_t>& v, std::vector<real_t>& rv, SymType s) const 234 {lowerD1MatrixVector<>(m,v,rv);} lowerD1MatrixVector(const std::vector<real_t> & m,const std::vector<complex_t> & v,std::vector<complex_t> & rv,SymType s) const235 virtual void lowerD1MatrixVector(const std::vector<real_t>& m, const std::vector<complex_t>& v, std::vector<complex_t>& rv, SymType s) const 236 {lowerD1MatrixVector<>(m,v,rv);} lowerD1MatrixVector(const std::vector<complex_t> & m,const std::vector<real_t> & v,std::vector<complex_t> & rv,SymType s) const237 virtual void lowerD1MatrixVector(const std::vector<complex_t>& m, const std::vector<real_t>& v, std::vector<complex_t>& rv, SymType s) const 238 {lowerD1MatrixVector<>(m,v,rv);} lowerD1MatrixVector(const std::vector<complex_t> & m,const std::vector<complex_t> & v,std::vector<complex_t> & rv,SymType s) const239 virtual void lowerD1MatrixVector(const std::vector<complex_t>& m, const std::vector<complex_t>& v, std::vector<complex_t>& rv, SymType s) const 240 {lowerD1MatrixVector<>(m,v,rv);} 241 //@} 242 template<typename M, typename V, typename R> 243 void upperMatrixVector(const std::vector<M>&, const std::vector<V>&, std::vector<R>&) 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);} 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);} 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);} 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);} 254 //@} 255 template<typename M, typename V, typename R> 256 void upperD1MatrixVector(const std::vector<M>&, const std::vector<V>&, std::vector<R>&) const; 257 //@{ 258 //! 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) const259 virtual void upperD1MatrixVector(const std::vector<real_t>&m, const std::vector<real_t>&v, std::vector<real_t>&rv, SymType s) const 260 {upperD1MatrixVector<>(m,v,rv);} upperD1MatrixVector(const std::vector<real_t> & m,const std::vector<complex_t> & v,std::vector<complex_t> & rv,SymType s) const261 virtual void upperD1MatrixVector(const std::vector<real_t>&m, const std::vector<complex_t>&v, std::vector<complex_t>&rv, SymType s) const 262 {upperD1MatrixVector<>(m,v,rv);} upperD1MatrixVector(const std::vector<complex_t> & m,const std::vector<real_t> & v,std::vector<complex_t> & rv,SymType s) const263 virtual void upperD1MatrixVector(const std::vector<complex_t>&m, const std::vector<real_t>&v, std::vector<complex_t>&rv, SymType s) const 264 {upperD1MatrixVector<>(m,v,rv);} upperD1MatrixVector(const std::vector<complex_t> & m,const std::vector<complex_t> & v,std::vector<complex_t> & rv,SymType s) const265 virtual void upperD1MatrixVector(const std::vector<complex_t>&m, const std::vector<complex_t>&v, std::vector<complex_t>&rv, SymType s) const 266 {upperD1MatrixVector<>(m,v,rv);} 267 //@} 268 template<typename M, typename V, typename R> 269 void diagonalMatrixVector(const std::vector<M>&, const std::vector<V>&, std::vector<R>&) const; 270 //@{ 271 //! diag Matrix * Vector (Scalar) diagonalMatrixVector(const std::vector<real_t> & m,const std::vector<real_t> & v,std::vector<real_t> & rv,SymType s) const272 virtual void diagonalMatrixVector(const std::vector<real_t>&m, const std::vector<real_t>&v, std::vector<real_t>&rv, SymType s) const 273 {diagonalMatrixVector<>(m,v,rv);} diagonalMatrixVector(const std::vector<real_t> & m,const std::vector<complex_t> & v,std::vector<complex_t> & rv,SymType s) const274 virtual void diagonalMatrixVector(const std::vector<real_t>&m, const std::vector<complex_t>&v, std::vector<complex_t>&rv, SymType s) const 275 {diagonalMatrixVector<>(m,v,rv);} diagonalMatrixVector(const std::vector<complex_t> & m,const std::vector<real_t> & v,std::vector<complex_t> & rv,SymType s) const276 virtual void diagonalMatrixVector(const std::vector<complex_t>&m, const std::vector<real_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<complex_t> & v,std::vector<complex_t> & rv,SymType s) const278 virtual void diagonalMatrixVector(const std::vector<complex_t>&m, const std::vector<complex_t>&v, std::vector<complex_t>&rv, SymType s) const 279 {diagonalMatrixVector<>(m,v,rv);} 280 //@} 281 282 /* ------------------------------------------------------------------------------------------------------------------------- 283 SOR and SSOR stuff 284 ----------------------------------------------------------------------------------------------------------------------------*/ 285 /*! special template partial Matrix x Vector multiplications [w*D] v used in SSOR algorithm 286 for D + L + U matrix splitting (D diagonal, L & U strict lower & upper trangular parts) */ 287 template<typename M, typename V, typename R> sorDiagonalMatrixVector(const std::vector<M> & m,const std::vector<V> & v,std::vector<R> & r,const real_t w) const288 void sorDiagonalMatrixVector(const std::vector<M>& m, const std::vector<V>& v, std::vector<R>& r, const real_t w) const 289 { 290 typename std::vector<M>::const_iterator it_m = m.begin() + 1; 291 typename std::vector<V>::const_iterator it_vb = v.begin(); 292 typename std::vector<R>::iterator it_rb = r.begin(), it_re=r.end();; 293 bzSorDiagonalMatrixVector(it_m, it_vb, it_rb, it_re, w); 294 } 295 //@{ 296 //! specializations of partial matrix std::vector multiplications sorDiagonalMatrixVector(const std::vector<real_t> & m,const std::vector<real_t> & v,std::vector<real_t> & rv,const real_t w) const297 void sorDiagonalMatrixVector(const std::vector<real_t>& m, const std::vector<real_t>& v, std::vector<real_t>& rv, const real_t w) const 298 { sorDiagonalMatrixVector<>(m, v, rv, w); } sorDiagonalMatrixVector(const std::vector<real_t> & m,const std::vector<complex_t> & v,std::vector<complex_t> & rv,const real_t w) const299 void sorDiagonalMatrixVector(const std::vector<real_t>& m, const std::vector<complex_t>& v, std::vector<complex_t>& rv, const real_t w) const 300 { sorDiagonalMatrixVector<>(m, v, rv, w); } sorDiagonalMatrixVector(const std::vector<complex_t> & m,const std::vector<real_t> & v,std::vector<complex_t> & rv,const real_t w) const301 void sorDiagonalMatrixVector(const std::vector<complex_t>& m, const std::vector<real_t>& v, std::vector<complex_t>& rv, const real_t w) const 302 { sorDiagonalMatrixVector<>(m, v, rv, w); } sorDiagonalMatrixVector(const std::vector<complex_t> & m,const std::vector<complex_t> & v,std::vector<complex_t> & rv,const real_t w) const303 void sorDiagonalMatrixVector(const std::vector<complex_t>& m, const std::vector<complex_t>& v, std::vector<complex_t>& rv, const real_t w) const 304 { sorDiagonalMatrixVector<>(m, v, rv, w); } 305 //@} 306 307 /*! special template partial Matrix x Vector multiplications [w*D+L] v used in SSOR algorithm 308 for D + L + U matrix splitting (D diagonal, L & U strict lower & upper trangular parts) */ 309 template<typename M, typename V, typename R> sorLowerMatrixVector(const std::vector<M> & m,const std::vector<V> & v,std::vector<R> & r,const real_t w,const SymType sym) const310 void sorLowerMatrixVector(const std::vector<M>& m, const std::vector<V>& v, 311 std::vector<R>& r, const real_t w, const SymType sym) const 312 { 313 typename std::vector<M>::const_iterator it_m = m.begin() + 1, it_l=it_m+diagonalSize(); 314 typename std::vector<V>::const_iterator it_vb = v.begin(), it_ve=v.end(); 315 typename std::vector<R>::iterator it_rb = r.begin(), it_re=r.end(); 316 bzSorDiagonalMatrixVector(it_m, it_vb, it_rb, it_re, w); 317 DenseStorage::lowerMatrixVector(it_l, it_vb, it_ve, it_rb, it_re, sym); 318 } 319 //@{ 320 //! specializations of partial matrix vector multiplications sorLowerMatrixVector(const std::vector<real_t> & m,const std::vector<real_t> & v,std::vector<real_t> & rv,const real_t w,const SymType sym) const321 void sorLowerMatrixVector(const std::vector<real_t>& m, const std::vector<real_t>& v, std::vector<real_t>& rv, const real_t w, const SymType sym) const 322 { sorLowerMatrixVector<>(m, v, rv, w, sym); } sorLowerMatrixVector(const std::vector<real_t> & m,const std::vector<complex_t> & v,std::vector<complex_t> & rv,const real_t w,const SymType sym) const323 void sorLowerMatrixVector(const std::vector<real_t>& m, const std::vector<complex_t>& v, std::vector<complex_t>& rv, const real_t w, const SymType sym) const 324 { sorLowerMatrixVector<>(m, v, rv, w, sym); } sorLowerMatrixVector(const std::vector<complex_t> & m,const std::vector<real_t> & v,std::vector<complex_t> & rv,const real_t w,const SymType sym) const325 void sorLowerMatrixVector(const std::vector<complex_t>& m, const std::vector<real_t>& v, std::vector<complex_t>& rv, const real_t w, const SymType sym) const 326 { sorLowerMatrixVector<>(m, v, rv, w, sym); } sorLowerMatrixVector(const std::vector<complex_t> & m,const std::vector<complex_t> & v,std::vector<complex_t> & rv,const real_t w,const SymType sym) const327 void sorLowerMatrixVector(const std::vector<complex_t>& m, const std::vector<complex_t>& v, std::vector<complex_t>& rv, const real_t w, const SymType sym) const 328 { sorLowerMatrixVector<>(m, v, rv, w, sym); } 329 //@} 330 331 /*! special template partial Matrix x Vector multiplications [w*D+U] v used in SSOR algorithm 332 for D + L + U matrix splitting (D diagonal, L & U strict lower & upper trangular parts) */ 333 template<typename M, typename V, typename R> sorUpperMatrixVector(const std::vector<M> & m,const std::vector<V> & v,std::vector<R> & r,const real_t w,const SymType sym) const334 void sorUpperMatrixVector(const std::vector<M>& m, const std::vector<V>& v, 335 std::vector<R>& r, const real_t w, const SymType sym) const 336 { 337 typename std::vector<M>::const_iterator it_m = m.begin() + 1, it_u=it_m+diagonalSize()+lowerPartSize(); 338 typename std::vector<V>::const_iterator it_vb = v.begin(), it_ve=v.end(); 339 typename std::vector<R>::iterator it_rb = r.begin(), it_re=r.end(); 340 bzSorDiagonalMatrixVector(it_m, it_vb, it_rb, it_re, w); 341 DenseStorage::upperMatrixVector(it_u, it_vb, it_ve, it_rb, it_re, sym); 342 } 343 //@{ 344 //! specializations of partial matrix vector multiplications sorUpperMatrixVector(const std::vector<real_t> & m,const std::vector<real_t> & v,std::vector<real_t> & rv,const real_t w,const SymType sym) const345 void sorUpperMatrixVector(const std::vector<real_t>& m, const std::vector<real_t>& v, std::vector<real_t>& rv, const real_t w, const SymType sym) const 346 { sorUpperMatrixVector<>(m, v, rv, w, sym); } sorUpperMatrixVector(const std::vector<real_t> & m,const std::vector<complex_t> & v,std::vector<complex_t> & rv,const real_t w,const SymType sym) const347 void sorUpperMatrixVector(const std::vector<real_t>& m, const std::vector<complex_t>& v, std::vector<complex_t>& rv, const real_t w, const SymType sym) const 348 { sorUpperMatrixVector<>(m, v, rv, w, sym); } sorUpperMatrixVector(const std::vector<complex_t> & m,const std::vector<real_t> & v,std::vector<complex_t> & rv,const real_t w,const SymType sym) const349 void sorUpperMatrixVector(const std::vector<complex_t>& m, const std::vector<real_t>& v, std::vector<complex_t>& rv, const real_t w, const SymType sym) const 350 { sorUpperMatrixVector<>(m, v, rv, w, sym); } sorUpperMatrixVector(const std::vector<complex_t> & m,const std::vector<complex_t> & v,std::vector<complex_t> & rv,const real_t w,const SymType sym) const351 void sorUpperMatrixVector(const std::vector<complex_t>& m, const std::vector<complex_t>& v, std::vector<complex_t>& rv, const real_t w, const SymType sym) const 352 { sorUpperMatrixVector<>(m, v, rv, w, sym); } 353 //@} 354 355 /*! special template diagonal and triangular solvers D/w x = b used in SSOR algorithm 356 for D + L + U matrix splitting (D diagonal, L & U strict lower & upper triangular parts) */ 357 template<typename M, typename V, typename R> sorDiagonalSolver(const std::vector<M> & m,const std::vector<R> & b,std::vector<V> & x,const real_t w) const358 void sorDiagonalSolver(const std::vector<M>& m, const std::vector<R>& b, 359 std::vector<V>& x, const real_t w) const 360 { 361 typename std::vector<M>::const_iterator it_m = m.begin() + 1; 362 typename std::vector<R>::const_iterator it_b = b.begin(); 363 typename std::vector<V>::iterator it_xb = x.begin(), it_xe = x.end(); 364 bzSorDiagonalSolver(it_m, it_b, it_xb, it_xe, w); 365 } 366 //@{ 367 //! specializations of diagonal solvers sorDiagonalSolver(const std::vector<real_t> & m,const std::vector<real_t> & b,std::vector<real_t> & x,const real_t w) const368 void sorDiagonalSolver(const std::vector<real_t>& m, const std::vector<real_t>& b, std::vector<real_t>& x, const real_t w) const 369 { sorDiagonalSolver<>(m, b, x, w); } sorDiagonalSolver(const std::vector<real_t> & m,const std::vector<complex_t> & b,std::vector<complex_t> & x,const real_t w) const370 void sorDiagonalSolver(const std::vector<real_t>& m, const std::vector<complex_t>& b, std::vector<complex_t>& x, const real_t w) const 371 { sorDiagonalSolver<>(m, b, x, w); } sorDiagonalSolver(const std::vector<complex_t> & m,const std::vector<real_t> & b,std::vector<complex_t> & x,const real_t w) const372 void sorDiagonalSolver(const std::vector<complex_t>& m, const std::vector<real_t>& b, std::vector<complex_t>& x, const real_t w) const 373 { sorDiagonalSolver<>(m, b, x, w); } sorDiagonalSolver(const std::vector<complex_t> & m,const std::vector<complex_t> & b,std::vector<complex_t> & x,const real_t w) const374 void sorDiagonalSolver(const std::vector<complex_t>& m, const std::vector<complex_t>& b, std::vector<complex_t>& x, const real_t w) const 375 { sorDiagonalSolver<>(m, b, x, w); } 376 //@} 377 378 /*! special template diagonal and triangular solvers (D/w+L) x = b used in SSOR algorithm 379 for D + L + U matrix splitting (D diagonal, L & U strict lower & upper triangular parts) */ 380 template<typename M, typename V, typename R> sorLowerSolver(const std::vector<M> & m,const std::vector<R> & b,std::vector<V> & x,const real_t w) const381 void sorLowerSolver(const std::vector<M>& m, const std::vector<R>& b, 382 std::vector<V>& x, const real_t w) const 383 { 384 typename std::vector<M>::const_iterator it_d = m.begin() + 1, it_l = it_d + diagonalSize(); 385 typename std::vector<R>::const_iterator it_b = b.begin(); 386 typename std::vector<V>::iterator it_xb = x.begin(), it_xe = x.end(); 387 bzSorLowerSolver(it_d, it_l, it_b, it_xb, it_xe, w); 388 } 389 //@{ 390 //! specializations of lower triangular part solvers sorLowerSolver(const std::vector<real_t> & m,const std::vector<real_t> & b,std::vector<real_t> & x,const real_t w) const391 void sorLowerSolver(const std::vector<real_t>& m, const std::vector<real_t>& b, std::vector<real_t>& x, const real_t w) const 392 { sorLowerSolver<>(m, b, x, w); } sorLowerSolver(const std::vector<real_t> & m,const std::vector<complex_t> & b,std::vector<complex_t> & x,const real_t w) const393 void sorLowerSolver(const std::vector<real_t>& m, const std::vector<complex_t>& b, std::vector<complex_t>& x, const real_t w) const 394 { sorLowerSolver<>(m, b, x, w); } sorLowerSolver(const std::vector<complex_t> & m,const std::vector<real_t> & b,std::vector<complex_t> & x,const real_t w) const395 void sorLowerSolver(const std::vector<complex_t>& m, const std::vector<real_t>& b, std::vector<complex_t>& x, const real_t w) const 396 { sorLowerSolver<>(m, b, x, w); } sorLowerSolver(const std::vector<complex_t> & m,const std::vector<complex_t> & b,std::vector<complex_t> & x,const real_t w) const397 void sorLowerSolver(const std::vector<complex_t>& m, const std::vector<complex_t>& b, std::vector<complex_t>& x, const real_t w) const 398 { sorLowerSolver<>(m, b, x, w); } 399 //@} 400 401 /*! special template diagonal and triangular solvers (D/w+U) x = b used in SSOR algorithm 402 for D + L + U matrix splitting (D diagonal, L & U strict lower & upper triangular parts) */ 403 template<typename M, typename V, typename R> sorUpperSolver(const std::vector<M> & m,const std::vector<R> & b,std::vector<V> & x,const real_t w,const SymType sym) const404 void sorUpperSolver(const std::vector<M>& m, const std::vector<R>& b, 405 std::vector<V>& x, const real_t w, const SymType sym) const 406 { 407 // *** LO AND BEHOLD we use reverse_iterators here 408 typename std::vector<M>::const_reverse_iterator it_ru = m.rbegin(), it_rd = it_ru + upperPartSize() + lowerPartSize(); 409 typename std::vector<R>::const_reverse_iterator it_rbb = b.rbegin(); 410 typename std::vector<V>::reverse_iterator it_rxb = x.rbegin(), it_rxe = x.rend(); 411 bzSorUpperSolver(it_rd, it_ru, it_rbb, it_rxb, it_rxe, w); 412 } 413 //@{ 414 //! specializations of upper triangular part solvers sorUpperSolver(const std::vector<real_t> & m,const std::vector<real_t> & b,std::vector<real_t> & x,const real_t w,const SymType sym) const415 void sorUpperSolver(const std::vector<real_t>& m, const std::vector<real_t>& b, std::vector<real_t>& x, const real_t w, const SymType sym) const 416 { sorUpperSolver<>(m, b, x, w, sym); } sorUpperSolver(const std::vector<real_t> & m,const std::vector<complex_t> & b,std::vector<complex_t> & x,const real_t w,const SymType sym) const417 void sorUpperSolver(const std::vector<real_t>& m, const std::vector<complex_t>& b, std::vector<complex_t>& x, const real_t w, const SymType sym) const 418 { sorUpperSolver<>(m, b, x, w, sym); } sorUpperSolver(const std::vector<complex_t> & m,const std::vector<real_t> & b,std::vector<complex_t> & x,const real_t w,const SymType sym) const419 void sorUpperSolver(const std::vector<complex_t>& m, const std::vector<real_t>& b, std::vector<complex_t>& x, const real_t w, const SymType sym) const 420 { sorUpperSolver<>(m, b, x, w, sym); } sorUpperSolver(const std::vector<complex_t> & m,const std::vector<complex_t> & b,std::vector<complex_t> & x,const real_t w,const SymType sym) const421 void sorUpperSolver(const std::vector<complex_t>& m, const std::vector<complex_t>& b, std::vector<complex_t>& x, const real_t w, const SymType sym) const 422 { sorUpperSolver<>(m, b, x, w, sym); } 423 //@} 424 425 /* ---------------------------------------------------------------------------------------------------------------------------- 426 UMFPACK stuff 427 ----------------------------------------------------------------------------------------------------------------------------*/ 428 //! conversion to umfpack format 429 template<typename M, typename OrdinalType> 430 void toUmfPack(const std::vector<M>& values, std::vector<OrdinalType>& colPointer, std::vector<OrdinalType>& rowIndex, std::vector<M>& mat) const; 431 //@{ 432 //! specialization of 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) const433 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 434 { toUmfPack<>(values,colPointer,rowIndex,mat);} 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) const435 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 436 { toUmfPack<>(values,colPointer,rowIndex,mat);} 437 //@} 438 }; 439 440 } // end of namespace xlifepp 441 442 #endif // DUAL_DENSE_STORAGE_HPP 443