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
dbleCHOL.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 "dRowVector.h"
31 #include "dbleCHOL.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 (dpotrf, DPOTRF) (F77_CONST_CHAR_ARG_DECL,
44  const octave_idx_type&, double*,
45  const octave_idx_type&, octave_idx_type&
47 
48  F77_RET_T
49  F77_FUNC (dpotri, DPOTRI) (F77_CONST_CHAR_ARG_DECL,
50  const octave_idx_type&, double*,
51  const octave_idx_type&, octave_idx_type&
53 
54  F77_RET_T
55  F77_FUNC (dpocon, DPOCON) (F77_CONST_CHAR_ARG_DECL,
56  const octave_idx_type&, double*,
57  const octave_idx_type&, const double&,
58  double&, double*, octave_idx_type*,
59  octave_idx_type&
61 #ifdef HAVE_QRUPDATE
62 
63  F77_RET_T
64  F77_FUNC (dch1up, DCH1UP) (const octave_idx_type&, double*,
65  const octave_idx_type&, double*, double*);
66 
67  F77_RET_T
68  F77_FUNC (dch1dn, DCH1DN) (const octave_idx_type&, double*,
69  const octave_idx_type&, double*, double*,
70  octave_idx_type&);
71 
72  F77_RET_T
73  F77_FUNC (dchinx, DCHINX) (const octave_idx_type&, double*,
74  const octave_idx_type&, const octave_idx_type&,
75  double*, double*, octave_idx_type&);
76 
77  F77_RET_T
78  F77_FUNC (dchdex, DCHDEX) (const octave_idx_type&, double*,
79  const octave_idx_type&, const octave_idx_type&,
80  double*);
81 
82  F77_RET_T
83  F77_FUNC (dchshx, DCHSHX) (const octave_idx_type&, double*,
84  const octave_idx_type&, const octave_idx_type&,
85  const octave_idx_type&, double*);
86 #endif
87 }
88 
90 CHOL::init (const Matrix& 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) ("CHOL 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.0;
111  }
112  double *h = chol_mat.fortran_vec ();
113 
114  // Calculate the norm of the matrix, for later use.
115  double anorm = 0;
116  if (calc_cond)
117  anorm = xnorm (a, 1);
118 
119  F77_XFCN (dpotrf, DPOTRF, (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 dpocon_info = 0;
129 
130  // Now calculate the condition number for non-singular matrix.
131  Array<double> z (dim_vector (3*n, 1));
132  double *pz = z.fortran_vec ();
134  octave_idx_type *piz = iz.fortran_vec ();
135  F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h,
136  n, anorm, xrcond, pz, piz, dpocon_info
137  F77_CHAR_ARG_LEN (1)));
138 
139  if (dpocon_info != 0)
140  info = -1;
141  }
142 
143  return info;
144 }
145 
146 static Matrix
148 {
149  Matrix 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  Matrix tmp = r;
160  double *v = tmp.fortran_vec ();
161 
162  if (info == 0)
163  {
164  F77_XFCN (dpotri, DPOTRI, (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.
186 Matrix
187 CHOL::inverse (void) const
188 {
189  return chol2inv_internal (chol_mat);
190 }
191 
192 void
193 CHOL::set (const Matrix& R)
194 {
195  if (R.is_square ())
196  chol_mat = R;
197  else
198  (*current_liboctave_error_handler) ("CHOL requires square matrix");
199 }
200 
201 #ifdef HAVE_QRUPDATE
202 
203 void
205 {
207 
208  if (u.length () == n)
209  {
210  ColumnVector utmp = u;
211 
212  OCTAVE_LOCAL_BUFFER (double, w, n);
213 
214  F77_XFCN (dch1up, DCH1UP, (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  ColumnVector utmp = u;
231 
232  OCTAVE_LOCAL_BUFFER (double, w, n);
233 
234  F77_XFCN (dch1dn, DCH1DN, (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  ColumnVector utmp = u;
257 
258  OCTAVE_LOCAL_BUFFER (double, w, n);
259 
260  chol_mat.resize (n+1, n+1);
261 
262  F77_XFCN (dchinx, DCHINX, (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 (double, w, n);
279 
280  F77_XFCN (dchdex, DCHDEX, (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 (double, w, 2*n);
297 
298  F77_XFCN (dchshx, DCHSHX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
299  i + 1, j + 1, w));
300  }
301 }
302 
303 #else
304 
305 void
306 CHOL::update (const ColumnVector& u)
307 {
308  warn_qrupdate_once ();
309 
311 
312  if (u.length () == n)
313  {
315  + Matrix (u) * Matrix (u).transpose (), false);
316  }
317  else
318  (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
319 }
320 
321 static bool
322 singular (const Matrix& a)
323 {
324  for (octave_idx_type i = 0; i < a.rows (); i++)
325  if (a(i,i) == 0.0) return true;
326  return false;
327 }
328 
330 CHOL::downdate (const ColumnVector& u)
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  - Matrix (u) * Matrix (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  Matrix 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 
443 Matrix
444 chol2inv (const Matrix& r)
445 {
446  return chol2inv_internal (r);
447 }
F77_RET_T const octave_idx_type double const octave_idx_type octave_idx_type & F77_CHAR_ARG_LEN_DECL
Definition: dbleCHOL.cc:44
void shift_sym(octave_idx_type i, octave_idx_type j)
Definition: dbleCHOL.cc:288
#define F77_CHAR_ARG_LEN(l)
Definition: f77-fcn.h:253
octave_idx_type downdate(const ColumnVector &u)
Definition: dbleCHOL.cc:222
void resize(octave_idx_type nr, octave_idx_type nc, double rfv=0)
Definition: dMatrix.h:130
void delete_elements(const idx_vector &i)
Deleting elements.
Definition: Array.cc:1393
void update(const ColumnVector &u)
Definition: dbleCHOL.cc:204
Matrix chol2inv(const Matrix &r)
Definition: dbleCHOL.cc:444
#define F77_XFCN(f, F, args)
Definition: f77-fcn.h:51
octave_idx_type rows(void) const
Definition: Array.h:313
void delete_sym(octave_idx_type j)
Definition: dbleCHOL.cc:270
octave_idx_type insert_sym(const ColumnVector &u, octave_idx_type j)
Definition: dbleCHOL.cc:244
#define F77_CONST_CHAR_ARG2(x, l)
Definition: f77-fcn.h:251
liboctave_error_handler current_liboctave_error_handler
Definition: lo-error.c:38
double xrcond
Definition: dbleCHOL.h:89
F77_RET_T F77_FUNC(dpotrf, DPOTRF)(F77_CONST_CHAR_ARG_DECL
octave_idx_type init(const Matrix &a, bool calc_cond)
Definition: dbleCHOL.cc:90
std::complex< double > w(std::complex< double > z, double relerr=0)
Matrix inverse(void) const
Definition: dbleCHOL.cc:187
Matrix transpose(void) const
Definition: dMatrix.h:114
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
Definition: dMatrix.h:35
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
Matrix chol_mat
Definition: dbleCHOL.h:87
#define F77_CONST_CHAR_ARG_DECL
Definition: f77-fcn.h:255
#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
void set(const Matrix &R)
Definition: dbleCHOL.cc:193
Array< T > index(const idx_vector &i) const
Indexing without resizing.
Definition: Array.cc:716
static Matrix chol2inv_internal(const Matrix &r)
Definition: dbleCHOL.cc:147