GNU Octave  4.0.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
floatCHOL.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1994-2015 John W. Eaton
4 Copyright (C) 2008-2009 Jaroslav Hajek
5 
6 This file is part of Octave.
7 
8 Octave is free software; you can redistribute it and/or modify it
9 under the terms of the GNU General Public License as published by the
10 Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12 
13 Octave is distributed in the hope that it will be useful, but WITHOUT
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
16 for more details.
17 
18 You should have received a copy of the GNU General Public License
19 along with Octave; see the file COPYING. If not, see
20 <http://www.gnu.org/licenses/>.
21 
22 */
23 
24 #ifdef HAVE_CONFIG_H
25 #include <config.h>
26 #endif
27 
28 #include <vector>
29 
30 #include "fRowVector.h"
31 #include "floatCHOL.h"
32 #include "f77-fcn.h"
33 #include "lo-error.h"
34 #include "oct-locbuf.h"
35 #include "oct-norm.h"
36 #ifndef HAVE_QRUPDATE
37 #include "dbleQR.h"
38 #endif
39 
40 extern "C"
41 {
42  F77_RET_T
43  F77_FUNC (spotrf, SPOTRF) (F77_CONST_CHAR_ARG_DECL,
44  const octave_idx_type&, float*,
45  const octave_idx_type&, octave_idx_type&
47 
48  F77_RET_T
49  F77_FUNC (spotri, SPOTRI) (F77_CONST_CHAR_ARG_DECL,
50  const octave_idx_type&, float*,
51  const octave_idx_type&, octave_idx_type&
53 
54  F77_RET_T
55  F77_FUNC (spocon, SPOCON) (F77_CONST_CHAR_ARG_DECL,
56  const octave_idx_type&, float*,
57  const octave_idx_type&, const float&,
58  float&, float*, octave_idx_type*,
59  octave_idx_type&
61 #ifdef HAVE_QRUPDATE
62 
63  F77_RET_T
64  F77_FUNC (sch1up, SCH1UP) (const octave_idx_type&, float*,
65  const octave_idx_type&, float*, float*);
66 
67  F77_RET_T
68  F77_FUNC (sch1dn, SCH1DN) (const octave_idx_type&, float*,
69  const octave_idx_type&, float*, float*,
70  octave_idx_type&);
71 
72  F77_RET_T
73  F77_FUNC (schinx, SCHINX) (const octave_idx_type&, float*,
74  const octave_idx_type&, const octave_idx_type&,
75  float*, float*, octave_idx_type&);
76 
77  F77_RET_T
78  F77_FUNC (schdex, SCHDEX) (const octave_idx_type&, float*,
79  const octave_idx_type&, const octave_idx_type&,
80  float*);
81 
82  F77_RET_T
83  F77_FUNC (schshx, SCHSHX) (const octave_idx_type&, float*,
84  const octave_idx_type&, const octave_idx_type&,
85  const octave_idx_type&, float*);
86 #endif
87 }
88 
90 FloatCHOL::init (const FloatMatrix& a, bool calc_cond)
91 {
92  octave_idx_type a_nr = a.rows ();
93  octave_idx_type a_nc = a.cols ();
94 
95  if (a_nr != a_nc)
96  {
97  (*current_liboctave_error_handler) ("FloatCHOL requires square matrix");
98  return -1;
99  }
100 
101  octave_idx_type n = a_nc;
102  octave_idx_type info;
103 
104  chol_mat.clear (n, n);
105  for (octave_idx_type j = 0; j < n; j++)
106  {
107  for (octave_idx_type i = 0; i <= j; i++)
108  chol_mat.xelem (i, j) = a(i, j);
109  for (octave_idx_type i = j+1; i < n; i++)
110  chol_mat.xelem (i, j) = 0.0f;
111  }
112  float *h = chol_mat.fortran_vec ();
113 
114  // Calculate the norm of the matrix, for later use.
115  float anorm = 0;
116  if (calc_cond)
117  anorm = xnorm (a, 1);
118 
119  F77_XFCN (spotrf, SPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1),
120  n, h, n, info
121  F77_CHAR_ARG_LEN (1)));
122 
123  xrcond = 0.0;
124  if (info > 0)
125  chol_mat.resize (info - 1, info - 1);
126  else if (calc_cond)
127  {
128  octave_idx_type spocon_info = 0;
129 
130  // Now calculate the condition number for non-singular matrix.
131  Array<float> z (dim_vector (3*n, 1));
132  float *pz = z.fortran_vec ();
134  octave_idx_type *piz = iz.fortran_vec ();
135  F77_XFCN (spocon, SPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h,
136  n, anorm, xrcond, pz, piz, spocon_info
137  F77_CHAR_ARG_LEN (1)));
138 
139  if (spocon_info != 0)
140  info = -1;
141  }
142 
143  return info;
144 }
145 
146 static FloatMatrix
148 {
149  FloatMatrix retval;
150 
151  octave_idx_type r_nr = r.rows ();
152  octave_idx_type r_nc = r.cols ();
153 
154  if (r_nr == r_nc)
155  {
156  octave_idx_type n = r_nc;
157  octave_idx_type info = 0;
158 
159  FloatMatrix tmp = r;
160  float *v = tmp.fortran_vec ();
161 
162  if (info == 0)
163  {
164  F77_XFCN (spotri, SPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n,
165  v, n, info
166  F77_CHAR_ARG_LEN (1)));
167 
168  // If someone thinks of a more graceful way of doing this (or
169  // faster for that matter :-)), please let me know!
170 
171  if (n > 1)
172  for (octave_idx_type j = 0; j < r_nc; j++)
173  for (octave_idx_type i = j+1; i < r_nr; i++)
174  tmp.xelem (i, j) = tmp.xelem (j, i);
175 
176  retval = tmp;
177  }
178  }
179  else
180  (*current_liboctave_error_handler) ("chol2inv requires square matrix");
181 
182  return retval;
183 }
184 
185 // Compute the inverse of a matrix using the Cholesky factorization.
187 FloatCHOL::inverse (void) const
188 {
189  return chol2inv_internal (chol_mat);
190 }
191 
192 void
194 {
195  if (R.is_square ())
196  chol_mat = R;
197  else
198  (*current_liboctave_error_handler) ("FloatCHOL requires square matrix");
199 }
200 
201 #ifdef HAVE_QRUPDATE
202 
203 void
205 {
207 
208  if (u.length () == n)
209  {
210  FloatColumnVector utmp = u;
211 
212  OCTAVE_LOCAL_BUFFER (float, w, n);
213 
214  F77_XFCN (sch1up, SCH1UP, (n, chol_mat.fortran_vec (), chol_mat.rows (),
215  utmp.fortran_vec (), w));
216  }
217  else
218  (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
219 }
220 
223 {
224  octave_idx_type info = -1;
225 
227 
228  if (u.length () == n)
229  {
230  FloatColumnVector utmp = u;
231 
232  OCTAVE_LOCAL_BUFFER (float, w, n);
233 
234  F77_XFCN (sch1dn, SCH1DN, (n, chol_mat.fortran_vec (), chol_mat.rows (),
235  utmp.fortran_vec (), w, info));
236  }
237  else
238  (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
239 
240  return info;
241 }
242 
245 {
246  octave_idx_type info = -1;
247 
249 
250  if (u.length () != n + 1)
251  (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
252  else if (j < 0 || j > n)
253  (*current_liboctave_error_handler) ("cholinsert: index out of range");
254  else
255  {
256  FloatColumnVector utmp = u;
257 
258  OCTAVE_LOCAL_BUFFER (float, w, n);
259 
260  chol_mat.resize (n+1, n+1);
261 
262  F77_XFCN (schinx, SCHINX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
263  j + 1, utmp.fortran_vec (), w, info));
264  }
265 
266  return info;
267 }
268 
269 void
271 {
273 
274  if (j < 0 || j > n-1)
275  (*current_liboctave_error_handler) ("choldelete: index out of range");
276  else
277  {
278  OCTAVE_LOCAL_BUFFER (float, w, n);
279 
280  F77_XFCN (schdex, SCHDEX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
281  j + 1, w));
282 
283  chol_mat.resize (n-1, n-1);
284  }
285 }
286 
287 void
289 {
291 
292  if (i < 0 || i > n-1 || j < 0 || j > n-1)
293  (*current_liboctave_error_handler) ("cholshift: index out of range");
294  else
295  {
296  OCTAVE_LOCAL_BUFFER (float, w, 2*n);
297 
298  F77_XFCN (schshx, SCHSHX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
299  i + 1, j + 1, w));
300  }
301 }
302 
303 #else
304 
305 void
307 {
308  warn_qrupdate_once ();
309 
311 
312  if (u.length () == n)
313  {
315  + FloatMatrix (u) * FloatMatrix (u).transpose (), false);
316  }
317  else
318  (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
319 }
320 
321 static bool
322 singular (const FloatMatrix& a)
323 {
324  for (octave_idx_type i = 0; i < a.rows (); i++)
325  if (a(i,i) == 0.0f) return true;
326  return false;
327 }
328 
331 {
332  warn_qrupdate_once ();
333 
334  octave_idx_type info = -1;
335 
337 
338  if (u.length () == n)
339  {
340  if (singular (chol_mat))
341  info = 2;
342  else
343  {
344  info = init (chol_mat.transpose () * chol_mat
345  - FloatMatrix (u) * FloatMatrix (u).transpose (), false);
346  if (info) info = 1;
347  }
348  }
349  else
350  (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
351 
352  return info;
353 }
354 
357 {
358  warn_qrupdate_once ();
359 
360  octave_idx_type info = -1;
361 
363 
364  if (u.length () != n + 1)
365  (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
366  else if (j < 0 || j > n)
367  (*current_liboctave_error_handler) ("cholinsert: index out of range");
368  else
369  {
370  if (singular (chol_mat))
371  info = 2;
372  else
373  {
375  FloatMatrix a1 (n+1, n+1);
376  for (octave_idx_type k = 0; k < n+1; k++)
377  for (octave_idx_type l = 0; l < n+1; l++)
378  {
379  if (l == j)
380  a1(k, l) = u(k);
381  else if (k == j)
382  a1(k, l) = u(l);
383  else
384  a1(k, l) = a(k < j ? k : k-1, l < j ? l : l-1);
385  }
386  info = init (a1, false);
387  if (info) info = 1;
388  }
389  }
390 
391  return info;
392 }
393 
394 void
396 {
397  warn_qrupdate_once ();
398 
400 
401  if (j < 0 || j > n-1)
402  (*current_liboctave_error_handler) ("choldelete: index out of range");
403  else
404  {
406  a.delete_elements (1, idx_vector (j));
407  a.delete_elements (0, idx_vector (j));
408  init (a, false);
409  }
410 }
411 
412 void
414 {
415  warn_qrupdate_once ();
416 
418 
419  if (i < 0 || i > n-1 || j < 0 || j > n-1)
420  (*current_liboctave_error_handler) ("cholshift: index out of range");
421  else
422  {
425  for (octave_idx_type k = 0; k < n; k++) p(k) = k;
426  if (i < j)
427  {
428  for (octave_idx_type k = i; k < j; k++) p(k) = k+1;
429  p(j) = i;
430  }
431  else if (j < i)
432  {
433  p(j) = i;
434  for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1;
435  }
436 
437  init (a.index (idx_vector (p), idx_vector (p)), false);
438  }
439 }
440 
441 #endif
442 
445 {
446  return chol2inv_internal (r);
447 }
octave_idx_type init(const FloatMatrix &a, bool calc_cond)
Definition: floatCHOL.cc:90
#define F77_CHAR_ARG_LEN(l)
Definition: f77-fcn.h:253
void resize(octave_idx_type nr, octave_idx_type nc, float rfv=0)
Definition: fMatrix.h:134
void delete_elements(const idx_vector &i)
Deleting elements.
Definition: Array.cc:1393
FloatMatrix transpose(void) const
Definition: fMatrix.h:118
void shift_sym(octave_idx_type i, octave_idx_type j)
Definition: floatCHOL.cc:288
float xrcond
Definition: floatCHOL.h:91
#define F77_XFCN(f, F, args)
Definition: f77-fcn.h:51
octave_idx_type rows(void) const
Definition: Array.h:313
F77_RET_T const octave_idx_type float const octave_idx_type octave_idx_type & F77_CHAR_ARG_LEN_DECL
Definition: floatCHOL.cc:44
static FloatMatrix chol2inv_internal(const FloatMatrix &r)
Definition: floatCHOL.cc:147
#define F77_CONST_CHAR_ARG2(x, l)
Definition: f77-fcn.h:251
F77_RET_T F77_FUNC(spotrf, SPOTRF)(F77_CONST_CHAR_ARG_DECL
liboctave_error_handler current_liboctave_error_handler
Definition: lo-error.c:38
void update(const FloatColumnVector &u)
Definition: floatCHOL.cc:204
octave_idx_type downdate(const FloatColumnVector &u)
Definition: floatCHOL.cc:222
FloatMatrix chol2inv(const FloatMatrix &r)
Definition: floatCHOL.cc:444
std::complex< double > w(std::complex< double > z, double relerr=0)
octave_idx_type insert_sym(const FloatColumnVector &u, octave_idx_type j)
Definition: floatCHOL.cc:244
bool is_square(void) const
Definition: Array.h:470
OCTAVE_API double xnorm(const ColumnVector &x, double p)
Definition: oct-norm.cc:536
#define F77_RET_T
Definition: f77-fcn.h:264
T & xelem(octave_idx_type n)
Definition: Array.h:353
octave_idx_type length(void) const
Number of elements in the array.
Definition: Array.h:267
void clear(void)
Definition: Array.cc:84
#define F77_CONST_CHAR_ARG_DECL
Definition: f77-fcn.h:255
FloatMatrix inverse(void) const
Definition: floatCHOL.cc:187
void delete_sym(octave_idx_type j)
Definition: floatCHOL.cc:270
void set(const FloatMatrix &R)
Definition: floatCHOL.cc:193
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition: oct-locbuf.h:197
const T * fortran_vec(void) const
Definition: Array.h:481
octave_idx_type cols(void) const
Definition: Array.h:321
FloatMatrix chol_mat
Definition: floatCHOL.h:89
Array< T > index(const idx_vector &i) const
Indexing without resizing.
Definition: Array.cc:716