1 // SPDX-License-Identifier: Apache-2.0
2 //
3 // Copyright 2008-2016 Conrad Sanderson (http://conradsanderson.id.au)
4 // Copyright 2008-2016 National ICT Australia (NICTA)
5 //
6 // Licensed under the Apache License, Version 2.0 (the "License");
7 // you may not use this file except in compliance with the License.
8 // You may obtain a copy of the License at
9 // http://www.apache.org/licenses/LICENSE-2.0
10 //
11 // Unless required by applicable law or agreed to in writing, software
12 // distributed under the License is distributed on an "AS IS" BASIS,
13 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 // See the License for the specific language governing permissions and
15 // limitations under the License.
16 // ------------------------------------------------------------------------
17
18
19 //! \addtogroup fn_qr
20 //! @{
21
22
23
24 //! QR decomposition
25 template<typename T1>
26 inline
27 bool
qr(Mat<typename T1::elem_type> & Q,Mat<typename T1::elem_type> & R,const Base<typename T1::elem_type,T1> & X,const typename arma_blas_type_only<typename T1::elem_type>::result * junk=nullptr)28 qr
29 (
30 Mat<typename T1::elem_type>& Q,
31 Mat<typename T1::elem_type>& R,
32 const Base<typename T1::elem_type,T1>& X,
33 const typename arma_blas_type_only<typename T1::elem_type>::result* junk = nullptr
34 )
35 {
36 arma_extra_debug_sigprint();
37 arma_ignore(junk);
38
39 arma_debug_check( (&Q == &R), "qr(): Q and R are the same object" );
40
41 const bool status = auxlib::qr(Q, R, X);
42
43 if(status == false)
44 {
45 Q.soft_reset();
46 R.soft_reset();
47 arma_debug_warn_level(3, "qr(): decomposition failed");
48 }
49
50 return status;
51 }
52
53
54
55 //! economical QR decomposition
56 template<typename T1>
57 inline
58 bool
qr_econ(Mat<typename T1::elem_type> & Q,Mat<typename T1::elem_type> & R,const Base<typename T1::elem_type,T1> & X,const typename arma_blas_type_only<typename T1::elem_type>::result * junk=nullptr)59 qr_econ
60 (
61 Mat<typename T1::elem_type>& Q,
62 Mat<typename T1::elem_type>& R,
63 const Base<typename T1::elem_type,T1>& X,
64 const typename arma_blas_type_only<typename T1::elem_type>::result* junk = nullptr
65 )
66 {
67 arma_extra_debug_sigprint();
68 arma_ignore(junk);
69
70 arma_debug_check( (&Q == &R), "qr_econ(): Q and R are the same object" );
71
72 const bool status = auxlib::qr_econ(Q, R, X);
73
74 if(status == false)
75 {
76 Q.soft_reset();
77 R.soft_reset();
78 arma_debug_warn_level(3, "qr_econ(): decomposition failed");
79 }
80
81 return status;
82 }
83
84
85
86 //! QR decomposition with pivoting
87 template<typename T1>
88 inline
89 typename enable_if2< is_supported_blas_type<typename T1::elem_type>::value, bool >::result
qr(Mat<typename T1::elem_type> & Q,Mat<typename T1::elem_type> & R,Mat<uword> & P,const Base<typename T1::elem_type,T1> & X,const char * P_mode="matrix")90 qr
91 (
92 Mat<typename T1::elem_type>& Q,
93 Mat<typename T1::elem_type>& R,
94 Mat<uword>& P,
95 const Base<typename T1::elem_type,T1>& X,
96 const char* P_mode = "matrix"
97 )
98 {
99 arma_extra_debug_sigprint();
100
101 arma_debug_check( (&Q == &R), "qr(): Q and R are the same object" );
102
103 const char sig = (P_mode != nullptr) ? P_mode[0] : char(0);
104
105 arma_debug_check( ((sig != 'm') && (sig != 'v')), "qr(): argument 'P_mode' must be \"vector\" or \"matrix\"" );
106
107 bool status = false;
108
109 if(sig == 'v')
110 {
111 status = auxlib::qr_pivot(Q, R, P, X);
112 }
113 else
114 if(sig == 'm')
115 {
116 Mat<uword> P_vec;
117
118 status = auxlib::qr_pivot(Q, R, P_vec, X);
119
120 if(status)
121 {
122 // construct P
123
124 const uword N = P_vec.n_rows;
125
126 P.zeros(N,N);
127
128 for(uword row=0; row < N; ++row) { P.at(P_vec[row], row) = uword(1); }
129 }
130 }
131
132 if(status == false)
133 {
134 Q.soft_reset();
135 R.soft_reset();
136 P.soft_reset();
137 arma_debug_warn_level(3, "qr(): decomposition failed");
138 }
139
140 return status;
141 }
142
143
144
145 //! @}
146