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 //#define XDEBUG
24 
25 
26 #include "TMV_Blas.h"
27 #include "tmv/TMV_SymMatrix.h"
28 #include "tmv/TMV_Matrix.h"
29 #include "tmv/TMV_Vector.h"
30 #include "tmv/TMV_SymLDLD.h"
31 #include "tmv/TMV_SymCHD.h"
32 #include "tmv/TMV_SymSVD.h"
33 #include "tmv/TMV_VIt.h"
34 #include "tmv/TMV_SymMatrixArith.h"
35 #include "tmv/TMV_DiagMatrix.h"
36 #include "TMV_IntegerDet.h"
37 #include <iostream>
38 
39 #ifdef XDEBUG
40 using std::cerr;
41 using std::endl;
42 #endif
43 
44 namespace tmv {
45 
46 #define RT TMV_RealType(T)
47 
48     //
49     // Access
50     //
51 
52     template <class T>
cref(ptrdiff_t i,ptrdiff_t j) const53     T GenSymMatrix<T>::cref(ptrdiff_t i, ptrdiff_t j) const
54     {
55         if ((uplo() == Upper && i<=j) || (uplo() == Lower && i>=j)) {
56             const T* mi = cptr() + i*stepi() + j*stepj();
57             return isconj() ? TMV_CONJ(*mi) : *mi;
58         } else {
59             const T* mi = cptr() + j*stepi() + i*stepj();
60             return issym() != isconj() ? *mi : TMV_CONJ(*mi);
61         }
62     }
63 
64     template <class T, int A>
ref(ptrdiff_t i,ptrdiff_t j)65     typename SymMatrixView<T,A>::reference SymMatrixView<T,A>::ref(
66         ptrdiff_t i, ptrdiff_t j)
67     {
68         if ((uplo() == Upper && i<=j) || (uplo() == Lower && i>=j)) {
69             T* mi = ptr() + i*stepi() + j*stepj();
70             return RefHelper<T>::makeRef(mi,ct());
71         } else {
72             T* mi = ptr() + j*stepi() + i*stepj();
73             return RefHelper<T>::makeRef(
74                 mi, this->issym() != this->isconj() ? NonConj : Conj);
75         }
76     }
77 
78     template <class T>
setDiv() const79     void GenSymMatrix<T>::setDiv() const
80     {
81         if (!this->divIsSet()) {
82             DivType dt = this->getDivType();
83             TMVAssert(dt == tmv::LU || dt == tmv::CH || dt == tmv::SV);
84             TMVAssert(isherm() || dt != tmv::CH);
85             switch (dt) {
86               case LU :
87                    this->divider.reset(
88                        new SymLDLDiv<T>(*this,this->divIsInPlace()));
89                    break;
90               case CH :
91                    this->divider.reset(
92                        new HermCHDiv<T>(*this,this->divIsInPlace()));
93                    break;
94               case SV :
95                    if (isherm())
96                        this->divider.reset(
97                            new HermSVDiv<T>(*this,this->divIsInPlace()));
98                    else
99                        this->divider.reset(
100                            new SymSVDiv<T>(*this,this->divIsInPlace()));
101                    break;
102               default :
103                    // The above assert should have already failed
104                    // so go ahead and fall through.
105                    break;
106             }
107         }
108     }
109 
110 #ifdef INST_INT
111     template <>
setDiv() const112     void GenSymMatrix<int>::setDiv() const
113     { TMVAssert(TMV_FALSE); }
114     template <>
setDiv() const115     void GenSymMatrix<std::complex<int> >::setDiv() const
116     { TMVAssert(TMV_FALSE); }
117 #endif
118 
119     template <class T>
QInverse() const120     QuotXS<T,T> GenSymMatrix<T>::QInverse() const
121     { return QuotXS<T,T>(T(1),*this); }
122 
123     template <class T> template <class T1>
doMakeInverse(SymMatrixView<T1> sinv) const124     void GenSymMatrix<T>::doMakeInverse(SymMatrixView<T1> sinv) const
125     {
126         TMVAssert(issym() == sinv.issym());
127         TMVAssert(isherm() == sinv.isherm());
128 
129         this->setDiv();
130         const SymDivider<T>* sdiv = dynamic_cast<const SymDivider<T>*>(
131             this->getDiv());
132         TMVAssert(sdiv);
133         sdiv->makeInverse(sinv);
134     }
135 
136     template <class T>
divIsLUDiv() const137     bool GenSymMatrix<T>::divIsLUDiv() const
138     { return dynamic_cast<const SymLDLDiv<T>*>(this->getDiv()); }
139 
140     template <class T>
divIsCHDiv() const141     bool GenSymMatrix<T>::divIsCHDiv() const
142     { return dynamic_cast<const HermCHDiv<T>*>(this->getDiv()); }
143 
144     template <class T>
divIsHermSVDiv() const145     bool GenSymMatrix<T>::divIsHermSVDiv() const
146     { return dynamic_cast<const HermSVDiv<T>*>(this->getDiv()); }
147 
148     template <class T>
divIsSymSVDiv() const149     bool GenSymMatrix<T>::divIsSymSVDiv() const
150     { return dynamic_cast<const SymSVDiv<T>*>(this->getDiv()); }
151 
152 #ifdef INST_INT
153     template <>
divIsLUDiv() const154     bool GenSymMatrix<int>::divIsLUDiv() const
155     { return false; }
156     template <>
divIsCHDiv() const157     bool GenSymMatrix<int>::divIsCHDiv() const
158     { return false; }
159     template <>
divIsHermSVDiv() const160     bool GenSymMatrix<int>::divIsHermSVDiv() const
161     { return false; }
162     template <>
divIsSymSVDiv() const163     bool GenSymMatrix<int>::divIsSymSVDiv() const
164     { return false; }
165 
166     template <>
divIsLUDiv() const167     bool GenSymMatrix<std::complex<int> >::divIsLUDiv() const
168     { return false; }
169     template <>
divIsCHDiv() const170     bool GenSymMatrix<std::complex<int> >::divIsCHDiv() const
171     { return false; }
172     template <>
divIsHermSVDiv() const173     bool GenSymMatrix<std::complex<int> >::divIsHermSVDiv() const
174     { return false; }
175     template <>
divIsSymSVDiv() const176     bool GenSymMatrix<std::complex<int> >::divIsSymSVDiv() const
177     { return false; }
178 #endif
179 
180     //
181     // OK? (SubMatrix, etc.)
182     //
183 
184     template <class T>
hasSubMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t j1,ptrdiff_t j2,ptrdiff_t istep,ptrdiff_t jstep) const185     bool GenSymMatrix<T>::hasSubMatrix(
186         ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const
187     {
188         if (i1==i2 || j1==j2) return true; // no elements, so whatever...
189         bool ok = true;
190         ptrdiff_t i2x = i2-istep;
191         ptrdiff_t j2x = j2-jstep;
192         if (istep == 0) {
193             ok = false;
194             std::cerr<<"istep ("<<istep<<") can not be 0\n";
195         }
196         if (i1 < 0 || i1 >= size()) {
197             ok = false;
198             std::cerr<<"first col element ("<<i1<<") must be in 0 -- ";
199             std::cerr<<size()-1<<std::endl;
200         }
201         if (i2x < 0 || i2x >= size()) {
202             ok = false;
203             std::cerr<<"last col element ("<<i2x<<") must be in 0 -- ";
204             std::cerr<<size()-1<<std::endl;
205         }
206         if ((i2-i1)%istep != 0) {
207             ok = false;
208             std::cerr<<"col range ("<<i2-i1<<") must be multiple of istep (";
209             std::cerr<<istep<<")\n";
210         }
211         if ((i2-i1)/istep < 0) {
212             ok = false;
213             std::cerr<<"n col elements ("<<(i2-i1)/istep<<") must be nonnegative\n";
214         }
215         if (jstep == 0) {
216             ok = false;
217             std::cerr<<"jstep ("<<jstep<<") can not be 0\n";
218         }
219         if (j1 < 0 || j1 >= size()) {
220             ok = false;
221             std::cerr<<"first row element ("<<j1<<") must be in 0 -- ";
222             std::cerr<<size()-1<<std::endl;
223         }
224         if (j2-jstep < 0 || j2-jstep >= size()) {
225             ok = false;
226             std::cerr<<"last row element ("<<j2-jstep<<") must be in 0 -- ";
227             std::cerr<<size()-1<<std::endl;
228         }
229         if ((j2-j1)%jstep != 0) {
230             ok = false;
231             std::cerr<<"row range ("<<j2-j1<<") must be multiple of istep (";
232             std::cerr<<jstep<<")\n";
233         }
234         if ((j2-j1)/jstep < 0) {
235             ok = false;
236             std::cerr<<"n row elements ("<<(j2-j1)/jstep<<") must be nonnegative\n";
237         }
238         if ((i1<j1 && i2x>j2x) || (i1>j1 && i2x<j2x)) {
239             ok = false;
240             std::cerr<<"Upper left ("<<i1<<','<<j1<<") and lower right (";
241             std::cerr<<i2x<<','<<j2x<<") corners must be in same triangle\n";
242         }
243         if ((i2x<j1 && i1>j2x) || (i2x>j1 && i1<j2x)) {
244             ok = false;
245             std::cerr<<"Upper right ("<<i1<<','<<j2x<<") and lower left (";
246             std::cerr<<i2x<<','<<j1<<") corners must be in same triangle\n";
247         }
248         return ok;
249     }
250 
251     template <class T>
hasSubVector(ptrdiff_t i,ptrdiff_t j,ptrdiff_t istep,ptrdiff_t jstep,ptrdiff_t n) const252     bool GenSymMatrix<T>::hasSubVector(
253         ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t n) const
254     {
255         if (n==0) return true;
256         bool ok = true;
257         if (istep == 0 && jstep == 0) {
258             ok = false;
259             std::cerr<<"istep ("<<istep<<") and jstep ("<<jstep;
260             std::cerr<<") can not both be 0\n";
261         }
262         if (i<0 || i >= size()) {
263             ok = false;
264             std::cerr<<"i ("<<i<<") must be in 0 -- "<<size()-1<<std::endl;
265         }
266         if (j<0 || j >= size()) {
267             ok = false;
268             std::cerr<<"j ("<<j<<") must be in 0 -- "<<size()-1<<std::endl;
269         }
270         ptrdiff_t i2 = i+istep*(n-1);
271         ptrdiff_t j2 = j+jstep*(n-1);
272         if (i2 < 0 || i2 >= size()) {
273             ok = false;
274             std::cerr<<"last element's i ("<<i2<<") must be in 0 -- ";
275             std::cerr<<size()-1<<std::endl;
276         }
277         if (j2 < 0 || j2 >= size()) {
278             ok = false;
279             std::cerr<<"last element's j ("<<j2<<") must be in 0 -- ";
280             std::cerr<<size()-1<<std::endl;
281         }
282         if ((i<j && i2>j2) || (i>j && i2<j2)) {
283             ok = false;
284             std::cerr<<"First ("<<i<<','<<j<<") and last ("<<i2<<','<<j2;
285             std::cerr<<") elements must be in same triangle\n";
286         }
287         return ok;
288     }
289 
290     template <class T>
hasSubSymMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t istep) const291     bool GenSymMatrix<T>::hasSubSymMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const
292     {
293         if (i1==i2) return true;
294         bool ok=true;
295         if (istep == 0) {
296             ok = false;
297             std::cerr<<"istep ("<<istep<<") can not be 0\n";
298         }
299         if (i1<0 || i1 >= size()) {
300             ok = false;
301             std::cerr<<"first diag element ("<<i1<<") must be in 0 -- ";
302             std::cerr<<size()-1<<std::endl;
303         }
304         if (i2-istep<0 || i2-istep >= size()) {
305             ok = false;
306             std::cerr<<"last diag element ("<<i2-istep<<") must be in 0 -- ";
307             std::cerr<<size()-1<<std::endl;
308         }
309         if ((i2-i1)%istep != 0) {
310             ok = false;
311             std::cerr<<"range ("<<i2-i1<<") must be multiple of istep (";
312             std::cerr<<istep<<")\n";
313         }
314         if ((i2-i1)/istep < 0) {
315             ok = false;
316             std::cerr<<"n diag elements ("<<(i2-i1)/istep<<") must be nonnegative\n";
317         }
318         return ok;
319     }
320 
321     template <class T>
hasSubMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t j1,ptrdiff_t j2,ptrdiff_t istep,ptrdiff_t jstep) const322     bool ConstSymMatrixView<T,FortranStyle>::hasSubMatrix(
323         ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const
324     {
325         if (i1==i2 || j1==j2) return true; // no elements, so whatever...
326         bool ok = true;
327         if (istep == 0) {
328             ok = false;
329             std::cerr<<"istep ("<<istep<<") can not be 0\n";
330         }
331         if (i1 < 1 || i1 > this->size()) {
332             ok = false;
333             std::cerr<<"first col element ("<<i1<<") must be in 1 -- ";
334             std::cerr<<this->size()<<std::endl;
335         }
336         if (i2 < 1 || i2 > this->size()) {
337             ok = false;
338             std::cerr<<"last col element ("<<i2-istep<<") must be in 1 -- ";
339             std::cerr<<this->size()<<std::endl;
340         }
341         if ((i2-i1)%istep != 0) {
342             ok = false;
343             std::cerr<<"col range ("<<i2-i1<<") must be multiple of istep (";
344             std::cerr<<istep<<")\n";
345         }
346         if ((i2-i1)/istep < 0) {
347             ok = false;
348             std::cerr<<"n col elements ("<<(i2-i1)/istep+1<<") must be positive\n";
349         }
350         if (jstep == 0) {
351             ok = false;
352             std::cerr<<"jstep ("<<jstep<<") can not be 0\n";
353         }
354         if (j1 < 1 || j1 > this->size()) {
355             ok = false;
356             std::cerr<<"first row element ("<<j1<<") must be in 1 -- ";
357             std::cerr<<this->size()<<std::endl;
358         }
359         if (j2 < 1 || j2 > this->size()) {
360             ok = false;
361             std::cerr<<"last row element ("<<j2-jstep<<") must be in 1 -- ";
362             std::cerr<<this->size()<<std::endl;
363         }
364         if ((j2-j1)%jstep != 0) {
365             ok = false;
366             std::cerr<<"row range ("<<j2-j1<<") must be multiple of istep (";
367             std::cerr<<jstep<<")\n";
368         }
369         if ((j2-j1)/jstep < 0) {
370             ok = false;
371             std::cerr<<"n row elements ("<<(j2-j1)/jstep+1<<") must be positive\n";
372         }
373         if ((i1<j1 && i2>j2) || (i1>j1 && i2<j2)) {
374             ok = false;
375             std::cerr<<"Upper left ("<<i1<<','<<j1<<") and lower right (";
376             std::cerr<<i2<<','<<j2<<") corners must be in same triangle\n";
377         }
378         if ((i2<j1 && i1>j2) || (i2>j1 && i1<j2)) {
379             ok = false;
380             std::cerr<<"Upper right ("<<i1<<','<<j2<<") and lower left (";
381             std::cerr<<i2<<','<<j1<<") corners must be in same triangle\n";
382         }
383         return ok;
384     }
385 
386     template <class T>
hasSubVector(ptrdiff_t i,ptrdiff_t j,ptrdiff_t istep,ptrdiff_t jstep,ptrdiff_t n) const387     bool ConstSymMatrixView<T,FortranStyle>::hasSubVector(
388         ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t n) const
389     {
390         if (n==0) return true;
391         bool ok = true;
392         if (istep == 0 && jstep == 0) {
393             ok = false;
394             std::cerr<<"istep ("<<istep<<") and jstep ("<<jstep;
395             std::cerr<<") can not both be 0\n";
396         }
397         if (i < 1 || i > this->size()) {
398             ok = false;
399             std::cerr<<"i ("<<i<<") must be in 1 -- "<<this->size()<<std::endl;
400         }
401         if (i < 1 || j > this->size()) {
402             ok = false;
403             std::cerr<<"j ("<<j<<") must be in 1 -- "<<this->size()<<std::endl;
404         }
405         ptrdiff_t i2 = i+istep*(n-1);
406         ptrdiff_t j2 = j+jstep*(n-1);
407         if (i2 < 1 || i2 > this->size()) {
408             ok = false;
409             std::cerr<<"last element's i ("<<i2<<") must be in 1 -- ";
410             std::cerr<<this->size()<<std::endl;
411         }
412         if (j2 < 1 || j2 > this->size()) {
413             ok = false;
414             std::cerr<<"last element's j ("<<j2<<") must be in 1 -- ";
415             std::cerr<<this->size()<<std::endl;
416         }
417         if ((i<j && i2>j2) || (i>j && i2<j2)) {
418             ok = false;
419             std::cerr<<"First ("<<i<<','<<j<<") and last ("<<i2<<','<<j2;
420             std::cerr<<") elements must be in same triangle\n";
421         }
422         return ok;
423     }
424 
425     template <class T>
hasSubSymMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t istep) const426     bool ConstSymMatrixView<T,FortranStyle>::hasSubSymMatrix(
427         ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const
428     {
429         if (i1==i2) return true;
430         bool ok=true;
431         if (istep == 0) {
432             ok = false;
433             std::cerr<<"istep ("<<istep<<") can not be 0\n";
434         }
435         if (i1<1 || i1 > this->size()) {
436             ok = false;
437             std::cerr<<"first diag element ("<<i1<<") must be in 1 -- ";
438             std::cerr<<this->size()<<std::endl;
439         }
440         if (i2-istep<1 || i2-istep > this->size()) {
441             ok = false;
442             std::cerr<<"last diag element ("<<i2-istep<<") must be in 1 -- ";
443             std::cerr<<this->size()<<std::endl;
444         }
445         if ((i2-i1)%istep != 0) {
446             ok = false;
447             std::cerr<<"range ("<<i2-i1<<") must be multiple of istep ("<<istep;
448             std::cerr<<")\n";
449         }
450         if ((i2-i1)/istep < 0) {
451             ok = false;
452             std::cerr<<"n diag elements ("<<(i2-i1)/istep+1<<") must be positive\n";
453         }
454         return ok;
455     }
456 
457     //
458     // SwapRowsCols
459     //
460 
461     template <class T, int A>
swapRowsCols(ptrdiff_t i1,ptrdiff_t i2)462     SymMatrixView<T,A>& SymMatrixView<T,A>::swapRowsCols(ptrdiff_t i1, ptrdiff_t i2)
463     {
464         TMVAssert(i1<size());
465         TMVAssert(i2<size());
466         if (i1 == i2) return *this;
467         else {
468 #ifdef XDEBUG
469             Matrix<T> M(*this);
470             Matrix<T> M0(*this);
471             M.swapRows(i1,i2);
472             M.swapCols(i1,i2);
473 #endif
474             if (i1 > i2) TMV_SWAP(i1,i2);
475             // Now i1 < i2
476             if (uplo() == Upper) transpose().swapRowsCols(i1,i2);
477             else {
478                 Swap(row(i1,0,i1),row(i2,0,i1));
479                 Swap(row(i2,i1+1,i2),col(i1,i1+1,i2));
480                 if (!this->issym()) {
481                     row(i2,i1,i2).conjugateSelf(); // Conj m(i2,i1) as well
482                     col(i1,i1+1,i2).conjugateSelf();
483                 }
484                 Swap(col(i1,i2+1,size()),col(i2,i2+1,size()));
485                 diag().swap(i1,i2);
486             }
487 #ifdef XDEBUG
488             if (Norm(M-*this) > 1.e-5*TMV_MAX(RT(1),Norm(M))) {
489                 cerr<<"swapRowsCols: i1,i2 = "<<i1<<','<<i2<<endl;
490                 cerr<<"M0 = "<<TMV_Text(*this)<<"  "<<M0<<endl;
491                 cerr<<"M = "<<M<<endl;
492                 cerr<<"S = "<<*this<<endl;
493                 abort();
494             }
495 #endif
496             return *this;
497         }
498     }
499 
500     template <class T, int A>
permuteRowsCols(const ptrdiff_t * p,ptrdiff_t i1,ptrdiff_t i2)501     SymMatrixView<T,A>& SymMatrixView<T,A>::permuteRowsCols(
502         const ptrdiff_t* p, ptrdiff_t i1, ptrdiff_t i2)
503     {
504         const ptrdiff_t* pi = p+i1;
505         for(ptrdiff_t i=i1;i<i2;++i,++pi) {
506             TMVAssert(*pi < size());
507             swapRowsCols(i,*pi);
508         }
509         return *this;
510     }
511 
512     template <class T, int A>
reversePermuteRowsCols(const ptrdiff_t * p,ptrdiff_t i1,ptrdiff_t i2)513     SymMatrixView<T,A>& SymMatrixView<T,A>::reversePermuteRowsCols(
514         const ptrdiff_t* p, ptrdiff_t i1, ptrdiff_t i2)
515     {
516         const ptrdiff_t* pi = p+i2;
517         for(ptrdiff_t i=i2;i>i1;) {
518             --i; --pi;
519             TMVAssert(*pi < size());
520             swapRowsCols(i,*pi);
521         }
522         return *this;
523     }
524 
525 
526     //
527     // Norms
528     //
529 
530     template <class T>
det() const531     T GenSymMatrix<T>::det() const
532     { return DivHelper<T>::det(); }
533 
534     template <class T>
logDet(T * sign) const535     RT GenSymMatrix<T>::logDet(T* sign) const
536     { return DivHelper<T>::logDet(sign); }
537 
538     template <class T>
isSingular() const539     bool GenSymMatrix<T>::isSingular() const
540     { return DivHelper<T>::isSingular(); }
541 
542 #ifdef INST_INT
543     template <>
det() const544     int GenSymMatrix<int>::det() const
545     { return IntegerDet(*this); }
546 
547     template <>
det() const548     std::complex<int> GenSymMatrix<std::complex<int> >::det() const
549     { return IntegerDet(*this); }
550 
551     template <>
logDet(int *) const552     int GenSymMatrix<int>::logDet(int* ) const
553     { TMVAssert(TMV_FALSE); return 0; }
554 
555     template <>
logDet(std::complex<int> *) const556     int GenSymMatrix<std::complex<int> >::logDet(std::complex<int>* ) const
557     { TMVAssert(TMV_FALSE); return 0; }
558 
559     template <>
isSingular() const560     bool GenSymMatrix<int>::isSingular() const
561     { return det() == 0; }
562 
563     template <>
isSingular() const564     bool GenSymMatrix<std::complex<int> >::isSingular() const
565     { return det() == 0; }
566 #endif
567 
568     template <class T>
sumElements() const569     T GenSymMatrix<T>::sumElements() const
570     {
571         T sum = diag().sumElements();
572         if (size() > 1) {
573             T temp = upperTri().offDiag().sumElements();
574             if (issym()) {
575                 sum += RT(2) * temp;
576             } else {
577                 // temp + conj(temp) = 2*real(temp)
578                 sum += RT(2) * TMV_REAL(temp);
579             }
580         }
581         return sum;
582     }
583 
584     template <class T>
DoSumAbsElements(const GenSymMatrix<T> & m)585     static RT DoSumAbsElements(const GenSymMatrix<T>& m)
586     {
587         RT sum = m.diag().sumAbsElements();
588         if (m.size() > 1)
589             sum += RT(2) * m.upperTri().offDiag().sumAbsElements();
590         return sum;
591     }
592 
593     template <class T>
DoSumAbs2Elements(const GenSymMatrix<T> & m)594     static RT DoSumAbs2Elements(const GenSymMatrix<T>& m)
595     {
596         RT sum = m.diag().sumAbs2Elements();
597         if (m.size() > 1)
598             sum += RT(2) * m.upperTri().offDiag().sumAbs2Elements();
599         return sum;
600     }
601 
602 #ifdef INST_INT
DoSumAbsElements(const GenSymMatrix<std::complex<int>> &)603     static int DoSumAbsElements(const GenSymMatrix<std::complex<int> >& )
604     { TMVAssert(TMV_FALSE); return 0; }
605 #endif
606 
607     template <class T>
sumAbsElements() const608     RT GenSymMatrix<T>::sumAbsElements() const
609     { return DoSumAbsElements(*this); }
610 
611     template <class T>
sumAbs2Elements() const612     RT GenSymMatrix<T>::sumAbs2Elements() const
613     { return DoSumAbs2Elements(*this); }
614 
615     template <class T>
normSq(const RT scale) const616     RT GenSymMatrix<T>::normSq(const RT scale) const
617     {
618         RT sum = diag().normSq(scale);
619         if (size() > 1) sum += RT(2) * upperTri().offDiag().normSq(scale);
620         return sum;
621     }
622 
623     template <class T>
NonLapNorm1(const GenSymMatrix<T> & m)624     static RT NonLapNorm1(const GenSymMatrix<T>& m)
625     {
626         RT max(0);
627         for(ptrdiff_t j=0;j<m.size();++j) {
628             RT temp = m.col(j,0,j).norm1();
629             temp += m.col(j,j,m.size()).norm1();
630             if (temp > max) max = temp;
631         }
632         return max;
633     }
634 
635     template <class T>
NonLapNormF(const GenSymMatrix<T> & m)636     static RT NonLapNormF(const GenSymMatrix<T>& m)
637     {
638         const RT eps = TMV_Epsilon<T>();
639 
640         RT mmax = m.maxAbs2Element();
641         if (mmax == RT(0)) return RT(0);
642         else if (TMV_Underflow(mmax * mmax)) {
643             // Then we need to rescale, since underflow has caused
644             // rounding errors.
645             // Epsilon is a pure power of 2, so no rounding errors from
646             // rescaling.
647             const RT inveps = RT(1)/eps;
648             RT scale = inveps;
649             mmax *= scale;
650             const RT eps2 = eps*eps;
651             while (mmax < eps2) { scale *= inveps; mmax *= inveps; }
652             return TMV_SQRT(m.normSq(scale))/scale;
653         } else if (RT(1) / mmax == RT(0)) {
654             // Then mmax is already inf, so no hope of making it more accurate.
655             return mmax;
656         } else if (RT(1) / (mmax*mmax) == RT(0)) {
657             // Then we have overflow, so we need to rescale:
658             const RT inveps = RT(1)/eps;
659             RT scale = eps;
660             mmax *= scale;
661             while (mmax > inveps) { scale *= eps; mmax *= eps; }
662             return TMV_SQRT(m.normSq(scale))/scale;
663         } else {
664             return TMV_SQRT(m.normSq());
665         }
666     }
667 
668     template <class T>
NonLapMaxAbsElement(const GenSymMatrix<T> & m)669     static inline RT NonLapMaxAbsElement(const GenSymMatrix<T>& m)
670     { return m.upperTri().maxAbsElement(); }
671 
672     template <class T>
NonLapMaxAbs2Element(const GenSymMatrix<T> & m)673     static inline RT NonLapMaxAbs2Element(const GenSymMatrix<T>& m)
674     { return m.upperTri().maxAbs2Element(); }
675 
676 #ifdef XLAP
677     template <class T>
LapNorm(const char c,const GenSymMatrix<T> & m)678     static RT LapNorm(const char c, const GenSymMatrix<T>& m)
679     {
680         switch(c) {
681           case 'M' : return NonLapMaxAbsElement(m);
682           case '1' : return NonLapNorm1(m);
683           case 'F' : return NonLapNormF(m);
684           default : TMVAssert(TMV_FALSE);
685         }
686         return RT(0);
687     }
688 #ifdef INST_DOUBLE
689     template <>
LapNorm(const char c,const GenSymMatrix<double> & m)690     double LapNorm(const char c, const GenSymMatrix<double>& m)
691     {
692         TMVAssert(m.iscm() || m.isrm());
693         char cc = c;
694         int N = m.size();
695         int lda = m.iscm() ? m.stepj() : m.stepi();
696 #ifndef LAPNOWORK
697         int lwork = c=='1' ? N : 0;
698         AlignedArray<double> work(lwork);
699         VectorViewOf(work.get(),lwork).setZero();
700 #endif
701         double norm = LAPNAME(dlansy) (
702             LAPCM LAPV(cc),
703             (m.iscm() == (m.uplo()==Upper) ? LAPCH_UP : LAPCH_LO),
704             LAPV(N),LAPP(m.cptr()),LAPV(lda) LAPWK(work.get()) LAP1 LAP1);
705         return norm;
706     }
707     template <>
LapNorm(const char c,const GenSymMatrix<std::complex<double>> & m)708     double LapNorm(const char c, const GenSymMatrix<std::complex<double> >& m)
709     {
710         TMVAssert(m.iscm() || m.isrm());
711         char cc = c;
712         int N = m.size();
713         int lda = m.iscm() ? m.stepj() : m.stepi();
714 #ifndef LAPNOWORK
715         int lwork = c=='1' ? N : 0;
716         AlignedArray<double> work(lwork);
717         VectorViewOf(work.get(),lwork).setZero();
718 #endif
719         double norm = m.isherm() ?
720             LAPNAME(zlanhe) (
721                 LAPCM LAPV(cc),
722                 (m.iscm() == (m.uplo()==Upper) ? LAPCH_UP : LAPCH_LO),
723                 LAPV(N),LAPP(m.cptr()),LAPV(lda) LAPWK(work.get()) LAP1 LAP1) :
724             LAPNAME(zlansy) (
725                 LAPCM LAPV(cc),
726                 (m.iscm() == (m.uplo()==Upper) ? LAPCH_UP : LAPCH_LO),
727                 LAPV(N),LAPP(m.cptr()),LAPV(lda) LAPWK(work.get()) LAP1 LAP1);
728         return norm;
729     }
730 #endif
731 #ifdef INST_FLOAT
732     template <>
LapNorm(const char c,const GenSymMatrix<float> & m)733     float LapNorm(const char c, const GenSymMatrix<float>& m)
734     {
735         TMVAssert(m.iscm() || m.isrm());
736         char cc = c;
737         int N = m.size();
738         int lda = m.iscm() ? m.stepj() : m.stepi();
739 #ifndef LAPNOWORK
740         int lwork = c=='1' ? N : 0;
741         AlignedArray<float> work(lwork);
742         VectorViewOf(work.get(),lwork).setZero();
743 #endif
744         float norm = LAPNAME(slansy) (
745             LAPCM LAPV(cc),
746             (m.iscm() == (m.uplo()==Upper) ? LAPCH_UP : LAPCH_LO),
747             LAPV(N),LAPP(m.cptr()),LAPV(lda) LAPWK(work.get()) LAP1 LAP1);
748         return norm;
749     }
750     template <>
LapNorm(const char c,const GenSymMatrix<std::complex<float>> & m)751     float LapNorm(const char c, const GenSymMatrix<std::complex<float> >& m)
752     {
753         TMVAssert(m.iscm() || m.isrm());
754         char cc = c;
755         int N = m.size();
756         int lda = m.iscm() ? m.stepj() : m.stepi();
757 #ifndef LAPNOWORK
758         int lwork = c=='1' ? N : 0;
759         AlignedArray<float> work(lwork);
760         VectorViewOf(work.get(),lwork).setZero();
761 #endif
762         float norm = m.isherm() ?
763             LAPNAME(clanhe) (
764                 LAPCM LAPV(cc),
765                 (m.iscm() == (m.uplo()==Upper) ? LAPCH_UP : LAPCH_LO),
766                 LAPV(N),LAPP(m.cptr()),LAPV(lda) LAPWK(work.get()) LAP1 LAP1) :
767             LAPNAME(clansy) (
768                 LAPCM LAPV(cc),
769                 (m.iscm() == (m.uplo()==Upper) ? LAPCH_UP : LAPCH_LO),
770                 LAPV(N),LAPP(m.cptr()),LAPV(lda) LAPWK(work.get()) LAP1 LAP1);
771         return norm;
772     }
773 #endif
774 #endif // XLAP
775 
776 #ifdef INST_INT
NonLapNormF(const GenSymMatrix<int> &)777     static int NonLapNormF(const GenSymMatrix<int>& )
778     { TMVAssert(TMV_FALSE); return 0; }
NonLapNormF(const GenSymMatrix<std::complex<int>> &)779     static int NonLapNormF(const GenSymMatrix<std::complex<int> >& )
780     { TMVAssert(TMV_FALSE); return 0; }
NonLapNorm1(const GenSymMatrix<std::complex<int>> &)781     static int NonLapNorm1(const GenSymMatrix<std::complex<int> >& )
782     { TMVAssert(TMV_FALSE); return 0; }
NonLapMaxAbsElement(const GenSymMatrix<std::complex<int>> &)783     static int NonLapMaxAbsElement(const GenSymMatrix<std::complex<int> >& )
784     { TMVAssert(TMV_FALSE); return 0; }
785 #endif
786 
787     template <class T>
maxAbsElement() const788     RT GenSymMatrix<T>::maxAbsElement() const
789     {
790 #ifdef XLAP
791         if ((isrm() && stepi()>0) || (iscm() && stepj()>0))
792             return LapNorm('M',*this);
793         else
794 #endif
795             return NonLapMaxAbsElement(*this);
796     }
797     template <class T>
maxAbs2Element() const798     RT GenSymMatrix<T>::maxAbs2Element() const
799     {
800 #ifdef XLAP
801         if (Traits<T>::iscomplex) return NonLapMaxAbs2Element(*this);
802         else if ((isrm() && stepi()>0) || (iscm() && stepj()>0))
803             return LapNorm('M',*this);
804         else
805 #endif
806             return NonLapMaxAbs2Element(*this);
807     }
808     template <class T>
norm1() const809     RT GenSymMatrix<T>::norm1() const
810     {
811 #ifdef XLAP
812         if ((isrm() && stepi()>0) || (iscm() && stepj()>0))
813             return LapNorm('1',*this);
814         else
815 #endif
816             return NonLapNorm1(*this);
817     }
818     template <class T>
normF() const819     RT GenSymMatrix<T>::normF() const
820     {
821 #ifdef XLAP
822         if ((isrm() && stepi()>0) || (iscm() && stepj()>0))
823             return LapNorm('F',*this);
824         else
825 #endif
826             return NonLapNormF(*this);
827     }
828 
829     template <class T>
DoNorm2(const GenSymMatrix<T> & m)830     static RT DoNorm2(const GenSymMatrix<T>& m)
831     {
832         if (m.size() == 0) return RT(0);
833         DiagMatrix<RT> S(m.size());
834         if (m.isherm()) {
835             HermMatrix<T> m2(m);
836             SV_Decompose(m2.view(),S.view());
837         } else {
838             SymMatrix<T> m2(m);
839             SV_Decompose(m2.view(),S.view());
840         }
841         return S(0);
842     }
843 
844     template <class T>
DoCondition(const GenSymMatrix<T> & m)845     static RT DoCondition(const GenSymMatrix<T>& m)
846     {
847         if (m.size() == 0) return RT(1);
848         DiagMatrix<RT> S(m.size());
849         if (m.isherm()) {
850             HermMatrix<T> m2(m);
851             SV_Decompose(m2.view(),S.view());
852         } else {
853             SymMatrix<T> m2(m);
854             SV_Decompose(m2.view(),S.view());
855         }
856         return TMV_ABS(S(0)/S(S.size()-1));
857     }
858 
859 #ifdef INST_INT
DoNorm2(const GenSymMatrix<int> &)860     static int DoNorm2(const GenSymMatrix<int>& )
861     { TMVAssert(TMV_FALSE); return 0; }
DoCondition(const GenSymMatrix<int> &)862     static int DoCondition(const GenSymMatrix<int>& )
863     { TMVAssert(TMV_FALSE); return 0; }
DoNorm2(const GenSymMatrix<std::complex<int>> &)864     static int DoNorm2(const GenSymMatrix<std::complex<int> >& )
865     { TMVAssert(TMV_FALSE); return 0; }
DoCondition(const GenSymMatrix<std::complex<int>> &)866     static int DoCondition(const GenSymMatrix<std::complex<int> >& )
867     { TMVAssert(TMV_FALSE); return 0; }
868 #endif
869 
870     template <class T>
doNorm2() const871     RT GenSymMatrix<T>::doNorm2() const
872     { return tmv::DoNorm2(*this); }
873     template <class T>
doCondition() const874     RT GenSymMatrix<T>::doCondition() const
875     { return tmv::DoCondition(*this); }
876 
877 
878     //
879     // I/O
880     //
881 
882     template <class T>
write(const TMV_Writer & writer) const883     void GenSymMatrix<T>::write(const TMV_Writer& writer) const
884     {
885         const ptrdiff_t N = size();
886         writer.begin();
887         writer.writeCode(issym() ? "S" : "H");
888         writer.writeSize(N);
889         writer.writeSimpleSize(N);
890         writer.writeStart();
891         for(ptrdiff_t i=0;i<N;++i) {
892             writer.writeLParen();
893             for(ptrdiff_t j=0;j<i+1;++j) {
894                 if (j > 0) writer.writeSpace();
895                 writer.writeValue(cref(i,j));
896             }
897             if (!writer.isCompact()) {
898                 for(ptrdiff_t j=i+1;j<N;++j) {
899                     writer.writeSpace();
900                     writer.writeValue(cref(i,j));
901                 }
902             }
903             writer.writeRParen();
904             if (i < N-1) writer.writeRowEnd();
905         }
906         writer.writeFinal();
907         writer.end();
908     }
909 
910 #ifndef NOTHROW
911     template <class T>
912     class SymMatrixReadError : public ReadError
913     {
914     public :
915         SymMatrix<T> m;
916         ptrdiff_t i,j;
917         std::string exp,got;
918         ptrdiff_t s;
919         T v1, v2;
920         bool is, iseof, isbad;
921 
SymMatrixReadError(std::istream & _is)922         SymMatrixReadError(std::istream& _is) throw() :
923             ReadError("SymMatrix."),
924             i(0), j(0), exp(0), got(0), s(0), v1(0), v2(0),
925             is(_is), iseof(_is.eof()), isbad(_is.bad()) {}
SymMatrixReadError(std::istream & _is,const std::string & _e,const std::string & _g)926         SymMatrixReadError(
927             std::istream& _is,
928             const std::string& _e, const std::string& _g) throw() :
929             ReadError("SymMatrix."),
930             i(0), j(0), exp(_e), got(_g), s(0), v1(0), v2(0),
931             is(_is), iseof(_is.eof()), isbad(_is.bad()) {}
932 
SymMatrixReadError(ptrdiff_t _i,ptrdiff_t _j,const GenSymMatrix<T> & _m,std::istream & _is,const std::string & _e,const std::string & _g)933         SymMatrixReadError(
934             ptrdiff_t _i, ptrdiff_t _j, const GenSymMatrix<T>& _m,
935             std::istream& _is,
936             const std::string& _e, const std::string& _g) throw() :
937             ReadError("SymMatrix."),
938             m(_m), i(_i), j(_j), exp(_e), got(_g), s(m.size()), v1(0), v2(0),
939             is(_is), iseof(_is.eof()), isbad(_is.bad()) {}
SymMatrixReadError(ptrdiff_t _i,ptrdiff_t _j,const GenSymMatrix<T> & _m,std::istream & _is,T _v1=0,T _v2=0)940         SymMatrixReadError(
941             ptrdiff_t _i, ptrdiff_t _j, const GenSymMatrix<T>& _m,
942             std::istream& _is, T _v1=0, T _v2=0) throw() :
943             ReadError("SymMatrix."),
944             m(_m), i(_i), j(_j), s(m.size()), v1(_v1), v2(_v2),
945             is(_is), iseof(_is.eof()), isbad(_is.bad()) {}
SymMatrixReadError(const GenSymMatrix<T> & _m,std::istream & _is,ptrdiff_t _s)946         SymMatrixReadError(
947             const GenSymMatrix<T>& _m, std::istream& _is, ptrdiff_t _s) throw() :
948             ReadError("SymMatrix."),
949             m(_m), i(0), j(0), exp(0), got(0), s(_s), v1(0), v2(0),
950             is(_is), iseof(_is.eof()), isbad(_is.bad()) {}
951 
SymMatrixReadError(const SymMatrixReadError<T> & rhs)952         SymMatrixReadError(const SymMatrixReadError<T>& rhs) :
953             m(rhs.m), i(rhs.i), j(rhs.j), exp(rhs.exp), got(rhs.got),
954             s(rhs.s), v1(rhs.v1), v2(rhs.v2),
955             is(rhs.is), iseof(rhs.iseof), isbad(rhs.isbad) {}
~SymMatrixReadError()956         virtual ~SymMatrixReadError() throw() {}
957 
write(std::ostream & os) const958         virtual void write(std::ostream& os) const throw()
959         {
960             os<<"TMV Read Error: Reading istream input for SymMatrix\n";
961             if (exp != got) {
962                 os<<"Wrong format: expected '"<<exp<<"'";
963                 if (isReal(T()) && exp == "S") os<<" (or 'H')";
964                 os<<", got '"<<got<<"'.\n";
965             }
966             if (s != m.size()) {
967                 os<<"Wrong size: expected "<<m.size()<<", got "<<s<<".\n";
968             }
969             if (!is) {
970                 if (iseof) {
971                     os<<"Input stream reached end-of-file prematurely.\n";
972                 } else if (isbad) {
973                     os<<"Input stream is corrupted.\n";
974                 } else {
975                     os<<"Input stream cannot read next character.\n";
976                 }
977             }
978             if (v1 != v2) {
979                 os<<"Input matrix is not symmetric.\n";
980                 os<<"Lower triangle has the value "<<v1<<" at ("<<i<<","<<j<<")\n";
981                 os<<"Upper triangle has the value "<<v2<<" at ("<<j<<","<<i<<")\n";
982             }
983             if (m.size() > 0) {
984                 os<<"The portion of the SymMatrix which was successfully "
985                     "read is: \n";
986                 const ptrdiff_t N = m.size();
987                 for(ptrdiff_t ii=0;ii<i;++ii) {
988                     os<<"( ";
989                     for(ptrdiff_t jj=0;jj<N;++jj) os<<' '<<m.cref(ii,jj)<<' ';
990                     os<<" )\n";
991                 }
992                 os<<"( ";
993                 for(ptrdiff_t jj=0;jj<j;++jj) os<<' '<<m.cref(i,jj)<<' ';
994                 os<<" )\n";
995             }
996         }
997     };
998 
999     template <class T>
1000     class HermMatrixReadError : public ReadError
1001     {
1002     public :
1003         HermMatrix<T> m;
1004         ptrdiff_t i,j;
1005         std::string exp,got;
1006         ptrdiff_t s;
1007         T v1, v2;
1008         bool is, iseof, isbad;
1009 
HermMatrixReadError(std::istream & _is)1010         HermMatrixReadError(std::istream& _is) throw() :
1011             ReadError("HermMatrix."),
1012             i(0), j(0), exp(0), got(0), s(0), v1(0), v2(0),
1013             is(_is), iseof(_is.eof()), isbad(_is.bad()) {}
HermMatrixReadError(std::istream & _is,const std::string & _e,const std::string & _g)1014         HermMatrixReadError(
1015             std::istream& _is,
1016             const std::string& _e, const std::string& _g) throw() :
1017             ReadError("HermMatrix."),
1018             i(0), j(0), exp(_e), got(_g), s(0), v1(0), v2(0),
1019             is(_is), iseof(_is.eof()), isbad(_is.bad()) {}
1020 
HermMatrixReadError(ptrdiff_t _i,ptrdiff_t _j,const GenSymMatrix<T> & _m,std::istream & _is,const std::string & _e,const std::string & _g)1021         HermMatrixReadError(
1022             ptrdiff_t _i, ptrdiff_t _j, const GenSymMatrix<T>& _m,
1023             std::istream& _is,
1024             const std::string& _e, const std::string& _g) throw() :
1025             ReadError("HermMatrix."),
1026             m(_m), i(_i), j(_j), exp(_e), got(_g),
1027             s(m.size()), v1(0), v2(0),
1028             is(_is), iseof(_is.eof()), isbad(_is.bad()) {}
HermMatrixReadError(ptrdiff_t _i,ptrdiff_t _j,const GenSymMatrix<T> & _m,std::istream & _is,T _v1=0,T _v2=0)1029         HermMatrixReadError(
1030             ptrdiff_t _i, ptrdiff_t _j, const GenSymMatrix<T>& _m,
1031             std::istream& _is, T _v1=0, T _v2=0) throw() :
1032             ReadError("HermMatrix."),
1033             m(_m), i(_i), j(_j), exp(0), got(0),
1034             s(m.size()), v1(_v1), v2(_v2),
1035             is(_is), iseof(_is.eof()), isbad(_is.bad()) {}
HermMatrixReadError(const GenSymMatrix<T> & _m,std::istream & _is,ptrdiff_t _s)1036         HermMatrixReadError(
1037             const GenSymMatrix<T>& _m, std::istream& _is, ptrdiff_t _s) throw() :
1038             ReadError("HermMatrix."),
1039             m(_m), i(0), j(0), exp(0), got(0), s(_s), v1(0), v2(0),
1040             is(_is), iseof(_is.eof()), isbad(_is.bad()) {}
1041 
HermMatrixReadError(const HermMatrixReadError<T> & rhs)1042         HermMatrixReadError(const HermMatrixReadError<T>& rhs) :
1043             m(rhs.m), i(rhs.i), j(rhs.j), exp(rhs.exp), got(rhs.got),
1044             s(rhs.s), v1(rhs.v1), v2(rhs.v2),
1045             is(rhs.is), iseof(rhs.iseof), isbad(rhs.isbad) {}
~HermMatrixReadError()1046         virtual ~HermMatrixReadError() throw() {}
1047 
write(std::ostream & os) const1048         virtual void write(std::ostream& os) const throw()
1049         {
1050             os<<"TMV Read Error: Reading istream input for HermMatrix\n";
1051             if (exp != got) {
1052                 os<<"Wrong format: expected '"<<exp<<"'";
1053                 if (isReal(T()) && exp == "H") os<<" (or 'S')";
1054                 os<<", got '"<<got<<"'.\n";
1055             }
1056             if (s != m.size()) {
1057                 os<<"Wrong size: expected "<<m.size()<<", got "<<s<<".\n";
1058             }
1059             if (!is) {
1060                 if (iseof) {
1061                     os<<"Input stream reached end-of-file prematurely.\n";
1062                 } else if (isbad) {
1063                     os<<"Input stream is corrupted.\n";
1064                 } else {
1065                     os<<"Input stream cannot read next character.\n";
1066                 }
1067             }
1068             if (i==j && v1 != T(0)) {
1069                 os<<"Non-real value found on diagonal: "<<v1<<std::endl;
1070             }
1071             if (i!=j && v1 != v2) {
1072                 os<<"Input matrix is not Hermitian.\n";
1073                 os<<"Lower triangle has the value "<<v1<<" at ("<<i<<","<<j<<")\n";
1074                 os<<"Upper triangle has the value "<<v2<<" at ("<<j<<","<<i<<")\n";
1075             }
1076             if (m.size() > 0) {
1077                 os<<"The portion of the HermMatrix which was successfully "
1078                     "read is: \n";
1079                 const ptrdiff_t N = m.size();
1080                 for(ptrdiff_t ii=0;ii<i;++ii) {
1081                     os<<"( ";
1082                     for(ptrdiff_t jj=0;jj<N;++jj) os<<' '<<m.cref(ii,jj)<<' ';
1083                     os<<" )\n";
1084                 }
1085                 os<<"( ";
1086                 for(ptrdiff_t jj=0;jj<j;++jj) os<<' '<<m.cref(i,jj)<<' ';
1087                 os<<" )\n";
1088             }
1089         }
1090     };
1091 #endif
1092 
1093     template <class T>
FinishRead(const TMV_Reader & reader,SymMatrixView<T> m)1094     static void FinishRead(const TMV_Reader& reader, SymMatrixView<T> m)
1095     {
1096         const ptrdiff_t N = m.size();
1097         std::string exp, got;
1098         T temp;
1099         if (!reader.readStart(exp,got)) {
1100 #ifdef NOTHROW
1101             std::cerr<<"SymMatrix Read Error: "<<got<<" != "<<exp<<std::endl;
1102             exit(1);
1103 #else
1104             if (m.issym())
1105                 throw SymMatrixReadError<T>(0,0,m,reader.getis(),exp,got);
1106             else
1107                 throw HermMatrixReadError<T>(0,0,m,reader.getis(),exp,got);
1108 #endif
1109         }
1110         for(ptrdiff_t i=0;i<N;++i) {
1111             if (!reader.readLParen(exp,got)) {
1112 #ifdef NOTHROW
1113                 std::cerr<<"SymMatrix Read Error: "<<got<<" != "<<exp<<std::endl;
1114                 exit(1);
1115 #else
1116                 if (m.issym())
1117                     throw SymMatrixReadError<T>(i,0,m,reader.getis(),exp,got);
1118                 else
1119                     throw HermMatrixReadError<T>(i,0,m,reader.getis(),exp,got);
1120 #endif
1121             }
1122             for(ptrdiff_t j=0;j<i+1;++j) {
1123                 if (j>0 && !reader.readSpace(exp,got)) {
1124 #ifdef NOTHROW
1125                     std::cerr<<"SymMatrix Read Error: "<<got<<" != "<<exp<<std::endl;
1126                     exit(1);
1127 #else
1128                     if (m.issym())
1129                         throw SymMatrixReadError<T>(i,j,m,reader.getis(),exp,got);
1130                     else
1131                         throw HermMatrixReadError<T>(i,j,m,reader.getis(),exp,got);
1132 #endif
1133                 }
1134                 if (!reader.readValue(temp)) {
1135 #ifdef NOTHROW
1136                     std::cerr<<"SymMatrix Read Error: reading value\n";
1137                     exit(1);
1138 #else
1139                     if (m.issym())
1140                         throw SymMatrixReadError<T>(i,j,m,reader.getis());
1141                     else
1142                         throw HermMatrixReadError<T>(i,j,m,reader.getis());
1143 #endif
1144                 }
1145                 if (j == i && !m.issym() && TMV_IMAG(temp) != RT(0)) {
1146 #ifdef NOTHROW
1147                     std::cerr<<"SymMatrix Read Error: "<<temp<<" not real\n";
1148                     exit(1);
1149 #else
1150                     throw HermMatrixReadError<T>(i,j,m,reader.getis(),temp);
1151 #endif
1152                 }
1153                 if (j < i && !reader.isCompact()) {
1154                     T temp2 = m.issym() ? m.cref(j,i) : TMV_CONJ(m.cref(j,i));
1155                     if (temp != temp2) {
1156 #ifdef NOTHROW
1157                         std::cerr<<"SymMatrix Read Error: "<<temp<<" != "<<temp2<<"\n";
1158                         exit(1);
1159 #else
1160                         if (m.issym())
1161                             throw SymMatrixReadError<T>(i,j,m,reader.getis(),temp,m.cref(j,i));
1162                         else
1163                             throw HermMatrixReadError<T>(i,j,m,reader.getis(),temp,m.cref(j,i));
1164 #endif
1165                     }
1166                 } else {
1167                     m.ref(i,j) = temp;
1168                 }
1169             }
1170             if (!reader.isCompact()) {
1171                 for(ptrdiff_t j=i+1;j<N;++j) {
1172                     if (!reader.readSpace(exp,got)) {
1173 #ifdef NOTHROW
1174                         std::cerr<<"SymMatrix Read Error: "<<got<<" != "<<exp<<std::endl;
1175                         exit(1);
1176 #else
1177                         if (m.issym())
1178                             throw SymMatrixReadError<T>(i,j,m,reader.getis(),exp,got);
1179                         else
1180                             throw HermMatrixReadError<T>(i,j,m,reader.getis(),exp,got);
1181 #endif
1182                     }
1183                     if (!reader.readValue(temp)) {
1184 #ifdef NOTHROW
1185                         std::cerr<<"SymMatrix Read Error: reading value\n";
1186                         exit(1);
1187 #else
1188                         if (m.issym())
1189                             throw SymMatrixReadError<T>(i,j,m,reader.getis());
1190                         else
1191                             throw HermMatrixReadError<T>(i,j,m,reader.getis());
1192 #endif
1193                     }
1194                     m.ref(i,j) = temp;
1195                 }
1196             }
1197             if (!reader.readRParen(exp,got)) {
1198 #ifdef NOTHROW
1199                 std::cerr<<"SymMatrix Read Error: "<<got<<" != "<<exp<<std::endl;
1200                 exit(1);
1201 #else
1202                 if (m.issym())
1203                     throw SymMatrixReadError<T>(i,N,m,reader.getis(),exp,got);
1204                 else
1205                     throw HermMatrixReadError<T>(i,N,m,reader.getis(),exp,got);
1206 #endif
1207             }
1208             if (i < N-1 && !reader.readRowEnd(exp,got)) {
1209 #ifdef NOTHROW
1210                 std::cerr<<"SymMatrix Read Error: "<<got<<" != "<<exp<<std::endl;
1211                 exit(1);
1212 #else
1213                 if (m.issym())
1214                     throw SymMatrixReadError<T>(i,N,m,reader.getis(),exp,got);
1215                 else
1216                     throw HermMatrixReadError<T>(i,N,m,reader.getis(),exp,got);
1217 #endif
1218             }
1219         }
1220         if (!reader.readFinal(exp,got)) {
1221 #ifdef NOTHROW
1222             std::cerr<<"SymMatrix Read Error: "<<got<<" != "<<exp<<std::endl;
1223             exit(1);
1224 #else
1225             if (m.issym())
1226                 throw SymMatrixReadError<T>(N,0,m,reader.getis(),exp,got);
1227             else
1228                 throw HermMatrixReadError<T>(N,0,m,reader.getis(),exp,got);
1229 #endif
1230         }
1231     }
1232 
1233     template <class T, int A>
read(const TMV_Reader & reader)1234     void SymMatrix<T,A>::read(const TMV_Reader& reader)
1235     {
1236         std::string exp,got;
1237         bool ok = isReal(T()) ?
1238             reader.readCode("S","H",exp,got) : reader.readCode("S",exp,got);
1239         if (!ok) {
1240 #ifdef NOTHROW
1241             std::cerr<<"SymMatrix Read Error: "<<got<<" != "<<exp<<std::endl;
1242             exit(1);
1243 #else
1244             throw SymMatrixReadError<T>(reader.getis(),exp,got);
1245 #endif
1246         }
1247         ptrdiff_t s=size();
1248         if (!reader.readSize(s,exp,got)) {
1249 #ifdef NOTHROW
1250             std::cerr<<"SymMatrix Read Error: reading size\n";
1251             exit(1);
1252 #else
1253             throw SymMatrixReadError<T>(reader.getis(),exp,got);
1254 #endif
1255         }
1256         if (s != size()) resize(s);
1257         s=size();
1258         if (!reader.readSimpleSize(s,exp,got)) {
1259 #ifdef NOTHROW
1260             std::cerr<<"SymMatrix Read Error: reading size\n";
1261             exit(1);
1262 #else
1263             throw SymMatrixReadError<T>(reader.getis(),exp,got);
1264 #endif
1265         }
1266         if (s != size()) {
1267 #ifdef NOTHROW
1268             std::cerr<<"SymMatrix Read Error: Wrong size\n";
1269             exit(1);
1270 #else
1271             throw SymMatrixReadError<T>(*this,reader.getis(),s);
1272 #endif
1273         }
1274         SymMatrixView<T> v = view();
1275         FinishRead(reader,v);
1276     }
1277 
1278     template <class T, int A>
read(const TMV_Reader & reader)1279     void HermMatrix<T,A>::read(const TMV_Reader& reader)
1280     {
1281         std::string exp,got;
1282         bool ok = isReal(T()) ?
1283             reader.readCode("S","H",exp,got) : reader.readCode("H",exp,got);
1284         if (!ok) {
1285 #ifdef NOTHROW
1286             std::cerr<<"HermMatrix Read Error: "<<got<<" != "<<exp<<std::endl;
1287             exit(1);
1288 #else
1289             throw HermMatrixReadError<T>(reader.getis(),exp,got);
1290 #endif
1291         }
1292         ptrdiff_t s=size();
1293         if (!reader.readSize(s,exp,got)) {
1294 #ifdef NOTHROW
1295             std::cerr<<"HermMatrix Read Error: reading size\n";
1296             exit(1);
1297 #else
1298             throw HermMatrixReadError<T>(reader.getis(),exp,got);
1299 #endif
1300         }
1301         if (s != size()) resize(s);
1302         s=size();
1303         if (!reader.readSimpleSize(s,exp,got)) {
1304 #ifdef NOTHROW
1305             std::cerr<<"HermMatrix Read Error: reading size\n";
1306             exit(1);
1307 #else
1308             throw HermMatrixReadError<T>(reader.getis(),exp,got);
1309 #endif
1310         }
1311         if (s != size()) {
1312 #ifdef NOTHROW
1313             std::cerr<<"HermMatrix Read Error: Wrong size\n";
1314             exit(1);
1315 #else
1316             throw HermMatrixReadError<T>(*this,reader.getis(),s);
1317 #endif
1318         }
1319         SymMatrixView<T> v = view();
1320         FinishRead(reader,v);
1321     }
1322 
1323     template <class T, int A>
read(const TMV_Reader & reader)1324     void SymMatrixView<T,A>::read(const TMV_Reader& reader)
1325     {
1326         std::string exp,got;
1327         bool ok = isReal(T()) ?
1328             reader.readCode("S","H",exp,got) :
1329             reader.readCode(issym()?"S":"H",exp,got);
1330         if (!ok) {
1331 #ifdef NOTHROW
1332             std::cerr<<"SymMatrix Read Error: "<<got<<" != "<<exp<<std::endl;
1333             exit(1);
1334 #else
1335             if (issym())
1336                 throw SymMatrixReadError<T>(reader.getis(),exp,got);
1337             else
1338                 throw HermMatrixReadError<T>(reader.getis(),exp,got);
1339 #endif
1340         }
1341         ptrdiff_t s=size();
1342         if (!reader.readSize(s,exp,got)) {
1343 #ifdef NOTHROW
1344             std::cerr<<"SymMatrix Read Error: reading size\n";
1345             exit(1);
1346 #else
1347             if (issym())
1348                 throw SymMatrixReadError<T>(reader.getis(),exp,got);
1349             else
1350                 throw HermMatrixReadError<T>(reader.getis(),exp,got);
1351 #endif
1352         }
1353         if (s != size()) {
1354 #ifdef NOTHROW
1355             std::cerr<<"SymMatrix Read Error: Wrong size\n";
1356             exit(1);
1357 #else
1358             if (issym())
1359                 throw SymMatrixReadError<T>(*this,reader.getis(),s);
1360             else
1361                 throw HermMatrixReadError<T>(*this,reader.getis(),s);
1362 #endif
1363         }
1364         s=size();
1365         if (!reader.readSimpleSize(s,exp,got)) {
1366 #ifdef NOTHROW
1367             std::cerr<<"SymMatrix Read Error: reading size\n";
1368             exit(1);
1369 #else
1370             if (issym())
1371                 throw SymMatrixReadError<T>(reader.getis(),exp,got);
1372             else
1373                 throw HermMatrixReadError<T>(reader.getis(),exp,got);
1374 #endif
1375         }
1376         if (s != size()) {
1377 #ifdef NOTHROW
1378             std::cerr<<"SymMatrix Read Error: Wrong size\n";
1379             exit(1);
1380 #else
1381             if (issym())
1382                 throw SymMatrixReadError<T>(*this,reader.getis(),s);
1383             else
1384                 throw HermMatrixReadError<T>(*this,reader.getis(),s);
1385 #endif
1386         }
1387         FinishRead(reader,*this);
1388     }
1389 
1390 
1391 
1392 #define InstFile "TMV_SymMatrix.inst"
1393 #include "TMV_Inst.h"
1394 #undef InstFile
1395 
1396 } // namespace tmv
1397 
1398 
1399