auxlib_meat.hpp
Go to the documentation of this file.
1 // Copyright (C) 2008-2011 NICTA (www.nicta.com.au)
2 // Copyright (C) 2008-2011 Conrad Sanderson
3 // Copyright (C) 2009 Edmund Highcock
4 // Copyright (C) 2011 James Sanders
5 // Copyright (C) 2011 Stanislav Funiak
6 //
7 // This file is part of the Armadillo C++ library.
8 // It is provided without any warranty of fitness
9 // for any purpose. You can redistribute this file
10 // and/or modify it under the terms of the GNU
11 // Lesser General Public License (LGPL) as published
12 // by the Free Software Foundation, either version 3
13 // of the License or (at your option) any later version.
14 // (see http://www.opensource.org/licenses for more info)
15 
16 
19 
20 
21 
23 template<typename eT, typename T1>
24 inline
25 bool
26 auxlib::inv(Mat<eT>& out, const Base<eT,T1>& X, const bool slow)
27  {
28  arma_extra_debug_sigprint();
29 
30  out = X.get_ref();
31 
32  arma_debug_check( (out.is_square() == false), "inv(): given matrix is not square" );
33 
34  bool status = false;
35 
36  const uword N = out.n_rows;
37 
38  if( (N <= 4) && (slow == false) )
39  {
40  status = auxlib::inv_inplace_tinymat(out, N);
41  }
42 
43  if( (N > 4) || (status == false) )
44  {
45  status = auxlib::inv_inplace_lapack(out);
46  }
47 
48  return status;
49  }
50 
51 
52 
53 template<typename eT>
54 inline
55 bool
56 auxlib::inv(Mat<eT>& out, const Mat<eT>& X, const bool slow)
57  {
58  arma_extra_debug_sigprint();
59 
60  arma_debug_check( (X.is_square() == false), "inv(): given matrix is not square" );
61 
62  bool status = false;
63 
64  const uword N = X.n_rows;
65 
66  if( (N <= 4) && (slow == false) )
67  {
68  status = (&out != &X) ? auxlib::inv_noalias_tinymat(out, X, N) : auxlib::inv_inplace_tinymat(out, N);
69  }
70 
71  if( (N > 4) || (status == false) )
72  {
73  out = X;
74  status = auxlib::inv_inplace_lapack(out);
75  }
76 
77  return status;
78  }
79 
80 
81 
82 template<typename eT>
83 inline
84 bool
85 auxlib::inv_noalias_tinymat(Mat<eT>& out, const Mat<eT>& X, const uword N)
86  {
87  arma_extra_debug_sigprint();
88 
89  bool det_ok = true;
90 
91  out.set_size(N,N);
92 
93  switch(N)
94  {
95  case 1:
96  {
97  out[0] = eT(1) / X[0];
98  };
99  break;
100 
101  case 2:
102  {
103  const eT* Xm = X.memptr();
104 
105  const eT a = Xm[pos<0,0>::n2];
106  const eT b = Xm[pos<0,1>::n2];
107  const eT c = Xm[pos<1,0>::n2];
108  const eT d = Xm[pos<1,1>::n2];
109 
110  const eT tmp_det = (a*d - b*c);
111 
112  if(tmp_det != eT(0))
113  {
114  eT* outm = out.memptr();
115 
116  outm[pos<0,0>::n2] = d / tmp_det;
117  outm[pos<0,1>::n2] = -b / tmp_det;
118  outm[pos<1,0>::n2] = -c / tmp_det;
119  outm[pos<1,1>::n2] = a / tmp_det;
120  }
121  else
122  {
123  det_ok = false;
124  }
125  };
126  break;
127 
128  case 3:
129  {
130  const eT* X_col0 = X.colptr(0);
131  const eT a11 = X_col0[0];
132  const eT a21 = X_col0[1];
133  const eT a31 = X_col0[2];
134 
135  const eT* X_col1 = X.colptr(1);
136  const eT a12 = X_col1[0];
137  const eT a22 = X_col1[1];
138  const eT a32 = X_col1[2];
139 
140  const eT* X_col2 = X.colptr(2);
141  const eT a13 = X_col2[0];
142  const eT a23 = X_col2[1];
143  const eT a33 = X_col2[2];
144 
145  const eT tmp_det = a11*(a33*a22 - a32*a23) - a21*(a33*a12-a32*a13) + a31*(a23*a12 - a22*a13);
146 
147  if(tmp_det != eT(0))
148  {
149  eT* out_col0 = out.colptr(0);
150  out_col0[0] = (a33*a22 - a32*a23) / tmp_det;
151  out_col0[1] = -(a33*a21 - a31*a23) / tmp_det;
152  out_col0[2] = (a32*a21 - a31*a22) / tmp_det;
153 
154  eT* out_col1 = out.colptr(1);
155  out_col1[0] = -(a33*a12 - a32*a13) / tmp_det;
156  out_col1[1] = (a33*a11 - a31*a13) / tmp_det;
157  out_col1[2] = -(a32*a11 - a31*a12) / tmp_det;
158 
159  eT* out_col2 = out.colptr(2);
160  out_col2[0] = (a23*a12 - a22*a13) / tmp_det;
161  out_col2[1] = -(a23*a11 - a21*a13) / tmp_det;
162  out_col2[2] = (a22*a11 - a21*a12) / tmp_det;
163  }
164  else
165  {
166  det_ok = false;
167  }
168  };
169  break;
170 
171  case 4:
172  {
173  const eT tmp_det = det(X);
174 
175  if(tmp_det != eT(0))
176  {
177  const eT* Xm = X.memptr();
178  eT* outm = out.memptr();
179 
180  outm[pos<0,0>::n4] = ( Xm[pos<1,2>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,1>::n4] - Xm[pos<1,3>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,1>::n4] + Xm[pos<1,3>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,2>::n4] - Xm[pos<1,1>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,2>::n4] - Xm[pos<1,2>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,3>::n4] + Xm[pos<1,1>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,3>::n4] ) / tmp_det;
181  outm[pos<1,0>::n4] = ( Xm[pos<1,3>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,0>::n4] - Xm[pos<1,2>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,0>::n4] - Xm[pos<1,3>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,2>::n4] + Xm[pos<1,0>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,2>::n4] + Xm[pos<1,2>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,3>::n4] - Xm[pos<1,0>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,3>::n4] ) / tmp_det;
182  outm[pos<2,0>::n4] = ( Xm[pos<1,1>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,0>::n4] - Xm[pos<1,3>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,0>::n4] + Xm[pos<1,3>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,1>::n4] - Xm[pos<1,0>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,1>::n4] - Xm[pos<1,1>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,3>::n4] + Xm[pos<1,0>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,3>::n4] ) / tmp_det;
183  outm[pos<3,0>::n4] = ( Xm[pos<1,2>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,0>::n4] - Xm[pos<1,1>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,0>::n4] - Xm[pos<1,2>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,1>::n4] + Xm[pos<1,0>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,1>::n4] + Xm[pos<1,1>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,2>::n4] - Xm[pos<1,0>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,2>::n4] ) / tmp_det;
184 
185  outm[pos<0,1>::n4] = ( Xm[pos<0,3>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,1>::n4] - Xm[pos<0,2>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,1>::n4] - Xm[pos<0,3>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,2>::n4] + Xm[pos<0,1>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,2>::n4] + Xm[pos<0,2>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,3>::n4] - Xm[pos<0,1>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,3>::n4] ) / tmp_det;
186  outm[pos<1,1>::n4] = ( Xm[pos<0,2>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,0>::n4] - Xm[pos<0,3>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,0>::n4] + Xm[pos<0,3>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,2>::n4] - Xm[pos<0,0>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,2>::n4] - Xm[pos<0,2>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,3>::n4] + Xm[pos<0,0>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,3>::n4] ) / tmp_det;
187  outm[pos<2,1>::n4] = ( Xm[pos<0,3>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,0>::n4] - Xm[pos<0,1>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,0>::n4] - Xm[pos<0,3>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,1>::n4] + Xm[pos<0,0>::n4]*Xm[pos<2,3>::n4]*Xm[pos<3,1>::n4] + Xm[pos<0,1>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,3>::n4] - Xm[pos<0,0>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,3>::n4] ) / tmp_det;
188  outm[pos<3,1>::n4] = ( Xm[pos<0,1>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,0>::n4] - Xm[pos<0,2>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,0>::n4] + Xm[pos<0,2>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,1>::n4] - Xm[pos<0,0>::n4]*Xm[pos<2,2>::n4]*Xm[pos<3,1>::n4] - Xm[pos<0,1>::n4]*Xm[pos<2,0>::n4]*Xm[pos<3,2>::n4] + Xm[pos<0,0>::n4]*Xm[pos<2,1>::n4]*Xm[pos<3,2>::n4] ) / tmp_det;
189 
190  outm[pos<0,2>::n4] = ( Xm[pos<0,2>::n4]*Xm[pos<1,3>::n4]*Xm[pos<3,1>::n4] - Xm[pos<0,3>::n4]*Xm[pos<1,2>::n4]*Xm[pos<3,1>::n4] + Xm[pos<0,3>::n4]*Xm[pos<1,1>::n4]*Xm[pos<3,2>::n4] - Xm[pos<0,1>::n4]*Xm[pos<1,3>::n4]*Xm[pos<3,2>::n4] - Xm[pos<0,2>::n4]*Xm[pos<1,1>::n4]*Xm[pos<3,3>::n4] + Xm[pos<0,1>::n4]*Xm[pos<1,2>::n4]*Xm[pos<3,3>::n4] ) / tmp_det;
191  outm[pos<1,2>::n4] = ( Xm[pos<0,3>::n4]*Xm[pos<1,2>::n4]*Xm[pos<3,0>::n4] - Xm[pos<0,2>::n4]*Xm[pos<1,3>::n4]*Xm[pos<3,0>::n4] - Xm[pos<0,3>::n4]*Xm[pos<1,0>::n4]*Xm[pos<3,2>::n4] + Xm[pos<0,0>::n4]*Xm[pos<1,3>::n4]*Xm[pos<3,2>::n4] + Xm[pos<0,2>::n4]*Xm[pos<1,0>::n4]*Xm[pos<3,3>::n4] - Xm[pos<0,0>::n4]*Xm[pos<1,2>::n4]*Xm[pos<3,3>::n4] ) / tmp_det;
192  outm[pos<2,2>::n4] = ( Xm[pos<0,1>::n4]*Xm[pos<1,3>::n4]*Xm[pos<3,0>::n4] - Xm[pos<0,3>::n4]*Xm[pos<1,1>::n4]*Xm[pos<3,0>::n4] + Xm[pos<0,3>::n4]*Xm[pos<1,0>::n4]*Xm[pos<3,1>::n4] - Xm[pos<0,0>::n4]*Xm[pos<1,3>::n4]*Xm[pos<3,1>::n4] - Xm[pos<0,1>::n4]*Xm[pos<1,0>::n4]*Xm[pos<3,3>::n4] + Xm[pos<0,0>::n4]*Xm[pos<1,1>::n4]*Xm[pos<3,3>::n4] ) / tmp_det;
193  outm[pos<3,2>::n4] = ( Xm[pos<0,2>::n4]*Xm[pos<1,1>::n4]*Xm[pos<3,0>::n4] - Xm[pos<0,1>::n4]*Xm[pos<1,2>::n4]*Xm[pos<3,0>::n4] - Xm[pos<0,2>::n4]*Xm[pos<1,0>::n4]*Xm[pos<3,1>::n4] + Xm[pos<0,0>::n4]*Xm[pos<1,2>::n4]*Xm[pos<3,1>::n4] + Xm[pos<0,1>::n4]*Xm[pos<1,0>::n4]*Xm[pos<3,2>::n4] - Xm[pos<0,0>::n4]*Xm[pos<1,1>::n4]*Xm[pos<3,2>::n4] ) / tmp_det;
194 
195  outm[pos<0,3>::n4] = ( Xm[pos<0,3>::n4]*Xm[pos<1,2>::n4]*Xm[pos<2,1>::n4] - Xm[pos<0,2>::n4]*Xm[pos<1,3>::n4]*Xm[pos<2,1>::n4] - Xm[pos<0,3>::n4]*Xm[pos<1,1>::n4]*Xm[pos<2,2>::n4] + Xm[pos<0,1>::n4]*Xm[pos<1,3>::n4]*Xm[pos<2,2>::n4] + Xm[pos<0,2>::n4]*Xm[pos<1,1>::n4]*Xm[pos<2,3>::n4] - Xm[pos<0,1>::n4]*Xm[pos<1,2>::n4]*Xm[pos<2,3>::n4] ) / tmp_det;
196  outm[pos<1,3>::n4] = ( Xm[pos<0,2>::n4]*Xm[pos<1,3>::n4]*Xm[pos<2,0>::n4] - Xm[pos<0,3>::n4]*Xm[pos<1,2>::n4]*Xm[pos<2,0>::n4] + Xm[pos<0,3>::n4]*Xm[pos<1,0>::n4]*Xm[pos<2,2>::n4] - Xm[pos<0,0>::n4]*Xm[pos<1,3>::n4]*Xm[pos<2,2>::n4] - Xm[pos<0,2>::n4]*Xm[pos<1,0>::n4]*Xm[pos<2,3>::n4] + Xm[pos<0,0>::n4]*Xm[pos<1,2>::n4]*Xm[pos<2,3>::n4] ) / tmp_det;
197  outm[pos<2,3>::n4] = ( Xm[pos<0,3>::n4]*Xm[pos<1,1>::n4]*Xm[pos<2,0>::n4] - Xm[pos<0,1>::n4]*Xm[pos<1,3>::n4]*Xm[pos<2,0>::n4] - Xm[pos<0,3>::n4]*Xm[pos<1,0>::n4]*Xm[pos<2,1>::n4] + Xm[pos<0,0>::n4]*Xm[pos<1,3>::n4]*Xm[pos<2,1>::n4] + Xm[pos<0,1>::n4]*Xm[pos<1,0>::n4]*Xm[pos<2,3>::n4] - Xm[pos<0,0>::n4]*Xm[pos<1,1>::n4]*Xm[pos<2,3>::n4] ) / tmp_det;
198  outm[pos<3,3>::n4] = ( Xm[pos<0,1>::n4]*Xm[pos<1,2>::n4]*Xm[pos<2,0>::n4] - Xm[pos<0,2>::n4]*Xm[pos<1,1>::n4]*Xm[pos<2,0>::n4] + Xm[pos<0,2>::n4]*Xm[pos<1,0>::n4]*Xm[pos<2,1>::n4] - Xm[pos<0,0>::n4]*Xm[pos<1,2>::n4]*Xm[pos<2,1>::n4] - Xm[pos<0,1>::n4]*Xm[pos<1,0>::n4]*Xm[pos<2,2>::n4] + Xm[pos<0,0>::n4]*Xm[pos<1,1>::n4]*Xm[pos<2,2>::n4] ) / tmp_det;
199  }
200  else
201  {
202  det_ok = false;
203  }
204  };
205  break;
206 
207  default:
208  ;
209  }
210 
211  return det_ok;
212  }
213 
214 
215 
216 template<typename eT>
217 inline
218 bool
219 auxlib::inv_inplace_tinymat(Mat<eT>& X, const uword N)
220  {
221  arma_extra_debug_sigprint();
222 
223  bool det_ok = true;
224 
225  // for more info, see:
226  // http://www.dr-lex.34sp.com/random/matrix_inv.html
227  // http://www.cvl.iis.u-tokyo.ac.jp/~miyazaki/tech/teche23.html
228  // http://www.euclideanspace.com/maths/algebra/matrix/functions/inverse/fourD/index.htm
229  // http://www.geometrictools.com//LibFoundation/Mathematics/Wm4Matrix4.inl
230 
231  switch(N)
232  {
233  case 1:
234  {
235  X[0] = eT(1) / X[0];
236  };
237  break;
238 
239  case 2:
240  {
241  const eT a = X[pos<0,0>::n2];
242  const eT b = X[pos<0,1>::n2];
243  const eT c = X[pos<1,0>::n2];
244  const eT d = X[pos<1,1>::n2];
245 
246  const eT tmp_det = (a*d - b*c);
247 
248  if(tmp_det != eT(0))
249  {
250  X[pos<0,0>::n2] = d / tmp_det;
251  X[pos<0,1>::n2] = -b / tmp_det;
252  X[pos<1,0>::n2] = -c / tmp_det;
253  X[pos<1,1>::n2] = a / tmp_det;
254  }
255  else
256  {
257  det_ok = false;
258  }
259  };
260  break;
261 
262  case 3:
263  {
264  eT* X_col0 = X.colptr(0);
265  eT* X_col1 = X.colptr(1);
266  eT* X_col2 = X.colptr(2);
267 
268  const eT a11 = X_col0[0];
269  const eT a21 = X_col0[1];
270  const eT a31 = X_col0[2];
271 
272  const eT a12 = X_col1[0];
273  const eT a22 = X_col1[1];
274  const eT a32 = X_col1[2];
275 
276  const eT a13 = X_col2[0];
277  const eT a23 = X_col2[1];
278  const eT a33 = X_col2[2];
279 
280  const eT tmp_det = a11*(a33*a22 - a32*a23) - a21*(a33*a12-a32*a13) + a31*(a23*a12 - a22*a13);
281 
282  if(tmp_det != eT(0))
283  {
284  X_col0[0] = (a33*a22 - a32*a23) / tmp_det;
285  X_col0[1] = -(a33*a21 - a31*a23) / tmp_det;
286  X_col0[2] = (a32*a21 - a31*a22) / tmp_det;
287 
288  X_col1[0] = -(a33*a12 - a32*a13) / tmp_det;
289  X_col1[1] = (a33*a11 - a31*a13) / tmp_det;
290  X_col1[2] = -(a32*a11 - a31*a12) / tmp_det;
291 
292  X_col2[0] = (a23*a12 - a22*a13) / tmp_det;
293  X_col2[1] = -(a23*a11 - a21*a13) / tmp_det;
294  X_col2[2] = (a22*a11 - a21*a12) / tmp_det;
295  }
296  else
297  {
298  det_ok = false;
299  }
300  };
301  break;
302 
303  case 4:
304  {
305  const eT tmp_det = det(X);
306 
307  if(tmp_det != eT(0))
308  {
309  const Mat<eT> A(X);
310 
311  const eT* Am = A.memptr();
312  eT* Xm = X.memptr();
313 
314  Xm[pos<0,0>::n4] = ( Am[pos<1,2>::n4]*Am[pos<2,3>::n4]*Am[pos<3,1>::n4] - Am[pos<1,3>::n4]*Am[pos<2,2>::n4]*Am[pos<3,1>::n4] + Am[pos<1,3>::n4]*Am[pos<2,1>::n4]*Am[pos<3,2>::n4] - Am[pos<1,1>::n4]*Am[pos<2,3>::n4]*Am[pos<3,2>::n4] - Am[pos<1,2>::n4]*Am[pos<2,1>::n4]*Am[pos<3,3>::n4] + Am[pos<1,1>::n4]*Am[pos<2,2>::n4]*Am[pos<3,3>::n4] ) / tmp_det;
315  Xm[pos<1,0>::n4] = ( Am[pos<1,3>::n4]*Am[pos<2,2>::n4]*Am[pos<3,0>::n4] - Am[pos<1,2>::n4]*Am[pos<2,3>::n4]*Am[pos<3,0>::n4] - Am[pos<1,3>::n4]*Am[pos<2,0>::n4]*Am[pos<3,2>::n4] + Am[pos<1,0>::n4]*Am[pos<2,3>::n4]*Am[pos<3,2>::n4] + Am[pos<1,2>::n4]*Am[pos<2,0>::n4]*Am[pos<3,3>::n4] - Am[pos<1,0>::n4]*Am[pos<2,2>::n4]*Am[pos<3,3>::n4] ) / tmp_det;
316  Xm[pos<2,0>::n4] = ( Am[pos<1,1>::n4]*Am[pos<2,3>::n4]*Am[pos<3,0>::n4] - Am[pos<1,3>::n4]*Am[pos<2,1>::n4]*Am[pos<3,0>::n4] + Am[pos<1,3>::n4]*Am[pos<2,0>::n4]*Am[pos<3,1>::n4] - Am[pos<1,0>::n4]*Am[pos<2,3>::n4]*Am[pos<3,1>::n4] - Am[pos<1,1>::n4]*Am[pos<2,0>::n4]*Am[pos<3,3>::n4] + Am[pos<1,0>::n4]*Am[pos<2,1>::n4]*Am[pos<3,3>::n4] ) / tmp_det;
317  Xm[pos<3,0>::n4] = ( Am[pos<1,2>::n4]*Am[pos<2,1>::n4]*Am[pos<3,0>::n4] - Am[pos<1,1>::n4]*Am[pos<2,2>::n4]*Am[pos<3,0>::n4] - Am[pos<1,2>::n4]*Am[pos<2,0>::n4]*Am[pos<3,1>::n4] + Am[pos<1,0>::n4]*Am[pos<2,2>::n4]*Am[pos<3,1>::n4] + Am[pos<1,1>::n4]*Am[pos<2,0>::n4]*Am[pos<3,2>::n4] - Am[pos<1,0>::n4]*Am[pos<2,1>::n4]*Am[pos<3,2>::n4] ) / tmp_det;
318 
319  Xm[pos<0,1>::n4] = ( Am[pos<0,3>::n4]*Am[pos<2,2>::n4]*Am[pos<3,1>::n4] - Am[pos<0,2>::n4]*Am[pos<2,3>::n4]*Am[pos<3,1>::n4] - Am[pos<0,3>::n4]*Am[pos<2,1>::n4]*Am[pos<3,2>::n4] + Am[pos<0,1>::n4]*Am[pos<2,3>::n4]*Am[pos<3,2>::n4] + Am[pos<0,2>::n4]*Am[pos<2,1>::n4]*Am[pos<3,3>::n4] - Am[pos<0,1>::n4]*Am[pos<2,2>::n4]*Am[pos<3,3>::n4] ) / tmp_det;
320  Xm[pos<1,1>::n4] = ( Am[pos<0,2>::n4]*Am[pos<2,3>::n4]*Am[pos<3,0>::n4] - Am[pos<0,3>::n4]*Am[pos<2,2>::n4]*Am[pos<3,0>::n4] + Am[pos<0,3>::n4]*Am[pos<2,0>::n4]*Am[pos<3,2>::n4] - Am[pos<0,0>::n4]*Am[pos<2,3>::n4]*Am[pos<3,2>::n4] - Am[pos<0,2>::n4]*Am[pos<2,0>::n4]*Am[pos<3,3>::n4] + Am[pos<0,0>::n4]*Am[pos<2,2>::n4]*Am[pos<3,3>::n4] ) / tmp_det;
321  Xm[pos<2,1>::n4] = ( Am[pos<0,3>::n4]*Am[pos<2,1>::n4]*Am[pos<3,0>::n4] - Am[pos<0,1>::n4]*Am[pos<2,3>::n4]*Am[pos<3,0>::n4] - Am[pos<0,3>::n4]*Am[pos<2,0>::n4]*Am[pos<3,1>::n4] + Am[pos<0,0>::n4]*Am[pos<2,3>::n4]*Am[pos<3,1>::n4] + Am[pos<0,1>::n4]*Am[pos<2,0>::n4]*Am[pos<3,3>::n4] - Am[pos<0,0>::n4]*Am[pos<2,1>::n4]*Am[pos<3,3>::n4] ) / tmp_det;
322  Xm[pos<3,1>::n4] = ( Am[pos<0,1>::n4]*Am[pos<2,2>::n4]*Am[pos<3,0>::n4] - Am[pos<0,2>::n4]*Am[pos<2,1>::n4]*Am[pos<3,0>::n4] + Am[pos<0,2>::n4]*Am[pos<2,0>::n4]*Am[pos<3,1>::n4] - Am[pos<0,0>::n4]*Am[pos<2,2>::n4]*Am[pos<3,1>::n4] - Am[pos<0,1>::n4]*Am[pos<2,0>::n4]*Am[pos<3,2>::n4] + Am[pos<0,0>::n4]*Am[pos<2,1>::n4]*Am[pos<3,2>::n4] ) / tmp_det;
323 
324  Xm[pos<0,2>::n4] = ( Am[pos<0,2>::n4]*Am[pos<1,3>::n4]*Am[pos<3,1>::n4] - Am[pos<0,3>::n4]*Am[pos<1,2>::n4]*Am[pos<3,1>::n4] + Am[pos<0,3>::n4]*Am[pos<1,1>::n4]*Am[pos<3,2>::n4] - Am[pos<0,1>::n4]*Am[pos<1,3>::n4]*Am[pos<3,2>::n4] - Am[pos<0,2>::n4]*Am[pos<1,1>::n4]*Am[pos<3,3>::n4] + Am[pos<0,1>::n4]*Am[pos<1,2>::n4]*Am[pos<3,3>::n4] ) / tmp_det;
325  Xm[pos<1,2>::n4] = ( Am[pos<0,3>::n4]*Am[pos<1,2>::n4]*Am[pos<3,0>::n4] - Am[pos<0,2>::n4]*Am[pos<1,3>::n4]*Am[pos<3,0>::n4] - Am[pos<0,3>::n4]*Am[pos<1,0>::n4]*Am[pos<3,2>::n4] + Am[pos<0,0>::n4]*Am[pos<1,3>::n4]*Am[pos<3,2>::n4] + Am[pos<0,2>::n4]*Am[pos<1,0>::n4]*Am[pos<3,3>::n4] - Am[pos<0,0>::n4]*Am[pos<1,2>::n4]*Am[pos<3,3>::n4] ) / tmp_det;
326  Xm[pos<2,2>::n4] = ( Am[pos<0,1>::n4]*Am[pos<1,3>::n4]*Am[pos<3,0>::n4] - Am[pos<0,3>::n4]*Am[pos<1,1>::n4]*Am[pos<3,0>::n4] + Am[pos<0,3>::n4]*Am[pos<1,0>::n4]*Am[pos<3,1>::n4] - Am[pos<0,0>::n4]*Am[pos<1,3>::n4]*Am[pos<3,1>::n4] - Am[pos<0,1>::n4]*Am[pos<1,0>::n4]*Am[pos<3,3>::n4] + Am[pos<0,0>::n4]*Am[pos<1,1>::n4]*Am[pos<3,3>::n4] ) / tmp_det;
327  Xm[pos<3,2>::n4] = ( Am[pos<0,2>::n4]*Am[pos<1,1>::n4]*Am[pos<3,0>::n4] - Am[pos<0,1>::n4]*Am[pos<1,2>::n4]*Am[pos<3,0>::n4] - Am[pos<0,2>::n4]*Am[pos<1,0>::n4]*Am[pos<3,1>::n4] + Am[pos<0,0>::n4]*Am[pos<1,2>::n4]*Am[pos<3,1>::n4] + Am[pos<0,1>::n4]*Am[pos<1,0>::n4]*Am[pos<3,2>::n4] - Am[pos<0,0>::n4]*Am[pos<1,1>::n4]*Am[pos<3,2>::n4] ) / tmp_det;
328 
329  Xm[pos<0,3>::n4] = ( Am[pos<0,3>::n4]*Am[pos<1,2>::n4]*Am[pos<2,1>::n4] - Am[pos<0,2>::n4]*Am[pos<1,3>::n4]*Am[pos<2,1>::n4] - Am[pos<0,3>::n4]*Am[pos<1,1>::n4]*Am[pos<2,2>::n4] + Am[pos<0,1>::n4]*Am[pos<1,3>::n4]*Am[pos<2,2>::n4] + Am[pos<0,2>::n4]*Am[pos<1,1>::n4]*Am[pos<2,3>::n4] - Am[pos<0,1>::n4]*Am[pos<1,2>::n4]*Am[pos<2,3>::n4] ) / tmp_det;
330  Xm[pos<1,3>::n4] = ( Am[pos<0,2>::n4]*Am[pos<1,3>::n4]*Am[pos<2,0>::n4] - Am[pos<0,3>::n4]*Am[pos<1,2>::n4]*Am[pos<2,0>::n4] + Am[pos<0,3>::n4]*Am[pos<1,0>::n4]*Am[pos<2,2>::n4] - Am[pos<0,0>::n4]*Am[pos<1,3>::n4]*Am[pos<2,2>::n4] - Am[pos<0,2>::n4]*Am[pos<1,0>::n4]*Am[pos<2,3>::n4] + Am[pos<0,0>::n4]*Am[pos<1,2>::n4]*Am[pos<2,3>::n4] ) / tmp_det;
331  Xm[pos<2,3>::n4] = ( Am[pos<0,3>::n4]*Am[pos<1,1>::n4]*Am[pos<2,0>::n4] - Am[pos<0,1>::n4]*Am[pos<1,3>::n4]*Am[pos<2,0>::n4] - Am[pos<0,3>::n4]*Am[pos<1,0>::n4]*Am[pos<2,1>::n4] + Am[pos<0,0>::n4]*Am[pos<1,3>::n4]*Am[pos<2,1>::n4] + Am[pos<0,1>::n4]*Am[pos<1,0>::n4]*Am[pos<2,3>::n4] - Am[pos<0,0>::n4]*Am[pos<1,1>::n4]*Am[pos<2,3>::n4] ) / tmp_det;
332  Xm[pos<3,3>::n4] = ( Am[pos<0,1>::n4]*Am[pos<1,2>::n4]*Am[pos<2,0>::n4] - Am[pos<0,2>::n4]*Am[pos<1,1>::n4]*Am[pos<2,0>::n4] + Am[pos<0,2>::n4]*Am[pos<1,0>::n4]*Am[pos<2,1>::n4] - Am[pos<0,0>::n4]*Am[pos<1,2>::n4]*Am[pos<2,1>::n4] - Am[pos<0,1>::n4]*Am[pos<1,0>::n4]*Am[pos<2,2>::n4] + Am[pos<0,0>::n4]*Am[pos<1,1>::n4]*Am[pos<2,2>::n4] ) / tmp_det;
333  }
334  else
335  {
336  det_ok = false;
337  }
338  };
339  break;
340 
341  default:
342  ;
343  }
344 
345  return det_ok;
346  }
347 
348 
349 
350 template<typename eT>
351 inline
352 bool
353 auxlib::inv_inplace_lapack(Mat<eT>& out)
354  {
355  arma_extra_debug_sigprint();
356 
357  if(out.is_empty())
358  {
359  return true;
360  }
361 
362  #if defined(ARMA_USE_ATLAS)
363  {
364  podarray<int> ipiv(out.n_rows);
365 
366  int info = atlas::clapack_getrf(atlas::CblasColMajor, out.n_rows, out.n_cols, out.memptr(), out.n_rows, ipiv.memptr());
367 
368  if(info == 0)
369  {
370  info = atlas::clapack_getri(atlas::CblasColMajor, out.n_rows, out.memptr(), out.n_rows, ipiv.memptr());
371  }
372 
373  return (info == 0);
374  }
375  #elif defined(ARMA_USE_LAPACK)
376  {
377  blas_int n_rows = out.n_rows;
378  blas_int n_cols = out.n_cols;
379  blas_int info = 0;
380 
381  podarray<blas_int> ipiv(out.n_rows);
382 
383  // 84 was empirically found -- it is the maximum value suggested by LAPACK (as provided by ATLAS v3.6)
384  // based on tests with various matrix types on 32-bit and 64-bit machines
385  //
386  // the "work" array is deliberately long so that a secondary (time-consuming)
387  // memory allocation is avoided, if possible
388 
389  blas_int work_len = (std::max)(blas_int(1), n_rows*84);
390  podarray<eT> work( static_cast<uword>(work_len) );
391 
392  lapack::getrf(&n_rows, &n_cols, out.memptr(), &n_rows, ipiv.memptr(), &info);
393 
394  if(info == 0)
395  {
396  // query for optimum size of work_len
397 
398  blas_int work_len_tmp = -1;
399  lapack::getri(&n_rows, out.memptr(), &n_rows, ipiv.memptr(), work.memptr(), &work_len_tmp, &info);
400 
401  if(info == 0)
402  {
403  blas_int proposed_work_len = static_cast<blas_int>(access::tmp_real(work[0]));
404 
405  // if necessary, allocate more memory
406  if(work_len < proposed_work_len)
407  {
408  work_len = proposed_work_len;
409  work.set_size( static_cast<uword>(work_len) );
410  }
411  }
412 
413  lapack::getri(&n_rows, out.memptr(), &n_rows, ipiv.memptr(), work.memptr(), &work_len, &info);
414  }
415 
416  return (info == 0);
417  }
418  #else
419  {
420  arma_ignore(out);
421  arma_stop("inv(): use of ATLAS or LAPACK needs to be enabled");
422  return false;
423  }
424  #endif
425  }
426 
427 
428 
429 template<typename eT, typename T1>
430 inline
431 bool
432 auxlib::inv_tr(Mat<eT>& out, const Base<eT,T1>& X, const uword layout)
433  {
434  arma_extra_debug_sigprint();
435 
436  out = X.get_ref();
437 
438  arma_debug_check( (out.is_square() == false), "inv(): given matrix is not square" );
439 
440  if(out.is_empty())
441  {
442  return true;
443  }
444 
445  bool status;
446 
447  #if defined(ARMA_USE_LAPACK)
448  {
449  char uplo = (layout == 0) ? 'U' : 'L';
450  char diag = 'N';
451  blas_int n = blas_int(out.n_rows);
452  blas_int info = 0;
453 
454  lapack::trtri(&uplo, &diag, &n, out.memptr(), &n, &info);
455 
456  status = (info == 0);
457  }
458  #else
459  {
460  arma_ignore(layout);
461  arma_stop("inv(): use of LAPACK needs to be enabled");
462  status = false;
463  }
464  #endif
465 
466 
467  if(status == true)
468  {
469  if(layout == 0)
470  {
471  // upper triangular
472  out = trimatu(out);
473  }
474  else
475  {
476  // lower triangular
477  out = trimatl(out);
478  }
479  }
480 
481  return status;
482  }
483 
484 
485 
486 template<typename eT, typename T1>
487 inline
488 bool
489 auxlib::inv_sym(Mat<eT>& out, const Base<eT,T1>& X, const uword layout)
490  {
491  arma_extra_debug_sigprint();
492 
493  out = X.get_ref();
494 
495  arma_debug_check( (out.is_square() == false), "inv(): given matrix is not square" );
496 
497  if(out.is_empty())
498  {
499  return true;
500  }
501 
502  bool status;
503 
504  #if defined(ARMA_USE_LAPACK)
505  {
506  char uplo = (layout == 0) ? 'U' : 'L';
507  blas_int n = blas_int(out.n_rows);
508  blas_int lwork = n*n; // TODO: use lwork = -1 to determine optimal size
509  blas_int info = 0;
510 
511  podarray<blas_int> ipiv;
512  ipiv.set_size(out.n_rows);
513 
514  podarray<eT> work;
515  work.set_size( uword(lwork) );
516 
517  lapack::sytrf(&uplo, &n, out.memptr(), &n, ipiv.memptr(), work.memptr(), &lwork, &info);
518 
519  status = (info == 0);
520 
521  if(status == true)
522  {
523  lapack::sytri(&uplo, &n, out.memptr(), &n, ipiv.memptr(), work.memptr(), &info);
524 
525  out = (layout == 0) ? symmatu(out) : symmatl(out);
526 
527  status = (info == 0);
528  }
529  }
530  #else
531  {
532  arma_ignore(layout);
533  arma_stop("inv(): use of LAPACK needs to be enabled");
534  status = false;
535  }
536  #endif
537 
538  return status;
539  }
540 
541 
542 
543 template<typename eT, typename T1>
544 inline
545 bool
546 auxlib::inv_sympd(Mat<eT>& out, const Base<eT,T1>& X, const uword layout)
547  {
548  arma_extra_debug_sigprint();
549 
550  out = X.get_ref();
551 
552  arma_debug_check( (out.is_square() == false), "inv(): given matrix is not square" );
553 
554  if(out.is_empty())
555  {
556  return true;
557  }
558 
559  bool status;
560 
561  #if defined(ARMA_USE_LAPACK)
562  {
563  char uplo = (layout == 0) ? 'U' : 'L';
564  blas_int n = blas_int(out.n_rows);
565  blas_int info = 0;
566 
567  lapack::potrf(&uplo, &n, out.memptr(), &n, &info);
568 
569  status = (info == 0);
570 
571  if(status == true)
572  {
573  lapack::potri(&uplo, &n, out.memptr(), &n, &info);
574 
575  out = (layout == 0) ? symmatu(out) : symmatl(out);
576 
577  status = (info == 0);
578  }
579  }
580  #else
581  {
582  arma_ignore(layout);
583  arma_stop("inv(): use of LAPACK needs to be enabled");
584  status = false;
585  }
586  #endif
587 
588  return status;
589  }
590 
591 
592 
593 template<typename eT, typename T1>
594 inline
595 eT
596 auxlib::det(const Base<eT,T1>& X, const bool slow)
597  {
598  const unwrap<T1> tmp(X.get_ref());
599  const Mat<eT>& A = tmp.M;
600 
601  arma_debug_check( (A.is_square() == false), "det(): matrix is not square" );
602 
603  const bool make_copy = (is_Mat<T1>::value == true) ? true : false;
604 
605  if(slow == false)
606  {
607  const uword N = A.n_rows;
608 
609  switch(N)
610  {
611  case 0:
612  case 1:
613  case 2:
614  return auxlib::det_tinymat(A, N);
615  break;
616 
617  case 3:
618  case 4:
619  {
620  const eT tmp_det = auxlib::det_tinymat(A, N);
621  return (tmp_det != eT(0)) ? tmp_det : auxlib::det_lapack(A, make_copy);
622  }
623  break;
624 
625  default:
626  return auxlib::det_lapack(A, make_copy);
627  }
628  }
629  else
630  {
631  return auxlib::det_lapack(A, make_copy);
632  }
633  }
634 
635 
636 
637 template<typename eT>
638 inline
639 eT
640 auxlib::det_tinymat(const Mat<eT>& X, const uword N)
641  {
642  arma_extra_debug_sigprint();
643 
644  switch(N)
645  {
646  case 0:
647  return eT(1);
648  break;
649 
650  case 1:
651  return X[0];
652  break;
653 
654  case 2:
655  {
656  const eT* Xm = X.memptr();
657 
658  return ( Xm[pos<0,0>::n2]*Xm[pos<1,1>::n2] - Xm[pos<0,1>::n2]*Xm[pos<1,0>::n2] );
659  }
660  break;
661 
662  case 3:
663  {
664  // const double tmp1 = X.at(0,0) * X.at(1,1) * X.at(2,2);
665  // const double tmp2 = X.at(0,1) * X.at(1,2) * X.at(2,0);
666  // const double tmp3 = X.at(0,2) * X.at(1,0) * X.at(2,1);
667  // const double tmp4 = X.at(2,0) * X.at(1,1) * X.at(0,2);
668  // const double tmp5 = X.at(2,1) * X.at(1,2) * X.at(0,0);
669  // const double tmp6 = X.at(2,2) * X.at(1,0) * X.at(0,1);
670  // return (tmp1+tmp2+tmp3) - (tmp4+tmp5+tmp6);
671 
672  const eT* a_col0 = X.colptr(0);
673  const eT a11 = a_col0[0];
674  const eT a21 = a_col0[1];
675  const eT a31 = a_col0[2];
676 
677  const eT* a_col1 = X.colptr(1);
678  const eT a12 = a_col1[0];
679  const eT a22 = a_col1[1];
680  const eT a32 = a_col1[2];
681 
682  const eT* a_col2 = X.colptr(2);
683  const eT a13 = a_col2[0];
684  const eT a23 = a_col2[1];
685  const eT a33 = a_col2[2];
686 
687  return ( a11*(a33*a22 - a32*a23) - a21*(a33*a12-a32*a13) + a31*(a23*a12 - a22*a13) );
688  }
689  break;
690 
691  case 4:
692  {
693  const eT* Xm = X.memptr();
694 
695  const eT val = \
696  Xm[pos<0,3>::n4] * Xm[pos<1,2>::n4] * Xm[pos<2,1>::n4] * Xm[pos<3,0>::n4] \
697  - Xm[pos<0,2>::n4] * Xm[pos<1,3>::n4] * Xm[pos<2,1>::n4] * Xm[pos<3,0>::n4] \
698  - Xm[pos<0,3>::n4] * Xm[pos<1,1>::n4] * Xm[pos<2,2>::n4] * Xm[pos<3,0>::n4] \
699  + Xm[pos<0,1>::n4] * Xm[pos<1,3>::n4] * Xm[pos<2,2>::n4] * Xm[pos<3,0>::n4] \
700  + Xm[pos<0,2>::n4] * Xm[pos<1,1>::n4] * Xm[pos<2,3>::n4] * Xm[pos<3,0>::n4] \
701  - Xm[pos<0,1>::n4] * Xm[pos<1,2>::n4] * Xm[pos<2,3>::n4] * Xm[pos<3,0>::n4] \
702  - Xm[pos<0,3>::n4] * Xm[pos<1,2>::n4] * Xm[pos<2,0>::n4] * Xm[pos<3,1>::n4] \
703  + Xm[pos<0,2>::n4] * Xm[pos<1,3>::n4] * Xm[pos<2,0>::n4] * Xm[pos<3,1>::n4] \
704  + Xm[pos<0,3>::n4] * Xm[pos<1,0>::n4] * Xm[pos<2,2>::n4] * Xm[pos<3,1>::n4] \
705  - Xm[pos<0,0>::n4] * Xm[pos<1,3>::n4] * Xm[pos<2,2>::n4] * Xm[pos<3,1>::n4] \
706  - Xm[pos<0,2>::n4] * Xm[pos<1,0>::n4] * Xm[pos<2,3>::n4] * Xm[pos<3,1>::n4] \
707  + Xm[pos<0,0>::n4] * Xm[pos<1,2>::n4] * Xm[pos<2,3>::n4] * Xm[pos<3,1>::n4] \
708  + Xm[pos<0,3>::n4] * Xm[pos<1,1>::n4] * Xm[pos<2,0>::n4] * Xm[pos<3,2>::n4] \
709  - Xm[pos<0,1>::n4] * Xm[pos<1,3>::n4] * Xm[pos<2,0>::n4] * Xm[pos<3,2>::n4] \
710  - Xm[pos<0,3>::n4] * Xm[pos<1,0>::n4] * Xm[pos<2,1>::n4] * Xm[pos<3,2>::n4] \
711  + Xm[pos<0,0>::n4] * Xm[pos<1,3>::n4] * Xm[pos<2,1>::n4] * Xm[pos<3,2>::n4] \
712  + Xm[pos<0,1>::n4] * Xm[pos<1,0>::n4] * Xm[pos<2,3>::n4] * Xm[pos<3,2>::n4] \
713  - Xm[pos<0,0>::n4] * Xm[pos<1,1>::n4] * Xm[pos<2,3>::n4] * Xm[pos<3,2>::n4] \
714  - Xm[pos<0,2>::n4] * Xm[pos<1,1>::n4] * Xm[pos<2,0>::n4] * Xm[pos<3,3>::n4] \
715  + Xm[pos<0,1>::n4] * Xm[pos<1,2>::n4] * Xm[pos<2,0>::n4] * Xm[pos<3,3>::n4] \
716  + Xm[pos<0,2>::n4] * Xm[pos<1,0>::n4] * Xm[pos<2,1>::n4] * Xm[pos<3,3>::n4] \
717  - Xm[pos<0,0>::n4] * Xm[pos<1,2>::n4] * Xm[pos<2,1>::n4] * Xm[pos<3,3>::n4] \
718  - Xm[pos<0,1>::n4] * Xm[pos<1,0>::n4] * Xm[pos<2,2>::n4] * Xm[pos<3,3>::n4] \
719  + Xm[pos<0,0>::n4] * Xm[pos<1,1>::n4] * Xm[pos<2,2>::n4] * Xm[pos<3,3>::n4] \
720  ;
721 
722  return val;
723  }
724  break;
725 
726  default:
727  return eT(0);
728  ;
729  }
730  }
731 
732 
733 
735 template<typename eT>
736 inline
737 eT
738 auxlib::det_lapack(const Mat<eT>& X, const bool make_copy)
739  {
740  arma_extra_debug_sigprint();
741 
742  Mat<eT> X_copy;
743 
744  if(make_copy == true)
745  {
746  X_copy = X;
747  }
748 
749  Mat<eT>& tmp = (make_copy == true) ? X_copy : const_cast< Mat<eT>& >(X);
750 
751  if(tmp.is_empty())
752  {
753  return eT(1);
754  }
755 
756 
757  #if defined(ARMA_USE_ATLAS)
758  {
759  podarray<int> ipiv(tmp.n_rows);
760 
761  //const int info =
762  atlas::clapack_getrf(atlas::CblasColMajor, tmp.n_rows, tmp.n_cols, tmp.memptr(), tmp.n_rows, ipiv.memptr());
763 
764  // on output tmp appears to be L+U_alt, where U_alt is U with the main diagonal set to zero
765  eT val = tmp.at(0,0);
766  for(uword i=1; i < tmp.n_rows; ++i)
767  {
768  val *= tmp.at(i,i);
769  }
770 
771  int sign = +1;
772  for(uword i=0; i < tmp.n_rows; ++i)
773  {
774  if( int(i) != ipiv.mem[i] ) // NOTE: no adjustment required, as the clapack version of getrf() assumes counting from 0
775  {
776  sign *= -1;
777  }
778  }
779 
780  return ( (sign < 0) ? -val : val );
781  }
782  #elif defined(ARMA_USE_LAPACK)
783  {
784  podarray<blas_int> ipiv(tmp.n_rows);
785 
786  blas_int info = 0;
787  blas_int n_rows = blas_int(tmp.n_rows);
788  blas_int n_cols = blas_int(tmp.n_cols);
789 
790  lapack::getrf(&n_rows, &n_cols, tmp.memptr(), &n_rows, ipiv.memptr(), &info);
791 
792  // on output tmp appears to be L+U_alt, where U_alt is U with the main diagonal set to zero
793  eT val = tmp.at(0,0);
794  for(uword i=1; i < tmp.n_rows; ++i)
795  {
796  val *= tmp.at(i,i);
797  }
798 
799  blas_int sign = +1;
800  for(uword i=0; i < tmp.n_rows; ++i)
801  {
802  if( blas_int(i) != (ipiv.mem[i] - 1) ) // NOTE: adjustment of -1 is required as Fortran counts from 1
803  {
804  sign *= -1;
805  }
806  }
807 
808  return ( (sign < 0) ? -val : val );
809  }
810  #else
811  {
812  arma_ignore(X);
813  arma_ignore(make_copy);
814  arma_ignore(tmp);
815  arma_stop("det(): use of ATLAS or LAPACK needs to be enabled");
816  return eT(0);
817  }
818  #endif
819  }
820 
821 
822 
824 template<typename eT, typename T1>
825 inline
826 bool
827 auxlib::log_det(eT& out_val, typename get_pod_type<eT>::result& out_sign, const Base<eT,T1>& X)
828  {
829  arma_extra_debug_sigprint();
830 
831  typedef typename get_pod_type<eT>::result T;
832 
833  #if defined(ARMA_USE_ATLAS)
834  {
835  Mat<eT> tmp(X.get_ref());
836  arma_debug_check( (tmp.is_square() == false), "log_det(): given matrix is not square" );
837 
838  if(tmp.is_empty())
839  {
840  out_val = eT(0);
841  out_sign = T(1);
842  return true;
843  }
844 
845  podarray<int> ipiv(tmp.n_rows);
846 
847  const int info = atlas::clapack_getrf(atlas::CblasColMajor, tmp.n_rows, tmp.n_cols, tmp.memptr(), tmp.n_rows, ipiv.memptr());
848 
849  // on output tmp appears to be L+U_alt, where U_alt is U with the main diagonal set to zero
850 
851  sword sign = (is_complex<eT>::value == false) ? ( (access::tmp_real( tmp.at(0,0) ) < T(0)) ? -1 : +1 ) : +1;
852  eT val = (is_complex<eT>::value == false) ? std::log( (access::tmp_real( tmp.at(0,0) ) < T(0)) ? tmp.at(0,0)*T(-1) : tmp.at(0,0) ) : std::log( tmp.at(0,0) );
853 
854  for(uword i=1; i < tmp.n_rows; ++i)
855  {
856  const eT x = tmp.at(i,i);
857 
858  sign *= (is_complex<eT>::value == false) ? ( (access::tmp_real(x) < T(0)) ? -1 : +1 ) : +1;
859  val += (is_complex<eT>::value == false) ? std::log( (access::tmp_real(x) < T(0)) ? x*T(-1) : x ) : std::log(x);
860  }
861 
862  for(uword i=0; i < tmp.n_rows; ++i)
863  {
864  if( int(i) != ipiv.mem[i] ) // NOTE: no adjustment required, as the clapack version of getrf() assumes counting from 0
865  {
866  sign *= -1;
867  }
868  }
869 
870  out_val = val;
871  out_sign = T(sign);
872 
873  return (info == 0);
874  }
875  #elif defined(ARMA_USE_LAPACK)
876  {
877  Mat<eT> tmp(X.get_ref());
878  arma_debug_check( (tmp.is_square() == false), "log_det(): given matrix is not square" );
879 
880  if(tmp.is_empty())
881  {
882  out_val = eT(0);
883  out_sign = T(1);
884  return true;
885  }
886 
887  podarray<blas_int> ipiv(tmp.n_rows);
888 
889  blas_int info = 0;
890  blas_int n_rows = blas_int(tmp.n_rows);
891  blas_int n_cols = blas_int(tmp.n_cols);
892 
893  lapack::getrf(&n_rows, &n_cols, tmp.memptr(), &n_rows, ipiv.memptr(), &info);
894 
895  // on output tmp appears to be L+U_alt, where U_alt is U with the main diagonal set to zero
896 
897  sword sign = (is_complex<eT>::value == false) ? ( (access::tmp_real( tmp.at(0,0) ) < T(0)) ? -1 : +1 ) : +1;
898  eT val = (is_complex<eT>::value == false) ? std::log( (access::tmp_real( tmp.at(0,0) ) < T(0)) ? tmp.at(0,0)*T(-1) : tmp.at(0,0) ) : std::log( tmp.at(0,0) );
899 
900  for(uword i=1; i < tmp.n_rows; ++i)
901  {
902  const eT x = tmp.at(i,i);
903 
904  sign *= (is_complex<eT>::value == false) ? ( (access::tmp_real(x) < T(0)) ? -1 : +1 ) : +1;
905  val += (is_complex<eT>::value == false) ? std::log( (access::tmp_real(x) < T(0)) ? x*T(-1) : x ) : std::log(x);
906  }
907 
908  for(uword i=0; i < tmp.n_rows; ++i)
909  {
910  if( blas_int(i) != (ipiv.mem[i] - 1) ) // NOTE: adjustment of -1 is required as Fortran counts from 1
911  {
912  sign *= -1;
913  }
914  }
915 
916  out_val = val;
917  out_sign = T(sign);
918 
919  return (info == 0);
920  }
921  #else
922  {
923  out_val = eT(0);
924  out_sign = T(0);
925 
926  arma_stop("log_det(): use of ATLAS or LAPACK needs to be enabled");
927 
928  return false;
929  }
930  #endif
931  }
932 
933 
934 
936 template<typename eT, typename T1>
937 inline
938 bool
939 auxlib::lu(Mat<eT>& L, Mat<eT>& U, podarray<blas_int>& ipiv, const Base<eT,T1>& X)
940  {
941  arma_extra_debug_sigprint();
942 
943  U = X.get_ref();
944 
945  const uword U_n_rows = U.n_rows;
946  const uword U_n_cols = U.n_cols;
947 
948  if(U.is_empty())
949  {
950  L.set_size(U_n_rows, 0);
951  U.set_size(0, U_n_cols);
952  ipiv.reset();
953  return true;
954  }
955 
956  #if defined(ARMA_USE_ATLAS) || defined(ARMA_USE_LAPACK)
957  {
958  bool status;
959 
960  #if defined(ARMA_USE_ATLAS)
961  {
962  ipiv.set_size( (std::min)(U_n_rows, U_n_cols) );
963 
964  int info = atlas::clapack_getrf(atlas::CblasColMajor, U_n_rows, U_n_cols, U.memptr(), U_n_rows, ipiv.memptr());
965 
966  status = (info == 0);
967  }
968  #elif defined(ARMA_USE_LAPACK)
969  {
970  ipiv.set_size( (std::min)(U_n_rows, U_n_cols) );
971 
972  blas_int info = 0;
973 
974  blas_int n_rows = U_n_rows;
975  blas_int n_cols = U_n_cols;
976 
977 
978  lapack::getrf(&n_rows, &n_cols, U.memptr(), &n_rows, ipiv.memptr(), &info);
979 
980  // take into account that Fortran counts from 1
981  arrayops::inplace_minus(ipiv.memptr(), blas_int(1), ipiv.n_elem);
982 
983  status = (info == 0);
984  }
985  #endif
986 
987  L.copy_size(U);
988 
989  for(uword col=0; col < U_n_cols; ++col)
990  {
991  for(uword row=0; (row < col) && (row < U_n_rows); ++row)
992  {
993  L.at(row,col) = eT(0);
994  }
995 
996  if( L.in_range(col,col) == true )
997  {
998  L.at(col,col) = eT(1);
999  }
1000 
1001  for(uword row = (col+1); row < U_n_rows; ++row)
1002  {
1003  L.at(row,col) = U.at(row,col);
1004  U.at(row,col) = eT(0);
1005  }
1006  }
1007 
1008  return status;
1009  }
1010  #else
1011  {
1012  arma_stop("lu(): use of ATLAS or LAPACK needs to be enabled");
1013 
1014  return false;
1015  }
1016  #endif
1017  }
1018 
1019 
1020 
1021 template<typename eT, typename T1>
1022 inline
1023 bool
1024 auxlib::lu(Mat<eT>& L, Mat<eT>& U, Mat<eT>& P, const Base<eT,T1>& X)
1025  {
1026  arma_extra_debug_sigprint();
1027 
1028  podarray<blas_int> ipiv1;
1029  const bool status = auxlib::lu(L, U, ipiv1, X);
1030 
1031  if(status == true)
1032  {
1033  if(U.is_empty())
1034  {
1035  // L and U have been already set to the correct empty matrices
1036  P.eye(L.n_rows, L.n_rows);
1037  return true;
1038  }
1039 
1040  const uword n = ipiv1.n_elem;
1041  const uword P_rows = U.n_rows;
1042 
1043  podarray<blas_int> ipiv2(P_rows);
1044 
1045  const blas_int* ipiv1_mem = ipiv1.memptr();
1046  blas_int* ipiv2_mem = ipiv2.memptr();
1047 
1048  for(uword i=0; i<P_rows; ++i)
1049  {
1050  ipiv2_mem[i] = blas_int(i);
1051  }
1052 
1053  for(uword i=0; i<n; ++i)
1054  {
1055  const uword k = static_cast<uword>(ipiv1_mem[i]);
1056 
1057  if( ipiv2_mem[i] != ipiv2_mem[k] )
1058  {
1059  std::swap( ipiv2_mem[i], ipiv2_mem[k] );
1060  }
1061  }
1062 
1063  P.zeros(P_rows, P_rows);
1064 
1065  for(uword row=0; row<P_rows; ++row)
1066  {
1067  P.at(row, static_cast<uword>(ipiv2_mem[row])) = eT(1);
1068  }
1069 
1070  if(L.n_cols > U.n_rows)
1071  {
1072  L.shed_cols(U.n_rows, L.n_cols-1);
1073  }
1074 
1075  if(U.n_rows > L.n_cols)
1076  {
1077  U.shed_rows(L.n_cols, U.n_rows-1);
1078  }
1079  }
1080 
1081  return status;
1082  }
1083 
1084 
1085 
1086 template<typename eT, typename T1>
1087 inline
1088 bool
1089 auxlib::lu(Mat<eT>& L, Mat<eT>& U, const Base<eT,T1>& X)
1090  {
1091  arma_extra_debug_sigprint();
1092 
1093  podarray<blas_int> ipiv1;
1094  const bool status = auxlib::lu(L, U, ipiv1, X);
1095 
1096  if(status == true)
1097  {
1098  if(U.is_empty())
1099  {
1100  // L and U have been already set to the correct empty matrices
1101  return true;
1102  }
1103 
1104  const uword n = ipiv1.n_elem;
1105  const uword P_rows = U.n_rows;
1106 
1107  podarray<blas_int> ipiv2(P_rows);
1108 
1109  const blas_int* ipiv1_mem = ipiv1.memptr();
1110  blas_int* ipiv2_mem = ipiv2.memptr();
1111 
1112  for(uword i=0; i<P_rows; ++i)
1113  {
1114  ipiv2_mem[i] = blas_int(i);
1115  }
1116 
1117  for(uword i=0; i<n; ++i)
1118  {
1119  const uword k = static_cast<uword>(ipiv1_mem[i]);
1120 
1121  if( ipiv2_mem[i] != ipiv2_mem[k] )
1122  {
1123  std::swap( ipiv2_mem[i], ipiv2_mem[k] );
1124  L.swap_rows( static_cast<uword>(ipiv2_mem[i]), static_cast<uword>(ipiv2_mem[k]) );
1125  }
1126  }
1127 
1128  if(L.n_cols > U.n_rows)
1129  {
1130  L.shed_cols(U.n_rows, L.n_cols-1);
1131  }
1132 
1133  if(U.n_rows > L.n_cols)
1134  {
1135  U.shed_rows(L.n_cols, U.n_rows-1);
1136  }
1137  }
1138 
1139  return status;
1140  }
1141 
1142 
1143 
1145 template<typename eT, typename T1>
1146 inline
1147 bool
1148 auxlib::eig_sym(Col<eT>& eigval, const Base<eT,T1>& X)
1149  {
1150  arma_extra_debug_sigprint();
1151 
1152  #if defined(ARMA_USE_LAPACK)
1153  {
1154  Mat<eT> A(X.get_ref());
1155 
1156  arma_debug_check( (A.is_square() == false), "eig_sym(): given matrix is not square");
1157 
1158  if(A.is_empty())
1159  {
1160  eigval.reset();
1161  return true;
1162  }
1163 
1164  // rudimentary "better-than-nothing" test for symmetry
1165  //arma_debug_check( (A.at(A.n_rows-1, 0) != A.at(0, A.n_cols-1)), "auxlib::eig(): given matrix is not symmetric" );
1166 
1167  char jobz = 'N';
1168  char uplo = 'U';
1169 
1170  blas_int n_rows = A.n_rows;
1171  blas_int lwork = (std::max)(blas_int(1), 3*n_rows-1);
1172 
1173  eigval.set_size( static_cast<uword>(n_rows) );
1174  podarray<eT> work( static_cast<uword>(lwork) );
1175 
1176  blas_int info;
1177 
1178  arma_extra_debug_print("lapack::syev()");
1179  lapack::syev(&jobz, &uplo, &n_rows, A.memptr(), &n_rows, eigval.memptr(), work.memptr(), &lwork, &info);
1180 
1181  return (info == 0);
1182  }
1183  #else
1184  {
1185  arma_ignore(eigval);
1186  arma_ignore(X);
1187  arma_stop("eig_sym(): use of LAPACK needs to be enabled");
1188  return false;
1189  }
1190  #endif
1191  }
1192 
1193 
1194 
1196 template<typename T, typename T1>
1197 inline
1198 bool
1199 auxlib::eig_sym(Col<T>& eigval, const Base<std::complex<T>,T1>& X)
1200  {
1201  arma_extra_debug_sigprint();
1202 
1203  typedef typename std::complex<T> eT;
1204 
1205  #if defined(ARMA_USE_LAPACK)
1206  {
1207  Mat<eT> A(X.get_ref());
1208  arma_debug_check( (A.is_square() == false), "eig_sym(): given matrix is not hermitian");
1209 
1210  if(A.is_empty())
1211  {
1212  eigval.reset();
1213  return true;
1214  }
1215 
1216  char jobz = 'N';
1217  char uplo = 'U';
1218 
1219  blas_int n_rows = A.n_rows;
1220  blas_int lda = A.n_rows;
1221  blas_int lwork = (std::max)(blas_int(1), 2*n_rows - 1); // TODO: automatically find best size of lwork
1222 
1223  eigval.set_size( static_cast<uword>(n_rows) );
1224 
1225  podarray<eT> work( static_cast<uword>(lwork) );
1226  podarray<T> rwork( static_cast<uword>((std::max)(blas_int(1), 3*n_rows - 2)) );
1227 
1228  blas_int info;
1229 
1230  arma_extra_debug_print("lapack::heev()");
1231  lapack::heev(&jobz, &uplo, &n_rows, A.memptr(), &lda, eigval.memptr(), work.memptr(), &lwork, rwork.memptr(), &info);
1232 
1233  return (info == 0);
1234  }
1235  #else
1236  {
1237  arma_ignore(eigval);
1238  arma_ignore(X);
1239  arma_stop("eig_sym(): use of LAPACK needs to be enabled");
1240  return false;
1241  }
1242  #endif
1243  }
1244 
1245 
1246 
1248 template<typename eT, typename T1>
1249 inline
1250 bool
1251 auxlib::eig_sym(Col<eT>& eigval, Mat<eT>& eigvec, const Base<eT,T1>& X)
1252  {
1253  arma_extra_debug_sigprint();
1254 
1255  #if defined(ARMA_USE_LAPACK)
1256  {
1257  eigvec = X.get_ref();
1258 
1259  arma_debug_check( (eigvec.is_square() == false), "eig_sym(): given matrix is not square" );
1260 
1261  if(eigvec.is_empty())
1262  {
1263  eigval.reset();
1264  eigvec.reset();
1265  return true;
1266  }
1267 
1268  // rudimentary "better-than-nothing" test for symmetry
1269  //arma_debug_check( (A.at(A.n_rows-1, 0) != A.at(0, A.n_cols-1)), "auxlib::eig(): given matrix is not symmetric" );
1270 
1271  char jobz = 'V';
1272  char uplo = 'U';
1273 
1274  blas_int n_rows = eigvec.n_rows;
1275  blas_int lwork = (std::max)(blas_int(1), 3*n_rows-1);
1276 
1277  eigval.set_size( static_cast<uword>(n_rows) );
1278  podarray<eT> work( static_cast<uword>(lwork) );
1279 
1280  blas_int info;
1281 
1282  arma_extra_debug_print("lapack::syev()");
1283  lapack::syev(&jobz, &uplo, &n_rows, eigvec.memptr(), &n_rows, eigval.memptr(), work.memptr(), &lwork, &info);
1284 
1285  return (info == 0);
1286  }
1287  #else
1288  {
1289  arma_ignore(eigval);
1290  arma_ignore(eigvec);
1291  arma_stop("eig_sym(): use of LAPACK needs to be enabled");
1292 
1293  return false;
1294  }
1295  #endif
1296  }
1297 
1298 
1299 
1301 template<typename T, typename T1>
1302 inline
1303 bool
1304 auxlib::eig_sym(Col<T>& eigval, Mat< std::complex<T> >& eigvec, const Base<std::complex<T>,T1>& X)
1305  {
1306  arma_extra_debug_sigprint();
1307 
1308  typedef typename std::complex<T> eT;
1309 
1310  #if defined(ARMA_USE_LAPACK)
1311  {
1312  eigvec = X.get_ref();
1313 
1314  arma_debug_check( (eigvec.is_square() == false), "eig_sym(): given matrix is not hermitian" );
1315 
1316  if(eigvec.is_empty())
1317  {
1318  eigval.reset();
1319  eigvec.reset();
1320  return true;
1321  }
1322 
1323  char jobz = 'V';
1324  char uplo = 'U';
1325 
1326  blas_int n_rows = eigvec.n_rows;
1327  blas_int lda = eigvec.n_rows;
1328  blas_int lwork = (std::max)(blas_int(1), 2*n_rows - 1); // TODO: automatically find best size of lwork
1329 
1330  eigval.set_size( static_cast<uword>(n_rows) );
1331 
1332  podarray<eT> work( static_cast<uword>(lwork) );
1333  podarray<T> rwork( static_cast<uword>((std::max)(blas_int(1), 3*n_rows - 2)) );
1334 
1335  blas_int info;
1336 
1337  arma_extra_debug_print("lapack::heev()");
1338  lapack::heev(&jobz, &uplo, &n_rows, eigvec.memptr(), &lda, eigval.memptr(), work.memptr(), &lwork, rwork.memptr(), &info);
1339 
1340  return (info == 0);
1341  }
1342  #else
1343  {
1344  arma_ignore(eigval);
1345  arma_ignore(eigvec);
1346  arma_ignore(X);
1347  arma_stop("eig_sym(): use of LAPACK needs to be enabled");
1348  return false;
1349  }
1350  #endif
1351  }
1352 
1353 
1354 
1358 template<typename T, typename T1>
1359 inline
1360 bool
1361 auxlib::eig_gen
1362  (
1363  Col< std::complex<T> >& eigval,
1364  Mat<T>& l_eigvec,
1365  Mat<T>& r_eigvec,
1366  const Base<T,T1>& X,
1367  const char side
1368  )
1369  {
1370  arma_extra_debug_sigprint();
1371 
1372  #if defined(ARMA_USE_LAPACK)
1373  {
1374  char jobvl;
1375  char jobvr;
1376 
1377  switch(side)
1378  {
1379  case 'l': // left
1380  jobvl = 'V';
1381  jobvr = 'N';
1382  break;
1383 
1384  case 'r': // right
1385  jobvl = 'N';
1386  jobvr = 'V';
1387  break;
1388 
1389  case 'b': // both
1390  jobvl = 'V';
1391  jobvr = 'V';
1392  break;
1393 
1394  case 'n': // neither
1395  jobvl = 'N';
1396  jobvr = 'N';
1397  break;
1398 
1399  default:
1400  arma_stop("eig_gen(): parameter 'side' is invalid");
1401  return false;
1402  }
1403 
1404  Mat<T> A(X.get_ref());
1405  arma_debug_check( (A.is_square() == false), "eig_gen(): given matrix is not square" );
1406 
1407  if(A.is_empty())
1408  {
1409  eigval.reset();
1410  l_eigvec.reset();
1411  r_eigvec.reset();
1412  return true;
1413  }
1414 
1415  uword A_n_rows = A.n_rows;
1416 
1417  blas_int n_rows = A_n_rows;
1418  blas_int lda = A_n_rows;
1419  blas_int lwork = (std::max)(blas_int(1), 4*n_rows); // TODO: automatically find best size of lwork
1420 
1421  eigval.set_size(A_n_rows);
1422  l_eigvec.set_size(A_n_rows, A_n_rows);
1423  r_eigvec.set_size(A_n_rows, A_n_rows);
1424 
1425  podarray<T> work( static_cast<uword>(lwork) );
1426  podarray<T> rwork( static_cast<uword>((std::max)(blas_int(1), 3*n_rows)) );
1427 
1428  podarray<T> wr(A_n_rows);
1429  podarray<T> wi(A_n_rows);
1430 
1431  Mat<T> A_copy = A;
1432  blas_int info;
1433 
1434  arma_extra_debug_print("lapack::geev()");
1435  lapack::geev(&jobvl, &jobvr, &n_rows, A_copy.memptr(), &lda, wr.memptr(), wi.memptr(), l_eigvec.memptr(), &n_rows, r_eigvec.memptr(), &n_rows, work.memptr(), &lwork, &info);
1436 
1437 
1438  eigval.set_size(A_n_rows);
1439  for(uword i=0; i<A_n_rows; ++i)
1440  {
1441  eigval[i] = std::complex<T>(wr[i], wi[i]);
1442  }
1443 
1444  return (info == 0);
1445  }
1446  #else
1447  {
1448  arma_ignore(eigval);
1449  arma_ignore(l_eigvec);
1450  arma_ignore(r_eigvec);
1451  arma_ignore(X);
1452  arma_ignore(side);
1453  arma_stop("eig_gen(): use of LAPACK needs to be enabled");
1454  return false;
1455  }
1456  #endif
1457  }
1458 
1459 
1460 
1461 
1462 
1466 template<typename T, typename T1>
1467 inline
1468 bool
1469 auxlib::eig_gen
1470  (
1471  Col< std::complex<T> >& eigval,
1472  Mat< std::complex<T> >& l_eigvec,
1473  Mat< std::complex<T> >& r_eigvec,
1474  const Base< std::complex<T>, T1 >& X,
1475  const char side
1476  )
1477  {
1478  arma_extra_debug_sigprint();
1479 
1480  typedef typename std::complex<T> eT;
1481 
1482  #if defined(ARMA_USE_LAPACK)
1483  {
1484  char jobvl;
1485  char jobvr;
1486 
1487  switch(side)
1488  {
1489  case 'l': // left
1490  jobvl = 'V';
1491  jobvr = 'N';
1492  break;
1493 
1494  case 'r': // right
1495  jobvl = 'N';
1496  jobvr = 'V';
1497  break;
1498 
1499  case 'b': // both
1500  jobvl = 'V';
1501  jobvr = 'V';
1502  break;
1503 
1504  case 'n': // neither
1505  jobvl = 'N';
1506  jobvr = 'N';
1507  break;
1508 
1509  default:
1510  arma_stop("eig_gen(): parameter 'side' is invalid");
1511  return false;
1512  }
1513 
1514  Mat<eT> A(X.get_ref());
1515  arma_debug_check( (A.is_square() == false), "eig_gen(): given matrix is not square" );
1516 
1517  if(A.is_empty())
1518  {
1519  eigval.reset();
1520  l_eigvec.reset();
1521  r_eigvec.reset();
1522  return true;
1523  }
1524 
1525  uword A_n_rows = A.n_rows;
1526 
1527  blas_int n_rows = A_n_rows;
1528  blas_int lda = A_n_rows;
1529  blas_int lwork = (std::max)(blas_int(1), 4*n_rows); // TODO: automatically find best size of lwork
1530 
1531  eigval.set_size(A_n_rows);
1532  l_eigvec.set_size(A_n_rows, A_n_rows);
1533  r_eigvec.set_size(A_n_rows, A_n_rows);
1534 
1535  podarray<eT> work( static_cast<uword>(lwork) );
1536  podarray<T> rwork( static_cast<uword>((std::max)(blas_int(1), 3*n_rows)) ); // was 2,3
1537 
1538  blas_int info;
1539 
1540  arma_extra_debug_print("lapack::cx_geev()");
1541  lapack::cx_geev(&jobvl, &jobvr, &n_rows, A.memptr(), &lda, eigval.memptr(), l_eigvec.memptr(), &n_rows, r_eigvec.memptr(), &n_rows, work.memptr(), &lwork, rwork.memptr(), &info);
1542 
1543  return (info == 0);
1544  }
1545  #else
1546  {
1547  arma_ignore(eigval);
1548  arma_ignore(l_eigvec);
1549  arma_ignore(r_eigvec);
1550  arma_ignore(X);
1551  arma_ignore(side);
1552  arma_stop("eig_gen(): use of LAPACK needs to be enabled");
1553  return false;
1554  }
1555  #endif
1556  }
1557 
1558 
1559 
1560 template<typename eT, typename T1>
1561 inline
1562 bool
1563 auxlib::chol(Mat<eT>& out, const Base<eT,T1>& X)
1564  {
1565  arma_extra_debug_sigprint();
1566 
1567  #if defined(ARMA_USE_LAPACK)
1568  {
1569  out = X.get_ref();
1570 
1571  arma_debug_check( (out.is_square() == false), "chol(): given matrix is not square" );
1572 
1573  if(out.is_empty())
1574  {
1575  return true;
1576  }
1577 
1578  const uword out_n_rows = out.n_rows;
1579 
1580  char uplo = 'U';
1581  blas_int n = out_n_rows;
1582  blas_int info;
1583 
1584  lapack::potrf(&uplo, &n, out.memptr(), &n, &info);
1585 
1586  for(uword col=0; col<out_n_rows; ++col)
1587  {
1588  eT* colptr = out.colptr(col);
1589 
1590  for(uword row=(col+1); row < out_n_rows; ++row)
1591  {
1592  colptr[row] = eT(0);
1593  }
1594  }
1595 
1596  return (info == 0);
1597  }
1598  #else
1599  {
1600  arma_ignore(out);
1601  arma_stop("chol(): use of LAPACK needs to be enabled");
1602  return false;
1603  }
1604  #endif
1605  }
1606 
1607 
1608 
1609 template<typename eT, typename T1>
1610 inline
1611 bool
1612 auxlib::qr(Mat<eT>& Q, Mat<eT>& R, const Base<eT,T1>& X)
1613  {
1614  arma_extra_debug_sigprint();
1615 
1616  #if defined(ARMA_USE_LAPACK)
1617  {
1618  R = X.get_ref();
1619 
1620  const uword R_n_rows = R.n_rows;
1621  const uword R_n_cols = R.n_cols;
1622 
1623  if(R.is_empty())
1624  {
1625  Q.eye(R_n_rows, R_n_rows);
1626  return true;
1627  }
1628 
1629  blas_int m = static_cast<blas_int>(R_n_rows);
1630  blas_int n = static_cast<blas_int>(R_n_cols);
1631  blas_int work_len = (std::max)(blas_int(1),n);
1632  blas_int work_len_tmp;
1633  blas_int k = (std::min)(m,n);
1634  blas_int info;
1635 
1636  podarray<eT> tau( static_cast<uword>(k) );
1637  podarray<eT> work( static_cast<uword>(work_len) );
1638 
1639  // query for the optimum value of work_len
1640  work_len_tmp = -1;
1641  lapack::geqrf(&m, &n, R.memptr(), &m, tau.memptr(), work.memptr(), &work_len_tmp, &info);
1642 
1643  if(info == 0)
1644  {
1645  work_len = static_cast<blas_int>(access::tmp_real(work[0]));
1646  work.set_size( static_cast<uword>(work_len) );
1647  }
1648 
1649  lapack::geqrf(&m, &n, R.memptr(), &m, tau.memptr(), work.memptr(), &work_len, &info);
1650 
1651  Q.set_size(R_n_rows, R_n_rows);
1652 
1653  arrayops::copy( Q.memptr(), R.memptr(), (std::min)(Q.n_elem, R.n_elem) );
1654 
1655  //
1656  // construct R
1657 
1658  for(uword col=0; col < R_n_cols; ++col)
1659  {
1660  for(uword row=(col+1); row < R_n_rows; ++row)
1661  {
1662  R.at(row,col) = eT(0);
1663  }
1664  }
1665 
1666 
1667  if( (is_float<eT>::value == true) || (is_double<eT>::value == true) )
1668  {
1669  // query for the optimum value of work_len
1670  work_len_tmp = -1;
1671  lapack::orgqr(&m, &m, &k, Q.memptr(), &m, tau.memptr(), work.memptr(), &work_len_tmp, &info);
1672 
1673  if(info == 0)
1674  {
1675  work_len = static_cast<blas_int>(access::tmp_real(work[0]));
1676  work.set_size( static_cast<uword>(work_len) );
1677  }
1678 
1679  lapack::orgqr(&m, &m, &k, Q.memptr(), &m, tau.memptr(), work.memptr(), &work_len, &info);
1680  }
1681  else
1682  if( (is_supported_complex_float<eT>::value == true) || (is_supported_complex_double<eT>::value == true) )
1683  {
1684  // query for the optimum value of work_len
1685  work_len_tmp = -1;
1686  lapack::ungqr(&m, &m, &k, Q.memptr(), &m, tau.memptr(), work.memptr(), &work_len_tmp, &info);
1687 
1688  if(info == 0)
1689  {
1690  work_len = static_cast<blas_int>(access::tmp_real(work[0]));
1691  work.set_size( static_cast<uword>(work_len) );
1692  }
1693 
1694  lapack::ungqr(&m, &m, &k, Q.memptr(), &m, tau.memptr(), work.memptr(), &work_len, &info);
1695  }
1696 
1697  return (info == 0);
1698  }
1699  #else
1700  {
1701  arma_ignore(Q);
1702  arma_ignore(R);
1703  arma_ignore(X);
1704  arma_stop("qr(): use of LAPACK needs to be enabled");
1705  return false;
1706  }
1707  #endif
1708  }
1709 
1710 
1711 
1712 template<typename eT, typename T1>
1713 inline
1714 bool
1715 auxlib::svd(Col<eT>& S, const Base<eT,T1>& X, uword& X_n_rows, uword& X_n_cols)
1716  {
1717  arma_extra_debug_sigprint();
1718 
1719  #if defined(ARMA_USE_LAPACK)
1720  {
1721  Mat<eT> A(X.get_ref());
1722 
1723  X_n_rows = A.n_rows;
1724  X_n_cols = A.n_cols;
1725 
1726  if(A.is_empty())
1727  {
1728  S.reset();
1729  return true;
1730  }
1731 
1732  Mat<eT> U(1, 1);
1733  Mat<eT> V(1, A.n_cols);
1734 
1735  char jobu = 'N';
1736  char jobvt = 'N';
1737 
1738  blas_int m = A.n_rows;
1739  blas_int n = A.n_cols;
1740  blas_int lda = A.n_rows;
1741  blas_int ldu = U.n_rows;
1742  blas_int ldvt = V.n_rows;
1743  blas_int lwork = 2 * (std::max)(blas_int(1), (std::max)( (3*(std::min)(m,n) + (std::max)(m,n)), 5*(std::min)(m,n) ) );
1744  blas_int info;
1745 
1746  S.set_size( static_cast<uword>((std::min)(m, n)) );
1747 
1748  podarray<eT> work( static_cast<uword>(lwork) );
1749 
1750 
1751  // let gesvd_() calculate the optimum size of the workspace
1752  blas_int lwork_tmp = -1;
1753 
1754  lapack::gesvd<eT>
1755  (
1756  &jobu, &jobvt,
1757  &m,&n,
1758  A.memptr(), &lda,
1759  S.memptr(),
1760  U.memptr(), &ldu,
1761  V.memptr(), &ldvt,
1762  work.memptr(), &lwork_tmp,
1763  &info
1764  );
1765 
1766  if(info == 0)
1767  {
1768  blas_int proposed_lwork = static_cast<blas_int>(work[0]);
1769 
1770  if(proposed_lwork > lwork)
1771  {
1772  lwork = proposed_lwork;
1773  work.set_size( static_cast<uword>(lwork) );
1774  }
1775 
1776  lapack::gesvd<eT>
1777  (
1778  &jobu, &jobvt,
1779  &m, &n,
1780  A.memptr(), &lda,
1781  S.memptr(),
1782  U.memptr(), &ldu,
1783  V.memptr(), &ldvt,
1784  work.memptr(), &lwork,
1785  &info
1786  );
1787  }
1788 
1789  return (info == 0);
1790  }
1791  #else
1792  {
1793  arma_ignore(S);
1794  arma_ignore(X);
1795  arma_ignore(X_n_rows);
1796  arma_ignore(X_n_cols);
1797  arma_stop("svd(): use of LAPACK needs to be enabled");
1798  return false;
1799  }
1800  #endif
1801  }
1802 
1803 
1804 
1805 template<typename T, typename T1>
1806 inline
1807 bool
1808 auxlib::svd(Col<T>& S, const Base<std::complex<T>, T1>& X, uword& X_n_rows, uword& X_n_cols)
1809  {
1810  arma_extra_debug_sigprint();
1811 
1812  typedef std::complex<T> eT;
1813 
1814  #if defined(ARMA_USE_LAPACK)
1815  {
1816  Mat<eT> A(X.get_ref());
1817 
1818  X_n_rows = A.n_rows;
1819  X_n_cols = A.n_cols;
1820 
1821  if(A.is_empty())
1822  {
1823  S.reset();
1824  return true;
1825  }
1826 
1827  Mat<eT> U(1, 1);
1828  Mat<eT> V(1, A.n_cols);
1829 
1830  char jobu = 'N';
1831  char jobvt = 'N';
1832 
1833  blas_int m = A.n_rows;
1834  blas_int n = A.n_cols;
1835  blas_int lda = A.n_rows;
1836  blas_int ldu = U.n_rows;
1837  blas_int ldvt = V.n_rows;
1838  blas_int lwork = 2 * (std::max)(blas_int(1), 2*(std::min)(m,n)+(std::max)(m,n) );
1839  blas_int info;
1840 
1841  S.set_size( static_cast<uword>((std::min)(m,n)) );
1842 
1843  podarray<eT> work( static_cast<uword>(lwork) );
1844  podarray<T> rwork( static_cast<uword>(5*(std::min)(m,n)) );
1845 
1846  // let gesvd_() calculate the optimum size of the workspace
1847  blas_int lwork_tmp = -1;
1848 
1849  lapack::cx_gesvd<T>
1850  (
1851  &jobu, &jobvt,
1852  &m, &n,
1853  A.memptr(), &lda,
1854  S.memptr(),
1855  U.memptr(), &ldu,
1856  V.memptr(), &ldvt,
1857  work.memptr(), &lwork_tmp,
1858  rwork.memptr(),
1859  &info
1860  );
1861 
1862  if(info == 0)
1863  {
1864  blas_int proposed_lwork = static_cast<blas_int>(real(work[0]));
1865  if(proposed_lwork > lwork)
1866  {
1867  lwork = proposed_lwork;
1868  work.set_size( static_cast<uword>(lwork) );
1869  }
1870 
1871  lapack::cx_gesvd<T>
1872  (
1873  &jobu, &jobvt,
1874  &m, &n,
1875  A.memptr(), &lda,
1876  S.memptr(),
1877  U.memptr(), &ldu,
1878  V.memptr(), &ldvt,
1879  work.memptr(), &lwork,
1880  rwork.memptr(),
1881  &info
1882  );
1883  }
1884 
1885  return (info == 0);
1886  }
1887  #else
1888  {
1889  arma_ignore(S);
1890  arma_ignore(X);
1891  arma_ignore(X_n_rows);
1892  arma_ignore(X_n_cols);
1893 
1894  arma_stop("svd(): use of LAPACK needs to be enabled");
1895  return false;
1896  }
1897  #endif
1898  }
1899 
1900 
1901 
1902 template<typename eT, typename T1>
1903 inline
1904 bool
1905 auxlib::svd(Col<eT>& S, const Base<eT,T1>& X)
1906  {
1907  arma_extra_debug_sigprint();
1908 
1909  uword junk;
1910  return auxlib::svd(S, X, junk, junk);
1911  }
1912 
1913 
1914 
1915 template<typename T, typename T1>
1916 inline
1917 bool
1918 auxlib::svd(Col<T>& S, const Base<std::complex<T>, T1>& X)
1919  {
1920  arma_extra_debug_sigprint();
1921 
1922  uword junk;
1923  return auxlib::svd(S, X, junk, junk);
1924  }
1925 
1926 
1927 
1928 template<typename eT, typename T1>
1929 inline
1930 bool
1931 auxlib::svd(Mat<eT>& U, Col<eT>& S, Mat<eT>& V, const Base<eT,T1>& X)
1932  {
1933  arma_extra_debug_sigprint();
1934 
1935  #if defined(ARMA_USE_LAPACK)
1936  {
1937  Mat<eT> A(X.get_ref());
1938 
1939  if(A.is_empty())
1940  {
1941  U.eye(A.n_rows, A.n_rows);
1942  S.reset();
1943  V.eye(A.n_cols, A.n_cols);
1944  return true;
1945  }
1946 
1947  U.set_size(A.n_rows, A.n_rows);
1948  V.set_size(A.n_cols, A.n_cols);
1949 
1950  char jobu = 'A';
1951  char jobvt = 'A';
1952 
1953  blas_int m = A.n_rows;
1954  blas_int n = A.n_cols;
1955  blas_int lda = A.n_rows;
1956  blas_int ldu = U.n_rows;
1957  blas_int ldvt = V.n_rows;
1958  blas_int lwork = 2 * (std::max)(blas_int(1), (std::max)( (3*(std::min)(m,n) + (std::max)(m,n)), 5*(std::min)(m,n) ) );
1959  blas_int info;
1960 
1961 
1962  S.set_size( static_cast<uword>((std::min)(m,n)) );
1963  podarray<eT> work( static_cast<uword>(lwork) );
1964 
1965  // let gesvd_() calculate the optimum size of the workspace
1966  blas_int lwork_tmp = -1;
1967 
1968  lapack::gesvd<eT>
1969  (
1970  &jobu, &jobvt,
1971  &m, &n,
1972  A.memptr(), &lda,
1973  S.memptr(),
1974  U.memptr(), &ldu,
1975  V.memptr(), &ldvt,
1976  work.memptr(), &lwork_tmp,
1977  &info
1978  );
1979 
1980  if(info == 0)
1981  {
1982  blas_int proposed_lwork = static_cast<blas_int>(work[0]);
1983  if(proposed_lwork > lwork)
1984  {
1985  lwork = proposed_lwork;
1986  work.set_size( static_cast<uword>(lwork) );
1987  }
1988 
1989  lapack::gesvd<eT>
1990  (
1991  &jobu, &jobvt,
1992  &m, &n,
1993  A.memptr(), &lda,
1994  S.memptr(),
1995  U.memptr(), &ldu,
1996  V.memptr(), &ldvt,
1997  work.memptr(), &lwork,
1998  &info
1999  );
2000 
2001  op_strans::apply(V,V); // op_strans will work out that an in-place transpose can be done
2002  }
2003 
2004  return (info == 0);
2005  }
2006  #else
2007  {
2008  arma_ignore(U);
2009  arma_ignore(S);
2010  arma_ignore(V);
2011  arma_ignore(X);
2012  arma_stop("svd(): use of LAPACK needs to be enabled");
2013  return false;
2014  }
2015  #endif
2016  }
2017 
2018 
2019 
2020 template<typename T, typename T1>
2021 inline
2022 bool
2023 auxlib::svd(Mat< std::complex<T> >& U, Col<T>& S, Mat< std::complex<T> >& V, const Base< std::complex<T>, T1>& X)
2024  {
2025  arma_extra_debug_sigprint();
2026 
2027  typedef std::complex<T> eT;
2028 
2029  #if defined(ARMA_USE_LAPACK)
2030  {
2031  Mat<eT> A(X.get_ref());
2032 
2033  if(A.is_empty())
2034  {
2035  U.eye(A.n_rows, A.n_rows);
2036  S.reset();
2037  V.eye(A.n_cols, A.n_cols);
2038  return true;
2039  }
2040 
2041  U.set_size(A.n_rows, A.n_rows);
2042  V.set_size(A.n_cols, A.n_cols);
2043 
2044  char jobu = 'A';
2045  char jobvt = 'A';
2046 
2047  blas_int m = A.n_rows;
2048  blas_int n = A.n_cols;
2049  blas_int lda = A.n_rows;
2050  blas_int ldu = U.n_rows;
2051  blas_int ldvt = V.n_rows;
2052  blas_int lwork = 2 * (std::max)(blas_int(1), 2*(std::min)(m,n)+(std::max)(m,n) );
2053  blas_int info;
2054 
2055  S.set_size( static_cast<uword>((std::min)(m,n)) );
2056 
2057  podarray<eT> work( static_cast<uword>(lwork) );
2058  podarray<T> rwork( static_cast<uword>(5*(std::min)(m,n)) );
2059 
2060  // let gesvd_() calculate the optimum size of the workspace
2061  blas_int lwork_tmp = -1;
2062  lapack::cx_gesvd<T>
2063  (
2064  &jobu, &jobvt,
2065  &m, &n,
2066  A.memptr(), &lda,
2067  S.memptr(),
2068  U.memptr(), &ldu,
2069  V.memptr(), &ldvt,
2070  work.memptr(), &lwork_tmp,
2071  rwork.memptr(),
2072  &info
2073  );
2074 
2075  if(info == 0)
2076  {
2077  blas_int proposed_lwork = static_cast<blas_int>(real(work[0]));
2078  if(proposed_lwork > lwork)
2079  {
2080  lwork = proposed_lwork;
2081  work.set_size( static_cast<uword>(lwork) );
2082  }
2083 
2084  lapack::cx_gesvd<T>
2085  (
2086  &jobu, &jobvt,
2087  &m, &n,
2088  A.memptr(), &lda,
2089  S.memptr(),
2090  U.memptr(), &ldu,
2091  V.memptr(), &ldvt,
2092  work.memptr(), &lwork,
2093  rwork.memptr(),
2094  &info
2095  );
2096 
2097  op_htrans::apply(V,V); // op_htrans will work out that an in-place transpose can be done
2098  }
2099 
2100  return (info == 0);
2101  }
2102  #else
2103  {
2104  arma_ignore(U);
2105  arma_ignore(S);
2106  arma_ignore(V);
2107  arma_ignore(X);
2108  arma_stop("svd(): use of LAPACK needs to be enabled");
2109  return false;
2110  }
2111  #endif
2112 
2113  }
2114 
2115 
2116 
2117 template<typename eT, typename T1>
2118 inline
2119 bool
2120 auxlib::svd_econ(Mat<eT>& U, Col<eT>& S, Mat<eT>& V, const Base<eT,T1>& X, const char mode)
2121  {
2122  arma_extra_debug_sigprint();
2123 
2124  #if defined(ARMA_USE_LAPACK)
2125  {
2126  Mat<eT> A(X.get_ref());
2127 
2128  blas_int m = A.n_rows;
2129  blas_int n = A.n_cols;
2130  blas_int lda = A.n_rows;
2131 
2132  S.set_size( static_cast<uword>((std::min)(m,n)) );
2133 
2134  blas_int ldu = 0;
2135  blas_int ldvt = 0;
2136 
2137  char jobu;
2138  char jobvt;
2139 
2140  switch(mode)
2141  {
2142  case 'l':
2143  jobu = 'S';
2144  jobvt = 'N';
2145 
2146  ldu = m;
2147  ldvt = 1;
2148 
2149  U.set_size( static_cast<uword>(ldu), static_cast<uword>((std::min)(m,n)) );
2150  V.reset();
2151 
2152  break;
2153 
2154 
2155  case 'r':
2156  jobu = 'N';
2157  jobvt = 'S';
2158 
2159  ldu = 1;
2160  ldvt = (std::min)(m,n);
2161 
2162  U.reset();
2163  V.set_size( static_cast<uword>(ldvt), static_cast<uword>(n) );
2164 
2165  break;
2166 
2167 
2168  case 'b':
2169  jobu = 'S';
2170  jobvt = 'S';
2171 
2172  ldu = m;
2173  ldvt = (std::min)(m,n);
2174 
2175  U.set_size( static_cast<uword>(ldu), static_cast<uword>((std::min)(m,n)) );
2176  V.set_size( static_cast<uword>(ldvt), static_cast<uword>(n) );
2177 
2178  break;
2179 
2180 
2181  default:
2182  U.reset();
2183  S.reset();
2184  V.reset();
2185  return false;
2186  }
2187 
2188 
2189  if(A.is_empty())
2190  {
2191  U.eye();
2192  S.reset();
2193  V.eye();
2194  return true;
2195  }
2196 
2197 
2198  blas_int lwork = 2 * (std::max)(blas_int(1), (std::max)( (3*(std::min)(m,n) + (std::max)(m,n)), 5*(std::min)(m,n) ) );
2199  blas_int info = 0;
2200 
2201 
2202  podarray<eT> work( static_cast<uword>(lwork) );
2203 
2204  // let gesvd_() calculate the optimum size of the workspace
2205  blas_int lwork_tmp = -1;
2206 
2207  lapack::gesvd<eT>
2208  (
2209  &jobu, &jobvt,
2210  &m, &n,
2211  A.memptr(), &lda,
2212  S.memptr(),
2213  U.memptr(), &ldu,
2214  V.memptr(), &ldvt,
2215  work.memptr(), &lwork_tmp,
2216  &info
2217  );
2218 
2219  if(info == 0)
2220  {
2221  blas_int proposed_lwork = static_cast<blas_int>(work[0]);
2222  if(proposed_lwork > lwork)
2223  {
2224  lwork = proposed_lwork;
2225  work.set_size( static_cast<uword>(lwork) );
2226  }
2227 
2228  lapack::gesvd<eT>
2229  (
2230  &jobu, &jobvt,
2231  &m, &n,
2232  A.memptr(), &lda,
2233  S.memptr(),
2234  U.memptr(), &ldu,
2235  V.memptr(), &ldvt,
2236  work.memptr(), &lwork,
2237  &info
2238  );
2239 
2240  op_strans::apply(V,V); // op_strans will work out that an in-place transpose can be done
2241  }
2242 
2243  return (info == 0);
2244  }
2245  #else
2246  {
2247  arma_ignore(U);
2248  arma_ignore(S);
2249  arma_ignore(V);
2250  arma_ignore(X);
2251  arma_ignore(mode);
2252  arma_stop("svd(): use of LAPACK needs to be enabled");
2253  return false;
2254  }
2255  #endif
2256  }
2257 
2258 
2259 
2260 template<typename T, typename T1>
2261 inline
2262 bool
2263 auxlib::svd_econ(Mat< std::complex<T> >& U, Col<T>& S, Mat< std::complex<T> >& V, const Base< std::complex<T>, T1>& X, const char mode)
2264  {
2265  arma_extra_debug_sigprint();
2266 
2267  typedef std::complex<T> eT;
2268 
2269  #if defined(ARMA_USE_LAPACK)
2270  {
2271  Mat<eT> A(X.get_ref());
2272 
2273  blas_int m = A.n_rows;
2274  blas_int n = A.n_cols;
2275  blas_int lda = A.n_rows;
2276 
2277  S.set_size( static_cast<uword>((std::min)(m,n)) );
2278 
2279  blas_int ldu = 0;
2280  blas_int ldvt = 0;
2281 
2282  char jobu;
2283  char jobvt;
2284 
2285  switch(mode)
2286  {
2287  case 'l':
2288  jobu = 'S';
2289  jobvt = 'N';
2290 
2291  ldu = m;
2292  ldvt = 1;
2293 
2294  U.set_size( static_cast<uword>(ldu), static_cast<uword>((std::min)(m,n)) );
2295  V.reset();
2296 
2297  break;
2298 
2299 
2300  case 'r':
2301  jobu = 'N';
2302  jobvt = 'S';
2303 
2304  ldu = 1;
2305  ldvt = (std::min)(m,n);
2306 
2307  U.reset();
2308  V.set_size( static_cast<uword>(ldvt), static_cast<uword>(n) );
2309 
2310  break;
2311 
2312 
2313  case 'b':
2314  jobu = 'S';
2315  jobvt = 'S';
2316 
2317  ldu = m;
2318  ldvt = (std::min)(m,n);
2319 
2320  U.set_size( static_cast<uword>(ldu), static_cast<uword>((std::min)(m,n)) );
2321  V.set_size( static_cast<uword>(ldvt), static_cast<uword>(n) );
2322 
2323  break;
2324 
2325 
2326  default:
2327  U.reset();
2328  S.reset();
2329  V.reset();
2330  return false;
2331  }
2332 
2333 
2334  if(A.is_empty())
2335  {
2336  U.eye();
2337  S.reset();
2338  V.eye();
2339  return true;
2340  }
2341 
2342 
2343  blas_int lwork = 2 * (std::max)(blas_int(1), (std::max)( (3*(std::min)(m,n) + (std::max)(m,n)), 5*(std::min)(m,n) ) );
2344  blas_int info = 0;
2345 
2346 
2347  podarray<eT> work( static_cast<uword>(lwork) );
2348  podarray<T> rwork( static_cast<uword>(5*(std::min)(m,n)) );
2349 
2350  // let gesvd_() calculate the optimum size of the workspace
2351  blas_int lwork_tmp = -1;
2352 
2353  lapack::cx_gesvd<T>
2354  (
2355  &jobu, &jobvt,
2356  &m, &n,
2357  A.memptr(), &lda,
2358  S.memptr(),
2359  U.memptr(), &ldu,
2360  V.memptr(), &ldvt,
2361  work.memptr(), &lwork_tmp,
2362  rwork.memptr(),
2363  &info
2364  );
2365 
2366  if(info == 0)
2367  {
2368  blas_int proposed_lwork = static_cast<blas_int>(real(work[0]));
2369  if(proposed_lwork > lwork)
2370  {
2371  lwork = proposed_lwork;
2372  work.set_size( static_cast<uword>(lwork) );
2373  }
2374 
2375  lapack::cx_gesvd<T>
2376  (
2377  &jobu, &jobvt,
2378  &m, &n,
2379  A.memptr(), &lda,
2380  S.memptr(),
2381  U.memptr(), &ldu,
2382  V.memptr(), &ldvt,
2383  work.memptr(), &lwork,
2384  rwork.memptr(),
2385  &info
2386  );
2387 
2388  op_htrans::apply(V,V); // op_strans will work out that an in-place transpose can be done
2389  }
2390 
2391  return (info == 0);
2392  }
2393  #else
2394  {
2395  arma_ignore(U);
2396  arma_ignore(S);
2397  arma_ignore(V);
2398  arma_ignore(X);
2399  arma_ignore(mode);
2400  arma_stop("svd(): use of LAPACK needs to be enabled");
2401  return false;
2402  }
2403  #endif
2404  }
2405 
2406 
2407 
2410 template<typename eT>
2411 inline
2412 bool
2413 auxlib::solve(Mat<eT>& out, Mat<eT>& A, const Mat<eT>& B, const bool slow)
2414  {
2415  arma_extra_debug_sigprint();
2416 
2417  if(A.is_empty() || B.is_empty())
2418  {
2419  out.zeros(A.n_cols, B.n_cols);
2420  return true;
2421  }
2422  else
2423  {
2424  const uword A_n_rows = A.n_rows;
2425 
2426  bool status = false;
2427 
2428  if( (A_n_rows <= 4) && (slow == false) )
2429  {
2430  Mat<eT> A_inv;
2431 
2432  status = auxlib::inv_noalias_tinymat(A_inv, A, A_n_rows);
2433 
2434  if(status == true)
2435  {
2436  out.set_size(A_n_rows, B.n_cols);
2437 
2438  gemm_emul<false,false,false,false>::apply(out, A_inv, B);
2439 
2440  return true;
2441  }
2442  }
2443 
2444  if( (A_n_rows > 4) || (status == false) )
2445  {
2446  #if defined(ARMA_USE_ATLAS)
2447  {
2448  podarray<int> ipiv(A_n_rows);
2449 
2450  out = B;
2451 
2452  int info = atlas::clapack_gesv<eT>(atlas::CblasColMajor, A_n_rows, B.n_cols, A.memptr(), A_n_rows, ipiv.memptr(), out.memptr(), A_n_rows);
2453 
2454  return (info == 0);
2455  }
2456  #elif defined(ARMA_USE_LAPACK)
2457  {
2458  blas_int n = A_n_rows;
2459  blas_int lda = A_n_rows;
2460  blas_int ldb = A_n_rows;
2461  blas_int nrhs = B.n_cols;
2462  blas_int info;
2463 
2464  podarray<blas_int> ipiv(A_n_rows);
2465 
2466  out = B;
2467 
2468  lapack::gesv<eT>(&n, &nrhs, A.memptr(), &lda, ipiv.memptr(), out.memptr(), &ldb, &info);
2469 
2470  return (info == 0);
2471  }
2472  #else
2473  {
2474  arma_stop("solve(): use of ATLAS or LAPACK needs to be enabled");
2475  return false;
2476  }
2477  #endif
2478  }
2479  }
2480 
2481  return true;
2482  }
2483 
2484 
2485 
2488 template<typename eT>
2489 inline
2490 bool
2491 auxlib::solve_od(Mat<eT>& out, Mat<eT>& A, const Mat<eT>& B)
2492  {
2493  arma_extra_debug_sigprint();
2494 
2495  #if defined(ARMA_USE_LAPACK)
2496  {
2497  if(A.is_empty() || B.is_empty())
2498  {
2499  out.zeros(A.n_cols, B.n_cols);
2500  return true;
2501  }
2502 
2503  char trans = 'N';
2504 
2505  blas_int m = A.n_rows;
2506  blas_int n = A.n_cols;
2507  blas_int lda = A.n_rows;
2508  blas_int ldb = A.n_rows;
2509  blas_int nrhs = B.n_cols;
2510  blas_int lwork = n + (std::max)(n, nrhs);
2511  blas_int info;
2512 
2513  Mat<eT> tmp = B;
2514 
2515  podarray<eT> work( static_cast<uword>(lwork) );
2516 
2517  arma_extra_debug_print("lapack::gels()");
2518 
2519  // NOTE: the dgels() function in the lapack library supplied by ATLAS 3.6 seems to have problems
2520 
2521  lapack::gels<eT>
2522  (
2523  &trans, &m, &n, &nrhs,
2524  A.memptr(), &lda,
2525  tmp.memptr(), &ldb,
2526  work.memptr(), &lwork,
2527  &info
2528  );
2529 
2530  arma_extra_debug_print("lapack::gels() -- finished");
2531 
2532  out.set_size(A.n_cols, B.n_cols);
2533 
2534  for(uword col=0; col<B.n_cols; ++col)
2535  {
2536  arrayops::copy( out.colptr(col), tmp.colptr(col), A.n_cols );
2537  }
2538 
2539  return (info == 0);
2540  }
2541  #else
2542  {
2543  arma_ignore(out);
2544  arma_ignore(A);
2545  arma_ignore(B);
2546  arma_stop("solve(): use of LAPACK needs to be enabled");
2547  return false;
2548  }
2549  #endif
2550  }
2551 
2552 
2553 
2556 template<typename eT>
2557 inline
2558 bool
2559 auxlib::solve_ud(Mat<eT>& out, Mat<eT>& A, const Mat<eT>& B)
2560  {
2561  arma_extra_debug_sigprint();
2562 
2563  #if defined(ARMA_USE_LAPACK)
2564  {
2565  if(A.is_empty() || B.is_empty())
2566  {
2567  out.zeros(A.n_cols, B.n_cols);
2568  return true;
2569  }
2570 
2571  char trans = 'N';
2572 
2573  blas_int m = A.n_rows;
2574  blas_int n = A.n_cols;
2575  blas_int lda = A.n_rows;
2576  blas_int ldb = A.n_cols;
2577  blas_int nrhs = B.n_cols;
2578  blas_int lwork = m + (std::max)(m,nrhs);
2579  blas_int info;
2580 
2581 
2582  Mat<eT> tmp;
2583  tmp.zeros(A.n_cols, B.n_cols);
2584 
2585  for(uword col=0; col<B.n_cols; ++col)
2586  {
2587  eT* tmp_colmem = tmp.colptr(col);
2588 
2589  arrayops::copy( tmp_colmem, B.colptr(col), B.n_rows );
2590 
2591  for(uword row=B.n_rows; row<A.n_cols; ++row)
2592  {
2593  tmp_colmem[row] = eT(0);
2594  }
2595  }
2596 
2597  podarray<eT> work( static_cast<uword>(lwork) );
2598 
2599  arma_extra_debug_print("lapack::gels()");
2600 
2601  // NOTE: the dgels() function in the lapack library supplied by ATLAS 3.6 seems to have problems
2602 
2603  lapack::gels<eT>
2604  (
2605  &trans, &m, &n, &nrhs,
2606  A.memptr(), &lda,
2607  tmp.memptr(), &ldb,
2608  work.memptr(), &lwork,
2609  &info
2610  );
2611 
2612  arma_extra_debug_print("lapack::gels() -- finished");
2613 
2614  out.set_size(A.n_cols, B.n_cols);
2615 
2616  for(uword col=0; col<B.n_cols; ++col)
2617  {
2618  arrayops::copy( out.colptr(col), tmp.colptr(col), A.n_cols );
2619  }
2620 
2621  return (info == 0);
2622  }
2623  #else
2624  {
2625  arma_ignore(out);
2626  arma_ignore(A);
2627  arma_ignore(B);
2628  arma_stop("solve(): use of LAPACK needs to be enabled");
2629  return false;
2630  }
2631  #endif
2632  }
2633 
2634 
2635 
2636 //
2637 // solve_tr
2638 
2639 template<typename eT>
2640 inline
2641 bool
2642 auxlib::solve_tr(Mat<eT>& out, const Mat<eT>& A, const Mat<eT>& B, const uword layout)
2643  {
2644  arma_extra_debug_sigprint();
2645 
2646  #if defined(ARMA_USE_LAPACK)
2647  {
2648  if(A.is_empty() || B.is_empty())
2649  {
2650  out.zeros(A.n_cols, B.n_cols);
2651  return true;
2652  }
2653 
2654  out = B;
2655 
2656  char uplo = (layout == 0) ? 'U' : 'L';
2657  char trans = 'N';
2658  char diag = 'N';
2659  blas_int n = blas_int(A.n_rows);
2660  blas_int nrhs = blas_int(B.n_cols);
2661  blas_int info = 0;
2662 
2663  lapack::trtrs<eT>(&uplo, &trans, &diag, &n, &nrhs, A.memptr(), &n, out.memptr(), &n, &info);
2664 
2665  return (info == 0);
2666  }
2667  #else
2668  {
2669  arma_ignore(out);
2670  arma_ignore(A);
2671  arma_ignore(B);
2672  arma_ignore(layout);
2673  arma_stop("solve(): use of LAPACK needs to be enabled");
2674  return false;
2675  }
2676  #endif
2677  }
2678 
2679 
2680 
2681 //
2682 // Schur decomposition
2683 
2684 template<typename eT>
2685 inline
2686 bool
2687 auxlib::schur_dec(Mat<eT>& Z, Mat<eT>& T, const Mat<eT>& A)
2688  {
2689  arma_extra_debug_sigprint();
2690 
2691  #if defined(ARMA_USE_LAPACK)
2692  {
2693  arma_debug_check( (A.is_square() == false), "schur_dec(): matrix A is not square" );
2694 
2695  if(A.is_empty())
2696  {
2697  Z.reset();
2698  T.reset();
2699  return true;
2700  }
2701 
2702  const uword A_n_rows = A.n_rows;
2703 
2704  char jobvs = 'V'; // get Schur vectors (Z)
2705  char sort = 'N'; // do not sort eigenvalues/vectors
2706  blas_int* select = 0; // pointer to sorting function
2707  blas_int n = blas_int(A_n_rows);
2708  blas_int sdim = 0; // output for sorting
2709 
2710  blas_int lwork = 3 * n; // workspace must be at least 3 * n (if set to -1, optimal size is output in work(0) and nothing else is done
2711 
2712  podarray<eT> work( static_cast<uword>(lwork) );
2713  podarray<blas_int> bwork(A_n_rows);
2714 
2715  blas_int info = 0;
2716 
2717  Z.set_size(A_n_rows, A_n_rows);
2718  T = A;
2719 
2720  podarray<eT> wr(A_n_rows); // output for eigenvalues
2721  podarray<eT> wi(A_n_rows); // output for eigenvalues
2722 
2723  lapack::gees(&jobvs, &sort, select, &n, T.memptr(), &n, &sdim, wr.memptr(), wi.memptr(), Z.memptr(), &n, work.memptr(), &lwork, bwork.memptr(), &info);
2724 
2725  return (info == 0);
2726  }
2727  #else
2728  {
2729  arma_ignore(Z);
2730  arma_ignore(T);
2731  arma_stop("schur_dec(): use of LAPACK needs to be enabled");
2732  return false;
2733  }
2734  #endif
2735  }
2736 
2737 
2738 
2739 template<typename cT>
2740 inline
2741 bool
2742 auxlib::schur_dec(Mat<std::complex<cT> >& Z, Mat<std::complex<cT> >& T, const Mat<std::complex<cT> >& A)
2743  {
2744  arma_extra_debug_sigprint();
2745 
2746  #if defined(ARMA_USE_LAPACK)
2747  {
2748  arma_debug_check( (A.is_square() == false), "schur_dec(): matrix A is not square" );
2749 
2750  if(A.is_empty())
2751  {
2752  Z.reset();
2753  T.reset();
2754  return true;
2755  }
2756 
2757  typedef std::complex<cT> eT;
2758 
2759  const uword A_n_rows = A.n_rows;
2760 
2761  char jobvs = 'V'; // get Schur vectors (Z)
2762  char sort = 'N'; // do not sort eigenvalues/vectors
2763  blas_int* select = 0; // pointer to sorting function
2764  blas_int n = blas_int(A_n_rows);
2765  blas_int sdim = 0; // output for sorting
2766 
2767  blas_int lwork = 3 * n; // workspace must be at least 3 * n (if set to -1, optimal size is output in work(0) and nothing else is done
2768 
2769  podarray<eT> work( static_cast<uword>(lwork) );
2770  podarray<blas_int> bwork(A_n_rows);
2771 
2772  blas_int info = 0;
2773 
2774  Z.set_size(A_n_rows, A_n_rows);
2775  T = A;
2776 
2777  podarray<eT> w(A_n_rows); // output for eigenvalues
2778  podarray<cT> rwork(A_n_rows);
2779 
2780  lapack::cx_gees(&jobvs, &sort, select, &n, T.memptr(), &n, &sdim, w.memptr(), Z.memptr(), &n, work.memptr(), &lwork, rwork.memptr(), bwork.memptr(), &info);
2781 
2782  return (info == 0);
2783  }
2784  #else
2785  {
2786  arma_ignore(Z);
2787  arma_ignore(T);
2788  arma_stop("schur_dec(): use of LAPACK needs to be enabled");
2789  return false;
2790  }
2791  #endif
2792  }
2793 
2794 
2795 
2796 //
2797 // syl (solution of the Sylvester equation AX + XB = C)
2798 
2799 template<typename eT>
2800 inline
2801 bool
2802 auxlib::syl(Mat<eT>& X, const Mat<eT>& A, const Mat<eT>& B, const Mat<eT>& C)
2803  {
2804  arma_extra_debug_sigprint();
2805 
2806  arma_debug_check( (A.is_square() == false), "syl(): matrix A is not square" );
2807  arma_debug_check( (B.is_square() == false), "syl(): matrix B is not square" );
2808 
2809  arma_debug_check( (C.n_rows != A.n_rows) || (C.n_cols != B.n_cols), "syl(): matrices are not conformant" );
2810 
2811  if(A.is_empty() || B.is_empty() || C.is_empty())
2812  {
2813  X.reset();
2814  return true;
2815  }
2816 
2817  bool status;
2818 
2819  #if defined(ARMA_USE_LAPACK)
2820  {
2821  Mat<eT> Z1, Z2, T1, T2;
2822 
2823  status = auxlib::schur_dec(Z1, T1, A);
2824  if(status == false)
2825  {
2826  return false;
2827  }
2828 
2829  status = auxlib::schur_dec(Z2, T2, B);
2830  if(status == false)
2831  {
2832  return false;
2833  }
2834 
2835  char trana = 'N';
2836  char tranb = 'N';
2837  blas_int isgn = +1;
2838  blas_int m = blas_int(T1.n_rows);
2839  blas_int n = blas_int(T2.n_cols);
2840 
2841  eT scale = eT(0);
2842  blas_int info = 0;
2843 
2844  Mat<eT> Y = trans(Z1) * C * Z2;
2845 
2846  lapack::trsyl<eT>(&trana, &tranb, &isgn, &m, &n, T1.memptr(), &m, T2.memptr(), &n, Y.memptr(), &m, &scale, &info);
2847 
2848  //Y /= scale;
2849  Y /= (-scale);
2850 
2851  X = Z1 * Y * trans(Z2);
2852 
2853  status = (info == 0);
2854  }
2855  #else
2856  {
2857  arma_stop("syl(): use of LAPACK needs to be enabled");
2858  return false;
2859  }
2860  #endif
2861 
2862 
2863  return status;
2864  }
2865 
2866 
2867 
2868 //
2869 // lyap (solution of the continuous Lyapunov equation AX + XA^H + Q = 0)
2870 
2871 template<typename eT>
2872 inline
2873 bool
2874 auxlib::lyap(Mat<eT>& X, const Mat<eT>& A, const Mat<eT>& Q)
2875  {
2876  arma_extra_debug_sigprint();
2877 
2878  arma_debug_check( (A.is_square() == false), "lyap(): matrix A is not square");
2879  arma_debug_check( (Q.is_square() == false), "lyap(): matrix Q is not square");
2880  arma_debug_check( (A.n_rows != Q.n_rows), "lyap(): matrices A and Q have different dimensions");
2881 
2882  Mat<eT> htransA;
2883  op_htrans::apply_noalias(htransA, A);
2884 
2885  const Mat<eT> mQ = -Q;
2886 
2887  return auxlib::syl(X, A, htransA, mQ);
2888  }
2889 
2890 
2891 
2892 //
2893 // dlyap (solution of the discrete Lyapunov equation AXA^H - X + Q = 0)
2894 
2895 template<typename eT>
2896 inline
2897 bool
2898 auxlib::dlyap(Mat<eT>& X, const Mat<eT>& A, const Mat<eT>& Q)
2899  {
2900  arma_extra_debug_sigprint();
2901 
2902  arma_debug_check( (A.is_square() == false), "dlyap(): matrix A is not square");
2903  arma_debug_check( (Q.is_square() == false), "dlyap(): matrix Q is not square");
2904  arma_debug_check( (A.n_rows != Q.n_rows), "dlyap(): matrices A and Q have different dimensions");
2905 
2906  const Col<eT> vecQ = reshape(Q, Q.n_elem, 1);
2907 
2908  const Mat<eT> M = eye< Mat<eT> >(Q.n_elem, Q.n_elem) - kron(conj(A), A);
2909 
2910  Col<eT> vecX;
2911 
2912  const bool status = solve(vecX, M, vecQ);
2913 
2914  if(status == true)
2915  {
2916  X = reshape(vecX, Q.n_rows, Q.n_cols);
2917  return true;
2918  }
2919  else
2920  {
2921  X.reset();
2922  return false;
2923  }
2924  }
2925 
2926 
2927 
podarray::n_elem
arma_aligned const uword n_elem
number of elements held
Definition: podarray_bones.hpp:33
access::tmp_real
static arma_inline const eT & tmp_real(const eT &X)
internal function to obtain the real part of either a plain number or a complex number ...
Definition: access.hpp:26
kron
arma_inline const Glue< T1, T2, glue_kron > kron(const Base< typename T1::elem_type, T1 > &A, const Base< typename T1::elem_type, T2 > &B)
kronecker product of two matrices, with the matrices having the same element type ...
Definition: fn_kron.hpp:26
Mat::memptr
arma_inline arma_warn_unused eT * memptr()
returns a pointer to array of eTs used by the matrix
Definition: Mat_meat.hpp:4024
blas_int
int blas_int
Definition: typedef_blas_int.hpp:23
podarray
A lightweight array for POD types. If the amount of memory requested is small, the stack is used...
Definition: podarray_bones.hpp:29
auxlib::solve_tr
static bool solve_tr(Mat< eT > &out, const Mat< eT > &A, const Mat< eT > &B, const uword layout)
Definition: auxlib_meat.hpp:2642
Mat::set_size
void set_size(const uword in_elem)
change the matrix to have user specified dimensions (data is not preserved)
Definition: Mat_meat.hpp:4211
symmatl
arma_inline const Op< T1, op_symmat > symmatl(const Base< typename T1::elem_type, T1 > &X)
Definition: fn_symmat.hpp:33
Mat::eye
const Mat & eye()
Definition: Mat_meat.hpp:4518
arrayops::copy
arma_hot static arma_inline void copy(eT *dest, const eT *src, const uword n_elem)
auxlib::det_tinymat
static eT det_tinymat(const Mat< eT > &X, const uword N)
Definition: auxlib_meat.hpp:640
Base::get_ref
arma_inline const derived & get_ref() const
Definition: Base_meat.hpp:22
arrayops::inplace_minus
static arma_hot void inplace_minus(eT *dest, const eT *src, const uword n_elem)
Definition: arrayops_meat.hpp:255
Mat::shed_rows
void shed_rows(const uword in_row1, const uword in_row2)
remove specified rows
Definition: Mat_meat.hpp:2419
auxlib::lu
static bool lu(Mat< eT > &L, Mat< eT > &U, podarray< blas_int > &ipiv, const Base< eT, T1 > &X)
immediate LU decomposition of a matrix using ATLAS or LAPACK
Definition: auxlib_meat.hpp:939
op_htrans::apply
static arma_inline void apply(Mat< eT > &out, const Mat< eT > &A, const typename arma_not_cx< eT >::result *junk=0)
Definition: op_htrans_meat.hpp:63
op_strans::apply
static void apply(Mat< eT > &out, const Mat< eT > &A)
Definition: op_strans_meat.hpp:148
auxlib::eig_gen
static bool eig_gen(Col< std::complex< T > > &eigval, Mat< T > &l_eigvec, Mat< T > &r_eigvec, const Base< T, T1 > &X, const char side)
Definition: auxlib_meat.hpp:1362
auxlib::lyap
static bool lyap(Mat< eT > &X, const Mat< eT > &A, const Mat< eT > &Q)
Definition: auxlib_meat.hpp:2874
Mat::n_cols
const uword n_cols
number of columns in the matrix (read-only)
Definition: Mat_bones.hpp:30
gemm_emul::apply
static arma_hot void apply(Mat< eT > &C, const Mat< eT > &A, const Mat< eT > &B, const eT alpha=eT(1), const eT beta=eT(0), const typename arma_not_cx< eT >::result *junk=0)
Definition: gemm.hpp:235
conj
arma_inline const T1 & conj(const Base< typename T1::pod_type, T1 > &A)
Definition: fn_elem.hpp:430
Mat::n_elem
const uword n_elem
number of elements in the matrix (read-only)
Definition: Mat_bones.hpp:31
max
arma_inline const Op< T1, op_max > max(const Base< typename T1::elem_type, T1 > &X, const uword dim=0)
Delayed &#39;maximum values&#39; operation. The dimension, along which the maxima are found, is set via &#39;dim&#39;. For dim = 0, the maximum value of each column is found (i.e. searches by traversing across rows). For dim = 1, the maximum value of each row is found (i.e. searches by traversing across columns). The default is dim = 0.
Definition: fn_max.hpp:28
Mat::n_rows
const uword n_rows
number of rows in the matrix (read-only)
Definition: Mat_bones.hpp:29
auxlib::eig_sym
static bool eig_sym(Col< eT > &eigval, const Base< eT, T1 > &X)
immediate eigenvalues of a symmetric real matrix using LAPACK
Definition: auxlib_meat.hpp:1148
arma_extra_debug_print
#define arma_extra_debug_print
Definition: debug.hpp:1118
podarray::reset
void reset()
Definition: podarray_meat.hpp:166
auxlib::solve
static bool solve(Mat< eT > &out, Mat< eT > &A, const Mat< eT > &B, const bool slow=false)
Definition: auxlib_meat.hpp:2413
auxlib::det_lapack
static eT det_lapack(const Mat< eT > &X, const bool make_copy)
immediate determinant of a matrix using ATLAS or LAPACK
Definition: auxlib_meat.hpp:738
min
arma_inline const Op< T1, op_min > min(const Base< typename T1::elem_type, T1 > &X, const uword dim=0)
Delayed &#39;minimum values&#39; operation. The dimension, along which the minima are found, is set via &#39;dim&#39;. For dim = 0, the minimum value of each column is found (i.e. searches by traversing across rows). For dim = 1, the minimum value of each row is found (i.e. searches by traversing across columns). The default is dim = 0.
Definition: fn_min.hpp:27
Mat::colptr
arma_inline arma_warn_unused eT * colptr(const uword in_col)
returns a pointer to array of eTs for a specified column; no bounds check
Definition: Mat_meat.hpp:4000
uword
u32 uword
Definition: typedef.hpp:85
auxlib::inv_inplace_lapack
static bool inv_inplace_lapack(Mat< eT > &out)
Definition: auxlib_meat.hpp:353
trimatu
arma_inline const Op< T1, op_trimat > trimatu(const Base< typename T1::elem_type, T1 > &X)
Definition: fn_trimat.hpp:21
is_float
Definition: traits.hpp:504
Col
Class for column vectors (matrices with only one column)
Definition: Col_bones.hpp:20
trans
arma_inline const Op< T1, op_htrans > trans(const Base< typename T1::elem_type, T1 > &X)
Definition: fn_trans.hpp:21
reshape
const Op< T1, op_reshape > reshape(const Base< typename T1::elem_type, T1 > &X, const uword in_n_rows, const uword in_n_cols, const uword dim=0)
Definition: fn_reshape.hpp:22
arma_debug_check
#define arma_debug_check
Definition: debug.hpp:1084
auxlib::schur_dec
static bool schur_dec(Mat< eT > &Z, Mat< eT > &T, const Mat< eT > &A)
Definition: auxlib_meat.hpp:2687
auxlib::inv_inplace_tinymat
static bool inv_inplace_tinymat(Mat< eT > &out, const uword N)
Definition: auxlib_meat.hpp:219
auxlib::inv_sym
static bool inv_sym(Mat< eT > &out, const Base< eT, T1 > &X, const uword layout)
Definition: auxlib_meat.hpp:489
symmatu
arma_inline const Op< T1, op_symmat > symmatu(const Base< typename T1::elem_type, T1 > &X)
Definition: fn_symmat.hpp:21
is_Mat
Definition: traits.hpp:39
Mat::at
arma_inline arma_warn_unused eT & at(const uword i)
linear element accessor (treats the matrix as a vector); no bounds check.
Definition: Mat_meat.hpp:3692
Mat::is_square
arma_inline arma_warn_unused bool is_square() const
returns true if the object has the same number of non-zero rows and columnns
Definition: Mat_meat.hpp:3860
arma_ignore
#define arma_ignore(variable)
Definition: compiler_setup.hpp:23
real
arma_inline const T1 & real(const Base< typename T1::pod_type, T1 > &X)
Definition: fn_elem.hpp:50
auxlib::chol
static bool chol(Mat< eT > &out, const Base< eT, T1 > &X)
Definition: auxlib_meat.hpp:1563
auxlib::inv
static bool inv(Mat< eT > &out, const Base< eT, T1 > &X, const bool slow=false)
immediate matrix inverse
Definition: auxlib_meat.hpp:26
auxlib::inv_noalias_tinymat
static bool inv_noalias_tinymat(Mat< eT > &out, const Mat< eT > &X, const uword N)
Definition: auxlib_meat.hpp:85
is_complex
Definition: traits.hpp:524
op_htrans::apply_noalias
static arma_inline void apply_noalias(Mat< eT > &out, const Mat< eT > &A, const typename arma_not_cx< eT >::result *junk=0)
Definition: op_htrans_meat.hpp:22
Mat::in_range
arma_inline arma_warn_unused bool in_range(const uword i) const
returns true if the given index is currently in range
Definition: Mat_meat.hpp:3884
Mat::reset
void reset()
Definition: Mat_meat.hpp:4553
arma_stop
void arma_cold arma_stop(const T1 &x)
print a message to get_stream_err1() and/or throw a logic_error exception
Definition: debug.hpp:98
Mat::copy_size
void copy_size(const Mat< eT2 > &m)
change the matrix (without preserving data) to have the same dimensions as the given matrix ...
Definition: Mat_meat.hpp:4303
trimatl
arma_inline const Op< T1, op_trimat > trimatl(const Base< typename T1::elem_type, T1 > &X)
Definition: fn_trimat.hpp:33
auxlib::dlyap
static bool dlyap(Mat< eT > &X, const Mat< eT > &A, const Mat< eT > &Q)
Definition: auxlib_meat.hpp:2898
arma_extra_debug_sigprint
#define arma_extra_debug_sigprint
Definition: debug.hpp:1116
auxlib::log_det
static bool log_det(eT &out_val, typename get_pod_type< eT >::result &out_sign, const Base< eT, T1 > &X)
immediate log determinant of a matrix using ATLAS or LAPACK
Definition: auxlib_meat.hpp:827
auxlib::solve_ud
static bool solve_ud(Mat< eT > &out, Mat< eT > &A, const Mat< eT > &B)
Definition: auxlib_meat.hpp:2559
sort
arma_inline const Op< T1, op_sort > sort(const Base< typename T1::elem_type, T1 > &X, const uword sort_type=0, const uword dim=0)
Definition: fn_sort.hpp:21
podarray::memptr
arma_inline eT * memptr()
Definition: podarray_meat.hpp:215
Mat
Dense matrix class.
Definition: forward_bones.hpp:19
Mat::shed_cols
void shed_cols(const uword in_col1, const uword in_col2)
remove specified columns
Definition: Mat_meat.hpp:2453
log
arma_inline const eOp< T1, eop_log > log(const Base< typename T1::elem_type, T1 > &A)
Definition: fn_elem.hpp:156
Mat::zeros
const Mat & zeros()
Definition: Mat_meat.hpp:4331
podarray::set_size
void set_size(const uword new_n_elem)
Definition: podarray_meat.hpp:154
auxlib::det
static eT det(const Base< eT, T1 > &X, const bool slow=false)
Definition: auxlib_meat.hpp:596
auxlib::inv_sympd
static bool inv_sympd(Mat< eT > &out, const Base< eT, T1 > &X, const uword layout)
Definition: auxlib_meat.hpp:546
get_pod_type::result
T1 result
Definition: traits.hpp:20
podarray::mem
arma_aligned const eT *const mem
pointer to memory used by the object
Definition: podarray_bones.hpp:34
is_double
Definition: traits.hpp:514
auxlib::solve_od
static bool solve_od(Mat< eT > &out, Mat< eT > &A, const Mat< eT > &B)
Definition: auxlib_meat.hpp:2491
auxlib::syl
static bool syl(Mat< eT > &X, const Mat< eT > &A, const Mat< eT > &B, const Mat< eT > &C)
Definition: auxlib_meat.hpp:2802
auxlib::inv_tr
static bool inv_tr(Mat< eT > &out, const Base< eT, T1 > &X, const uword layout)
Definition: auxlib_meat.hpp:432
sword
s32 sword
Definition: typedef.hpp:86
auxlib::svd
static bool svd(Col< eT > &S, const Base< eT, T1 > &X, uword &n_rows, uword &n_cols)
Definition: auxlib_meat.hpp:1715
auxlib::svd_econ
static bool svd_econ(Mat< eT > &U, Col< eT > &S, Mat< eT > &V, const Base< eT, T1 > &X, const char mode)
Definition: auxlib_meat.hpp:2120
is_supported_complex_float
Definition: traits.hpp:568
auxlib::qr
static bool qr(Mat< eT > &Q, Mat< eT > &R, const Base< eT, T1 > &X)
Definition: auxlib_meat.hpp:1612
Mat::swap_rows
void swap_rows(const uword in_row1, const uword in_row2)
Definition: Mat_meat.hpp:2331
Mat::is_empty
arma_inline arma_warn_unused bool is_empty() const
returns true if the matrix has no elements
Definition: Mat_meat.hpp:3812
is_supported_complex_double
Definition: traits.hpp:578

armadillo_matrix
Author(s):
autogenerated on Fri Apr 16 2021 02:31:56