1 /////////////////////////////////////////////////////////////////////////////// 2 // // 3 // The Template Matrix/Vector Library for C++ was created by Mike Jarvis // 4 // Copyright (C) 1998 - 2016 // 5 // All rights reserved // 6 // // 7 // The project is hosted at https://code.google.com/p/tmv-cpp/ // 8 // where you can find the current version and current documention. // 9 // // 10 // For concerns or problems with the software, Mike may be contacted at // 11 // mike_jarvis17 [at] gmail. // 12 // // 13 // This software is licensed under a FreeBSD license. The file // 14 // TMV_LICENSE should have bee included with this distribution. // 15 // It not, you can get a copy from https://code.google.com/p/tmv-cpp/. // 16 // // 17 // Essentially, you can use this software however you want provided that // 18 // you include the TMV_LICENSE file in any distribution that uses it. // 19 // // 20 /////////////////////////////////////////////////////////////////////////////// 21 22 23 #ifndef TMV_SymMatrixArithFunc_H 24 #define TMV_SymMatrixArithFunc_H 25 26 #include "tmv/TMV_BaseSymMatrix.h" 27 #include "tmv/TMV_VectorArithFunc.h" 28 #include "tmv/TMV_MatrixArithFunc.h" 29 #include "tmv/TMV_TriMatrixArithFunc.h" 30 #include "tmv/TMV_Array.h" 31 32 #define CT std::complex<T> 33 34 namespace tmv { 35 36 // y (+)= alpha * A * x 37 template <bool add, typename T, typename Ta, typename Tx> 38 void MultMV( 39 const T alpha, const GenSymMatrix<Ta>& A, 40 const GenVector<Tx>& x, VectorView<T> y); 41 42 // A = alpha * A 43 template <typename T> MultXM(const T alpha,SymMatrixView<T> A)44 inline void MultXM(const T alpha, SymMatrixView<T> A) 45 { MultXM(alpha,A.upperTri()); } 46 47 // B += alpha * A 48 template <typename T, typename Ta> AddMM(const T alpha,const GenSymMatrix<Ta> & A,SymMatrixView<T> B)49 inline void AddMM( 50 const T alpha, const GenSymMatrix<Ta>& A, 51 SymMatrixView<T> B) 52 { AddMM(alpha,A.upperTri(),B.upperTri()); } 53 template <typename T, typename Ta> 54 void AddMM( 55 const T alpha, const GenSymMatrix<Ta>& A, MatrixView<T> B); 56 57 // C = alpha * A + beta * B 58 template <typename T, typename Ta, typename Tb> AddMM(const T alpha,const GenSymMatrix<Ta> & A,const T beta,const GenSymMatrix<Tb> & B,SymMatrixView<T> C)59 inline void AddMM( 60 const T alpha, const GenSymMatrix<Ta>& A, 61 const T beta, const GenSymMatrix<Tb>& B, SymMatrixView<T> C) 62 { AddMM(alpha,A.upperTri(),beta,B.upperTri(),C.upperTri()); } 63 template <typename T, typename Ta, typename Tb> 64 void AddMM( 65 const T alpha, const GenSymMatrix<Ta>& A, 66 const T beta, const GenSymMatrix<Tb>& B, MatrixView<T> C); 67 template <typename T, typename Ta, typename Tb> 68 void AddMM( 69 const T alpha, const GenSymMatrix<Ta>& A, 70 const T beta, const GenMatrix<Tb>& B, MatrixView<T> C); 71 template <typename T, typename Ta, typename Tb> AddMM(const T alpha,const GenMatrix<Ta> & A,const T beta,const GenSymMatrix<Tb> & B,MatrixView<T> C)72 inline void AddMM( 73 const T alpha, const GenMatrix<Ta>& A, 74 const T beta, const GenSymMatrix<Tb>& B, MatrixView<T> C) 75 { AddMM(beta,B,alpha,A,C); } 76 77 // C (+)= alpha * A * B 78 template <bool add, typename T, typename Ta, typename Tb> 79 void MultMM( 80 const T alpha, const GenSymMatrix<Ta>& A, 81 const GenMatrix<Tb>& B, MatrixView<T> C); 82 template <bool add, typename T, typename Ta, typename Tb> MultMM(const T alpha,const GenMatrix<Ta> & A,const GenSymMatrix<Tb> & B,MatrixView<T> C)83 inline void MultMM( 84 const T alpha, const GenMatrix<Ta>& A, 85 const GenSymMatrix<Tb>& B, MatrixView<T> C) 86 { MultMM<add>(alpha,B.transpose(),A.transpose(),C.transpose()); } 87 template <bool add, typename T, typename Ta, typename Tb> 88 void MultMM( 89 const T alpha, const GenSymMatrix<Ta>& A, 90 const GenSymMatrix<Tb>& B, MatrixView<T> ); 91 92 // A (+)= alpha * (x^xT) (or x* if A is Herm) 93 template <bool add, typename T, typename Tx> 94 void Rank1Update( 95 const T alpha, const GenVector<Tx>& x, SymMatrixView<T> A); 96 97 // B (+)= alpha * (A * AT) (or At if B is Herm) 98 template <bool add, typename T, typename Ta> 99 void RankKUpdate( 100 const T alpha, const GenMatrix<Ta>& A, SymMatrixView<T> B); 101 102 template <bool add, typename T, typename Ta> 103 void RankKUpdate( 104 const T alpha, const GenLowerTriMatrix<Ta>& A, 105 SymMatrixView<T> B); 106 template <bool add, typename T, typename Ta> 107 void RankKUpdate( 108 const T alpha, const GenUpperTriMatrix<Ta>& A, 109 SymMatrixView<T> B); 110 111 // These two don't have += forms: they must called explicitly 112 // A (+)= alpha * (x^y + y^x) (or x^y* + y^x* is A is Herm) 113 template <bool add, typename T, typename Tx, typename Ty> 114 void Rank2Update( 115 const T alpha, const GenVector<Tx>& x, 116 const GenVector<Ty>& y, SymMatrixView<T> A); 117 template <bool add, typename T, typename Tx, typename Ty, int A3> Rank2Update(const T alpha,const GenVector<Tx> & x,const GenVector<Ty> & y,SymMatrix<T,A3> & A)118 inline void Rank2Update( 119 const T alpha, const GenVector<Tx>& x, 120 const GenVector<Ty>& y, SymMatrix<T,A3>& A) 121 { Rank2Update(alpha,x,y,A.view()); } 122 template <bool add, typename T, typename Tx, typename Ty, int A3> Rank2Update(const T alpha,const GenVector<Tx> & x,const GenVector<Ty> & y,HermMatrix<T,A3> & A)123 inline void Rank2Update( 124 const T alpha, const GenVector<Tx>& x, 125 const GenVector<Ty>& y, HermMatrix<T,A3>& A) 126 { Rank2Update(alpha,x,y,A.view()); } 127 128 // C (+)= alpha * (A * BT + B*AT) (or A*Bt + B*At if C is Herm) 129 template <bool add, typename T, typename Ta, typename Tb> 130 void Rank2KUpdate( 131 const T alpha, const GenMatrix<Ta>& A, 132 const GenMatrix<Tb>& B, SymMatrixView<T> C); 133 template <bool add, typename T, typename Ta, typename Tb, int A3> Rank2KUpdate(const T alpha,const GenMatrix<Ta> & A,const GenMatrix<Tb> & B,SymMatrix<T,A3> & C)134 inline void Rank2KUpdate( 135 const T alpha, const GenMatrix<Ta>& A, 136 const GenMatrix<Tb>& B, SymMatrix<T,A3>& C) 137 { Rank2KUpdate(alpha,A,B,C.view()); } 138 template <bool add, typename T, typename Ta, typename Tb, int A3> Rank2KUpdate(const T alpha,const GenMatrix<Ta> & A,const GenMatrix<Tb> & B,HermMatrix<T,A3> & C)139 inline void Rank2KUpdate( 140 const T alpha, const GenMatrix<Ta>& A, 141 const GenMatrix<Tb>& B, HermMatrix<T,A3>& C) 142 { Rank2KUpdate(alpha,A,B,C.view()); } 143 144 // C (+)= alpha * A * B 145 // This also needs to be called explicitly. 146 // This is to prevent the programmer from doing this accidentally 147 // when alpha * A * B is not necessarily symmetrix/hermitian. 148 template <bool add, typename T, typename Ta, typename Tb> 149 void SymMultMM( 150 const T alpha, const GenMatrix<Ta>& A, 151 const GenMatrix<Tb>& B, SymMatrixView<T> C); 152 template <bool add, typename T, typename Ta, typename Tb, int A3> SymMultMM(const T alpha,const GenMatrix<Ta> & A,const GenMatrix<Tb> & B,SymMatrix<T,A3> & C)153 inline void SymMultMM( 154 const T alpha, const GenMatrix<Ta>& A, 155 const GenMatrix<Tb>& B, SymMatrix<T,A3>& C) 156 { SymMultMM(alpha,A,B,C.view()); } 157 template <bool add, typename T, typename Ta, typename Tb, int A3> SymMultMM(const T alpha,const GenMatrix<Ta> & A,const GenMatrix<Tb> & B,HermMatrix<T,A3> & C)158 inline void SymMultMM( 159 const T alpha, const GenMatrix<Ta>& A, 160 const GenMatrix<Tb>& B, HermMatrix<T,A3>& C) 161 { SymMultMM(alpha,A,B,C.view()); } 162 163 template <typename T> 164 class SymMatrixComposite : public GenSymMatrix<T> 165 { 166 public : SymMatrixComposite()167 inline SymMatrixComposite() {} SymMatrixComposite(const SymMatrixComposite<T> &)168 inline SymMatrixComposite(const SymMatrixComposite<T>&) {} ~SymMatrixComposite()169 virtual inline ~SymMatrixComposite() {} 170 171 // Definitions are in TMV_MultSV.cpp 172 const T* cptr() const; 173 ptrdiff_t stepi() const; 174 ptrdiff_t stepj() const; 175 sym()176 inline SymType sym() const { return Sym; } uplo()177 inline UpLoType uplo() const { return Lower; } ct()178 inline ConjType ct() const { return NonConj; } 179 180 private : 181 mutable AlignedArray<T> itsm; 182 183 SymMatrixComposite<T>& operator=(const SymMatrixComposite<T>&); 184 }; 185 186 template <typename T> 187 class SymMatrixComposite<CT> : 188 public MatrixComposite<CT>, 189 virtual public AssignableToSymMatrix<CT> 190 { 191 public : 192 colsize()193 inline ptrdiff_t colsize() const { return this->size(); } rowsize()194 inline ptrdiff_t rowsize() const { return this->size(); } 195 assignToM(MatrixView<TMV_RealType (T)> m0)196 inline void assignToM(MatrixView<TMV_RealType(T)> m0) const 197 { MatrixComposite<CT>::assignToM(m0); } assignToM(MatrixView<TMV_ComplexType (T)> m0)198 inline void assignToM(MatrixView<TMV_ComplexType(T)> m0) const 199 { MatrixComposite<CT>::assignToM(m0); } 200 }; 201 202 // Specialize allowed complex combinations: 203 template <bool add, typename T, typename Ta, typename Tx> MultMV(const T alpha,const GenSymMatrix<Ta> & A,const GenVector<Tx> & x,VectorView<CT> y)204 inline void MultMV( 205 const T alpha, const GenSymMatrix<Ta>& A, 206 const GenVector<Tx>& x, VectorView<CT> y) 207 { MultMV<add>(CT(alpha),A,x,y); } 208 209 template <typename T> MultXM(const T alpha,SymMatrixView<CT> A)210 inline void MultXM(const T alpha, SymMatrixView<CT> A) 211 { MultXM(CT(alpha),A); } 212 213 template <typename T, typename Ta> AddMM(const T alpha,const GenSymMatrix<Ta> & A,SymMatrixView<CT> B)214 inline void AddMM( 215 const T alpha, const GenSymMatrix<Ta>& A, SymMatrixView<CT> B) 216 { AddMM(CT(alpha),A,B); } 217 template <typename T, typename Ta> AddMM(const T alpha,const GenSymMatrix<Ta> & A,MatrixView<CT> B)218 inline void AddMM( 219 const T alpha, const GenSymMatrix<Ta>& A, MatrixView<CT> B) 220 { AddMM(CT(alpha),A,B); } 221 template <typename T, typename Ta, typename Tb> AddMM(const T alpha,const GenSymMatrix<Ta> & A,const T beta,const GenSymMatrix<Tb> & B,SymMatrixView<CT> C)222 inline void AddMM( 223 const T alpha, const GenSymMatrix<Ta>& A, 224 const T beta, const GenSymMatrix<Tb>& B, SymMatrixView<CT> C) 225 { AddMM(CT(alpha),A,CT(beta),B,C); } 226 template <typename T, typename Ta, typename Tb> AddMM(const T alpha,const GenSymMatrix<Ta> & A,const T beta,const GenSymMatrix<Tb> & B,MatrixView<CT> C)227 inline void AddMM( 228 const T alpha, const GenSymMatrix<Ta>& A, 229 const T beta, const GenSymMatrix<Tb>& B, MatrixView<CT> C) 230 { AddMM(CT(alpha),A,CT(beta),B,C); } 231 template <typename T, typename Ta, typename Tb> AddMM(const T alpha,const GenSymMatrix<Ta> & A,const T beta,const GenMatrix<Tb> & B,MatrixView<CT> C)232 inline void AddMM( 233 const T alpha, const GenSymMatrix<Ta>& A, 234 const T beta, const GenMatrix<Tb>& B, MatrixView<CT> C) 235 { AddMM(CT(alpha),A,CT(beta),B,C); } 236 template <typename T, typename Ta, typename Tb> AddMM(const T alpha,const GenMatrix<Ta> & A,const T beta,const GenSymMatrix<Tb> & B,MatrixView<CT> C)237 inline void AddMM( 238 const T alpha, const GenMatrix<Ta>& A, 239 const T beta, const GenSymMatrix<Tb>& B, MatrixView<CT> C) 240 { AddMM(beta,B,alpha,A,C); } 241 template <typename T> AddMM(const CT alpha,const GenSymMatrix<CT> & A,const CT beta,const GenSymMatrix<T> & B,MatrixView<CT> C)242 inline void AddMM( 243 const CT alpha, const GenSymMatrix<CT>& A, 244 const CT beta, const GenSymMatrix<T>& B, MatrixView<CT> C) 245 { AddMM(beta,B,alpha,A,C); } 246 247 template <bool add, typename T, typename Ta, typename Tb> MultMM(const T alpha,const GenSymMatrix<Ta> & A,const GenMatrix<Tb> & B,MatrixView<CT> C)248 inline void MultMM( 249 const T alpha, const GenSymMatrix<Ta>& A, 250 const GenMatrix<Tb>& B, MatrixView<CT> C) 251 { MultMM<add>(CT(alpha),A,B,C); } 252 template <bool add, typename T, typename Ta, typename Tb> MultMM(const T alpha,const GenMatrix<Ta> & A,const GenSymMatrix<Tb> & B,MatrixView<CT> C)253 inline void MultMM( 254 const T alpha, const GenMatrix<Ta>& A, 255 const GenSymMatrix<Tb>& B, MatrixView<CT> C) 256 { MultMM<add>(CT(alpha),A,B,C); } 257 template <bool add, typename T, typename Ta, typename Tb> MultMM(const T alpha,const GenSymMatrix<Ta> & A,const GenSymMatrix<Tb> & B,MatrixView<CT> C)258 inline void MultMM( 259 const T alpha, const GenSymMatrix<Ta>& A, 260 const GenSymMatrix<Tb>& B, MatrixView<CT> C) 261 { MultMM<add>(CT(alpha),A,B,C); } 262 263 template <bool add, typename T, typename Tx> Rank1Update(const T alpha,const GenVector<Tx> & x,SymMatrixView<CT> A)264 inline void Rank1Update( 265 const T alpha, const GenVector<Tx>& x, SymMatrixView<CT> A) 266 { Rank1Update<add>(CT(alpha),x,A); } 267 template <bool add, typename T, typename Ta> RankKUpdate(const T alpha,const GenMatrix<Ta> & A,SymMatrixView<CT> B)268 inline void RankKUpdate( 269 const T alpha, const GenMatrix<Ta>& A, SymMatrixView<CT> B) 270 { RankKUpdate<add>(CT(alpha),A,B); } 271 template <bool add, typename T, typename Ta> RankKUpdate(const T alpha,const GenLowerTriMatrix<Ta> & A,SymMatrixView<CT> B)272 inline void RankKUpdate( 273 const T alpha, const GenLowerTriMatrix<Ta>& A, 274 SymMatrixView<CT> B) 275 { RankKUpdate<add>(CT(alpha),A,B); } 276 template <bool add, typename T, typename Ta> RankKUpdate(const T alpha,const GenUpperTriMatrix<Ta> & A,SymMatrixView<CT> B)277 inline void RankKUpdate( 278 const T alpha, const GenUpperTriMatrix<Ta>& A, 279 SymMatrixView<CT> B) 280 { RankKUpdate<add>(CT(alpha),A,B); } 281 282 template <bool add, typename T, typename Tx, typename Ty> Rank2Update(const T alpha,const GenVector<Tx> & x,const GenVector<Ty> & y,SymMatrixView<CT> A)283 inline void Rank2Update( 284 const T alpha, const GenVector<Tx>& x, 285 const GenVector<Ty>& y, SymMatrixView<CT> A) 286 { Rank2Update<add>(CT(alpha),x,y,A); } 287 template <bool add, typename T, typename Ta, typename Tb> Rank2KUpdate(const T alpha,const GenMatrix<Ta> & A,const GenMatrix<Tb> & B,SymMatrixView<CT> C)288 inline void Rank2KUpdate( 289 const T alpha, const GenMatrix<Ta>& A, 290 const GenMatrix<Tb>& B, SymMatrixView<CT> C) 291 { Rank2KUpdate<add>(CT(alpha),A,B,C); } 292 293 template <bool add, typename T, typename Ta, typename Tb> SymMultMM(const T alpha,const GenMatrix<Ta> & A,const GenMatrix<Tb> & B,SymMatrixView<CT> C)294 inline void SymMultMM( 295 const T alpha, const GenMatrix<Ta>& A, 296 const GenMatrix<Tb>& B, SymMatrixView<CT> C) 297 { SymMultMM<add>(CT(alpha),A,B,C); } 298 299 // Specialize disallowed complex combinations: 300 template <bool add, typename T, typename Ta, typename Tb> MultMV(const CT,const GenSymMatrix<Ta> &,const GenVector<Tb> &,VectorView<T>)301 inline void MultMV( 302 const CT , const GenSymMatrix<Ta>& , 303 const GenVector<Tb>& , VectorView<T> ) 304 { TMVAssert(TMV_FALSE); } 305 306 template <typename T> MultXM(const CT,SymMatrixView<T>)307 inline void MultXM(const CT , SymMatrixView<T> ) 308 { TMVAssert(TMV_FALSE); } 309 310 template <typename T, typename Ta, typename Tb> AddMM(const CT,const GenSymMatrix<Ta> &,const CT,const GenSymMatrix<Tb> &,SymMatrixView<T>)311 inline void AddMM( 312 const CT , const GenSymMatrix<Ta>& , 313 const CT , const GenSymMatrix<Tb>& , SymMatrixView<T> ) 314 { TMVAssert(TMV_FALSE); } 315 template <typename T, typename Ta, typename Tb> AddMM(const CT,const GenSymMatrix<Ta> &,const CT,const GenSymMatrix<Tb> &,MatrixView<T>)316 inline void AddMM( 317 const CT , const GenSymMatrix<Ta>& , 318 const CT , const GenSymMatrix<Tb>& , MatrixView<T> ) 319 { TMVAssert(TMV_FALSE); } 320 template <typename T, typename Ta, typename Tb> AddMM(const CT,const GenSymMatrix<Ta> &,const CT,const GenMatrix<Tb> &,MatrixView<T>)321 inline void AddMM( 322 const CT , const GenSymMatrix<Ta>& , 323 const CT , const GenMatrix<Tb>& , MatrixView<T> ) 324 { TMVAssert(TMV_FALSE); } 325 template <typename T, typename Ta, typename Tb> AddMM(const CT,const GenMatrix<Ta> &,const CT,const GenSymMatrix<Tb> &,MatrixView<T>)326 inline void AddMM( 327 const CT , const GenMatrix<Ta>& , 328 const CT , const GenSymMatrix<Tb>& , MatrixView<T> ) 329 { TMVAssert(TMV_FALSE); } 330 331 template <bool add, typename T, typename Ta, typename Tb> MultMM(const CT,const GenSymMatrix<Ta> &,const GenMatrix<Tb> &,MatrixView<T>)332 inline void MultMM( 333 const CT , const GenSymMatrix<Ta>& , 334 const GenMatrix<Tb>& , MatrixView<T> ) 335 { TMVAssert(TMV_FALSE); } 336 template <bool add, typename T, typename Ta, typename Tb> MultMM(const CT,const GenMatrix<Ta> &,const GenSymMatrix<Tb> &,MatrixView<T>)337 inline void MultMM( 338 const CT , const GenMatrix<Ta>& , 339 const GenSymMatrix<Tb>& , MatrixView<T> ) 340 { TMVAssert(TMV_FALSE); } 341 template <bool add, typename T, typename Ta, typename Tb> MultMM(const CT,const GenSymMatrix<Ta> &,const GenSymMatrix<Tb> &,MatrixView<T>)342 inline void MultMM( 343 const CT , const GenSymMatrix<Ta>& , 344 const GenSymMatrix<Tb>& , MatrixView<T> ) 345 { TMVAssert(TMV_FALSE); } 346 347 348 } // namespace tmv 349 350 #undef CT 351 352 #endif 353