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