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_QRD.h"
26 #include "TMV_QRDiv.h"
27 #include "tmv/TMV_Matrix.h"
28 #include "tmv/TMV_TriMatrix.h"
29 #include "tmv/TMV_DiagMatrix.h"
30 #include "tmv/TMV_TriMatrixArith.h"
31 #include "tmv/TMV_MatrixArith.h"
32 #include "tmv/TMV_PackedQ.h"
33 #include <ostream>
34 
35 namespace tmv {
36 
37 #define RT TMV_RealType(T)
38 
39     template <class T>
40     struct QRDiv<T>::QRDiv_Impl
41     {
42     public :
43         QRDiv_Impl(const GenMatrix<T>& A, bool inplace);
44 
45         const bool istrans;
46         const bool inplace;
47         AlignedArray<T> Aptr1;
48         T* Aptr;
49         MatrixView<T> QRx;
50         Vector<T> beta;
51         mutable RT logdet;
52         mutable T signdet;
53         mutable bool donedet;
54     };
55 
56 #define APTR1_SIZE (A.colsize()*A.rowsize())
57 #define APTR1 (inplace ? 0 : APTR1_SIZE)
58 #define APTR (inplace ? A.nonConst().ptr() : Aptr1.get())
59 #define QRX (istrans ? \
60              (inplace ? A.nonConst().transpose() : \
61               MatrixViewOf(Aptr,A.rowsize(),A.colsize(),ColMajor)) : \
62              (inplace ? A.nonConst().view() : \
63               MatrixViewOf(Aptr,A.colsize(),A.rowsize(),ColMajor)))
64 
65     template <class T>
QRDiv_Impl(const GenMatrix<T> & A,bool _inplace)66     QRDiv<T>::QRDiv_Impl::QRDiv_Impl(const GenMatrix<T>& A, bool _inplace) :
67         istrans(A.colsize()<A.rowsize()),
68         inplace(_inplace && (A.iscm() || A.isrm())),
69         Aptr1(APTR1), Aptr(APTR), QRx(QRX),
70         beta(QRx.rowsize()), logdet(0), signdet(1), donedet(false)
71     {
72 #ifdef LAP
73         // Otherwise I get "Conditional jump or move" errors in valgrind,
74         // although they don't seem to cause any problems somehow.
75         if (!inplace) VectorViewOf(Aptr,APTR1_SIZE).setZero();
76 #endif
77     }
78 
79 #undef QRX
80 #undef APTR
81 
82     template <class T>
QRDiv(const GenMatrix<T> & A,bool inplace)83     QRDiv<T>::QRDiv(const GenMatrix<T>& A, bool inplace) :
84         pimpl(new QRDiv_Impl(A,inplace))
85     {
86         if (pimpl->istrans) {
87             if (inplace) TMVAssert(A.transpose() == pimpl->QRx);
88             else pimpl->QRx = A.transpose();
89         }
90         else {
91             if (inplace) TMVAssert(A == pimpl->QRx);
92             else pimpl->QRx = A;
93         }
94         QR_Decompose(pimpl->QRx,pimpl->beta.view(),pimpl->signdet);
95     }
96 
~QRDiv()97     template <class T> QRDiv<T>::~QRDiv() {}
98 
99     template <class T> template <class T1>
doLDivEq(MatrixView<T1> m) const100     void QRDiv<T>::doLDivEq(MatrixView<T1> m) const
101     {
102         if (pimpl->istrans)
103             QR_LDivEq(pimpl->QRx,pimpl->beta,0,m.transpose(),
104                      pimpl->QRx.rowsize());
105         else
106             QR_LDivEq(pimpl->QRx,pimpl->beta,0,m,pimpl->QRx.rowsize());
107     }
108 
109     template <class T> template <class T1>
doRDivEq(MatrixView<T1> m) const110     void QRDiv<T>::doRDivEq(MatrixView<T1> m) const
111     {
112         if (pimpl->istrans)
113             QR_RDivEq(pimpl->QRx,pimpl->beta,0,m.transpose(),
114                      pimpl->QRx.rowsize());
115         else
116             QR_RDivEq(pimpl->QRx,pimpl->beta,0,m,pimpl->QRx.rowsize());
117     }
118 
119     template <class T> template <class T1, class T2>
doLDiv(const GenMatrix<T1> & m,MatrixView<T2> x) const120     void QRDiv<T>::doLDiv(const GenMatrix<T1>& m, MatrixView<T2> x) const
121     {
122         TMVAssert(m.colsize() == (pimpl->istrans ?  pimpl->QRx.rowsize() :
123                                   pimpl->QRx.colsize()));
124         TMVAssert(x.colsize() == (pimpl->istrans ?  pimpl->QRx.colsize() :
125                                   pimpl->QRx.rowsize()));
126         if (pimpl->istrans)
127             QR_RDiv(pimpl->QRx,pimpl->beta,0,m.transpose(),x.transpose(),
128                     pimpl->QRx.rowsize());
129         else QR_LDiv(pimpl->QRx,pimpl->beta,0,m,x,pimpl->QRx.rowsize());
130     }
131 
132     template <class T> template <class T1, class T2>
doRDiv(const GenMatrix<T1> & m,MatrixView<T2> x) const133     void QRDiv<T>::doRDiv(const GenMatrix<T1>& m, MatrixView<T2> x) const
134     {
135         TMVAssert(m.colsize() == x.colsize());
136         TMVAssert(m.rowsize() == (pimpl->istrans ? pimpl->QRx.colsize() :
137                                   pimpl->QRx.rowsize()));
138         TMVAssert(x.rowsize() == (pimpl->istrans ? pimpl->QRx.rowsize() :
139                                   pimpl->QRx.colsize()));
140         if (pimpl->istrans)
141             QR_LDiv(pimpl->QRx,pimpl->beta,0,m.transpose(),x.transpose(),
142                     pimpl->QRx.rowsize());
143         else
144             QR_RDiv(pimpl->QRx,pimpl->beta,0,m,x,pimpl->QRx.rowsize());
145     }
146 
147     template <class T>
det() const148     T QRDiv<T>::det() const
149     {
150         if (!pimpl->donedet) {
151             T s;
152             pimpl->logdet = DiagMatrixViewOf(pimpl->QRx.diag()).logDet(&s);
153             pimpl->signdet *= s;
154             pimpl->donedet = true;
155         }
156         if (pimpl->signdet == T(0)) return T(0);
157         else return pimpl->signdet * TMV_EXP(pimpl->logdet);
158     }
159 
160     template <class T>
logDet(T * sign) const161     RT QRDiv<T>::logDet(T* sign) const
162     {
163         if (!pimpl->donedet) {
164             T s;
165             pimpl->logdet = DiagMatrixViewOf(pimpl->QRx.diag()).logDet(&s);
166             pimpl->signdet *= s;
167             pimpl->donedet = true;
168         }
169         if (sign) *sign = pimpl->signdet;
170         return pimpl->logdet;
171     }
172 
173     template <class T> template <class T1>
doMakeInverse(MatrixView<T1> minv) const174     void QRDiv<T>::doMakeInverse(MatrixView<T1> minv) const
175     {
176         if (pimpl->istrans)
177             QR_Inverse(pimpl->QRx,pimpl->beta,0,minv.transpose(),
178                        pimpl->QRx.rowsize());
179         else
180             QR_Inverse(pimpl->QRx,pimpl->beta,0,minv,pimpl->QRx.rowsize());
181     }
182 
183     template <class T>
doMakeInverseATA(MatrixView<T> ata) const184     void QRDiv<T>::doMakeInverseATA(MatrixView<T> ata) const
185     {
186         // At A = Rt Qt Q R = Rt R
187         // (At A)^-1 = (Rt R)^-1 = R^-1 * Rt^-1
188 
189         UpperTriMatrixView<T> rinv = ata.upperTri();
190         rinv = pimpl->QRx.upperTri().inverse();
191         ata = rinv * rinv.adjoint();
192     }
193 
194     template <class T>
isSingular() const195     bool QRDiv<T>::isSingular() const
196     {
197         return pimpl->QRx.diag().minAbs2Element() <=
198             TMV_Epsilon<T>() * pimpl->QRx.diag().maxAbs2Element();
199     }
200 
201     template <class T>
isTrans() const202     bool QRDiv<T>::isTrans() const
203     { return pimpl->istrans; }
204 
205     template <class T>
getQ() const206     PackedQ<T> QRDiv<T>::getQ() const
207     { return PackedQ<T>(pimpl->QRx,pimpl->beta); }
208 
209     template <class T>
getR() const210     ConstUpperTriMatrixView<T> QRDiv<T>::getR() const
211     { return pimpl->QRx.upperTri(); }
212 
213     template <class T>
getQR() const214     const GenMatrix<T>& QRDiv<T>::getQR() const
215     { return pimpl->QRx; }
216 
getBeta() const217     template <class T> const GenVector<T>& QRDiv<T>::getBeta() const
218     { return pimpl->beta; }
219 
220     template <class T>
checkDecomp(const BaseMatrix<T> & m,std::ostream * fout) const221     bool QRDiv<T>::checkDecomp(const BaseMatrix<T>& m, std::ostream* fout) const
222     {
223         Matrix<T> mm = m;
224         bool printmat = fout && m.colsize() < 100 && m.rowsize() < 100;
225         if (printmat) {
226             *fout << "QRDiv:\n";
227             *fout << "M = "<<
228                 (pimpl->istrans?mm.transpose():mm.view())<<std::endl;
229             *fout << "Q = "<<getQ()<<std::endl;
230             *fout << "R = "<<getR()<<std::endl;
231         }
232         Matrix<T> qr = getQ()*getR();
233         RT nm = Norm(qr-(pimpl->istrans ? mm.transpose() : mm.view()));
234         nm /= Norm(getQ())*Norm(getR());
235         if (printmat) {
236             *fout << "QR = "<<qr<<std::endl;
237         }
238         RT kappa = mm.doCondition();
239         if (fout) {
240             *fout << "Norm(M-QR)/Norm(QR) = "<<nm<<" <? ";
241             *fout << kappa<<" * "<<RT(mm.colsize())<<" * "<<TMV_Epsilon<T>();
242             *fout <<" = "<<kappa*RT(mm.colsize())*TMV_Epsilon<T>()<<std::endl;
243 
244         }
245         return nm < kappa*RT(mm.colsize())*TMV_Epsilon<T>();
246     }
247 
248     template <class T>
colsize() const249     ptrdiff_t QRDiv<T>::colsize() const
250     { return pimpl->istrans ? pimpl->QRx.rowsize() : pimpl->QRx.colsize(); }
251 
252     template <class T>
rowsize() const253     ptrdiff_t QRDiv<T>::rowsize() const
254     { return pimpl->istrans ? pimpl->QRx.colsize() : pimpl->QRx.rowsize(); }
255 
256 #ifdef INST_INT
257 #undef INST_INT
258 #endif
259 
260 #define InstFile "TMV_QRD.inst"
261 #include "TMV_Inst.h"
262 #undef InstFile
263 
264 } // namespace tmv
265 
266 
267