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 #include "TMV_Blas.h"
25 #include "tmv/TMV_Matrix.h"
26 #include "tmv/TMV_Vector.h"
27 #include "tmv/TMV_DiagMatrix.h"
28 #include "tmv/TMV_Divider.h"
29 #include "tmv/TMV_LUD.h"
30 #include "tmv/TMV_QRD.h"
31 #include "tmv/TMV_QRPD.h"
32 #include "tmv/TMV_SVD.h"
33 #include "tmv/TMV_MatrixArith.h"
34 #include "tmv/TMV_VIt.h"
35 #include "TMV_IntegerDet.h"
36 #include <iostream>
37 
38 namespace tmv {
39 
40 #define RT TMV_RealType(T)
41 
42 #ifdef TMV_BLOCKSIZE
43 #define PERM_BLOCKSIZE TMV_BLOCKSIZE/2
44 #else
45 #define PERM_BLOCKSIZE 32
46 #endif
47 
48     //
49     // Access
50     //
51 
52     template <class T>
cref(ptrdiff_t i,ptrdiff_t j) const53     T GenMatrix<T>::cref(ptrdiff_t i, ptrdiff_t j) const
54     {
55         const T* mi = cptr() + i*stepi() + j*stepj();
56         return isconj() ? TMV_CONJ(*mi) : *mi;
57     }
58 
59     template <class T, int A>
ref(ptrdiff_t i,ptrdiff_t j)60     typename MatrixView<T,A>::reference MatrixView<T,A>::ref(ptrdiff_t i, ptrdiff_t j)
61     {
62         T* mi = ptr() + i*itssi + j*stepj();
63         return RefHelper<T>::makeRef(mi,ct());
64     }
65 
66     template <class T>
setDiv() const67     void GenMatrix<T>::setDiv() const
68     {
69         if (!this->divIsSet()) {
70             DivType dt = this->getDivType();
71             TMVAssert(dt == tmv::LU || dt == tmv::QR ||
72                       dt == tmv::QRP || dt == tmv::SV);
73             switch (dt) {
74               case LU :
75                    this->divider.reset(
76                        new LUDiv<T>(*this,this->divIsInPlace()));
77                    break;
78               case QR :
79                    this->divider.reset(
80                        new QRDiv<T>(*this,this->divIsInPlace()));
81                    break;
82               case QRP :
83                    this->divider.reset(
84                        new QRPDiv<T>(*this,this->divIsInPlace()));
85                    break;
86               case SV :
87                    this->divider.reset(
88                        new SVDiv<T>(*this,this->divIsInPlace()));
89                    break;
90               default :
91                    // The above assert should have already failed
92                    // so go ahead and fall through.
93                    break;
94             }
95         }
96     }
97 
98 #ifdef INST_INT
99     template <>
setDiv() const100     void GenMatrix<int>::setDiv() const
101     { TMVAssert(TMV_FALSE); }
102     template <>
setDiv() const103     void GenMatrix<std::complex<int> >::setDiv() const
104     { TMVAssert(TMV_FALSE); }
105 #endif
106 
107     // Note: These need to be in the .cpp file, not the .h file for
108     // dynamic libraries.  Apparently the typeinfo used for dynamic_cast
109     // doesn't get shared correctly by different modules, so the
110     // dynamic_cast fails when called in one module for an object that
111     // was created in a different module.
112     //
113     // So putting these functions here puts the dynamic cast in the shared
114     // library, which is also where it is created (by setDiv() above).
115     template <class T>
divIsLUDiv() const116     bool GenMatrix<T>::divIsLUDiv() const
117     { return dynamic_cast<const LUDiv<T>*>(this->getDiv()); }
118 
119     template <class T>
divIsQRDiv() const120     bool GenMatrix<T>::divIsQRDiv() const
121     { return dynamic_cast<const QRDiv<T>*>(this->getDiv()); }
122 
123     template <class T>
divIsQRPDiv() const124     bool GenMatrix<T>::divIsQRPDiv() const
125     { return dynamic_cast<const QRPDiv<T>*>(this->getDiv()); }
126 
127     template <class T>
divIsSVDiv() const128     bool GenMatrix<T>::divIsSVDiv() const
129     { return dynamic_cast<const SVDiv<T>*>(this->getDiv()); }
130 
131 
132 #ifdef INST_INT
133     template <>
divIsLUDiv() const134     bool GenMatrix<int>::divIsLUDiv() const
135     { return false; }
136     template <>
divIsQRDiv() const137     bool GenMatrix<int>::divIsQRDiv() const
138     { return false; }
139     template <>
divIsQRPDiv() const140     bool GenMatrix<int>::divIsQRPDiv() const
141     { return false; }
142     template <>
divIsSVDiv() const143     bool GenMatrix<int>::divIsSVDiv() const
144     { return false; }
145 
146     template <>
divIsLUDiv() const147     bool GenMatrix<std::complex<int> >::divIsLUDiv() const
148     { return false; }
149     template <>
divIsQRDiv() const150     bool GenMatrix<std::complex<int> >::divIsQRDiv() const
151     { return false; }
152     template <>
divIsQRPDiv() const153     bool GenMatrix<std::complex<int> >::divIsQRPDiv() const
154     { return false; }
155     template <>
divIsSVDiv() const156     bool GenMatrix<std::complex<int> >::divIsSVDiv() const
157     { return false; }
158 #endif
159 
160     //
161     // OK? (SubMatrix, SubVector)
162     //
163 
164     template <class T>
hasSubMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t j1,ptrdiff_t j2,ptrdiff_t istep,ptrdiff_t jstep) const165     bool GenMatrix<T>::hasSubMatrix(
166         ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const
167     {
168         if (i1==i2 || j1==j2) return true; // no elements, so whatever...
169         bool ok = true;
170         if (istep == 0) {
171             ok = false;
172             std::cerr<<"istep ("<<istep<<") can not be 0\n";
173         }
174         if (i1 < 0 || i1 >= colsize()) {
175             ok = false;
176             std::cerr<<"first col element ("<<i1<<") must be in 0 -- ";
177             std::cerr<<colsize()-1<<std::endl;
178         }
179         if (i2-istep < 0 || i2-istep >= colsize()) {
180             ok = false;
181             std::cerr<<"last col element ("<<i2-istep<<") must be in 0 -- ";
182             std::cerr<<colsize()-1<<std::endl;
183         }
184         if ((i2-i1)%istep != 0) {
185             ok = false;
186             std::cerr<<"col range ("<<i2-i1<<") must be multiple of istep (";
187             std::cerr<<istep<<")\n";
188         }
189         if ((i2-i1)/istep < 0) {
190             ok = false;
191             std::cerr<<"n col elements ("<<(i2-i1)/istep<<") must be nonnegative\n";
192         }
193         if (jstep == 0) {
194             ok = false;
195             std::cerr<<"jstep ("<<jstep<<") can not be 0\n";
196         }
197         if (j1 < 0 || j1 >= rowsize()) {
198             ok = false;
199             std::cerr<<"first row element ("<<j1<<") must be in 0 -- ";
200             std::cerr<<rowsize()-1<<std::endl;
201         }
202         if (j2-jstep < 0 || j2-jstep >= rowsize()) {
203             ok = false;
204             std::cerr<<"last row element ("<<j2-jstep<<") must be in 0 -- ";
205             std::cerr<<rowsize()-1<<std::endl;
206         }
207         if ((j2-j1)%jstep != 0) {
208             ok = false;
209             std::cerr<<"row range ("<<j2-j1<<") must be multiple of istep (";
210             std::cerr<<jstep<<")\n";
211         }
212         if ((j2-j1)/jstep < 0) {
213             ok = false;
214             std::cerr<<"n row elements ("<<(j2-j1)/jstep<<") must be nonnegative\n";
215         }
216         return ok;
217     }
218 
219     template <class T>
hasSubVector(ptrdiff_t i,ptrdiff_t j,ptrdiff_t istep,ptrdiff_t jstep,ptrdiff_t size) const220     bool GenMatrix<T>::hasSubVector(
221         ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t size) const
222     {
223         if (size==0) return true;
224         bool ok = true;
225         if (istep == 0 && jstep == 0) {
226             ok = false;
227             std::cerr<<"istep ("<<istep<<") and jstep ("<<jstep;
228             std::cerr<<") can not both be 0\n";
229         }
230         if (i < 0 || i >= colsize()) {
231             ok = false;
232             std::cerr<<"i ("<<i<<") must be in 0 -- "<<colsize()-1<<std::endl;
233         }
234         if (j < 0 || j >= rowsize()) {
235             ok = false;
236             std::cerr<<"j ("<<j<<") must be in 0 -- "<<rowsize()-1<<std::endl;
237         }
238         ptrdiff_t i2 = i+istep*(size-1);
239         ptrdiff_t j2 = j+jstep*(size-1);
240         if (i2 < 0 || i2 >= colsize()) {
241             ok = false;
242             std::cerr<<"last element's i ("<<i2<<") must be in 0 -- ";
243             std::cerr<<colsize()-1<<std::endl;
244         }
245         if (j2 < 0 || j2 >= rowsize()) {
246             ok = false;
247             std::cerr<<"last element's j ("<<j2<<") must be in 0 -- ";
248             std::cerr<<rowsize()-1<<std::endl;
249         }
250         return ok;
251     }
252 
253     template <class T>
hasSubMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t j1,ptrdiff_t j2,ptrdiff_t istep,ptrdiff_t jstep) const254     bool ConstMatrixView<T,FortranStyle>::hasSubMatrix(
255         ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const
256     {
257         if (i1==i2 || j1==j2) return true; // no elements, so whatever...
258         bool ok = true;
259         if (istep == 0) {
260             ok = false;
261             std::cerr<<"istep ("<<istep<<") can not be 0\n";
262         }
263         if (i1 < 1 || i1 > this->colsize()) {
264             ok = false;
265             std::cerr<<"first col element ("<<i1<<") must be in 1 -- ";
266             std::cerr<<this->colsize()<<std::endl;
267         }
268         if (i2 < 1 || i2 > this->colsize()) {
269             ok = false;
270             std::cerr<<"last col element ("<<i2<<") must be in 1 -- ";
271             std::cerr<<this->colsize()<<std::endl;
272         }
273         if ((i2-i1)%istep != 0) {
274             ok = false;
275             std::cerr<<"col range ("<<i2-i1<<") must be multiple of istep (";
276             std::cerr<<istep<<")\n";
277         }
278         if ((i2-i1)/istep < 0) {
279             ok = false;
280             std::cerr<<"n col elements ("<<(i2-i1)/istep+1<<") must be positive\n";
281         }
282         if (jstep == 0) {
283             ok = false;
284             std::cerr<<"jstep ("<<jstep<<") can not be 0\n";
285         }
286         if (j1 < 1 || j1 > this->rowsize()) {
287             ok = false;
288             std::cerr<<"first row element ("<<j1<<") must be in 1 -- ";
289             std::cerr<<this->rowsize()<<std::endl;
290         }
291         if (j2 < 1 || j2 > this->rowsize()) {
292             ok = false;
293             std::cerr<<"last row element ("<<j2<<") must be in 1 -- ";
294             std::cerr<<this->rowsize()<<std::endl;
295         }
296         if ((j2-j1)%jstep != 0) {
297             ok = false;
298             std::cerr<<"row range ("<<j2-j1<<") must be multiple of istep (";
299             std::cerr<<jstep<<")\n";
300         }
301         if ((j2-j1)/jstep < 0) {
302             ok = false;
303             std::cerr<<"n row elements ("<<(j2-j1)/jstep+1<<") must be positive\n";
304         }
305         return ok;
306     }
307 
308     template <class T>
hasSubVector(ptrdiff_t i,ptrdiff_t j,ptrdiff_t istep,ptrdiff_t jstep,ptrdiff_t size) const309     bool ConstMatrixView<T,FortranStyle>::hasSubVector(
310         ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t size) const
311     {
312         if (size==0) return true;
313         bool ok = true;
314         if (istep == 0 && jstep == 0) {
315             ok = false;
316             std::cerr<<"istep ("<<istep<<") and jstep ("<<jstep;
317             std::cerr<<") can not both be 0\n";
318         }
319         if (i < 1 || i > this->colsize()) {
320             ok = false;
321             std::cerr<<"i ("<<i<<") must be in 1 -- "<<this->colsize()<<std::endl;
322         }
323         if (j < 1 || j > this->rowsize()) {
324             ok = false;
325             std::cerr<<"j ("<<j<<") must be in 1 -- "<<this->rowsize()<<std::endl;
326         }
327         ptrdiff_t i2 = i+istep*(size-1);
328         ptrdiff_t j2 = j+jstep*(size-1);
329         if (i2 < 1 || i2 > this->colsize()) {
330             ok = false;
331             std::cerr<<"last element's i ("<<i2<<") must be in 1 -- ";
332             std::cerr<<this->colsize()<<std::endl;
333         }
334         if (j2 < 1 || j2 > this->rowsize()) {
335             ok = false;
336             std::cerr<<"last element's j ("<<j2<<") must be in 1 -- ";
337             std::cerr<<this->rowsize()<<std::endl;
338         }
339         return ok;
340     }
341 
342 
343     //
344     // Norms, det, etc.
345     //
346 
347     template <class T>
det() const348     T GenMatrix<T>::det() const
349     { return DivHelper<T>::det(); }
350 
351     template <class T>
logDet(T * sign) const352     RT GenMatrix<T>::logDet(T* sign) const
353     { return DivHelper<T>::logDet(sign); }
354 
355     template <class T>
isSingular() const356     bool GenMatrix<T>::isSingular() const
357     { return DivHelper<T>::isSingular(); }
358 
359 #ifdef INST_INT
360     template <>
det() const361     int GenMatrix<int>::det() const
362     { return IntegerDet(*this); }
363 
364     template <>
det() const365     std::complex<int> GenMatrix<std::complex<int> >::det() const
366     { return IntegerDet(*this); }
367 
368     template <>
logDet(int *) const369     int GenMatrix<int>::logDet(int* ) const
370     { TMVAssert(TMV_FALSE); return 0; }
371 
372     template <>
logDet(std::complex<int> *) const373     int GenMatrix<std::complex<int> >::logDet(std::complex<int>* ) const
374     { TMVAssert(TMV_FALSE); return 0; }
375 
376     template <>
isSingular() const377     bool GenMatrix<int>::isSingular() const
378     { return det() == 0; }
379 
380     template <>
isSingular() const381     bool GenMatrix<std::complex<int> >::isSingular() const
382     { return det() == 0; }
383 #endif
384 
385     template <class T>
sumElements() const386     T GenMatrix<T>::sumElements() const
387     {
388         if (canLinearize()) return constLinearView().sumElements();
389         else {
390             T sum(0);
391             if (iscm()) {
392                 const ptrdiff_t N = rowsize();
393                 for(ptrdiff_t j=0;j<N;++j) sum += col(j).sumElements();
394             } else {
395                 const ptrdiff_t M = colsize();
396                 for(ptrdiff_t i=0;i<M;++i) sum += row(i).sumElements();
397             }
398             return sum;
399         }
400     }
401 
402     template <class T>
DoSumAbsElements(const GenMatrix<T> & m)403     static RT DoSumAbsElements(const GenMatrix<T>& m)
404     {
405         if (m.canLinearize()) return m.constLinearView().sumAbsElements();
406         else {
407             RT sum(0);
408             if (m.iscm()) {
409                 const ptrdiff_t N = m.rowsize();
410                 for(ptrdiff_t j=0;j<N;++j) sum += m.col(j).sumAbsElements();
411             } else {
412                 const ptrdiff_t M = m.colsize();
413                 for(ptrdiff_t i=0;i<M;++i) sum += m.row(i).sumAbsElements();
414             }
415             return sum;
416         }
417     }
418 
419     template <class T>
DoSumAbs2Elements(const GenMatrix<T> & m)420     static RT DoSumAbs2Elements(const GenMatrix<T>& m)
421     {
422         if (m.canLinearize()) return m.constLinearView().sumAbs2Elements();
423         else {
424             RT sum(0);
425             if (m.iscm()) {
426                 const ptrdiff_t N = m.rowsize();
427                 for(ptrdiff_t j=0;j<N;++j) sum += m.col(j).sumAbs2Elements();
428             } else {
429                 const ptrdiff_t M = m.colsize();
430                 for(ptrdiff_t i=0;i<M;++i) sum += m.row(i).sumAbs2Elements();
431             }
432             return sum;
433         }
434     }
435 
436 #ifdef INST_INT
DoSumAbsElements(const GenMatrix<std::complex<int>> &)437     static int DoSumAbsElements(const GenMatrix<std::complex<int> >& )
438     { TMVAssert(TMV_FALSE); return 0; }
439 #endif
440 
441     template <class T>
sumAbsElements() const442     RT GenMatrix<T>::sumAbsElements() const
443     { return DoSumAbsElements(*this); }
444 
445     template <class T>
sumAbs2Elements() const446     RT GenMatrix<T>::sumAbs2Elements() const
447     { return DoSumAbs2Elements(*this); }
448 
449     template <class T>
normSq(const RT scale) const450     RT GenMatrix<T>::normSq(const RT scale) const
451     {
452         if (canLinearize()) return constLinearView().normSq(scale);
453         else {
454             RT sum(0);
455             if (isrm()) {
456                 const ptrdiff_t M = colsize();
457                 for(ptrdiff_t i=0;i<M;++i) sum += row(i).normSq(scale);
458             } else {
459                 const ptrdiff_t N = rowsize();
460                 for(ptrdiff_t j=0;j<N;++j) sum += col(j).normSq(scale);
461             }
462             return sum;
463         }
464     }
465 
466     template <class T>
NonLapMaxAbsElement(const GenMatrix<T> & m)467     static RT NonLapMaxAbsElement(const GenMatrix<T>& m)
468     {
469         if (m.canLinearize()) return m.constLinearView().maxAbsElement();
470         else {
471             RT max(0);
472             if (m.iscm()) {
473                 const ptrdiff_t N = m.rowsize();
474                 for(ptrdiff_t j=0;j<N;++j) {
475                     RT temp = m.col(j).normInf();
476                     if (temp > max) max = temp;
477                 }
478             } else {
479                 const ptrdiff_t M = m.colsize();
480                 for(ptrdiff_t i=0;i<M;++i) {
481                     RT temp = m.row(i).normInf();
482                     if (temp > max) max = temp;
483                 }
484             }
485             return max;
486         }
487     }
488 
489     template <class T>
NonLapMaxAbs2Element(const GenMatrix<T> & m)490     static RT NonLapMaxAbs2Element(const GenMatrix<T>& m)
491     {
492         if (m.canLinearize()) return m.constLinearView().maxAbs2Element();
493         else {
494             RT max(0);
495             if (m.iscm()) {
496                 const ptrdiff_t N = m.rowsize();
497                 for(ptrdiff_t j=0;j<N;++j) {
498                     RT temp = m.col(j).maxAbs2Element();
499                     if (temp > max) max = temp;
500                 }
501             } else {
502                 const ptrdiff_t M = m.colsize();
503                 for(ptrdiff_t i=0;i<M;++i) {
504                     RT temp = m.row(i).maxAbs2Element();
505                     if (temp > max) max = temp;
506                 }
507             }
508             return max;
509         }
510     }
511 
512     template <class T>
NonLapNorm1(const GenMatrix<T> & m)513     static RT NonLapNorm1(const GenMatrix<T>& m)
514     {
515         RT max(0);
516         const ptrdiff_t N = m.rowsize();
517         for(ptrdiff_t j=0;j<N;++j) {
518             RT temp = m.col(j).norm1();
519             if (temp > max) max = temp;
520         }
521         return max;
522     }
523 
524     template <class T>
NonLapNormF(const GenMatrix<T> & m)525     static RT NonLapNormF(const GenMatrix<T>& m)
526     {
527         const RT eps = TMV_Epsilon<T>();
528 
529         RT mmax = m.maxAbs2Element();
530         RT norm;
531         if (mmax == RT(0)) {
532             // Then norm is also 0
533             norm = RT(0);
534         } else if (TMV_Underflow(mmax * mmax)) {
535             // Then we need to rescale, since underflow has caused
536             // rounding errors.
537             // Epsilon is a pure power of 2, so no rounding errors from
538             // rescaling.
539             const RT inveps = RT(1)/eps;
540             RT scale = inveps;
541             mmax *= scale;
542             const RT eps2 = eps*eps;
543             while (mmax < eps2) { scale *= inveps; mmax *= inveps; }
544             norm = TMV_SQRT(m.normSq(scale))/scale;
545         } else if (RT(1) / mmax == RT(0)) {
546             // Then mmax is already inf, so no hope of making it more accurate.
547             norm = mmax;
548         } else if (RT(1) / (mmax*mmax) == RT(0)) {
549             // Then we have overflow, so we need to rescale:
550             const RT inveps = RT(1)/eps;
551             RT scale = eps;
552             mmax *= scale;
553             while (mmax > inveps) { scale *= eps; mmax *= eps; }
554             norm = TMV_SQRT(m.normSq(scale))/scale;
555         } else {
556             norm = TMV_SQRT(m.normSq());
557         }
558         return norm;
559     }
560 
561     template <class T>
NonLapNormInf(const GenMatrix<T> & m)562     static inline RT NonLapNormInf(const GenMatrix<T>& m)
563     { return NonLapNorm1(m.transpose()); }
564 
565 #ifdef XLAP
566     template <class T>
LapNorm(const char c,const GenMatrix<T> & m)567     static RT LapNorm(const char c, const GenMatrix<T>& m)
568     {
569         switch(c) {
570           case 'M' : return NonLapMaxAbsElement(m);
571           case '1' : return NonLapNorm1(m);
572           case 'F' : return NonLapNormF(m);
573           case 'I' : return NonLapNormInf(m);
574           default : TMVAssert(TMV_FALSE);
575         }
576         return RT(0);
577     }
578 #ifdef INST_DOUBLE
579     template <>
LapNorm(const char c,const GenMatrix<double> & m)580     double LapNorm(const char c, const GenMatrix<double>& m)
581     {
582         TMVAssert(m.iscm());
583         char cc = c;
584         int M = m.colsize();
585         int N = m.rowsize();
586         int lda = m.stepj();
587 #ifndef LAPNOWORK
588         int lwork = c=='I' ? M : 0;
589         AlignedArray<double> work(lwork);
590         VectorViewOf(work.get(),lwork).setZero();
591 #endif
592         double norm = LAPNAME(dlange) (
593             LAPCM LAPV(cc),LAPV(M),LAPV(N),
594             LAPP(m.cptr()),LAPV(lda) LAPWK(work.get()));
595         return norm;
596     }
597     template <>
LapNorm(const char c,const GenMatrix<std::complex<double>> & m)598     double LapNorm(const char c, const GenMatrix<std::complex<double> >& m)
599     {
600         TMVAssert(m.iscm());
601         char cc = c;
602         int M = m.colsize();
603         int N = m.rowsize();
604         int lda = m.stepj();
605 #ifndef LAPNOWORK
606         int lwork = c=='I' ? M : 0;
607         AlignedArray<double> work(lwork);
608         VectorViewOf(work.get(),lwork).setZero();
609 #endif
610         double norm = LAPNAME(zlange) (
611             LAPCM LAPV(cc),LAPV(M),LAPV(N),
612             LAPP(m.cptr()),LAPV(lda) LAPWK(work.get()));
613         return norm;
614     }
615 #endif
616 #ifdef INST_FLOAT
617     template <>
LapNorm(const char c,const GenMatrix<float> & m)618     float LapNorm(const char c, const GenMatrix<float>& m)
619     {
620         TMVAssert(m.iscm());
621         char cc = c;
622         int M = m.colsize();
623         int N = m.rowsize();
624         int lda = m.stepj();
625 #ifndef LAPNOWORK
626         int lwork = c=='I' ? M : 0;
627         AlignedArray<float> work(lwork);
628         VectorViewOf(work.get(),lwork).setZero();
629 #endif
630         double norm = LAPNAME(slange) (
631             LAPCM LAPV(cc),LAPV(M),LAPV(N),
632             LAPP(m.cptr()),LAPV(lda) LAPWK(work.get()));
633         return norm;
634     }
635     template <>
LapNorm(const char c,const GenMatrix<std::complex<float>> & m)636     float LapNorm(const char c, const GenMatrix<std::complex<float> >& m)
637     {
638         TMVAssert(m.iscm());
639         char cc = c;
640         int M = m.colsize();
641         int N = m.rowsize();
642         int lda = m.stepj();
643 #ifndef LAPNOWORK
644         int lwork = c=='I' ? M : 0;
645         AlignedArray<float> work(lwork);
646         VectorViewOf(work.get(),lwork).setZero();
647 #endif
648         double norm = LAPNAME(clange) (
649             LAPCM LAPV(cc),LAPV(M),LAPV(N),
650             LAPP(m.cptr()),LAPV(lda) LAPWK(work.get()));
651         return norm;
652     }
653 #endif
654 #endif // XLAP
655 
656 #ifdef INST_INT
NonLapNormF(const GenMatrix<int> &)657     static int NonLapNormF(const GenMatrix<int>& )
658     { TMVAssert(TMV_FALSE); return 0; }
NonLapNormF(const GenMatrix<std::complex<int>> &)659     static int NonLapNormF(const GenMatrix<std::complex<int> >& )
660     { TMVAssert(TMV_FALSE); return 0; }
NonLapNorm1(const GenMatrix<std::complex<int>> &)661     static int NonLapNorm1(const GenMatrix<std::complex<int> >& )
662     { TMVAssert(TMV_FALSE); return 0; }
NonLapMaxAbsElement(const GenMatrix<std::complex<int>> &)663     static int NonLapMaxAbsElement(const GenMatrix<std::complex<int> >& )
664     { TMVAssert(TMV_FALSE); return 0; }
665 #endif
666 
667     template <class T>
maxAbsElement() const668     RT GenMatrix<T>::maxAbsElement() const
669     {
670 #ifdef XLAP
671         if (isrm() && stepi() > 0) return LapNorm('M',transpose());
672         else if (iscm() && stepj() > 0) return LapNorm('M',*this);
673         else
674 #endif
675             return NonLapMaxAbsElement(*this);
676     }
677     template <class T>
maxAbs2Element() const678     RT GenMatrix<T>::maxAbs2Element() const
679     {
680 #ifdef XLAP
681         if (Traits<T>::iscomplex) return NonLapMaxAbs2Element(*this);
682         else if (isrm() && stepi() > 0) return LapNorm('M',transpose());
683         else if (iscm() && stepj() > 0) return LapNorm('M',*this);
684         else
685 #endif
686             return NonLapMaxAbs2Element(*this);
687     }
688     template <class T>
norm1() const689     RT GenMatrix<T>::norm1() const
690     {
691 #ifdef XLAP
692         if (isrm() && stepi() > 0) return LapNorm('I',transpose());
693         else if (iscm() && stepj() > 0) return LapNorm('1',*this);
694         else
695 #endif
696             return NonLapNorm1(*this);
697     }
698     template <class T>
normF() const699     RT GenMatrix<T>::normF() const
700     {
701 #ifdef XLAP
702         if (isrm() && stepi() > 0) return LapNorm('F',transpose());
703         else if (iscm() && stepj() > 0) return LapNorm('F',*this);
704         else
705 #endif
706             return NonLapNormF(*this);
707     }
708 
709     template <class T>
DoNorm2(const GenMatrix<T> & m)710     static RT DoNorm2(const GenMatrix<T>& m)
711     {
712         if (m.colsize() < m.rowsize()) return DoNorm2(m.transpose());
713         if (m.rowsize() == 0) return RT(0);
714         Matrix<T> m2(m);
715         DiagMatrix<RT> S(m.rowsize());
716         SV_Decompose(m2.view(),S.view(),false);
717         return S(0);
718     }
719 
720     template <class T>
DoCondition(const GenMatrix<T> & m)721     static RT DoCondition(const GenMatrix<T>& m)
722     {
723         if (m.colsize() < m.rowsize()) return DoCondition(m.transpose());
724         if (m.rowsize() == 0) return RT(1);
725         Matrix<T> m2(m);
726         DiagMatrix<RT> S(m.rowsize());
727         SV_Decompose(m2.view(),S.view(),false);
728         return S(0)/S(S.size()-1);
729     }
730 
731 #ifdef INST_INT
DoNorm2(const GenMatrix<int> &)732     static int DoNorm2(const GenMatrix<int>& )
733     { TMVAssert(TMV_FALSE); return 0; }
DoCondition(const GenMatrix<int> &)734     static int DoCondition(const GenMatrix<int>& )
735     { TMVAssert(TMV_FALSE); return 0; }
DoNorm2(const GenMatrix<std::complex<int>> &)736     static int DoNorm2(const GenMatrix<std::complex<int> >& )
737     { TMVAssert(TMV_FALSE); return 0; }
DoCondition(const GenMatrix<std::complex<int>> &)738     static int DoCondition(const GenMatrix<std::complex<int> >& )
739     { TMVAssert(TMV_FALSE); return 0; }
740 #endif
741 
742     template <class T>
doNorm2() const743     RT GenMatrix<T>::doNorm2() const
744     { return tmv::DoNorm2(*this); }
745 
746     template <class T>
doCondition() const747     RT GenMatrix<T>::doCondition() const
748     { return tmv::DoCondition(*this); }
749 
750     template <class T>
QInverse() const751     QuotXM<T,T> GenMatrix<T>::QInverse() const
752     { return QuotXM<T,T>(T(1),*this); }
753 
754     //
755     // Modifying Functions
756     //
757 
758     template <class T, int A>
clip(RT thresh)759     MatrixView<T,A>& MatrixView<T,A>::clip(RT thresh)
760     {
761         TMVAssert(A==CStyle);
762         if (this->canLinearize()) linearView().clip(thresh);
763         else {
764             if (this->isrm()) {
765                 const ptrdiff_t M = colsize();
766                 for(ptrdiff_t i=0;i<M;++i) row(i).clip(thresh);
767             } else {
768                 const ptrdiff_t N = rowsize();
769                 for(ptrdiff_t j=0;j<N;++j) col(j).clip(thresh);
770             }
771         }
772         return *this;
773     }
774 
775     template <class T, int A>
setZero()776     MatrixView<T,A>& MatrixView<T,A>::setZero()
777     {
778         TMVAssert(A==CStyle);
779         if (this->canLinearize()) linearView().setZero();
780         else {
781             if (this->isrm()) {
782                 const ptrdiff_t M = colsize();
783                 for(ptrdiff_t i=0;i<M;++i) row(i).setZero();
784             } else {
785                 const ptrdiff_t N = rowsize();
786                 for(ptrdiff_t j=0;j<N;++j) col(j).setZero();
787             }
788         }
789         return *this;
790     }
791 
792     template <class T, int A>
setAllTo(const T & x)793     MatrixView<T,A>& MatrixView<T,A>::setAllTo(const T& x)
794     {
795         TMVAssert(A==CStyle);
796         if (this->canLinearize()) linearView().setAllTo(x);
797         else {
798             if (this->isrm()) {
799                 const ptrdiff_t M = colsize();
800                 for(ptrdiff_t i=0;i<M;++i) row(i).setAllTo(x);
801             } else  {
802                 const ptrdiff_t N = rowsize();
803                 for(ptrdiff_t j=0;j<N;++j) col(j).setAllTo(x);
804             }
805         }
806         return *this;
807     }
808 
809     template <class T, int A>
addToAll(const T & x)810     MatrixView<T,A>& MatrixView<T,A>::addToAll(const T& x)
811     {
812         TMVAssert(A==CStyle);
813         if (this->canLinearize()) linearView().addToAll(x);
814         else {
815             if (this->isrm()) {
816                 const ptrdiff_t M = colsize();
817                 for(ptrdiff_t i=0;i<M;++i) row(i).addToAll(x);
818             } else  {
819                 const ptrdiff_t N = rowsize();
820                 for(ptrdiff_t j=0;j<N;++j) col(j).addToAll(x);
821             }
822         }
823         return *this;
824     }
825 
826     template <class T, int A>
transposeSelf()827     MatrixView<T,A>& MatrixView<T,A>::transposeSelf()
828     {
829         TMVAssert(A==CStyle);
830         TMVAssert(colsize() == rowsize());
831         const ptrdiff_t M = colsize();
832         for(ptrdiff_t i=1;i<M;++i) Swap(row(i,0,i),col(i,0,i));
833         return *this;
834     }
835 
836     template <class T, int A>
conjugateSelf()837     MatrixView<T,A>& MatrixView<T,A>::conjugateSelf()
838     {
839         TMVAssert(A==CStyle);
840         if (isComplex(T())) {
841             if (this->canLinearize()) linearView().conjugateSelf();
842             else {
843                 if (this->isrm()) {
844                     const ptrdiff_t M = colsize();
845                     for(ptrdiff_t i=0;i<M;++i) row(i).conjugateSelf();
846                 } else  {
847                     const ptrdiff_t N = rowsize();
848                     for(ptrdiff_t j=0;j<N;++j) col(j).conjugateSelf();
849                 }
850             }
851         }
852         return *this;
853     }
854 
855     template <class T, int A>
setToIdentity(const T & x)856     MatrixView<T,A>& MatrixView<T,A>::setToIdentity(const T& x)
857     {
858         TMVAssert(A==CStyle);
859         TMVAssert(colsize() == rowsize());
860         setZero();
861         diag().setAllTo(x);
862         return *this;
863     }
864 
865     template <class T, int A>
permuteRows(const ptrdiff_t * p,ptrdiff_t i1,ptrdiff_t i2)866     MatrixView<T,A>& MatrixView<T,A>::permuteRows(
867         const ptrdiff_t* p, ptrdiff_t i1, ptrdiff_t i2)
868     {
869         TMVAssert(A==CStyle);
870         TMVAssert(i2<=colsize());
871         TMVAssert(i1<=i2);
872         // This idea of doing the permutations a block at a time
873         // is cribbed from the LAPack code.  It does speed things up
874         // quite a bit for large matrices.  On my machine where BLOCKSIZE=64
875         // is optimal for most routines, blocks of 32 were optimal here,
876         // so I use BLOCKSIZE/2 in general.
877         const ptrdiff_t N = rowsize();
878         const ptrdiff_t Nx = N/PERM_BLOCKSIZE*PERM_BLOCKSIZE;
879         if (Nx != 0) {
880             for(ptrdiff_t j=0;j<Nx;) {
881                 ptrdiff_t j2 = j+PERM_BLOCKSIZE;
882                 const ptrdiff_t* pi = p+i1;
883                 for(ptrdiff_t i=i1;i<i2;++i,++pi) {
884                     TMVAssert(*pi < colsize());
885                     colRange(j,j2).swapRows(i,*pi);
886                 }
887                 j = j2;
888             }
889         }
890         if (Nx != N) {
891             const ptrdiff_t* pi = p+i1;
892             for(ptrdiff_t i=i1;i<i2;++i,++pi) {
893                 TMVAssert(*pi < colsize());
894                 colRange(Nx,N).swapRows(i,*pi);
895             }
896         }
897         return *this;
898     }
899 
900     template <class T, int A>
reversePermuteRows(const ptrdiff_t * p,ptrdiff_t i1,ptrdiff_t i2)901     MatrixView<T,A>& MatrixView<T,A>::reversePermuteRows(
902         const ptrdiff_t* p, ptrdiff_t i1, ptrdiff_t i2)
903     {
904         TMVAssert(A==CStyle);
905         TMVAssert(i2<=colsize());
906         TMVAssert(i1<=i2);
907         const ptrdiff_t N = rowsize();
908         const ptrdiff_t Nx = N/PERM_BLOCKSIZE*PERM_BLOCKSIZE;
909         if (Nx != 0) {
910             for(ptrdiff_t j=0;j<Nx;) {
911                 ptrdiff_t j2 = j+PERM_BLOCKSIZE;
912                 const ptrdiff_t* pi = p+i2;
913                 for(ptrdiff_t i=i2;i>i1;) {
914                     --i; --pi;
915                     TMVAssert(*pi < colsize());
916                     colRange(j,j2).swapRows(i,*pi);
917                 }
918                 j = j2;
919             }
920         }
921         if (Nx != N) {
922             const ptrdiff_t* pi = p+i2;
923             for(ptrdiff_t i=i2;i>i1;) {
924                 --i; --pi;
925                 TMVAssert(*pi < colsize());
926                 colRange(Nx,N).swapRows(i,*pi);
927             }
928         }
929         return *this;
930     }
931 
932     //
933     // Copy Matrices
934     //
935 
936     template <class T>
NonLapCopy(const GenMatrix<T> & m1,MatrixView<T> m2)937     static void NonLapCopy(const GenMatrix<T>& m1, MatrixView<T> m2)
938     {
939         TMVAssert(m2.rowsize() == m1.rowsize());
940         TMVAssert(m2.colsize() == m1.colsize());
941         TMVAssert(m1.ct()==NonConj);
942         TMVAssert(m2.ct()==NonConj);
943         const ptrdiff_t M = m2.colsize();
944         const ptrdiff_t N = m2.rowsize();
945 
946         if (m1.iscm() && m2.iscm()) {
947             const T* p1 = m1.cptr();
948             T* p2 = m2.ptr();
949             const ptrdiff_t s1 = m1.stepj();
950             const ptrdiff_t s2 = m2.stepj();
951             for(ptrdiff_t j=0;j<N;++j,p1+=s1,p2+=s2) {
952                 std::copy(p1,p1+M,p2);
953             }
954         } else if (M > N) {
955             if (shouldReverse(m1.stepi(),m2.stepi()))
956                 for(ptrdiff_t j=0;j<N;++j)
957                     DoCopySameType(m1.col(j).reverse(),m2.col(j).reverse());
958             else
959                 for(ptrdiff_t j=0;j<N;++j)
960                     DoCopySameType(m1.col(j),m2.col(j));
961         } else {
962             if (shouldReverse(m1.stepj(),m2.stepj()))
963                 for(ptrdiff_t i=0;i<M;++i)
964                     DoCopySameType(m1.row(i).reverse(),m2.row(i).reverse());
965             else
966                 for(ptrdiff_t i=0;i<M;++i)
967                     DoCopySameType(m1.row(i),m2.row(i));
968         }
969     }
970 #ifdef ELAP
971     template <class T>
LapCopy(const GenMatrix<T> & m1,MatrixView<T> m2)972     static inline void LapCopy(const GenMatrix<T>& m1, MatrixView<T> m2)
973     { NonLapCopy(m1,m2); }
974 #ifdef INST_DOUBLE
975     template <>
LapCopy(const GenMatrix<double> & m1,MatrixView<double> m2)976     void LapCopy(const GenMatrix<double>& m1, MatrixView<double> m2)
977     {
978         TMVAssert(m2.rowsize() == m1.rowsize());
979         TMVAssert(m2.colsize() == m1.colsize());
980         TMVAssert(m1.iscm());
981         TMVAssert(m2.iscm());
982         char c = 'A';
983         int m = m1.colsize();
984         int n = m1.rowsize();
985         int ld1 = m1.stepj();
986         int ld2 = m2.stepj();
987         LAPNAME(dlacpy) (
988             LAPCM LAPV(c),LAPV(m),LAPV(n),LAPP(m1.cptr()),LAPV(ld1),
989             LAPP(m2.ptr()),LAPV(ld2));
990     }
991     template <>
LapCopy(const GenMatrix<std::complex<double>> & m1,MatrixView<std::complex<double>> m2)992     void LapCopy(
993         const GenMatrix<std::complex<double> >& m1,
994         MatrixView<std::complex<double> > m2)
995     {
996         TMVAssert(m2.rowsize() == m1.rowsize());
997         TMVAssert(m2.colsize() == m1.colsize());
998         TMVAssert(m1.iscm());
999         TMVAssert(m2.iscm());
1000         TMVAssert(m1.ct() == NonConj);
1001         TMVAssert(m2.ct() == NonConj);
1002         char c = 'A';
1003         int m = m1.colsize();
1004         int n = m1.rowsize();
1005         int ld1 = m1.stepj();
1006         int ld2 = m2.stepj();
1007         LAPNAME(zlacpy) (
1008             LAPCM LAPV(c),LAPV(m),LAPV(n),LAPP(m1.cptr()),LAPV(ld1),
1009             LAPP(m2.ptr()),LAPV(ld2));
1010     }
1011 #endif
1012 #ifdef INST_FLOAT
1013     template <>
LapCopy(const GenMatrix<float> & m1,MatrixView<float> m2)1014     void LapCopy(const GenMatrix<float>& m1, MatrixView<float> m2)
1015     {
1016         TMVAssert(m2.rowsize() == m1.rowsize());
1017         TMVAssert(m2.colsize() == m1.colsize());
1018         TMVAssert(m1.iscm());
1019         TMVAssert(m2.iscm());
1020         char c = 'A';
1021         int m = m1.colsize();
1022         int n = m1.rowsize();
1023         int ld1 = m1.stepj();
1024         int ld2 = m2.stepj();
1025         LAPNAME(slacpy) (
1026             LAPCM LAPV(c),LAPV(m),LAPV(n),LAPP(m1.cptr()),LAPV(ld1),
1027             LAPP(m2.ptr()),LAPV(ld2));
1028     }
1029     template <>
LapCopy(const GenMatrix<std::complex<float>> & m1,MatrixView<std::complex<float>> m2)1030     void LapCopy(
1031         const GenMatrix<std::complex<float> >& m1,
1032         MatrixView<std::complex<float> > m2)
1033     {
1034         TMVAssert(m2.rowsize() == m1.rowsize());
1035         TMVAssert(m2.colsize() == m1.colsize());
1036         TMVAssert(m1.iscm());
1037         TMVAssert(m2.iscm());
1038         TMVAssert(m1.ct() == NonConj);
1039         TMVAssert(m2.ct() == NonConj);
1040         char c = 'A';
1041         int m = m1.colsize();
1042         int n = m1.rowsize();
1043         int ld1 = m1.stepj();
1044         int ld2 = m2.stepj();
1045         LAPNAME(clacpy) (
1046             LAPCM LAPV(c),LAPV(m),LAPV(n),LAPP(m1.cptr()),LAPV(ld1),
1047             LAPP(m2.ptr()),LAPV(ld2));
1048     }
1049 #endif
1050 #endif
1051     template <class T>
DoCopySameType(const GenMatrix<T> & m1,MatrixView<T> m2)1052     void DoCopySameType(const GenMatrix<T>& m1, MatrixView<T> m2)
1053     {
1054         TMVAssert(m2.rowsize() == m1.rowsize());
1055         TMVAssert(m2.colsize() == m1.colsize());
1056         TMVAssert(m2.rowsize() > 0);
1057         TMVAssert(m2.colsize() > 0);
1058         TMVAssert(m1.ct() == NonConj);
1059         TMVAssert(m2.ct() == NonConj);
1060         TMVAssert(!m2.isSameAs(m1));
1061 
1062 #ifdef ELAP
1063         if (m1.iscm() && m2.iscm() && m1.stepj()>0 && m2.stepj()>0)
1064             LapCopy(m1,m2);
1065         else
1066 #endif
1067             NonLapCopy(m1,m2);
1068     }
1069 
1070     //
1071     // Swap
1072     //
1073 
1074     template <class T>
Swap(MatrixView<T> m1,MatrixView<T> m2)1075     void Swap(MatrixView<T> m1, MatrixView<T> m2)
1076     {
1077         TMVAssert(m1.colsize() == m2.colsize());
1078         TMVAssert(m1.rowsize() == m2.rowsize());
1079         if (m1.canLinearize() && m2.canLinearize() &&
1080             m1.stepi() == m2.stepi() && m1.stepj() == m2.stepj()) {
1081             TMVAssert(m1.linearView().size() == m2.linearView().size());
1082             TMVAssert(m1.stepi()==m2.stepi() && m1.stepj()==m2.stepj());
1083             Swap(m1.linearView(),m2.linearView());
1084         } else {
1085             if (m1.isrm() && m2.isrm()) {
1086                 const ptrdiff_t M = m1.colsize();
1087                 for(ptrdiff_t i=0;i<M;++i) Swap(m1.row(i),m2.row(i));
1088             } else {
1089                 const ptrdiff_t N = m1.rowsize();
1090                 for(ptrdiff_t j=0;j<N;++j) Swap(m1.col(j),m2.col(j));
1091             }
1092         }
1093     }
1094 
1095     //
1096     // m1 == m2
1097     //
1098 
1099     template <class T1, class T2>
operator ==(const GenMatrix<T1> & m1,const GenMatrix<T2> & m2)1100     bool operator==(const GenMatrix<T1>& m1, const GenMatrix<T2>& m2)
1101     {
1102         if (m1.colsize() != m2.colsize()) return false;
1103         else if (m1.rowsize() != m2.rowsize()) return false;
1104         else if (m1.isSameAs(m2)) return true;
1105         else if (m1.stepi()==m2.stepi() && m1.stepj()==m2.stepj() &&
1106                  m1.canLinearize() && m2.canLinearize())
1107             return m1.constLinearView() == m2.constLinearView();
1108         else {
1109             const ptrdiff_t M = m1.colsize();
1110             for(ptrdiff_t i=0;i<M;++i)
1111                 if (m1.row(i) != m2.row(i)) return false;
1112             return true;
1113         }
1114     }
1115 
1116     //
1117     // I/O
1118     //
1119 
1120     template <class T>
write(const TMV_Writer & writer) const1121     void GenMatrix<T>::write(const TMV_Writer& writer) const
1122     {
1123         const ptrdiff_t M = colsize();
1124         const ptrdiff_t N = rowsize();
1125         writer.begin();
1126         writer.writeCode("M");
1127         writer.writeSize(M);
1128         writer.writeSize(N);
1129         writer.writeStart();
1130         for(ptrdiff_t i=0;i<M;++i) {
1131             writer.writeLParen();
1132             for(ptrdiff_t j=0;j<N;++j) {
1133                 if (j > 0) writer.writeSpace();
1134                 writer.writeValue(cref(i,j));
1135             }
1136             writer.writeRParen();
1137             if (i < M-1) writer.writeRowEnd();
1138         }
1139         writer.writeFinal();
1140         writer.end();
1141     }
1142 
1143 #ifndef NOTHROW
1144     template <class T>
1145     class MatrixReadError : public ReadError
1146     {
1147     public :
1148         Matrix<T> m;
1149         ptrdiff_t i,j;
1150         std::string exp,got;
1151         ptrdiff_t cs,rs;
1152         bool is,iseof,isbad;
1153 
MatrixReadError(std::istream & _is)1154         MatrixReadError(std::istream& _is) throw() :
1155             ReadError("Matrix."),
1156             i(0), j(0), cs(0), rs(0),
1157             is(_is), iseof(_is.eof()), isbad(_is.bad()) {}
MatrixReadError(std::istream & _is,const std::string & _e,const std::string & _g)1158         MatrixReadError(
1159             std::istream& _is,
1160             const std::string& _e, const std::string& _g) throw() :
1161             ReadError("Matrix."),
1162             i(0), j(0), exp(_e), got(_g), cs(0), rs(0),
1163             is(_is), iseof(_is.eof()), isbad(_is.bad()) {}
1164 
MatrixReadError(ptrdiff_t _i,ptrdiff_t _j,const GenMatrix<T> & _m,std::istream & _is)1165         MatrixReadError(
1166             ptrdiff_t _i, ptrdiff_t _j, const GenMatrix<T>& _m, std::istream& _is) throw() :
1167             ReadError("Matrix."),
1168             m(_m), i(_i), j(_j),
1169             cs(m.colsize()), rs(m.rowsize()),
1170             is(_is), iseof(_is.eof()), isbad(_is.bad()) {}
MatrixReadError(ptrdiff_t _i,ptrdiff_t _j,const GenMatrix<T> & _m,std::istream & _is,const std::string & _e,const std::string & _g)1171         MatrixReadError(
1172             ptrdiff_t _i, ptrdiff_t _j, const GenMatrix<T>& _m, std::istream& _is,
1173             const std::string& _e, const std::string& _g) throw() :
1174             ReadError("Matrix."),
1175             m(_m), i(_i), j(_j), exp(_e), got(_g),
1176             cs(m.colsize()), rs(m.rowsize()),
1177             is(_is), iseof(_is.eof()), isbad(_is.bad()) {}
MatrixReadError(const GenMatrix<T> & _m,std::istream & _is,ptrdiff_t _cs,ptrdiff_t _rs)1178         MatrixReadError(
1179             const GenMatrix<T>& _m,
1180             std::istream& _is, ptrdiff_t _cs, ptrdiff_t _rs) throw() :
1181             ReadError("Matrix."),
1182             m(_m), i(0), j(0), cs(_cs), rs(_rs),
1183             is(_is), iseof(_is.eof()), isbad(_is.bad()) {}
1184 
MatrixReadError(const MatrixReadError<T> & rhs)1185         MatrixReadError(const MatrixReadError<T>& rhs) :
1186             m(rhs.m), i(rhs.i), j(rhs.j), exp(rhs.exp), got(rhs.got),
1187             cs(rhs.cs), rs(rhs.rs),
1188             is(rhs.is), iseof(rhs.iseof), isbad(rhs.isbad) {}
~MatrixReadError()1189         ~MatrixReadError() throw() {}
1190 
write(std::ostream & os) const1191         void write(std::ostream& os) const throw()
1192         {
1193             os<<"TMV Read Error: Reading istream input for Matrix\n";
1194             if (exp != got) {
1195                 os<<"Wrong format: expected '"<<exp<<"', got '"<<got<<"'.\n";
1196             }
1197             if (cs != m.colsize()) {
1198                 os<<"Wrong column size: expected "<<m.colsize()<<
1199                     ", got "<<cs<<".\n";
1200             }
1201             if (rs != m.rowsize()) {
1202                 os<<"Wrong row size: expected "<<m.rowsize()<<
1203                     ", got "<<rs<<".\n";
1204             }
1205             if (!is) {
1206                 if (iseof) {
1207                     os<<"Input stream reached end-of-file prematurely.\n";
1208                 } else if (isbad) {
1209                     os<<"Input stream is corrupted.\n";
1210                 } else {
1211                     os<<"Input stream cannot read next character.\n";
1212                 }
1213             }
1214             if (m.colsize() > 0 || m.rowsize() > 0) {
1215                 const ptrdiff_t N = m.rowsize();
1216                 os<<"The portion of the Matrix which was successfully "
1217                     "read is: \n";
1218                 for(ptrdiff_t ii=0;ii<i;++ii) {
1219                     os<<"( ";
1220                     for(ptrdiff_t jj=0;jj<N;++jj) os<<' '<<m.cref(ii,jj)<<' ';
1221                     os<<" )\n";
1222                 }
1223                 os<<"( ";
1224                 for(ptrdiff_t jj=0;jj<j;++jj) os<<' '<<m.cref(i,jj)<<' ';
1225                 os<<" )\n";
1226             }
1227         }
1228     };
1229 #endif
1230 
1231     template <class T>
FinishRead(const TMV_Reader & reader,MatrixView<T> m)1232     static void FinishRead(const TMV_Reader& reader, MatrixView<T> m)
1233     {
1234         const ptrdiff_t M = m.colsize();
1235         const ptrdiff_t N = m.rowsize();
1236         std::string exp, got;
1237         T temp;
1238         if (!reader.readStart(exp,got)) {
1239 #ifdef NOTHROW
1240             std::cerr<<"Matrix Read Error: "<<got<<" != "<<exp<<std::endl;
1241             exit(1);
1242 #else
1243             throw MatrixReadError<T>(0,0,m,reader.getis(),exp,got);
1244 #endif
1245         }
1246         for(ptrdiff_t i=0;i<M;++i) {
1247             if (!reader.readLParen(exp,got)) {
1248 #ifdef NOTHROW
1249                 std::cerr<<"Matrix Read Error: "<<got<<" != "<<exp<<std::endl;
1250                 exit(1);
1251 #else
1252                 throw MatrixReadError<T>(i,0,m,reader.getis(),exp,got);
1253 #endif
1254             }
1255             for(ptrdiff_t j=0;j<N;++j) {
1256                 if (j>0) {
1257                     if (!reader.readSpace(exp,got)) {
1258 #ifdef NOTHROW
1259                         std::cerr<<"Matrix Read Error: "<<got<<" != "<<exp<<std::endl;
1260                         exit(1);
1261 #else
1262                         throw MatrixReadError<T>(i,j,m,reader.getis(),exp,got);
1263 #endif
1264                     }
1265                 }
1266                 if (!reader.readValue(temp)) {
1267 #ifdef NOTHROW
1268                     std::cerr<<"Matrix Read Error: reading value\n";
1269                     exit(1);
1270 #else
1271                     throw MatrixReadError<T>(i,j,m,reader.getis());
1272 #endif
1273                 }
1274                 m.ref(i,j) = temp;
1275             }
1276             if (!reader.readRParen(exp,got)) {
1277 #ifdef NOTHROW
1278                 std::cerr<<"Matrix Read Error: "<<got<<" != "<<exp<<std::endl;
1279                 exit(1);
1280 #else
1281                 throw MatrixReadError<T>(i,N,m,reader.getis(),exp,got);
1282 #endif
1283             }
1284             if (i < M-1 && !reader.readRowEnd(exp,got)) {
1285 #ifdef NOTHROW
1286                 std::cerr<<"Matrix Read Error: "<<got<<" != "<<exp<<std::endl;
1287                 exit(1);
1288 #else
1289                 throw MatrixReadError<T>(i,N,m,reader.getis(),exp,got);
1290 #endif
1291             }
1292         }
1293         if (!reader.readFinal(exp,got)) {
1294 #ifdef NOTHROW
1295             std::cerr<<"Matrix Read Error: "<<got<<" != "<<exp<<std::endl;
1296             exit(1);
1297 #else
1298             throw MatrixReadError<T>(M,0,m,reader.getis(),exp,got);
1299 #endif
1300         }
1301     }
1302 
1303     template <class T, int A>
read(const TMV_Reader & reader)1304     void Matrix<T,A>::read(const TMV_Reader& reader)
1305     {
1306         std::string exp,got;
1307         if (!reader.readCode("M",exp,got)) {
1308 #ifdef NOTHROW
1309             std::cerr<<"Matrix Read Error: "<<got<<" != "<<exp<<std::endl;
1310             exit(1);
1311 #else
1312             throw MatrixReadError<T>(reader.getis(),exp,got);
1313 #endif
1314         }
1315         ptrdiff_t cs=colsize(), rs=rowsize();
1316         if (!reader.readSize(cs,exp,got) ||
1317             !reader.readSize(rs,exp,got)) {
1318 #ifdef NOTHROW
1319             std::cerr<<"Matrix Read Error: reading size\n";
1320             exit(1);
1321 #else
1322             throw MatrixReadError<T>(reader.getis(),exp,got);
1323 #endif
1324         }
1325         if (cs != colsize() || rs != rowsize()) resize(cs,rs);
1326         MatrixView<T> v = view();
1327         FinishRead(reader,v);
1328     }
1329 
1330     template <class T, int A>
read(const TMV_Reader & reader)1331     void MatrixView<T,A>::read(const TMV_Reader& reader)
1332     {
1333         std::string exp,got;
1334         if (!reader.readCode("M",exp,got)) {
1335 #ifdef NOTHROW
1336             std::cerr<<"Matrix Read Error: "<<got<<" != "<<exp<<std::endl;
1337             exit(1);
1338 #else
1339             throw MatrixReadError<T>(reader.getis(),exp,got);
1340 #endif
1341         }
1342         ptrdiff_t cs=colsize(), rs=rowsize();
1343         if (!reader.readSize(cs,exp,got) ||
1344             !reader.readSize(rs,exp,got)) {
1345 #ifdef NOTHROW
1346             std::cerr<<"Matrix Read Error: reading size\n";
1347             exit(1);
1348 #else
1349             throw MatrixReadError<T>(reader.getis(),exp,got);
1350 #endif
1351         }
1352         if (cs != colsize() || rs != rowsize()) {
1353 #ifdef NOTHROW
1354             std::cerr<<"Matrix Read Error: wrong size\n";
1355             exit(1);
1356 #else
1357             throw MatrixReadError<T>(*this,reader.getis(),cs,rs);
1358 #endif
1359         }
1360         FinishRead(reader,*this);
1361     }
1362 
1363 #undef RT
1364 
1365 #define InstFile "TMV_Matrix.inst"
1366 #include "TMV_Inst.h"
1367 #undef InstFile
1368 
1369 } // namespace tmv
1370 
1371 
1372