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/TMV_BaseMatrix.h" 25 #include "tmv/TMV_Vector.h" 26 #include "tmv/TMV_Matrix.h" 27 #include "tmv/TMV_Divider.h" 28 29 namespace tmv { 30 31 template <class T> DivHelper()32 DivHelper<T>::DivHelper() : divider(), divtype(tmv::XX) {} 33 34 template <class T> ~DivHelper()35 DivHelper<T>::~DivHelper() {} 36 37 template <class T> doDet() const38 T DivHelper<T>::doDet() const 39 { 40 TMVAssert(colsize() == rowsize()); 41 setDiv(); 42 T det = divider->det(); 43 doneDiv(); 44 return det; 45 } 46 47 template <class T> TMV_RealType(T)48 TMV_RealType(T) DivHelper<T>::doLogDet(T* sign) const 49 { 50 TMVAssert(colsize() == rowsize()); 51 setDiv(); 52 TMV_RealType(T) logdet = divider->logDet(sign); 53 doneDiv(); 54 return logdet; 55 } 56 57 template <class T, class T1> DoMakeInverse1(const Divider<T> & div,MatrixView<T1> minv)58 static inline void DoMakeInverse1( 59 const Divider<T>& div, MatrixView<T1> minv) 60 { div.makeInverse(minv); } 61 62 template <class T> template <class T1> doMakeInverse(MatrixView<T1> minv) const63 void DivHelper<T>::doMakeInverse(MatrixView<T1> minv) const 64 { 65 TMVAssert(minv.colsize() == rowsize()); 66 TMVAssert(minv.rowsize() == colsize()); 67 setDiv(); 68 DoMakeInverse1(*divider,minv); 69 doneDiv(); 70 } 71 72 template <class T> doMakeInverseATA(MatrixView<T> minv) const73 void DivHelper<T>::doMakeInverseATA(MatrixView<T> minv) const 74 { 75 TMVAssert(minv.colsize() == TMV_MIN(rowsize(),colsize())); 76 TMVAssert(minv.rowsize() == TMV_MIN(rowsize(),colsize())); 77 setDiv(); 78 divider->makeInverseATA(minv); 79 doneDiv(); 80 } 81 82 template <class T> doIsSingular() const83 bool DivHelper<T>::doIsSingular() const 84 { 85 setDiv(); 86 bool s = divider->isSingular(); 87 doneDiv(); 88 return s; 89 } 90 91 template <class T> TMV_RealType(T)92 TMV_RealType(T) DivHelper<T>::doNorm2() const 93 { 94 TMVAssert(divIsSet()); 95 TMVAssert(divtype == SV); 96 return divider->norm2(); 97 } 98 99 template <class T> TMV_RealType(T)100 TMV_RealType(T) DivHelper<T>::doCondition() const 101 { 102 TMVAssert(divIsSet()); 103 TMVAssert(divtype == SV); 104 return divider->condition(); 105 } 106 107 template <class T, class T1> DoRDivEq1(const Divider<T> & div,MatrixView<T1> m0)108 static inline void DoRDivEq1( 109 const Divider<T>& div, MatrixView<T1> m0) 110 { div.RDivEq(m0); } 111 112 template <class T, class T1> DoLDivEq1(const Divider<T> & div,MatrixView<T1> m0)113 static inline void DoLDivEq1( 114 const Divider<T>& div, MatrixView<T1> m0) 115 { div.LDivEq(m0); } 116 117 template <class T, class T1, class T0> DoLDiv1(const Divider<T> & div,const GenMatrix<T1> & m1,MatrixView<T0> m0)118 static inline void DoLDiv1( 119 const Divider<T>& div, 120 const GenMatrix<T1>& m1, MatrixView<T0> m0) 121 { div.LDiv(m1,m0); } 122 123 template <class T, class T1, class T0> DoRDiv1(const Divider<T> & div,const GenMatrix<T1> & m1,MatrixView<T0> m0)124 static inline void DoRDiv1( 125 const Divider<T>& div, 126 const GenMatrix<T1>& m1, MatrixView<T0> m0) 127 { div.RDiv(m1,m0); } 128 129 template <class T> template <class T1> doLDivEq(VectorView<T1> v) const130 void DivHelper<T>::doLDivEq(VectorView<T1> v) const 131 { 132 TMVAssert(colsize() == rowsize()); 133 TMVAssert(colsize() == v.size()); 134 setDiv(); 135 DoLDivEq1(*divider,ColVectorViewOf(v)); 136 doneDiv(); 137 } 138 139 template <class T> template <class T1> doLDivEq(MatrixView<T1> m) const140 void DivHelper<T>::doLDivEq(MatrixView<T1> m) const 141 { 142 TMVAssert(colsize() == rowsize()); 143 TMVAssert(colsize() == m.colsize()); 144 setDiv(); 145 DoLDivEq1(*divider,m); 146 doneDiv(); 147 } 148 149 template <class T> template <class T1> doRDivEq(VectorView<T1> v) const150 void DivHelper<T>::doRDivEq(VectorView<T1> v) const 151 { 152 TMVAssert(colsize() == rowsize()); 153 TMVAssert(colsize() == v.size()); 154 setDiv(); 155 DoRDivEq1(*divider,RowVectorViewOf(v)); 156 doneDiv(); 157 } 158 159 template <class T> template <class T1> doRDivEq(MatrixView<T1> m) const160 void DivHelper<T>::doRDivEq(MatrixView<T1> m) const 161 { 162 TMVAssert(colsize() == rowsize()); 163 TMVAssert(colsize() == m.rowsize()); 164 setDiv(); 165 DoRDivEq1(*divider,m); 166 doneDiv(); 167 } 168 169 template <class T> template <class T1, class T0> doLDiv(const GenVector<T1> & v1,VectorView<T0> v0) const170 void DivHelper<T>::doLDiv( 171 const GenVector<T1>& v1, VectorView<T0> v0) const 172 { 173 TMVAssert(rowsize() == v0.size()); 174 TMVAssert(colsize() == v1.size()); 175 setDiv(); 176 DoLDiv1(*divider,ColVectorViewOf(v1),ColVectorViewOf(v0)); 177 doneDiv(); 178 } 179 180 template <class T> template <class T1, class T0> doLDiv(const GenMatrix<T1> & m1,MatrixView<T0> m0) const181 void DivHelper<T>::doLDiv( 182 const GenMatrix<T1>& m1, MatrixView<T0> m0) const 183 { 184 TMVAssert(rowsize() == m0.colsize()); 185 TMVAssert(colsize() == m1.colsize()); 186 TMVAssert(m1.rowsize() == m0.rowsize()); 187 setDiv(); 188 DoLDiv1(*divider,m1,m0); 189 doneDiv(); 190 } 191 192 template <class T> template <class T1, class T0> doRDiv(const GenVector<T1> & v1,VectorView<T0> v0) const193 void DivHelper<T>::doRDiv( 194 const GenVector<T1>& v1, VectorView<T0> v0) const 195 { 196 TMVAssert(rowsize() == v1.size()); 197 TMVAssert(colsize() == v0.size()); 198 setDiv(); 199 DoRDiv1(*divider,RowVectorViewOf(v1),RowVectorViewOf(v0)); 200 doneDiv(); 201 } 202 203 template <class T> template <class T1, class T0> doRDiv(const GenMatrix<T1> & m1,MatrixView<T0> m0) const204 void DivHelper<T>::doRDiv( 205 const GenMatrix<T1>& m1, MatrixView<T0> m0) const 206 { 207 TMVAssert(rowsize() == m1.rowsize()); 208 TMVAssert(colsize() == m0.rowsize()); 209 TMVAssert(m1.colsize() == m0.colsize()); 210 setDiv(); 211 DoRDiv1(*divider,m1,m0); 212 doneDiv(); 213 } 214 215 template <class T> divideUsing(DivType dt) const216 void DivHelper<T>::divideUsing(DivType dt) const 217 { 218 if (!(divtype & dt)) { 219 unsetDiv(); 220 divtype &= ~tmv::DivTypeFlags; 221 divtype |= dt; 222 } 223 } 224 225 template <class T> getDivType() const226 DivType DivHelper<T>::getDivType() const 227 { 228 if ((divtype & tmv::DivTypeFlags) == tmv::XX) resetDivType(); 229 return divtype & tmv::DivTypeFlags; 230 } 231 232 template <class T> resetDivType() const233 void DivHelper<T>::resetDivType() const 234 { divideUsing(getMatrix().isSquare() ? tmv::LU : tmv::QR); } 235 236 template <class T> divideInPlace() const237 void DivHelper<T>::divideInPlace() const 238 { divtype |= tmv::DivInPlaceFlag; saveDiv(); } 239 240 template <class T> dontDivideInPlace() const241 void DivHelper<T>::dontDivideInPlace() const 242 { divtype &= ~tmv::DivInPlaceFlag; } 243 244 template <class T> divIsInPlace() const245 bool DivHelper<T>::divIsInPlace() const 246 { return divtype & tmv::DivInPlaceFlag; } 247 248 template <class T> saveDiv() const249 void DivHelper<T>::saveDiv() const 250 { divtype |= tmv::SaveDivFlag; } 251 252 template <class T> dontSaveDiv() const253 void DivHelper<T>::dontSaveDiv() const 254 { divtype &= ~tmv::SaveDivFlag; } 255 256 template <class T> divIsSaved() const257 bool DivHelper<T>::divIsSaved() const 258 { return divtype & tmv::SaveDivFlag; } 259 260 template <class T> unsetDiv() const261 void DivHelper<T>::unsetDiv() const 262 { divider.reset(0); } 263 264 template <class T> resetDiv() const265 void DivHelper<T>::resetDiv() const 266 { unsetDiv(); setDiv(); } 267 268 template <class T> divIsSet() const269 bool DivHelper<T>::divIsSet() const 270 { return getDiv(); } 271 272 template <class T> doneDiv() const273 void DivHelper<T>::doneDiv() const 274 { if (!divIsSaved()) unsetDiv(); } 275 276 template <class T> getDiv() const277 const Divider<T>* DivHelper<T>::getDiv() const 278 { return divider.get(); } 279 280 template <class T> checkDecomp(std::ostream * fout) const281 bool DivHelper<T>::checkDecomp(std::ostream* fout) const 282 { 283 TMVAssert(divider.get()); 284 return divider->checkDecomp(getMatrix(),fout); 285 } 286 287 template <class T> checkDecomp(const BaseMatrix<T> & m2,std::ostream * fout) const288 bool DivHelper<T>::checkDecomp( 289 const BaseMatrix<T>& m2, std::ostream* fout) const 290 { 291 TMVAssert(divider.get()); 292 return divider->checkDecomp(m2,fout); 293 } 294 295 #define InstFile "TMV_BaseMatrix.inst" 296 #include "TMV_Inst.h" 297 #undef InstFile 298 299 } // namespace tmv 300 301 302