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