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_SymBandMatrixArithFunc_H 24 #define TMV_SymBandMatrixArithFunc_H 25 26 #include "tmv/TMV_BaseSymBandMatrix.h" 27 #include "tmv/TMV_BandMatrixArithFunc.h" 28 #include "tmv/TMV_Array.h" 29 30 #define CT std::complex<T> 31 32 namespace tmv { 33 34 // y (+)= alpha * A * x 35 template <bool add, typename T, typename Ta, typename Tx> 36 void MultMV( 37 const T alpha, const GenSymBandMatrix<Ta>& A, 38 const GenVector<Tx>& x, VectorView<T> y); 39 40 // A = alpha * A 41 template <typename T> MultXM(const T alpha,SymBandMatrixView<T> A)42 inline void MultXM(const T alpha, SymBandMatrixView<T> A) 43 { MultXM(alpha,A.upperBand()); } 44 45 // B += alpha * A 46 template <typename T, typename Ta> AddMM(const T alpha,const GenSymBandMatrix<Ta> & A,SymBandMatrixView<T> B)47 inline void AddMM( 48 const T alpha, const GenSymBandMatrix<Ta>& A, 49 SymBandMatrixView<T> B) 50 { AddMM(alpha,A.upperBand(),B.upperBand()); } 51 template <typename T, typename Ta> 52 void AddMM( 53 const T alpha, const GenSymBandMatrix<Ta>& A, 54 BandMatrixView<T> B); 55 template <typename T, typename Ta> AddMM(const T alpha,const GenSymBandMatrix<Ta> & A,MatrixView<T> B)56 inline void AddMM( 57 const T alpha, const GenSymBandMatrix<Ta>& A, 58 MatrixView<T> B) 59 { AddMM(alpha,A,BandMatrixViewOf(B,A.nlo(),A.nhi())); } 60 61 // C = alpha * A + beta * B 62 template <typename T, typename Ta, typename Tb> AddMM(const T alpha,const GenSymBandMatrix<Ta> & A,const T beta,const GenSymBandMatrix<Tb> & B,SymBandMatrixView<T> C)63 inline void AddMM( 64 const T alpha, const GenSymBandMatrix<Ta>& A, 65 const T beta, const GenSymBandMatrix<Tb>& B, 66 SymBandMatrixView<T> C) 67 { AddMM(alpha,A.upperBand(),beta,B.upperBand(),C.upperBand()); } 68 template <typename T, typename Ta, typename Tb> 69 void AddMM( 70 const T alpha, const GenSymBandMatrix<Ta>& A, 71 const T beta, const GenSymBandMatrix<Tb>& B, 72 BandMatrixView<T> C); 73 template <typename T, typename Ta, typename Tb> 74 void AddMM( 75 const T alpha, const GenSymBandMatrix<Ta>& A, 76 const T beta, const GenSymBandMatrix<Tb>& B, MatrixView<T> C); 77 template <typename T, typename Ta, typename Tb> 78 void AddMM( 79 const T alpha, const GenSymBandMatrix<Ta>& A, 80 const T beta, const GenBandMatrix<Tb>& B, BandMatrixView<T> C); 81 template <typename T, typename Ta, typename Tb> 82 void AddMM( 83 const T alpha, const GenSymBandMatrix<Ta>& A, 84 const T beta, const GenMatrix<Tb>& B, MatrixView<T> C); 85 template <typename T, typename Ta, typename Tb> AddMM(const T alpha,const GenBandMatrix<Ta> & A,const T beta,const GenSymBandMatrix<Tb> & B,BandMatrixView<T> C)86 inline void AddMM( 87 const T alpha, const GenBandMatrix<Ta>& A, 88 const T beta, const GenSymBandMatrix<Tb>& B, BandMatrixView<T> C) 89 { AddMM(beta,B,alpha,A,C); } 90 template <typename T, typename Ta, typename Tb> AddMM(const T alpha,const GenMatrix<Ta> & A,const T beta,const GenSymBandMatrix<Tb> & B,MatrixView<T> C)91 inline void AddMM( 92 const T alpha, const GenMatrix<Ta>& A, 93 const T beta, const GenSymBandMatrix<Tb>& B, MatrixView<T> C) 94 { AddMM(beta,B,alpha,A,C); } 95 96 97 // C (+)= alpha * A * B 98 template <bool add, typename T, typename Ta, typename Tb> 99 void MultMM( 100 const T alpha, const GenSymBandMatrix<Ta>& A, 101 const GenMatrix<Tb>& B, MatrixView<T> C); 102 template <bool add, typename T, typename Ta, typename Tb> MultMM(const T alpha,const GenMatrix<Ta> & A,const GenSymBandMatrix<Tb> & B,MatrixView<T> C)103 inline void MultMM( 104 const T alpha, const GenMatrix<Ta>& A, 105 const GenSymBandMatrix<Tb>& B, MatrixView<T> C) 106 { MultMM<add>(alpha,B.transpose(),A.transpose(),C.transpose()); } 107 template <bool add, typename T, typename Ta, typename Tb> 108 void MultMM( 109 const T alpha, const GenSymBandMatrix<Ta>& A, 110 const GenSymBandMatrix<Tb>& B, BandMatrixView<T> C); 111 template <bool add, typename T, typename Ta, typename Tb> 112 void MultMM( 113 const T alpha, const GenSymBandMatrix<Ta>& A, 114 const GenBandMatrix<Tb>& B, BandMatrixView<T> C); 115 template <bool add, typename T, typename Ta, typename Tb> MultMM(const T alpha,const GenBandMatrix<Ta> & A,const GenSymBandMatrix<Tb> & B,BandMatrixView<T> C)116 inline void MultMM( 117 const T alpha, const GenBandMatrix<Ta>& A, 118 const GenSymBandMatrix<Tb>& B, BandMatrixView<T> C) 119 { MultMM<add>(alpha,B.transpose(),A.transpose(),C.transpose()); } 120 121 template <typename T> 122 class SymBandMatrixComposite : public GenSymBandMatrix<T> 123 { 124 public : SymBandMatrixComposite()125 inline SymBandMatrixComposite() : itsm(0) {} SymBandMatrixComposite(const SymBandMatrixComposite<T> &)126 inline SymBandMatrixComposite(const SymBandMatrixComposite<T>&) : 127 itsm(0) {} ~SymBandMatrixComposite()128 virtual inline ~SymBandMatrixComposite() {} 129 130 // Definitions are in TMV_MultsBV.cpp 131 const T* cptr() const; 132 ptrdiff_t stepi() const; 133 ptrdiff_t stepj() const; 134 ptrdiff_t diagstep() const; 135 sym()136 inline SymType sym() const { return Sym; } uplo()137 inline UpLoType uplo() const { return Lower; } ct()138 inline ConjType ct() const { return NonConj; } 139 140 private: 141 mutable AlignedArray<T> itsm1; 142 mutable T* itsm; 143 144 SymBandMatrixComposite<T>& operator=(const SymBandMatrixComposite<T>&); 145 }; 146 147 148 template <typename T> 149 class SymBandMatrixComposite<CT> : 150 public BandMatrixComposite<CT>, 151 virtual public AssignableToSymBandMatrix<CT> 152 { 153 public : 154 // Need to respecify that size() is pure virtual here, 155 // since GenBandMatrix defines size for compaibility 156 // with AssignableToDiagMatrix, etc. 157 virtual ptrdiff_t size() const = 0; 158 159 using AssignableToSymMatrix<CT>::sym; 160 using AssignableToBandMatrix<CT>::nlo; nhi()161 inline ptrdiff_t nhi() const { return nlo(); } colsize()162 inline ptrdiff_t colsize() const { return size(); } rowsize()163 inline ptrdiff_t rowsize() const { return size(); } 164 assignToS(SymMatrixView<T> m0)165 inline void assignToS(SymMatrixView<T> m0) const 166 { 167 TMVAssert(m0.size() == size()); 168 TMVAssert(m0.sym() == sym()); 169 this->assignTosB(SymBandMatrixViewOf(m0,nlo())); 170 if (m0.size() > nlo()+1) 171 BandMatrixViewOf(m0.upperTri()).diagRange( 172 nlo()+1,m0.size()).setZero(); 173 } assignToS(SymMatrixView<CT> m0)174 inline void assignToS(SymMatrixView<CT> m0) const 175 { 176 TMVAssert(m0.size() == size()); 177 TMVAssert(m0.sym() == sym()); 178 this->assignTosB(SymBandMatrixViewOf(m0,nlo())); 179 if (m0.size() > nlo()+1) 180 BandMatrixViewOf(m0.upperTri()).diagRange( 181 nlo()+1,m0.size()).setZero(); 182 } 183 assignToM(MatrixView<T> m0)184 inline void assignToM(MatrixView<T> m0) const 185 { BandMatrixComposite<CT>::assignToM(m0); } assignToM(MatrixView<CT> m0)186 inline void assignToM(MatrixView<CT> m0) const 187 { BandMatrixComposite<CT>::assignToM(m0); } assignToB(BandMatrixView<T> m0)188 inline void assignToB(BandMatrixView<T> m0) const 189 { BandMatrixComposite<CT>::assignToB(m0); } assignToB(BandMatrixView<CT> m0)190 inline void assignToB(BandMatrixView<CT> m0) const 191 { BandMatrixComposite<CT>::assignToB(m0); } 192 }; 193 194 // Specialize allowed complex combinations: 195 template <bool add, typename T, typename Ta, typename Tx> MultMV(const T alpha,const GenSymBandMatrix<Ta> & A,const GenVector<Tx> & x,VectorView<CT> y)196 inline void MultMV( 197 const T alpha, const GenSymBandMatrix<Ta>& A, 198 const GenVector<Tx>& x, VectorView<CT> y) 199 { MultMV<add>(CT(alpha),A,x,y); } 200 201 template <typename T> MultXM(const T alpha,SymBandMatrixView<CT> A)202 inline void MultXM(const T alpha, SymBandMatrixView<CT> A) 203 { MultXM(CT(alpha),A); } 204 205 template <typename T, typename Ta> AddMM(const T alpha,const GenSymBandMatrix<Ta> & A,SymBandMatrixView<CT> B)206 inline void AddMM( 207 const T alpha, const GenSymBandMatrix<Ta>& A, 208 SymBandMatrixView<CT> B) 209 { AddMM(CT(alpha),A,B); } 210 template <typename T, typename Ta> AddMM(const T alpha,const GenSymBandMatrix<Ta> & A,BandMatrixView<CT> B)211 inline void AddMM( 212 const T alpha, const GenSymBandMatrix<Ta>& A, 213 BandMatrixView<CT> B) 214 { AddMM(CT(alpha),A,B); } 215 template <typename T, typename Ta> AddMM(const T alpha,const GenSymBandMatrix<Ta> & A,MatrixView<CT> B)216 inline void AddMM( 217 const T alpha, const GenSymBandMatrix<Ta>& A, MatrixView<CT> B) 218 { AddMM(CT(alpha),A,B); } 219 220 template <typename T, typename Ta, typename Tb> AddMM(const T alpha,const GenSymBandMatrix<Ta> & A,const T beta,const GenSymBandMatrix<Tb> & B,SymBandMatrixView<CT> C)221 inline void AddMM( 222 const T alpha, const GenSymBandMatrix<Ta>& A, 223 const T beta, const GenSymBandMatrix<Tb>& B, 224 SymBandMatrixView<CT> C) 225 { AddMM(CT(alpha),A,CT(beta),B,C); } 226 template <typename T, typename Ta, typename Tb> AddMM(const T alpha,const GenSymBandMatrix<Ta> & A,const T beta,const GenSymBandMatrix<Tb> & B,BandMatrixView<CT> C)227 inline void AddMM( 228 const T alpha, const GenSymBandMatrix<Ta>& A, 229 const T beta, const GenSymBandMatrix<Tb>& B, 230 BandMatrixView<CT> C) 231 { AddMM(CT(alpha),A,CT(beta),B,C); } 232 template <typename T, typename Ta, typename Tb> AddMM(const T alpha,const GenSymBandMatrix<Ta> & A,const T beta,const GenSymBandMatrix<Tb> & B,MatrixView<CT> C)233 inline void AddMM( 234 const T alpha, const GenSymBandMatrix<Ta>& A, 235 const T beta, const GenSymBandMatrix<Tb>& B, MatrixView<CT> C) 236 { AddMM(CT(alpha),A,CT(beta),B,C); } 237 template <typename T, typename Ta, typename Tb> AddMM(const T alpha,const GenSymBandMatrix<Ta> & A,const T beta,const GenBandMatrix<Tb> & B,BandMatrixView<CT> C)238 inline void AddMM( 239 const T alpha, const GenSymBandMatrix<Ta>& A, 240 const T beta, const GenBandMatrix<Tb>& B, BandMatrixView<CT> C) 241 { AddMM(CT(alpha),A,CT(beta),B,C); } 242 template <typename T, typename Ta, typename Tb> AddMM(const T alpha,const GenSymBandMatrix<Ta> & A,const T beta,const GenMatrix<Tb> & B,MatrixView<CT> C)243 inline void AddMM( 244 const T alpha, const GenSymBandMatrix<Ta>& A, 245 const T beta, const GenMatrix<Tb>& B, MatrixView<CT> C) 246 { AddMM(CT(alpha),A,CT(beta),B,C); } 247 template <typename T, typename Ta, typename Tb> AddMM(const T alpha,const GenBandMatrix<Ta> & A,const T beta,const GenSymBandMatrix<Tb> & B,BandMatrixView<CT> C)248 inline void AddMM( 249 const T alpha, const GenBandMatrix<Ta>& A, 250 const T beta, const GenSymBandMatrix<Tb>& B, 251 BandMatrixView<CT> C) 252 { AddMM(CT(alpha),A,CT(beta),B,C); } 253 template <typename T, typename Ta, typename Tb> AddMM(const T alpha,const GenMatrix<Ta> & A,const T beta,const GenSymBandMatrix<Tb> & B,MatrixView<CT> C)254 inline void AddMM( 255 const T alpha, const GenMatrix<Ta>& A, 256 const T beta, const GenSymBandMatrix<Tb>& B, MatrixView<CT> C) 257 { AddMM(CT(alpha),A,CT(beta),B,C); } 258 template <typename T> AddMM(const CT alpha,const GenSymBandMatrix<CT> & A,const CT beta,const GenSymBandMatrix<T> & B,MatrixView<CT> C)259 inline void AddMM( 260 const CT alpha, const GenSymBandMatrix<CT>& A, 261 const CT beta, const GenSymBandMatrix<T>& B, MatrixView<CT> C) 262 { AddMM(beta,B,alpha,A,C); } 263 template <typename T> AddMM(const CT alpha,const GenSymBandMatrix<CT> & A,const CT beta,const GenSymBandMatrix<T> & B,BandMatrixView<CT> C)264 inline void AddMM( 265 const CT alpha, const GenSymBandMatrix<CT>& A, 266 const CT beta, const GenSymBandMatrix<T>& B, 267 BandMatrixView<CT> C) 268 { AddMM(beta,B,alpha,A,C); } 269 270 template <bool add, typename T, typename Ta, typename Tb> MultMM(const T alpha,const GenSymBandMatrix<Ta> & A,const GenMatrix<Tb> & B,MatrixView<CT> C)271 inline void MultMM( 272 const T alpha, const GenSymBandMatrix<Ta>& A, 273 const GenMatrix<Tb>& B, MatrixView<CT> C) 274 { MultMM<add>(CT(alpha),A,B,C); } 275 template <bool add, typename T, typename Ta, typename Tb> MultMM(const T alpha,const GenMatrix<Ta> & A,const GenSymBandMatrix<Tb> & B,MatrixView<CT> C)276 inline void MultMM( 277 const T alpha, const GenMatrix<Ta>& A, 278 const GenSymBandMatrix<Tb>& B, MatrixView<CT> C) 279 { MultMM<add>(CT(alpha),A,B,C); } 280 template <bool add, typename T, typename Ta, typename Tb> MultMM(const T alpha,const GenSymBandMatrix<Ta> & A,const GenSymBandMatrix<Tb> & B,BandMatrixView<CT> C)281 inline void MultMM( 282 const T alpha, const GenSymBandMatrix<Ta>& A, 283 const GenSymBandMatrix<Tb>& B, BandMatrixView<CT> C) 284 { MultMM<add>(CT(alpha),A,B,C); } 285 template <bool add, typename T, typename Ta, typename Tb> MultMM(const T alpha,const GenSymBandMatrix<Ta> & A,const GenBandMatrix<Tb> & B,BandMatrixView<CT> C)286 inline void MultMM( 287 const T alpha, const GenSymBandMatrix<Ta>& A, 288 const GenBandMatrix<Tb>& B, BandMatrixView<CT> C) 289 { MultMM<add>(CT(alpha),A,B,C); } 290 template <bool add, typename T, typename Ta, typename Tb> MultMM(const T alpha,const GenBandMatrix<Ta> & A,const GenSymBandMatrix<Tb> & B,BandMatrixView<CT> C)291 inline void MultMM( 292 const T alpha, const GenBandMatrix<Ta>& A, 293 const GenSymBandMatrix<Tb>& B, BandMatrixView<CT> C) 294 { MultMM<add>(CT(alpha),A,B,C); } 295 296 // Specialize disallowed complex combinations: 297 template <bool add, typename T, typename Ta, typename Tb> MultMV(const CT,const GenSymBandMatrix<Ta> &,const GenVector<Tb> &,VectorView<T>)298 inline void MultMV( 299 const CT , const GenSymBandMatrix<Ta>& , 300 const GenVector<Tb>& , VectorView<T> ) 301 { TMVAssert(TMV_FALSE); } 302 303 template <typename T> MultXM(const CT,SymBandMatrixView<T>)304 inline void MultXM( 305 const CT , SymBandMatrixView<T> ) 306 { TMVAssert(TMV_FALSE); } 307 308 template <typename T, typename Ta> AddMM(const CT,const GenSymBandMatrix<Ta> &,SymBandMatrixView<T>)309 inline void AddMM( 310 const CT , const GenSymBandMatrix<Ta>& , SymBandMatrixView<T> ) 311 { TMVAssert(TMV_FALSE); } 312 template <typename T, typename Ta> AddMM(const CT,const GenSymBandMatrix<Ta> &,BandMatrixView<T>)313 inline void AddMM( 314 const CT , const GenSymBandMatrix<Ta>& , BandMatrixView<T> ) 315 { TMVAssert(TMV_FALSE); } 316 template <typename T, typename Ta> AddMM(const CT,const GenSymBandMatrix<Ta> &,MatrixView<T>)317 inline void AddMM( 318 const CT , const GenSymBandMatrix<Ta>& , MatrixView<T> ) 319 { TMVAssert(TMV_FALSE); } 320 template <typename T, typename Ta, typename Tb> AddMM(const CT,const GenSymBandMatrix<Ta> &,const CT,const GenSymBandMatrix<Tb> &,SymBandMatrixView<T>)321 inline void AddMM( 322 const CT , const GenSymBandMatrix<Ta>& , 323 const CT , const GenSymBandMatrix<Tb>& , SymBandMatrixView<T> ) 324 { TMVAssert(TMV_FALSE); } 325 template <typename T, typename Ta, typename Tb> AddMM(const CT,const GenSymBandMatrix<Ta> &,const CT,const GenSymBandMatrix<Tb> &,BandMatrixView<T>)326 inline void AddMM( 327 const CT , const GenSymBandMatrix<Ta>& , 328 const CT , const GenSymBandMatrix<Tb>& , BandMatrixView<T> ) 329 { TMVAssert(TMV_FALSE); } 330 template <typename T, typename Ta, typename Tb> AddMM(const CT,const GenSymBandMatrix<Ta> &,const CT,const GenSymBandMatrix<Tb> &,MatrixView<T>)331 inline void AddMM( 332 const CT , const GenSymBandMatrix<Ta>& , 333 const CT , const GenSymBandMatrix<Tb>& , MatrixView<T> ) 334 { TMVAssert(TMV_FALSE); } 335 template <typename T, typename Ta, typename Tb> AddMM(const CT,const GenSymBandMatrix<Ta> &,const CT,const GenBandMatrix<Tb> &,BandMatrixView<T>)336 inline void AddMM( 337 const CT , const GenSymBandMatrix<Ta>& , 338 const CT , const GenBandMatrix<Tb>& , BandMatrixView<T> ) 339 { TMVAssert(TMV_FALSE); } 340 template <typename T, typename Ta, typename Tb> AddMM(const CT,const GenSymBandMatrix<Ta> &,const CT,const GenMatrix<Tb> &,MatrixView<T>)341 inline void AddMM( 342 const CT , const GenSymBandMatrix<Ta>& , 343 const CT , const GenMatrix<Tb>& , MatrixView<T> ) 344 { TMVAssert(TMV_FALSE); } 345 template <typename T, typename Ta, typename Tb> AddMM(const CT,const GenBandMatrix<Ta> &,const CT,const GenSymBandMatrix<Tb> &,BandMatrixView<T>)346 inline void AddMM( 347 const CT , const GenBandMatrix<Ta>& , 348 const CT , const GenSymBandMatrix<Tb>& , BandMatrixView<T> ) 349 { TMVAssert(TMV_FALSE); } 350 template <typename T, typename Ta, typename Tb> AddMM(const CT,const GenMatrix<Ta> &,const CT,const GenSymBandMatrix<Tb> &,MatrixView<T>)351 inline void AddMM( 352 const CT , const GenMatrix<Ta>& , 353 const CT , const GenSymBandMatrix<Tb>& , MatrixView<T> ) 354 { TMVAssert(TMV_FALSE); } 355 356 template <bool add, typename T, typename Ta, typename Tb> MultMM(const CT,const GenSymBandMatrix<Ta> &,const GenMatrix<Tb> &,MatrixView<T>)357 inline void MultMM( 358 const CT , const GenSymBandMatrix<Ta>& , 359 const GenMatrix<Tb>& , MatrixView<T> ) 360 { TMVAssert(TMV_FALSE); } 361 template <bool add, typename T, typename Ta, typename Tb> MultMM(const CT,const GenMatrix<Ta> &,const GenSymBandMatrix<Tb> &,MatrixView<T>)362 inline void MultMM( 363 const CT , const GenMatrix<Ta>& , 364 const GenSymBandMatrix<Tb>& , MatrixView<T> ) 365 { TMVAssert(TMV_FALSE); } 366 template <bool add, typename T, typename Ta, typename Tb> MultMM(const CT,const GenSymBandMatrix<Ta> &,const GenSymBandMatrix<Tb> &,BandMatrixView<T>)367 inline void MultMM( 368 const CT , const GenSymBandMatrix<Ta>& , 369 const GenSymBandMatrix<Tb>& , BandMatrixView<T> ) 370 { TMVAssert(TMV_FALSE); } 371 template <bool add, typename T, typename Ta, typename Tb> MultMM(const CT,const GenSymBandMatrix<Ta> &,const GenBandMatrix<Tb> &,BandMatrixView<T>)372 inline void MultMM( 373 const CT , const GenSymBandMatrix<Ta>& , 374 const GenBandMatrix<Tb>& , BandMatrixView<T> ) 375 { TMVAssert(TMV_FALSE); } 376 template <bool add, typename T, typename Ta, typename Tb> MultMM(const CT,const GenBandMatrix<Ta> &,const GenSymBandMatrix<Tb> &,BandMatrixView<T>)377 inline void MultMM( 378 const CT , const GenBandMatrix<Ta>& , 379 const GenSymBandMatrix<Tb>& , BandMatrixView<T> ) 380 { TMVAssert(TMV_FALSE); } 381 382 } // namespace tmv 383 384 #undef CT 385 386 #endif 387