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 //---------------------------------------------------------------------------
24 //
25 // This file defines the TMV DiagMatrix class.
26 //
27 // The DiagMatrix class is provided for efficient storage of a diagonal
28 // matrix.  You can do most of the things that you can do with a
29 // regular Matrix, but it will do them more efficiently.
30 //
31 // Constructors:
32 //
33 //    DiagMatrix<T,A>(int size)
34 //        Makes a DiagMatrix with column size and row size = size
35 //        with _uninitialized_ values
36 //
37 //    DiagMatrix<T,A>(int size, T x)
38 //        Makes a DiagMatrix of size n with all values = x
39 //
40 //    DiagMatrix<T,A>(const Vector<T>& vv)
41 //        Make a DiagMatrix which copies the elements of vv.
42 //
43 //
44 // Access Functions
45 //
46 //    int colsize() const
47 //    int rowsize() const
48 //    int size() const
49 //        Return the dimensions of the DiagMatrix
50 //
51 //    T& operator()(int i)
52 //    T operator()(int i) const
53 //    T& operator()(int i, int j)
54 //    T operator()(int i, int j) const
55 //        Return the (i,j) element of the DiagMatrix
56 //        For the single paramter version, j=i
57 //
58 //    VectorView diag()
59 //    ConstVectorView diag() const
60 //        Return the diagonal of the DiagMatrix as a VectorView
61 //
62 //
63 // Modifying Functions - The same as the regular Matrix counterparts
64 //
65 //    DiagMatrix& setZero()
66 //    DiagMatrix& setAllTo(T x)
67 //    DiagMatrix& addToAll(T x)
68 //    DiagMatrix<T>& transposeSelf()
69 //        (Does nothing.)
70 //    DiagMatrix& conjugateSelf()
71 //    DiagMatrix& setToIdentity(x = 1)
72 //    void Swap(DiagMatrix& m1, DiagMatrix& m2)
73 //
74 //
75 // subDiagMatrix:
76 //
77 //    DiagMatrixView subDiagMatrix(int i1, int i2, int istep=1)
78 //        Returns a Sub-DiagMatrix which extends from i1 to i2 (step istep)
79 //        which refers to the same physical elements as the original.
80 //        As usual, i2 is the "one past the end" element.
81 //
82 //    DiagMatrixView realPart()
83 //    DiagMatrixView imagPart()
84 //        Returns a view to the real/imag elements of a complex DiagMatrix.
85 //
86 //    DiagMatrixView view()
87 //    DiagMatrixView transpose()
88 //        Returns a view of a DiagMatrix.
89 //
90 //    DiagMatrixView conjugate(m)
91 //    DiagMatrixView adjoint(m)
92 //        Returns a view to the conjugate of a DiagMatrix.
93 //
94 // Functions of DiagMatrices - Same as for regular Matrices:
95 //
96 //    m.det() or Det(m)
97 //    m.logDet() or m.logDet(T* sign) or LogDet(m)
98 //    m.trace() or Trace(m)
99 //    m.sumElements() or SumElements(m)
100 //    m.sumAbsElements() or SumAbsElements(m)
101 //    m.sumAbs2Elements() or SumAbs2Elements(m)
102 //    m.norm() or m.normF() or Norm(m) or NormF(m)
103 //    m.normSq() or NormSq(m)
104 //    m.norm1() or Norm1(m)
105 //    m.norm2() or Norm2(m)
106 //    m.normInf() or NormInf(m)
107 //    m.maxAbsElement() or MaxAbsElements(m)
108 //    m.maxAbs2Element() or MaxAbs2Elements(m)
109 //        (Note - for diagonal matrices,
110 //        norm1 = norm2 = normInf = maxAbsElement.)
111 //
112 //    m.inverse() or Inverse(m)
113 //    m.invertSelf()
114 //    m.makeInverse(minv) (takes either Matrix or DiagMatrix argument)
115 //    m.makeInverseATA(invata) (takes either Matrix or DiagMatrix argument)
116 //
117 // I/O:
118 //
119 //    os << d
120 //        Writes d to ostream os as a full matrix
121 //
122 //    os << CompactIO() << d
123 //        Writes only the diagonal Vector to os
124 //
125 //    is >> d
126 //    is >> CompactIO() >> d
127 //        Reads in d in either format
128 //
129 //
130 
131 #ifndef TMV_DiagMatrix_H
132 #define TMV_DiagMatrix_H
133 
134 #include "tmv/TMV_BaseDiagMatrix.h"
135 #include "tmv/TMV_BaseTriMatrix.h"
136 #include "tmv/TMV_Vector.h"
137 #include "tmv/TMV_TriMatrix.h"
138 #include "tmv/TMV_Matrix.h"
139 #include <vector>
140 
141 namespace tmv {
142 
143     template <typename T>
144     class GenDiagMatrix :
145         virtual public AssignableToDiagMatrix<T>,
146         virtual public AssignableToUpperTriMatrix<T>,
147         virtual public AssignableToLowerTriMatrix<T>,
148         public BaseMatrix<T>
149     {
150     public:
151 
152         typedef TMV_RealType(T) RT;
153         typedef TMV_ComplexType(T) CT;
154         typedef T value_type;
155         typedef RT real_type;
156         typedef CT complex_type;
157         typedef GenDiagMatrix<T> type;
158         typedef DiagMatrix<T> copy_type;
159         typedef ConstDiagMatrixView<T> const_view_type;
160         typedef const_view_type const_transpose_type;
161         typedef const_view_type const_conjugate_type;
162         typedef const_view_type const_adjoint_type;
163         typedef ConstDiagMatrixView<RT> const_realpart_type;
164         typedef ConstVectorView<T> const_vec_type;
165         typedef DiagMatrixView<T> nonconst_type;
166         typedef typename const_vec_type::const_iterator const_iterator;
167 
168         //
169         // Constructors
170         //
171 
172         inline GenDiagMatrix<T>() {}
GenDiagMatrix(const type &)173         inline GenDiagMatrix(const type&) {}
~GenDiagMatrix()174         virtual inline ~GenDiagMatrix() {}
175 
176         //
177         // Access Functions
178         //
179 
size()180         inline ptrdiff_t size() const { return diag().size(); }
colsize()181         inline ptrdiff_t colsize() const { return size(); }
rowsize()182         inline ptrdiff_t rowsize() const { return size(); }
dt()183         inline DiagType dt() const { return NonUnitDiag; }
184 
operator()185         inline T operator()(ptrdiff_t i, ptrdiff_t j) const
186         {
187             TMVAssert(i>=0 && i<size());
188             TMVAssert(j>=0 && j<size());
189             if (i==j) return diag()(i);
190             else return T(0);
191         }
192 
operator()193         inline T operator()(ptrdiff_t i) const
194         {
195             TMVAssert(i>=0 && i<size());
196             return diag()(i);
197         }
198 
diag()199         inline const_vec_type diag() const { return cdiag(); }
200 
201         template <typename T2>
isSameAs(const BaseMatrix<T2> &)202         inline bool isSameAs(const BaseMatrix<T2>& ) const
203         { return false; }
204 
isSameAs(const GenDiagMatrix<T> & m2)205         inline bool isSameAs(const GenDiagMatrix<T>& m2) const
206         {
207             if (this == &m2) return true;
208             else return (diag().isSameAs(m2.diag()));
209         }
210 
assignToM(MatrixView<RT> m2)211         inline void assignToM(MatrixView<RT> m2) const
212         {
213             TMVAssert(m2.colsize() == size());
214             TMVAssert(m2.rowsize() == size());
215             TMVAssert(isReal(T()));
216             m2.diag() = diag();
217             m2.upperTri().offDiag().setZero();
218             m2.lowerTri().offDiag().setZero();
219         }
assignToM(MatrixView<CT> m2)220         inline void assignToM(MatrixView<CT> m2) const
221         {
222             TMVAssert(m2.colsize() == size());
223             TMVAssert(m2.rowsize() == size());
224             m2.diag() = diag();
225             m2.upperTri().offDiag().setZero();
226             m2.lowerTri().offDiag().setZero();
227         }
228 
assignToU(UpperTriMatrixView<RT> m2)229         inline void assignToU(UpperTriMatrixView<RT> m2) const
230         {
231             TMVAssert(m2.size() == size());
232             TMVAssert(isReal(T()));
233             m2.diag() = diag();
234             m2.offDiag().setZero();
235         }
assignToU(UpperTriMatrixView<CT> m2)236         inline void assignToU(UpperTriMatrixView<CT> m2) const
237         {
238             TMVAssert(m2.size() == size());
239             m2.diag() = diag();
240             m2.offDiag().setZero();
241         }
242 
assignToL(LowerTriMatrixView<RT> m2)243         inline void assignToL(LowerTriMatrixView<RT> m2) const
244         {
245             TMVAssert(m2.size() == size());
246             TMVAssert(isReal(T()));
247             m2.diag() = diag();
248             m2.offDiag().setZero();
249         }
assignToL(LowerTriMatrixView<CT> m2)250         inline void assignToL(LowerTriMatrixView<CT> m2) const
251         {
252             TMVAssert(m2.size() == size());
253             m2.diag() = diag();
254             m2.offDiag().setZero();
255         }
256 
assignToD(DiagMatrixView<RT> m2)257         inline void assignToD(DiagMatrixView<RT> m2) const
258         {
259             TMVAssert(m2.size() == size());
260             TMVAssert(isReal(T()));
261             m2.diag() = diag();
262         }
assignToD(DiagMatrixView<CT> m2)263         inline void assignToD(DiagMatrixView<CT> m2) const
264         {
265             TMVAssert(m2.size() == size());
266             m2.diag() = diag();
267         }
268 
269         //
270         // subDiagMatrix
271         //
272 
cSubDiagMatrix(ptrdiff_t i1,ptrdiff_t i2)273         inline const_view_type cSubDiagMatrix(ptrdiff_t i1, ptrdiff_t i2) const
274         { return const_view_type(diag().cSubVector(i1,i2)); }
275 
subDiagMatrix(ptrdiff_t i1,ptrdiff_t i2)276         inline const_view_type subDiagMatrix(ptrdiff_t i1, ptrdiff_t i2) const
277         {
278             TMVAssert(diag().hasSubVector(i1,i2,1));
279             return cSubDiagMatrix(i1,i2);
280         }
281 
cSubDiagMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t istep)282         inline const_view_type cSubDiagMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const
283         { return const_view_type(diag().cSubVector(i1,i2,istep)); }
284 
subDiagMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t istep)285         inline const_view_type subDiagMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const
286         {
287             TMVAssert(diag().hasSubVector(i1,i2,istep));
288             return cSubDiagMatrix(i1,i2,istep);
289         }
290 
realPart()291         inline const_realpart_type realPart() const
292         { return const_realpart_type(diag().realPart()); }
293 
imagPart()294         inline const_realpart_type imagPart() const
295         { return const_realpart_type(diag().imagPart()); }
296 
view()297         inline const_view_type view() const
298         { return const_view_type(diag()); }
299 
transpose()300         inline const_view_type transpose() const
301         { return view(); }
302 
conjugate()303         inline const_view_type conjugate() const
304         { return const_view_type(diag().conjugate()); }
305 
adjoint()306         inline const_view_type adjoint() const
307         { return conjugate(); }
308 
nonConst()309         inline nonconst_type nonConst() const
310         { return nonconst_type(diag().nonConst()); }
311 
312 
313         //
314         // Functions of DiagMatrix
315         //
316 
317         T det() const;
318 
319         RT logDet(T* sign=0) const;
320 
trace()321         inline T trace() const
322         { return diag().sumElements(); }
323 
sumElements()324         inline T sumElements() const
325         { return diag().sumElements(); }
326 
sumAbsElements()327         inline RT sumAbsElements() const
328         { return diag().sumAbsElements(); }
329 
sumAbs2Elements()330         inline RT sumAbs2Elements() const
331         { return diag().sumAbs2Elements(); }
332 
norm()333         inline RT norm() const
334         { return normF(); }
335 
normF()336         inline RT normF() const
337         { return diag().norm(); }
338 
339         inline RT normSq(RT scale = RT(1)) const
340         { return diag().normSq(scale); }
341 
norm1()342         inline RT norm1() const
343         { return diag().maxAbsElement(); }
344 
doNorm2()345         inline RT doNorm2() const
346         { return diag().maxAbsElement(); }
347 
doCondition()348         inline RT doCondition() const
349         { return diag().maxAbsElement()/diag().minAbsElement(); }
350 
norm2()351         inline RT norm2() const
352         { return doNorm2(); }
353 
condition()354         inline RT condition() const
355         { return doCondition(); }
356 
normInf()357         inline RT normInf() const
358         { return diag().maxAbsElement(); }
359 
maxAbsElement()360         inline RT maxAbsElement() const
361         { return diag().maxAbsElement(); }
362 
maxAbs2Element()363         inline RT maxAbs2Element() const
364         { return diag().maxAbs2Element(); }
365 
isSingular()366         inline bool isSingular() const
367         { return diag().minAbsElement() == RT(0); }
368 
369         template <typename T1>
370         void doMakeInverse(MatrixView<T1> minv) const;
371 
372         template <typename T1>
373         void doMakeInverse(DiagMatrixView<T1> minv) const;
374 
375         void doMakeInverseATA(MatrixView<T> ata) const;
376         void doMakeInverseATA(DiagMatrixView<T> ata) const;
377 
makeInverse(MatrixView<T> minv)378         inline void makeInverse(MatrixView<T> minv) const
379         {
380             TMVAssert(minv.colsize() == size());
381             TMVAssert(minv.rowsize() == size());
382             doMakeInverse(minv);
383         }
384 
385         template <typename T1>
makeInverse(MatrixView<T1> minv)386         inline void makeInverse(MatrixView<T1> minv) const
387         {
388             TMVAssert(minv.colsize() == size());
389             TMVAssert(minv.rowsize() == size());
390             doMakeInverse(minv);
391         }
392 
393         template <typename T1>
makeInverse(DiagMatrixView<T1> minv)394         inline void makeInverse(DiagMatrixView<T1> minv) const
395         {
396             TMVAssert(minv.size() == size());
397             doMakeInverse(minv);
398         }
399 
400         QuotXD<T,T> QInverse() const;
inverse()401         inline QuotXD<T,T> inverse() const
402         { return QInverse(); }
403 
makeInverseATA(MatrixView<T> ata)404         inline void makeInverseATA(MatrixView<T> ata) const
405         {
406             TMVAssert(ata.colsize() == size());
407             TMVAssert(ata.rowsize() == size());
408             doMakeInverseATA(ata);
409         }
makeInverseATA(DiagMatrixView<T> ata)410         inline void makeInverseATA(DiagMatrixView<T> ata) const
411         {
412             TMVAssert(ata.size() == size());
413             doMakeInverseATA(ata);
414         }
415 
416         template <typename T1, int A>
makeInverse(DiagMatrix<T1,A> & minv)417         inline void makeInverse(DiagMatrix<T1,A>& minv) const
418         { makeInverse(minv.view()); }
419 
420         template <typename T1, int A>
makeInverse(Matrix<T1,A> & minv)421         inline void makeInverse(Matrix<T1,A>& minv) const
422         { makeInverse(minv.view()); }
423 
424         template <int A>
makeInverseATA(DiagMatrix<T,A> & minv)425         inline void makeInverseATA(DiagMatrix<T,A>& minv) const
426         { makeInverseATA(minv.view()); }
427 
428         template <int A>
makeInverseATA(Matrix<T,A> & minv)429         inline void makeInverseATA(Matrix<T,A>& minv) const
430         { makeInverseATA(minv.view()); }
431 
432 
433         //
434         // I/O
435         //
436 
437         void write(const TMV_Writer& writer) const;
438 
439         //
440         // Arithmetic Helpers
441         //
442 
443         template <typename T1>
444         void doLDivEq(VectorView<T1> v) const;
445 
446         template <typename T1, typename T0>
447         void doLDiv(const GenVector<T1>& v1, VectorView<T0> v0) const;
448 
449         template <typename T1>
450         void doLDivEq(MatrixView<T1> m) const;
451 
452         template <typename T1, typename T0>
453         void doLDiv(const GenMatrix<T1>& m1, MatrixView<T0> m0) const;
454 
455         template <typename T1>
LDivEq(VectorView<T1> v)456         inline void LDivEq(VectorView<T1> v) const
457         {
458             TMVAssert(v.size() == size());
459             doLDivEq(v);
460         }
461 
462         template <typename T1, typename T0>
LDiv(const GenVector<T1> & v1,VectorView<T0> v0)463         inline void LDiv(
464             const GenVector<T1>& v1, VectorView<T0> v0) const
465         {
466             TMVAssert(v1.size() == size());
467             TMVAssert(v0.size() == size());
468             doLDiv(v1,v0);
469         }
470 
471         template <typename T1>
RDivEq(VectorView<T1> v)472         inline void RDivEq(VectorView<T1> v) const
473         { LDivEq(v); }
474 
475         template <typename T1, typename T0>
RDiv(const GenVector<T1> & v1,VectorView<T0> v0)476         inline void RDiv(
477             const GenVector<T1>& v1, VectorView<T0> v0) const
478         { LDiv(v1,v0); }
479 
480         template <typename T1>
LDivEq(MatrixView<T1> m)481         inline void LDivEq(MatrixView<T1> m) const
482         {
483             TMVAssert(m.colsize() == size());
484             doLDivEq(m);
485         }
486         template <typename T1, typename T0>
LDiv(const GenMatrix<T1> & m1,MatrixView<T0> m0)487         inline void LDiv(
488             const GenMatrix<T1>& m1, MatrixView<T0> m0) const
489         {
490             TMVAssert(m1.colsize() == size());
491             TMVAssert(m0.colsize() == size());
492             TMVAssert(m1.rowsize() == m0.rowsize());
493             doLDiv(m1,m0);
494         }
495         template <typename T1>
RDivEq(MatrixView<T1> m)496         inline void RDivEq(MatrixView<T1> m) const
497         { LDivEq(m.transpose()); }
498 
499         template <typename T1, typename T0>
RDiv(const GenMatrix<T1> & m1,MatrixView<T0> m0)500         void RDiv(const GenMatrix<T1>& m1, MatrixView<T0> m0) const
501         { LDiv(m1.transpose(),m0.transpose()); }
502 
503         template <typename T1>
DivEq(DiagMatrixView<T1> m)504         void DivEq(DiagMatrixView<T1> m) const
505         { LDivEq(m.diag()); }
506 
507         template <typename T1, typename T0>
Div(const GenDiagMatrix<T1> & m1,DiagMatrixView<T0> m0)508         void Div(
509             const GenDiagMatrix<T1>& m1, DiagMatrixView<T0> m0) const
510         { LDiv(m1.diag(),m0.diag()); }
511 
512         // For easier compatibility with regular matrices:
divideInPlace()513         inline void divideInPlace() const {}
saveDiv()514         inline void saveDiv() const {}
setDiv()515         inline void setDiv() const {}
unsetDiv()516         inline void unsetDiv() const {}
resetDiv()517         inline void resetDiv() const {}
divIsSet()518         inline bool divIsSet() const { return true; }
divideUsing(DivType TMV_DEBUGPARAM (dt))519         inline void divideUsing(DivType TMV_DEBUGPARAM(dt)) const
520         { TMVAssert(dt == LU); }
521         inline bool checkDecomp(std::ostream* fout=0) const { return true; }
522         inline bool checkDecomp(
523             const BaseMatrix<T>& m2, std::ostream* fout=0) const
524         { return true; }
525 
cref(ptrdiff_t i,ptrdiff_t j)526         inline T cref(ptrdiff_t i, ptrdiff_t j) const
527         { return i==j ? cdiag()(i) : 0; }
cptr()528         inline const T* cptr() const
529         { return cdiag().cptr(); }
step()530         inline ptrdiff_t step() const
531         { return cdiag().step(); }
isconj()532         inline bool isconj() const
533         { return cdiag().isconj(); }
534 
begin()535         inline const_iterator begin() const
536         { return diag().begin(); }
end()537         inline const_iterator end() const
538         { return diag().end(); }
539 
540     protected :
541 
542         virtual ConstVectorView<T> cdiag() const = 0;
543 
544     private :
545 
546         type& operator=(const type&);
547 
548     }; // GenDiagMatrix
549 
550     template <typename T, int A>
551     class ConstDiagMatrixView : public GenDiagMatrix<T>
552     {
553     public :
554 
555         typedef GenDiagMatrix<T> base;
556         typedef ConstDiagMatrixView<T,A> type;
557         typedef ConstVectorView<T> const_vec_type;
558 
ConstDiagMatrixView(const type & rhs)559         inline ConstDiagMatrixView(const type& rhs) : itsdiag(rhs.cdiag())
560         { TMVAssert(Attrib<A>::viewok); }
561 
ConstDiagMatrixView(const base & rhs)562         inline ConstDiagMatrixView(const base& rhs) : itsdiag(rhs.diag())
563         { TMVAssert(Attrib<A>::viewok); }
564 
ConstDiagMatrixView(const GenVector<T> & v)565         explicit inline ConstDiagMatrixView(const GenVector<T>& v) : itsdiag(v)
566         { TMVAssert(Attrib<A>::viewok); }
567 
~ConstDiagMatrixView()568         virtual inline ~ConstDiagMatrixView() {}
569 
570     protected :
571 
572         ConstVectorView<T> itsdiag;
cdiag()573         inline ConstVectorView<T> cdiag() const { return itsdiag; }
574 
575     private :
576 
577         type& operator=(const type&);
578 
579     }; // ConstDiagMatrixView
580 
581     template <typename T>
582     class ConstDiagMatrixView<T,FortranStyle> :
583         public ConstDiagMatrixView<T,CStyle>
584     {
585     public :
586 
587         typedef TMV_RealType(T) RT;
588         typedef GenDiagMatrix<T> base;
589         typedef ConstDiagMatrixView<T,FortranStyle> type;
590         typedef ConstVectorView<T,FortranStyle> const_vec_type;
591         typedef ConstDiagMatrixView<T,FortranStyle> const_view_type;
592         typedef const_view_type const_transpose_type;
593         typedef const_view_type const_conjugate_type;
594         typedef const_view_type const_adjoint_type;
595         typedef ConstDiagMatrixView<RT,FortranStyle> const_realpart_type;
596         typedef ConstDiagMatrixView<T,CStyle> c_type;
597         typedef DiagMatrixView<T,FortranStyle> nonconst_type;
598 
ConstDiagMatrixView(const type & rhs)599         inline ConstDiagMatrixView(const type& rhs) : c_type(rhs) {}
600 
ConstDiagMatrixView(const base & rhs)601         inline ConstDiagMatrixView(const base& rhs) : c_type(rhs) {}
602 
ConstDiagMatrixView(const GenVector<T> & v)603         explicit inline ConstDiagMatrixView(const GenVector<T>& v) :
604             c_type(v) {}
605 
~ConstDiagMatrixView()606         virtual inline ~ConstDiagMatrixView() {}
607 
608         //
609         // Access Functions
610         //
611 
operator()612         inline T operator()(ptrdiff_t i, ptrdiff_t j) const
613         {
614             TMVAssert(i>0 && i<=c_type::size());
615             TMVAssert(j>0 && j<=c_type::size());
616             if (i==j) return diag()(i);
617             else return T(0);
618         }
619 
operator()620         inline T operator()(ptrdiff_t i) const
621         {
622             TMVAssert(i>0 && i<=c_type::size());
623             return diag()(i);
624         }
625 
diag()626         inline const_vec_type diag() const
627         { return const_vec_type(c_type::diag()); }
628 
629         //
630         // subDiagMatrix
631         //
632 
subDiagMatrix(ptrdiff_t i1,ptrdiff_t i2)633         inline const_view_type subDiagMatrix(ptrdiff_t i1, ptrdiff_t i2) const
634         {
635             TMVAssert(diag().hasSubVector(i1,i2,1));
636             return base::cSubDiagMatrix(i1-1,i2);
637         }
638 
subDiagMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t istep)639         inline const_view_type subDiagMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const
640         {
641             TMVAssert(diag().hasSubVector(i1,i2,istep));
642             return base::cSubDiagMatrix(i1-1,i2-1+istep,istep);
643         }
644 
realPart()645         inline const_realpart_type realPart() const
646         { return const_realpart_type(diag().realPart()); }
647 
imagPart()648         inline const_realpart_type imagPart() const
649         { return const_realpart_type(diag().imagPart()); }
650 
view()651         inline const_view_type view() const
652         { return const_view_type(diag()); }
653 
transpose()654         inline const_view_type transpose() const
655         { return view(); }
656 
conjugate()657         inline const_view_type conjugate() const
658         { return const_view_type(diag().conjugate()); }
659 
adjoint()660         inline const_view_type adjoint() const
661         { return conjugate(); }
662 
nonConst()663         inline nonconst_type nonConst() const
664         { return nonconst_type(diag().nonConst()); }
665 
666         type& operator=(const type&);
667 
668     }; // ConstDiagMatrixView - FortranStyle
669 
670     template <typename T, int A>
671     class DiagMatrixView : public GenDiagMatrix<T>
672     {
673     public:
674 
675         typedef TMV_RealType(T) RT;
676         typedef TMV_ComplexType(T) CT;
677         typedef GenDiagMatrix<T> base;
678         typedef DiagMatrixView<T,A> type;
679         typedef DiagMatrixView<T,A> view_type;
680         typedef view_type transpose_type;
681         typedef view_type conjugate_type;
682         typedef view_type adjoint_type;
683         typedef DiagMatrixView<RT,A> realpart_type;
684         typedef VectorView<T,A> vec_type;
685         typedef ConstDiagMatrixView<T,A> const_view_type;
686         typedef const_view_type const_transpose_type;
687         typedef const_view_type const_conjugate_type;
688         typedef const_view_type const_adjoint_type;
689         typedef ConstDiagMatrixView<RT,A> const_realpart_type;
690         typedef ConstVectorView<T,A> const_vec_type;
691         typedef typename RefHelper<T>::reference reference;
692         typedef typename vec_type::iterator iterator;
693         typedef typename const_vec_type::const_iterator const_iterator;
694 
695         //
696         // Constructors
697         //
698 
DiagMatrixView(const type & rhs)699         inline DiagMatrixView(const type& rhs) : itsdiag(rhs.itsdiag)
700         { TMVAssert(Attrib<A>::diagmatrixok); }
701 
DiagMatrixView(vec_type _diag)702         explicit inline DiagMatrixView(vec_type _diag) : itsdiag(_diag)
703         { TMVAssert(Attrib<A>::diagmatrixok); }
704 
~DiagMatrixView()705         virtual inline ~DiagMatrixView() {}
706 
707         //
708         // Op=
709         //
710 
711         inline type& operator=(const type& m2)
712         { m2.assignToD(*this); return *this; }
713 
714         inline type& operator=(const GenDiagMatrix<RT>& m2)
715         { m2.assignToD(*this); return *this; }
716 
717         inline type& operator=(const GenDiagMatrix<CT>& m2)
718         { m2.assignToD(*this); return *this; }
719 
720         template <typename T2>
721         inline type& operator=(const GenDiagMatrix<T2>& m2)
722         { itsdiag = m2.diag(); return *this; }
723 
724         inline type& operator=(const T& x)
725         { return setToIdentity(x); }
726 
727         inline type& operator=(const AssignableToDiagMatrix<RT>& m2)
728         {
729             TMVAssert(size() == m2.size());
730             m2.assignToD(*this);
731             return *this;
732         }
733 
734         inline type& operator=(const AssignableToDiagMatrix<CT>& m2)
735         {
736             TMVAssert(size() == m2.size());
737             m2.assignToD(*this);
738             return *this;
739         }
740 
741         typedef ListAssigner<T,iterator> MyListAssigner;
742         inline MyListAssigner operator<<(const T& x)
743         { return MyListAssigner(begin(),size(),x); }
744 
745 
746         //
747         // Access
748         //
749 
operator()750         inline reference operator()(ptrdiff_t i)
751         {
752             TMVAssert(i>=0 && i<size());
753             return diag()(i);
754         }
operator()755         inline reference operator()(ptrdiff_t i, ptrdiff_t TMV_DEBUGPARAM(j))
756         {
757             TMVAssert(i>=0 && i<size());
758             TMVAssert(i==j);
759             return diag()(i);
760         }
761 
diag()762         inline vec_type diag() { return itsdiag; }
763 
764         // Repeat const versions
operator()765         inline T operator()(ptrdiff_t i) const
766         { return base::operator()(i); }
operator()767         inline T operator()(ptrdiff_t i, ptrdiff_t j) const
768         { return base::operator()(i,j); }
diag()769         inline const_vec_type diag() const
770         { return base::diag(); }
771 
772         //
773         // Modifying Functions
774         //
775 
setZero()776         inline type& setZero()
777         { diag().setZero(); return *this; }
778 
setAllTo(const T & x)779         inline type& setAllTo(const T& x)
780         { diag().setAllTo(x); return *this; }
781 
addToAll(const T & x)782         inline type& addToAll(const T& x)
783         { diag().addToAll(x); return *this; }
784 
clip(RT thresh)785         inline type& clip(RT thresh)
786         { diag().clip(thresh); return *this; }
787 
transposeSelf()788         inline type& transposeSelf()
789         { return *this; }
790 
conjugateSelf()791         inline type& conjugateSelf()
792         { diag().conjugateSelf(); return *this; }
793 
794         type& invertSelf();
795 
796         inline type& setToIdentity(const T& x=T(1))
797         { return setAllTo(x); }
798 
799         //
800         // subDiagMatrix
801         //
802 
cSubDiagMatrix(ptrdiff_t i1,ptrdiff_t i2)803         inline view_type cSubDiagMatrix(ptrdiff_t i1, ptrdiff_t i2)
804         { return view_type(itsdiag.cSubVector(i1,i2)); }
805 
subDiagMatrix(ptrdiff_t i1,ptrdiff_t i2)806         inline view_type subDiagMatrix(ptrdiff_t i1, ptrdiff_t i2)
807         {
808             TMVAssert(diag().hasSubVector(i1,i2,1));
809             return cSubDiagMatrix(i1,i2);
810         }
811 
cSubDiagMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t istep)812         inline view_type cSubDiagMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep)
813         { return view_type(itsdiag.cSubVector(i1,i2,istep)); }
814 
subDiagMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t istep)815         inline view_type subDiagMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep)
816         {
817             TMVAssert(diag().hasSubVector(i1,i2,istep));
818             return cSubDiagMatrix(i1,i2,istep);
819         }
820 
realPart()821         inline realpart_type realPart()
822         { return realpart_type(diag().realPart()); }
823 
imagPart()824         inline realpart_type imagPart()
825         { return realpart_type(diag().imagPart()); }
826 
view()827         inline view_type view()
828         { return *this; }
829 
transpose()830         inline view_type transpose()
831         { return *this; }
832 
conjugate()833         inline view_type conjugate()
834         { return view_type(diag().conjugate()); }
835 
adjoint()836         inline view_type adjoint()
837         { return conjugate(); }
838 
cSubDiagMatrix(ptrdiff_t i1,ptrdiff_t i2)839         inline const_view_type cSubDiagMatrix(ptrdiff_t i1, ptrdiff_t i2) const
840         { return base::cSubDiagMatrix(i1,i2); }
subDiagMatrix(ptrdiff_t i1,ptrdiff_t i2)841         inline const_view_type subDiagMatrix(ptrdiff_t i1, ptrdiff_t i2) const
842         { return base::subDiagMatrix(i1,i2); }
cSubDiagMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t istep)843         inline const_view_type cSubDiagMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const
844         { return base::cSubDiagMatrix(i1,i2,istep); }
subDiagMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t istep)845         inline const_view_type subDiagMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const
846         { return base::subDiagMatrix(i1,i2,istep); }
realPart()847         inline const_realpart_type realPart() const
848         { return base::realPart(); }
imagPart()849         inline const_realpart_type imagPart() const
850         { return base::imagPart(); }
view()851         inline const_view_type view() const
852         { return base::view(); }
transpose()853         inline const_view_type transpose() const
854         { return base::transpose(); }
conjugate()855         inline const_view_type conjugate() const
856         { return base::conjugate(); }
adjoint()857         inline const_view_type adjoint() const
858         { return base::adjoint(); }
859 
860         //
861         // I/O
862         //
863 
864         void read(const TMV_Reader& reader);
865 
866         using GenDiagMatrix<T>::size;
867 
begin()868         inline iterator begin()
869         { return diag().begin(); }
end()870         inline iterator end()
871         { return diag().end(); }
872 
begin()873         inline const_iterator begin() const
874         { return diag().begin(); }
end()875         inline const_iterator end() const
876         { return diag().end(); }
877 
878     protected:
879 
880         VectorView<T> itsdiag;
cdiag()881         inline ConstVectorView<T> cdiag() const { return itsdiag; }
882 
883     }; // DiagMatrixView
884 
885     template <typename T>
886     class DiagMatrixView<T,FortranStyle> : public DiagMatrixView<T,CStyle>
887     {
888     public:
889 
890         typedef TMV_RealType(T) RT;
891         typedef TMV_ComplexType(T) CT;
892         typedef GenDiagMatrix<T> base;
893         typedef DiagMatrixView<T,FortranStyle> type;
894         typedef DiagMatrixView<T,CStyle> c_type;
895         typedef DiagMatrixView<T,FortranStyle> view_type;
896         typedef view_type transpose_type;
897         typedef view_type conjugate_type;
898         typedef view_type adjoint_type;
899         typedef DiagMatrixView<RT,FortranStyle> realpart_type;
900         typedef VectorView<T,FortranStyle> vec_type;
901         typedef ConstDiagMatrixView<T,FortranStyle> const_view_type;
902         typedef const_view_type const_transpose_type;
903         typedef const_view_type const_conjugate_type;
904         typedef const_view_type const_adjoint_type;
905         typedef ConstDiagMatrixView<RT,FortranStyle> const_realpart_type;
906         typedef ConstVectorView<T,FortranStyle> const_vec_type;
907         typedef typename RefHelper<T>::reference reference;
908 
909         //
910         // Constructors
911         //
912 
DiagMatrixView(const type & rhs)913         inline DiagMatrixView(const type& rhs) : c_type(rhs) {}
914 
DiagMatrixView(const c_type & rhs)915         inline DiagMatrixView(const c_type& rhs) : c_type(rhs) {}
916 
DiagMatrixView(VectorView<T> _diag)917         explicit inline DiagMatrixView(VectorView<T> _diag) :
918             c_type(_diag) {}
919 
~DiagMatrixView()920         virtual inline ~DiagMatrixView() {}
921 
922         //
923         // Op=
924         //
925 
926         inline type& operator=(const type& m2)
927         { c_type::operator=(m2); return *this; }
928 
929         inline type& operator=(const c_type& m2)
930         { c_type::operator=(m2); return *this; }
931 
932         inline type& operator=(const GenDiagMatrix<RT>& m2)
933         { c_type::operator=(m2); return *this; }
934 
935         inline type& operator=(const GenDiagMatrix<CT>& m2)
936         { c_type::operator=(m2); return *this; }
937 
938         template <typename T2>
939         inline type& operator=(const GenDiagMatrix<T2>& m2)
940         { c_type::operator=(m2); return *this; }
941 
942         inline type& operator=(const T& x)
943         { c_type::operator=(x); return *this; }
944 
945         inline type& operator=(const AssignableToDiagMatrix<RT>& m2)
946         { c_type::operator=(m2); return *this; }
947 
948         inline type& operator=(const AssignableToDiagMatrix<CT>& m2)
949         { c_type::operator=(m2); return *this; }
950 
951         typedef typename c_type::MyListAssigner MyListAssigner;
952         inline MyListAssigner operator<<(const T& x)
953         { return c_type::operator<<(x); }
954 
955 
956         //
957         // Access
958         //
959 
operator()960         inline reference operator()(ptrdiff_t i)
961         {
962             TMVAssert(i>0 && i<=c_type::size());
963             return diag()(i);
964         }
operator()965         inline reference operator()(ptrdiff_t i, ptrdiff_t TMV_DEBUGPARAM(j))
966         {
967             TMVAssert(i==j);
968             TMVAssert(i>0 && i<=c_type::size());
969             return diag()(i);
970         }
971 
diag()972         inline vec_type diag()
973         { return c_type::diag(); }
974 
operator()975         inline T operator()(ptrdiff_t i) const
976         {
977             TMVAssert(i>0 && i<=c_type::size());
978             return diag()(i);
979         }
operator()980         inline T operator()(ptrdiff_t i, ptrdiff_t TMV_DEBUGPARAM(j)) const
981         {
982             TMVAssert(i==j);
983             TMVAssert(i>0 && i<=c_type::size());
984             return diag()(i);
985         }
986 
diag()987         inline const_vec_type diag() const
988         { return c_type::diag(); }
989 
990         //
991         // Modifying Functions
992         //
993 
setZero()994         inline type& setZero()
995         { diag().setZero(); return *this; }
996 
setAllTo(const T & x)997         inline type& setAllTo(const T& x)
998         { diag().setAllTo(x); return *this; }
999 
addToAll(const T & x)1000         inline type& addToAll(const T& x)
1001         { diag().addToAll(x); return *this; }
1002 
clip(RT thresh)1003         inline type& clip(RT thresh)
1004         { diag().clip(thresh); return *this; }
1005 
transposeSelf()1006         inline type& transposeSelf()
1007         { return *this; }
1008 
conjugateSelf()1009         inline type& conjugateSelf()
1010         { diag().conjugateSelf(); return *this; }
1011 
invertSelf()1012         inline type& invertSelf()
1013         { return c_type::invertSelf(); }
1014 
1015         inline type& setToIdentity(const T& x=T(1))
1016         { diag().setAllTo(x); return *this; }
1017 
1018 
1019         //
1020         // subDiagMatrix
1021         //
1022 
subDiagMatrix(ptrdiff_t i1,ptrdiff_t i2)1023         inline view_type subDiagMatrix(ptrdiff_t i1, ptrdiff_t i2)
1024         {
1025             TMVAssert(diag().hasSubVector(i1,i2,1));
1026             return c_type::cSubDiagMatrix(i1-1,i2);
1027         }
1028 
cSubDiagMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t istep)1029         inline view_type cSubDiagMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep)
1030         {
1031             TMVAssert(diag().hasSubVector(i1,i2,istep));
1032             return c_type::cSubDiagMatrix(i1-1,i2-1+istep,istep);
1033         }
1034 
subDiagMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t istep)1035         inline view_type subDiagMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep)
1036         {
1037             TMVAssert(diag().hasSubVector(i1,i2,istep));
1038             return view_type(diag().subVector(i1,i2,istep));
1039         }
1040 
realPart()1041         inline realpart_type realPart()
1042         { return realpart_type(diag().realPart()); }
1043 
imagPart()1044         inline realpart_type imagPart()
1045         { return realpart_type(diag().imagPart()); }
1046 
view()1047         inline view_type view()
1048         { return *this; }
1049 
transpose()1050         inline view_type transpose()
1051         { return *this; }
1052 
conjugate()1053         inline view_type conjugate()
1054         { return view_type(diag().conjugate()); }
1055 
adjoint()1056         inline view_type adjoint()
1057         { return conjugate(); }
1058 
subDiagMatrix(ptrdiff_t i1,ptrdiff_t i2)1059         inline const_view_type subDiagMatrix(ptrdiff_t i1, ptrdiff_t i2) const
1060         {
1061             TMVAssert(diag().hasSubVector(i1,i2,1));
1062             return c_type::cSubDiagMatrix(i1-1,i2);
1063         }
1064 
cSubDiagMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t istep)1065         inline const_view_type cSubDiagMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const
1066         {
1067             TMVAssert(diag().hasSubVector(i1,i2,istep));
1068             return c_type::cSubDiagMatrix(i1-1,i2-1+istep,istep);
1069         }
1070 
subDiagMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t istep)1071         inline const_view_type subDiagMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const
1072         {
1073             TMVAssert(diag().hasSubVector(i1,i2,istep));
1074             return view_type(diag().subVector(i1,i2,istep));
1075         }
1076 
realPart()1077         inline const_realpart_type realPart() const
1078         { return realpart_type(diag().realPart()); }
1079 
imagPart()1080         inline const_realpart_type imagPart() const
1081         { return realpart_type(diag().imagPart()); }
1082 
view()1083         inline const_view_type view() const
1084         { return *this; }
1085 
transpose()1086         inline const_view_type transpose() const
1087         { return *this; }
1088 
conjugate()1089         inline const_view_type conjugate() const
1090         { return view_type(diag().conjugate()); }
1091 
adjoint()1092         inline const_view_type adjoint() const
1093         { return conjugate(); }
1094 
1095     }; // FortranStyle DiagMatrixView
1096 
1097 
1098     template <typename T, int A>
1099     class DiagMatrix : public GenDiagMatrix<T>
1100     {
1101     public:
1102 
1103         typedef TMV_RealType(T) RT;
1104         typedef TMV_ComplexType(T) CT;
1105         typedef GenDiagMatrix<T> base;
1106         typedef DiagMatrix<T,A> type;
1107         typedef DiagMatrixView<T,A> view_type;
1108         typedef view_type transpose_type;
1109         typedef view_type conjugate_type;
1110         typedef view_type adjoint_type;
1111         typedef DiagMatrixView<RT,A> realpart_type;
1112         typedef VectorView<T,A> vec_type;
1113         typedef ConstDiagMatrixView<T,A> const_view_type;
1114         typedef const_view_type const_transpose_type;
1115         typedef const_view_type const_conjugate_type;
1116         typedef const_view_type const_adjoint_type;
1117         typedef ConstDiagMatrixView<RT,A> const_realpart_type;
1118         typedef ConstVectorView<T,A> const_vec_type;
1119         typedef T& reference;
1120         typedef typename const_vec_type::const_iterator const_iterator;
1121         typedef typename vec_type::iterator iterator;
1122 
1123         //
1124         // Constructors
1125         //
1126 
itsdiag(_size)1127         inline explicit DiagMatrix(ptrdiff_t _size=0) : itsdiag(_size)
1128         { TMVAssert(_size >= 0); }
1129 
DiagMatrix(ptrdiff_t _size,const T & x)1130         inline DiagMatrix(ptrdiff_t _size, const T& x) : itsdiag(_size,x)
1131         { TMVAssert(_size >= 0); }
1132 
DiagMatrix(const GenVector<T> & rhs)1133         inline explicit DiagMatrix(const GenVector<T>& rhs) : itsdiag(rhs) {}
1134 
DiagMatrix(const GenMatrix<T> & m)1135         inline explicit DiagMatrix(const GenMatrix<T>& m) : itsdiag(m.diag()) {}
1136 
DiagMatrix(const type & rhs)1137         inline DiagMatrix(const type& rhs) : itsdiag(rhs.diag()) {}
1138 
DiagMatrix(const GenDiagMatrix<RT> & rhs)1139         inline DiagMatrix(const GenDiagMatrix<RT>& rhs) : itsdiag(rhs.size())
1140         { rhs.assignToD(view()); }
1141 
DiagMatrix(const GenDiagMatrix<CT> & rhs)1142         inline DiagMatrix(const GenDiagMatrix<CT>& rhs) : itsdiag(rhs.size())
1143         {
1144             TMVAssert(isComplex(T()));
1145             rhs.assignToD(view());
1146         }
1147 
1148         template <typename T2>
DiagMatrix(const GenDiagMatrix<T2> & rhs)1149         inline DiagMatrix(const GenDiagMatrix<T2>& rhs) : itsdiag(rhs.diag()) {}
1150 
DiagMatrix(const AssignableToDiagMatrix<RT> & m2)1151         inline DiagMatrix(const AssignableToDiagMatrix<RT>& m2) :
1152             itsdiag(m2.colsize())
1153         {
1154             TMVAssert(m2.colsize() == m2.rowsize());
1155             m2.assignToD(view());
1156         }
1157 
DiagMatrix(const AssignableToDiagMatrix<CT> & m2)1158         inline DiagMatrix(const AssignableToDiagMatrix<CT>& m2) :
1159             itsdiag(m2.colsize())
1160         {
1161             TMVAssert(m2.colsize() == m2.rowsize());
1162             TMVAssert(isComplex(T()));
1163             m2.assignToD(view());
1164         }
1165 
~DiagMatrix()1166         virtual inline ~DiagMatrix() {}
1167 
1168 
1169         //
1170         // Op=
1171         //
1172 
1173         inline type& operator=(const type& m2)
1174         {
1175             TMVAssert(m2.size() == size());
1176             m2.assignToD(view());
1177             return *this;
1178         }
1179 
1180         inline type& operator=(const GenDiagMatrix<RT>& m2)
1181         {
1182             TMVAssert(m2.size() == size());
1183             m2.assignToD(view());
1184             return *this;
1185         }
1186 
1187         inline type& operator=(const GenDiagMatrix<CT>& m2)
1188         {
1189             TMVAssert(m2.size() == size());
1190             TMVAssert(isComplex(T()));
1191             m2.assignToD(view());
1192             return *this;
1193         }
1194 
1195         template <typename T2>
1196         inline type& operator=(const GenDiagMatrix<T2>& m2)
1197         {
1198             TMVAssert(m2.size() == size());
1199             view() = m2;
1200             return *this;
1201         }
1202 
1203         inline type& operator=(const T& x) { view() = x; return *this; }
1204 
1205         inline type& operator=(const AssignableToDiagMatrix<RT>& m2)
1206         {
1207             TMVAssert(m2.size() == size());
1208             m2.assignToD(view());
1209             return *this;
1210         }
1211 
1212         inline type& operator=(const AssignableToDiagMatrix<CT>& m2)
1213         {
1214             TMVAssert(m2.size() == size());
1215             TMVAssert(isComplex(T()));
1216             m2.assignToD(view());
1217             return *this;
1218         }
1219 
1220         typedef ListAssigner<T,iterator> MyListAssigner;
1221         inline MyListAssigner operator<<(const T& x)
1222         { return MyListAssigner(begin(),size(),x); }
1223 
1224 
1225         //
1226         // Access
1227         //
1228 
operator()1229         inline T& operator()(ptrdiff_t i)
1230         {
1231             if (A==CStyle) {
1232                 TMVAssert(i>=0 && i<size());
1233                 return itsdiag(i);
1234             } else {
1235                 TMVAssert(i>0 && i<=size());
1236                 return itsdiag(i-1);
1237             }
1238         }
1239 
operator()1240         inline T& operator()(ptrdiff_t i, ptrdiff_t TMV_DEBUGPARAM(j))
1241         {
1242             if (A==CStyle) {
1243                 TMVAssert(i>=0 && i<size());
1244                 TMVAssert(j>=0 && j<size());
1245             } else {
1246                 TMVAssert(i>0 && i<=size());
1247                 TMVAssert(j>0 && j<=size());
1248             }
1249             TMVAssert(i==j);
1250             return operator()(i);
1251         }
1252 
operator()1253         inline T operator()(ptrdiff_t i) const
1254         {
1255             if (A==CStyle) {
1256                 TMVAssert(i>=0 && i<size());
1257                 return itsdiag(i);
1258             } else {
1259                 TMVAssert(i>0 && i<=size());
1260                 return itsdiag(i-1);
1261             }
1262         }
1263 
operator()1264         inline T operator()(ptrdiff_t i,ptrdiff_t j) const
1265         {
1266             if (A==CStyle) {
1267                 TMVAssert(i>=0 && i<size());
1268                 TMVAssert(j>=0 && j<size());
1269             } else {
1270                 TMVAssert(i>0 && i<=size());
1271                 TMVAssert(j>0 && j<=size());
1272             }
1273             if (i==j) return operator()(i);
1274             else return T(0);
1275         }
1276 
diag()1277         inline vec_type diag() { return itsdiag.view(); }
1278 
diag()1279         inline const_vec_type diag() const { return itsdiag.view(); }
1280 
1281         //
1282         // Modifying Functions
1283         //
1284 
setZero()1285         inline type& setZero()
1286         { itsdiag.setZero(); return *this; }
1287 
setAllTo(const T & x)1288         inline type& setAllTo(const T& x)
1289         { itsdiag.setAllTo(x); return *this; }
1290 
addToAll(const T & x)1291         inline type& addToAll(const T& x)
1292         { itsdiag.addToAll(x); return *this; }
1293 
clip(RT thresh)1294         inline type& clip(RT thresh)
1295         { diag().clip(thresh); return *this; }
1296 
transposeSelf()1297         inline type& transposeSelf()
1298         { return *this; }
1299 
conjugateSelf()1300         inline type& conjugateSelf()
1301         { itsdiag.conjugateSelf(); return *this; }
1302 
invertSelf()1303         inline type& invertSelf()
1304         { view().invertSelf(); return *this; }
1305 
1306         inline type& setToIdentity(const T& x=T(1))
1307         { itsdiag.setAllTo(x); return *this; }
1308 
1309         //
1310         // subDiagMatrix
1311         //
1312 
cSubDiagMatrix(ptrdiff_t i1,ptrdiff_t i2)1313         inline view_type cSubDiagMatrix(ptrdiff_t i1, ptrdiff_t i2)
1314         { return view_type(itsdiag.cSubVector(i1,i2)); }
1315 
subDiagMatrix(ptrdiff_t i1,ptrdiff_t i2)1316         inline view_type subDiagMatrix(ptrdiff_t i1, ptrdiff_t i2)
1317         {
1318             TMVAssert(diag().hasSubVector(i1,i2,1));
1319             if (A == int(FortranStyle)) { --i1; }
1320             return cSubDiagMatrix(i1,i2);
1321         }
1322 
cSubDiagMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t istep)1323         inline view_type cSubDiagMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep)
1324         { return view_type(itsdiag.cSubVector(i1,i2,istep)); }
1325 
subDiagMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t istep)1326         inline view_type subDiagMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep)
1327         {
1328             TMVAssert(diag().hasSubVector(i1,i2,istep));
1329             if (A == int(FortranStyle)) { --i1; i2 += istep-1; }
1330             return cSubDiagMatrix(i1,i2,istep);
1331         }
1332 
realPart()1333         inline realpart_type realPart()
1334         { return realpart_type(diag().realPart()); }
1335 
imagPart()1336         inline realpart_type imagPart()
1337         { return realpart_type(diag().imagPart()); }
1338 
cSubDiagMatrix(ptrdiff_t i1,ptrdiff_t i2)1339         inline const_view_type cSubDiagMatrix(ptrdiff_t i1, ptrdiff_t i2) const
1340         { return const_view_type(itsdiag.cSubVector(i1,i2)); }
1341 
subDiagMatrix(ptrdiff_t i1,ptrdiff_t i2)1342         inline const_view_type subDiagMatrix(ptrdiff_t i1, ptrdiff_t i2) const
1343         {
1344             TMVAssert(diag().hasSubVector(i1,i2,1));
1345             if (A == int(FortranStyle)) { --i1; }
1346             return cSubDiagMatrix(i1,i2);
1347         }
1348 
cSubDiagMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t istep)1349         inline const_view_type cSubDiagMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const
1350         { return const_view_type(itsdiag.cSubVector(i1,i2,istep)); }
1351 
subDiagMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t istep)1352         inline const_view_type subDiagMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const
1353         {
1354             TMVAssert(diag().hasSubVector(i1,i2,istep));
1355             if (A == int(FortranStyle)) { --i1; i2 += istep-1; }
1356             return cSubDiagMatrix(i1,i2,istep);
1357         }
1358 
realPart()1359         inline const_realpart_type realPart() const
1360         { return const_realpart_type(diag().realPart()); }
1361 
imagPart()1362         inline const_realpart_type imagPart() const
1363         { return const_realpart_type(diag().imagPart()); }
1364 
view()1365         inline const_view_type view() const
1366         { return const_view_type(itsdiag); }
1367 
transpose()1368         inline const_view_type transpose() const
1369         { return view(); }
1370 
conjugate()1371         inline const_view_type conjugate() const
1372         { return const_view_type(itsdiag.conjugate()); }
1373 
adjoint()1374         inline const_view_type adjoint() const
1375         { return conjugate(); }
1376 
view()1377         inline view_type view()
1378         { return view_type(itsdiag.view()); }
1379 
transpose()1380         inline view_type transpose()
1381         { return view(); }
1382 
conjugate()1383         inline view_type conjugate()
1384         { return view_type(itsdiag.conjugate()); }
1385 
adjoint()1386         inline view_type adjoint()
1387         { return conjugate(); }
1388 
1389         //
1390         // I/O
1391         //
1392 
1393         void read(const TMV_Reader& reader);
1394 
1395         using base::size;
1396 
resize(ptrdiff_t n)1397         inline void resize(ptrdiff_t n)
1398         {
1399             TMVAssert(n >= 0);
1400             itsdiag.resize(n);
1401         }
1402 
begin()1403         inline iterator begin()
1404         { return diag().begin(); }
end()1405         inline iterator end()
1406         { return diag().end(); }
1407 
begin()1408         inline const_iterator begin() const
1409         { return diag().begin(); }
end()1410         inline const_iterator end() const
1411         { return diag().end(); }
1412 
1413     protected :
1414 
1415         Vector<T> itsdiag;
cdiag()1416         inline ConstVectorView<T> cdiag() const { return itsdiag.view(); }
1417 
Swap(DiagMatrix<T,A> & m1,DiagMatrix<T,A> & m2)1418         friend void Swap(DiagMatrix<T,A>& m1, DiagMatrix<T,A>& m2)
1419         {
1420             TMVAssert(m1.size() == m2.size());
1421             Swap(m1.itsdiag,m2.itsdiag);
1422         }
1423 
1424     }; // DiagMatrix
1425 
1426     //---------------------------------------------------------------------------
1427 
1428     //
1429     // Special Creators:
1430     //   DiagMatrixViewOf(v)
1431     //   DiagMatrixViewOf(T* v, n)
1432     //   DiagMatrixViewOf(T* v, n, step)
1433     //
1434 
1435     template <typename T>
DiagMatrixViewOf(const GenVector<T> & v)1436     inline ConstDiagMatrixView<T> DiagMatrixViewOf(const GenVector<T>& v)
1437     { return ConstDiagMatrixView<T>(v); }
1438 
1439     template <typename T, int A>
DiagMatrixViewOf(const ConstVectorView<T,A> & v)1440     inline ConstDiagMatrixView<T,A> DiagMatrixViewOf(
1441         const ConstVectorView<T,A>& v)
1442     { return ConstDiagMatrixView<T,A>(v); }
1443 
1444     template <typename T, int A>
DiagMatrixViewOf(const Vector<T,A> & v)1445     inline ConstDiagMatrixView<T,A> DiagMatrixViewOf(const Vector<T,A>& v)
1446     { return ConstDiagMatrixView<T,A>(v.view()); }
1447 
1448     template <typename T, int A>
DiagMatrixViewOf(VectorView<T,A> v)1449     inline DiagMatrixView<T,A> DiagMatrixViewOf(VectorView<T,A> v)
1450     { return DiagMatrixView<T,A>(v); }
1451 
1452     template <typename T, int A>
DiagMatrixViewOf(Vector<T,A> & v)1453     inline DiagMatrixView<T,A> DiagMatrixViewOf(Vector<T,A>& v)
1454     { return DiagMatrixView<T,A>(v.view()); }
1455 
1456     template <typename T>
DiagMatrixViewOf(const T * m,ptrdiff_t size)1457     inline ConstDiagMatrixView<T> DiagMatrixViewOf(const T* m, ptrdiff_t size)
1458     {
1459         TMVAssert(size >= 0);
1460         return ConstDiagMatrixView<T>(VectorViewOf(m,size));
1461     }
1462 
1463     template <typename T>
DiagMatrixViewOf(T * m,ptrdiff_t size)1464     inline DiagMatrixView<T> DiagMatrixViewOf(T* m, ptrdiff_t size)
1465     {
1466         TMVAssert(size >= 0);
1467         return DiagMatrixView<T>(VectorViewOf(m,size));
1468     }
1469 
1470     template <typename T>
DiagMatrixViewOf(const T * m,ptrdiff_t size,ptrdiff_t step)1471     inline ConstDiagMatrixView<T> DiagMatrixViewOf(
1472         const T* m, ptrdiff_t size, ptrdiff_t step)
1473     {
1474         TMVAssert(size >= 0);
1475         return ConstDiagMatrixView<T>(VectorViewOf(m,size,step));
1476     }
1477 
1478     template <typename T>
DiagMatrixViewOf(T * m,ptrdiff_t size,ptrdiff_t step)1479     inline DiagMatrixView<T> DiagMatrixViewOf(T* m, ptrdiff_t size, ptrdiff_t step)
1480     {
1481         TMVAssert(size >= 0);
1482         return DiagMatrixView<T>(VectorViewOf(m,size,step));
1483     }
1484 
1485 
1486     //
1487     // Swap Matrices
1488     //
1489 
1490     template <typename T>
Swap(DiagMatrixView<T> m1,DiagMatrixView<T> m2)1491     inline void Swap(DiagMatrixView<T> m1, DiagMatrixView<T> m2)
1492     { Swap(m1.diag(),m2.diag()); }
1493 
1494     template <typename T, int A>
Swap(DiagMatrix<T,A> & m1,DiagMatrixView<T> m2)1495     inline void Swap(DiagMatrix<T,A>& m1, DiagMatrixView<T> m2)
1496     { Swap(m1.diag(),m2.diag()); }
1497 
1498     template <typename T, int A>
Swap(DiagMatrixView<T> m1,DiagMatrix<T,A> & m2)1499     inline void Swap(DiagMatrixView<T> m1, DiagMatrix<T,A>& m2)
1500     { Swap(m1.diag(),m2.diag()); }
1501 
1502     //
1503     // Views:
1504     //
1505 
1506     template <typename T>
Transpose(const GenDiagMatrix<T> & m)1507     inline ConstDiagMatrixView<T> Transpose(const GenDiagMatrix<T>& m)
1508     { return m.transpose(); }
1509 
1510     template <typename T, int A>
Transpose(const ConstDiagMatrixView<T,A> & m)1511     inline ConstDiagMatrixView<T,A> Transpose(const ConstDiagMatrixView<T,A>& m)
1512     { return m.transpose(); }
1513 
1514     template <typename T, int A>
Transpose(const DiagMatrix<T,A> & m)1515     inline ConstDiagMatrixView<T,A> Transpose(const DiagMatrix<T,A>& m)
1516     { return m.transpose(); }
1517 
1518     template <typename T, int A>
Transpose(DiagMatrixView<T,A> m)1519     inline DiagMatrixView<T,A> Transpose(DiagMatrixView<T,A> m)
1520     { return m.transpose(); }
1521 
1522     template <typename T, int A>
Transpose(DiagMatrix<T,A> & m)1523     inline DiagMatrixView<T,A> Transpose(DiagMatrix<T,A>& m)
1524     { return m.transpose(); }
1525 
1526     template <typename T>
Conjugate(const GenDiagMatrix<T> & m)1527     inline ConstDiagMatrixView<T> Conjugate(const GenDiagMatrix<T>& m)
1528     { return m.conjugate(); }
1529 
1530     template <typename T, int A>
Conjugate(const ConstDiagMatrixView<T,A> & m)1531     inline ConstDiagMatrixView<T,A> Conjugate(const ConstDiagMatrixView<T,A>& m)
1532     { return m.conjugate(); }
1533 
1534     template <typename T, int A>
Conjugate(const DiagMatrix<T,A> & m)1535     inline ConstDiagMatrixView<T,A> Conjugate(const DiagMatrix<T,A>& m)
1536     { return m.conjugate(); }
1537 
1538     template <typename T, int A>
Conjugate(DiagMatrixView<T,A> m)1539     inline DiagMatrixView<T,A> Conjugate(DiagMatrixView<T,A> m)
1540     { return m.conjugate(); }
1541 
1542     template <typename T, int A>
Conjugate(DiagMatrix<T,A> & m)1543     inline DiagMatrixView<T,A> Conjugate(DiagMatrix<T,A>& m)
1544     { return m.conjugate(); }
1545 
1546     template <typename T>
Adjoint(const GenDiagMatrix<T> & m)1547     inline ConstDiagMatrixView<T> Adjoint(const GenDiagMatrix<T>& m)
1548     { return m.adjoint(); }
1549 
1550     template <typename T, int A>
Adjoint(const ConstDiagMatrixView<T,A> & m)1551     inline ConstDiagMatrixView<T,A> Adjoint(const ConstDiagMatrixView<T,A>& m)
1552     { return m.adjoint(); }
1553 
1554     template <typename T, int A>
Adjoint(const DiagMatrix<T,A> & m)1555     inline ConstDiagMatrixView<T,A> Adjoint(const DiagMatrix<T,A>& m)
1556     { return m.adjoint(); }
1557 
1558     template <typename T, int A>
Adjoint(DiagMatrixView<T,A> m)1559     inline DiagMatrixView<T,A> Adjoint(DiagMatrixView<T,A> m)
1560     { return m.adjoint(); }
1561 
1562     template <typename T, int A>
Adjoint(DiagMatrix<T,A> & m)1563     inline DiagMatrixView<T,A> Adjoint(DiagMatrix<T,A>& m)
1564     { return m.adjoint(); }
1565 
1566     template <typename T>
Inverse(const GenDiagMatrix<T> & m)1567     inline QuotXD<T,T> Inverse(const GenDiagMatrix<T>& m)
1568     { return m.inverse(); }
1569 
1570 
1571 
1572     //
1573     // DiagMatrix ==, != DiagMatrix
1574     //
1575 
1576     template <typename T1, typename T2>
1577     inline bool operator==(
1578         const GenDiagMatrix<T1>& m1, const GenDiagMatrix<T2>& m2)
1579     { return m1.diag() == m2.diag(); }
1580 
1581     template <typename T1, typename T2>
1582     inline bool operator!=(
1583         const GenDiagMatrix<T1>& m1, const GenDiagMatrix<T2>& m2)
1584     { return !(m1 == m2); }
1585 
1586     template <typename T1, typename T2>
1587     inline bool operator==(
1588         const GenDiagMatrix<T1>& m1, const GenMatrix<T2>& m2)
1589     {
1590         return
1591             m1.diag() == m2.diag() &&
1592             m2.upperTri().offDiag().maxAbs2Element() == T2(0) &&
1593             m2.lowerTri().offDiag().maxAbs2Element() == T2(0);
1594     }
1595 
1596     template <typename T1, typename T2>
1597     inline bool operator==(
1598         const GenMatrix<T1>& m1, const GenDiagMatrix<T2>& m2)
1599     { return m2 == m1; }
1600 
1601     template <typename T1, typename T2>
1602     inline bool operator!=(
1603         const GenDiagMatrix<T1>& m1, const GenMatrix<T2>& m2)
1604     { return !(m1 == m2); }
1605 
1606     template <typename T1, typename T2>
1607     inline bool operator!=(
1608         const GenMatrix<T1>& m1, const GenDiagMatrix<T2>& m2)
1609     { return !(m1 == m2); }
1610 
1611     //
1612     // I/O
1613     //
1614 
1615     template <typename T>
1616     inline std::istream& operator>>(
1617         const TMV_Reader& reader, DiagMatrixView<T> m)
1618     { m.read(reader); return reader.getis(); }
1619 
1620     template <typename T, int A>
1621     inline std::istream& operator>>(
1622         const TMV_Reader& reader, DiagMatrix<T,A>& m)
1623     { m.read(reader); return reader.getis(); }
1624 
1625     template <typename T>
1626     inline std::istream& operator>>(
1627         std::istream& is, DiagMatrixView<T> m)
1628     { return is >> IOStyle() >> m; }
1629 
1630     template <typename T, int A>
1631     inline std::istream& operator>>(std::istream& is, DiagMatrix<T,A>& m)
1632     { return is >> IOStyle() >> m; }
1633 
1634 } // namespace tmv
1635 
1636 #endif
1637