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