GNU Octave  4.2.1
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
chol.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1994-2017 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 #if defined (HAVE_CONFIG_H)
25 # include "config.h"
26 #endif
27 
28 #include <vector>
29 
30 #include "CColVector.h"
31 #include "CMatrix.h"
32 #include "CRowVector.h"
33 #include "chol.h"
34 #include "dColVector.h"
35 #include "dMatrix.h"
36 #include "dRowVector.h"
37 #include "fCColVector.h"
38 #include "fCMatrix.h"
39 #include "fCRowVector.h"
40 #include "fColVector.h"
41 #include "fMatrix.h"
42 #include "fRowVector.h"
43 #include "lo-error.h"
44 #include "lo-lapack-proto.h"
45 #include "lo-qrupdate-proto.h"
46 #include "oct-locbuf.h"
47 #include "oct-norm.h"
48 
49 #if ! defined (HAVE_QRUPDATE)
50 # include "qr.h"
51 #endif
52 
53 static Matrix
54 chol2inv_internal (const Matrix& r, bool is_upper = true)
55 {
56  Matrix retval;
57 
58  octave_idx_type r_nr = r.rows ();
59  octave_idx_type r_nc = r.cols ();
60 
61  if (r_nr != r_nc)
62  (*current_liboctave_error_handler) ("chol2inv requires square matrix");
63 
64  octave_idx_type n = r_nc;
65  octave_idx_type info = 0;
66 
67  Matrix tmp = r;
68  double *v = tmp.fortran_vec ();
69 
70  if (info == 0)
71  {
72  if (is_upper)
73  F77_XFCN (dpotri, DPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n,
74  v, n, info
75  F77_CHAR_ARG_LEN (1)));
76  else
77  F77_XFCN (dpotri, DPOTRI, (F77_CONST_CHAR_ARG2 ("L", 1), n,
78  v, n, info
79  F77_CHAR_ARG_LEN (1)));
80 
81  // If someone thinks of a more graceful way of doing this (or
82  // faster for that matter :-)), please let me know!
83 
84  if (n > 1)
85  {
86  if (is_upper)
87  for (octave_idx_type j = 0; j < r_nc; j++)
88  for (octave_idx_type i = j+1; i < r_nr; i++)
89  tmp.xelem (i, j) = tmp.xelem (j, i);
90  else
91  for (octave_idx_type j = 0; j < r_nc; j++)
92  for (octave_idx_type i = j+1; i < r_nr; i++)
93  tmp.xelem (j, i) = tmp.xelem (i, j);
94  }
95 
96  retval = tmp;
97  }
98 
99  return retval;
100 }
101 
102 static FloatMatrix
103 chol2inv_internal (const FloatMatrix& r, bool is_upper = true)
104 {
105  FloatMatrix retval;
106 
107  octave_idx_type r_nr = r.rows ();
108  octave_idx_type r_nc = r.cols ();
109 
110  if (r_nr != r_nc)
111  (*current_liboctave_error_handler) ("chol2inv requires square matrix");
112 
113  octave_idx_type n = r_nc;
114  octave_idx_type info = 0;
115 
116  FloatMatrix tmp = r;
117  float *v = tmp.fortran_vec ();
118 
119  if (info == 0)
120  {
121  if (is_upper)
122  F77_XFCN (spotri, SPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n,
123  v, n, info
124  F77_CHAR_ARG_LEN (1)));
125  else
126  F77_XFCN (spotri, SPOTRI, (F77_CONST_CHAR_ARG2 ("L", 1), n,
127  v, n, info
128  F77_CHAR_ARG_LEN (1)));
129 
130  // If someone thinks of a more graceful way of doing this (or
131  // faster for that matter :-)), please let me know!
132 
133  if (n > 1)
134  {
135  if (is_upper)
136  for (octave_idx_type j = 0; j < r_nc; j++)
137  for (octave_idx_type i = j+1; i < r_nr; i++)
138  tmp.xelem (i, j) = tmp.xelem (j, i);
139  else
140  for (octave_idx_type j = 0; j < r_nc; j++)
141  for (octave_idx_type i = j+1; i < r_nr; i++)
142  tmp.xelem (j, i) = tmp.xelem (i, j);
143  }
144 
145  retval = tmp;
146  }
147 
148  return retval;
149 }
150 
151 static ComplexMatrix
152 chol2inv_internal (const ComplexMatrix& r, bool is_upper = true)
153 {
154  ComplexMatrix retval;
155 
156  octave_idx_type r_nr = r.rows ();
157  octave_idx_type r_nc = r.cols ();
158 
159  if (r_nr != r_nc)
160  (*current_liboctave_error_handler) ("chol2inv requires square matrix");
161 
162  octave_idx_type n = r_nc;
163  octave_idx_type info;
164 
165  ComplexMatrix tmp = r;
166 
167  if (is_upper)
168  F77_XFCN (zpotri, ZPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n,
169  F77_DBLE_CMPLX_ARG (tmp.fortran_vec ()), n, info
170  F77_CHAR_ARG_LEN (1)));
171  else
172  F77_XFCN (zpotri, ZPOTRI, (F77_CONST_CHAR_ARG2 ("L", 1), n,
173  F77_DBLE_CMPLX_ARG (tmp.fortran_vec ()), n, info
174  F77_CHAR_ARG_LEN (1)));
175 
176  // If someone thinks of a more graceful way of doing this (or
177  // faster for that matter :-)), please let me know!
178 
179  if (n > 1)
180  {
181  if (is_upper)
182  for (octave_idx_type j = 0; j < r_nc; j++)
183  for (octave_idx_type i = j+1; i < r_nr; i++)
184  tmp.xelem (i, j) = std::conj (tmp.xelem (j, i));
185  else
186  for (octave_idx_type j = 0; j < r_nc; j++)
187  for (octave_idx_type i = j+1; i < r_nr; i++)
188  tmp.xelem (j, i) = std::conj (tmp.xelem (i, j));
189  }
190 
191  retval = tmp;
192 
193  return retval;
194 }
195 
196 static FloatComplexMatrix
197 chol2inv_internal (const FloatComplexMatrix& r, bool is_upper = true)
198 {
199  FloatComplexMatrix retval;
200 
201  octave_idx_type r_nr = r.rows ();
202  octave_idx_type r_nc = r.cols ();
203 
204  if (r_nr != r_nc)
205  (*current_liboctave_error_handler) ("chol2inv requires square matrix");
206 
207  octave_idx_type n = r_nc;
208  octave_idx_type info;
209 
210  FloatComplexMatrix tmp = r;
211 
212  if (is_upper)
213  F77_XFCN (cpotri, CPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n,
214  F77_CMPLX_ARG (tmp.fortran_vec ()), n, info
215  F77_CHAR_ARG_LEN (1)));
216  else
217  F77_XFCN (cpotri, CPOTRI, (F77_CONST_CHAR_ARG2 ("L", 1), n,
218  F77_CMPLX_ARG (tmp.fortran_vec ()), n, info
219  F77_CHAR_ARG_LEN (1)));
220 
221  // If someone thinks of a more graceful way of doing this (or
222  // faster for that matter :-)), please let me know!
223 
224  if (n > 1)
225  {
226  if (is_upper)
227  for (octave_idx_type j = 0; j < r_nc; j++)
228  for (octave_idx_type i = j+1; i < r_nr; i++)
229  tmp.xelem (i, j) = std::conj (tmp.xelem (j, i));
230  else
231  for (octave_idx_type j = 0; j < r_nc; j++)
232  for (octave_idx_type i = j+1; i < r_nr; i++)
233  tmp.xelem (j, i) = std::conj (tmp.xelem (i, j));
234  }
235 
236  retval = tmp;
237 
238  return retval;
239 }
240 
241 namespace octave
242 {
243  namespace math
244  {
245  template <typename T>
246  T
247  chol2inv (const T& r)
248  {
249  return chol2inv_internal (r);
250  }
251 
252  // Compute the inverse of a matrix using the Cholesky factorization.
253  template <typename T>
254  T
255  chol<T>::inverse (void) const
256  {
257  return chol2inv_internal (chol_mat, is_upper);
258  }
259 
260  template <typename T>
261  void
262  chol<T>::set (const T& R)
263  {
264  if (! R.is_square ())
265  (*current_liboctave_error_handler) ("chol: requires square matrix");
266 
267  chol_mat = R;
268  }
269 
270 #if ! defined (HAVE_QRUPDATE)
271 
272  template <typename T>
273  void
274  chol<T>::update (const VT& u)
275  {
276  warn_qrupdate_once ();
277 
278  octave_idx_type n = chol_mat.rows ();
279 
280  if (u.numel () != n)
281  (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
282 
283  init (chol_mat.hermitian () * chol_mat + T (u) * T (u).hermitian (),
284  true, false);
285  }
286 
287  template <typename T>
288  bool
289  singular (const T& a)
290  {
291  static typename T::element_type zero (0);
292  for (octave_idx_type i = 0; i < a.rows (); i++)
293  if (a(i,i) == zero) return true;
294  return false;
295  }
296 
297  template <typename T>
298  octave_idx_type
299  chol<T>::downdate (const VT& u)
300  {
301  warn_qrupdate_once ();
302 
303  octave_idx_type info = -1;
304 
305  octave_idx_type n = chol_mat.rows ();
306 
307  if (u.numel () != n)
308  (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
309 
310  if (singular (chol_mat))
311  info = 2;
312  else
313  {
314  info = init (chol_mat.hermitian () * chol_mat
315  - T (u) * T (u).hermitian (), true, false);
316  if (info) info = 1;
317  }
318 
319  return info;
320  }
321 
322  template <typename T>
323  octave_idx_type
324  chol<T>::insert_sym (const VT& u, octave_idx_type j)
325  {
326  static typename T::element_type zero (0);
327 
328  warn_qrupdate_once ();
329 
330  octave_idx_type info = -1;
331 
332  octave_idx_type n = chol_mat.rows ();
333 
334  if (u.numel () != n + 1)
335  (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
336  if (j < 0 || j > n)
337  (*current_liboctave_error_handler) ("cholinsert: index out of range");
338 
339  if (singular (chol_mat))
340  info = 2;
341  else if (octave::math::imag (u(j)) != zero)
342  info = 3;
343  else
344  {
345  T a = chol_mat.hermitian () * chol_mat;
346  T a1 (n+1, n+1);
347  for (octave_idx_type k = 0; k < n+1; k++)
348  for (octave_idx_type l = 0; l < n+1; l++)
349  {
350  if (l == j)
351  a1(k, l) = u(k);
352  else if (k == j)
353  a1(k, l) = octave::math::conj (u(l));
354  else
355  a1(k, l) = a(k < j ? k : k-1, l < j ? l : l-1);
356  }
357  info = init (a1, true, false);
358  if (info) info = 1;
359  }
360 
361  return info;
362  }
363 
364  template <typename T>
365  void
366  chol<T>::delete_sym (octave_idx_type j)
367  {
368  warn_qrupdate_once ();
369 
370  octave_idx_type n = chol_mat.rows ();
371 
372  if (j < 0 || j > n-1)
373  (*current_liboctave_error_handler) ("choldelete: index out of range");
374 
375  T a = chol_mat.hermitian () * chol_mat;
376  a.delete_elements (1, idx_vector (j));
377  a.delete_elements (0, idx_vector (j));
378  init (a, true, false);
379  }
380 
381  template <typename T>
382  void
383  chol<T>::shift_sym (octave_idx_type i, octave_idx_type j)
384  {
385  warn_qrupdate_once ();
386 
387  octave_idx_type n = chol_mat.rows ();
388 
389  if (i < 0 || i > n-1 || j < 0 || j > n-1)
390  (*current_liboctave_error_handler) ("cholshift: index out of range");
391 
392  T a = chol_mat.hermitian () * chol_mat;
393  Array<octave_idx_type> p (dim_vector (n, 1));
394  for (octave_idx_type k = 0; k < n; k++) p(k) = k;
395  if (i < j)
396  {
397  for (octave_idx_type k = i; k < j; k++) p(k) = k+1;
398  p(j) = i;
399  }
400  else if (j < i)
401  {
402  p(j) = i;
403  for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1;
404  }
405 
406  init (a.index (idx_vector (p), idx_vector (p)), true, false);
407  }
408 
409 #endif
410 
411  // Specializations.
412 
413  template <>
414  octave_idx_type
415  chol<Matrix>::init (const Matrix& a, bool upper, bool calc_cond)
416  {
417  octave_idx_type a_nr = a.rows ();
418  octave_idx_type a_nc = a.cols ();
419 
420  if (a_nr != a_nc)
421  (*current_liboctave_error_handler) ("chol: requires square matrix");
422 
423  octave_idx_type n = a_nc;
424  octave_idx_type info;
425 
426  is_upper = upper;
427 
428  chol_mat.clear (n, n);
429  if (is_upper)
430  for (octave_idx_type j = 0; j < n; j++)
431  {
432  for (octave_idx_type i = 0; i <= j; i++)
433  chol_mat.xelem (i, j) = a(i, j);
434  for (octave_idx_type i = j+1; i < n; i++)
435  chol_mat.xelem (i, j) = 0.0;
436  }
437  else
438  for (octave_idx_type j = 0; j < n; j++)
439  {
440  for (octave_idx_type i = 0; i < j; i++)
441  chol_mat.xelem (i, j) = 0.0;
442  for (octave_idx_type i = j; i < n; i++)
443  chol_mat.xelem (i, j) = a(i, j);
444  }
445  double *h = chol_mat.fortran_vec ();
446 
447  // Calculate the norm of the matrix, for later use.
448  double anorm = 0;
449  if (calc_cond)
450  anorm = xnorm (a, 1);
451 
452  if (is_upper)
453  F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1), n, h, n, info
454  F77_CHAR_ARG_LEN (1)));
455  else
456  F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), n, h, n, info
457  F77_CHAR_ARG_LEN (1)));
458 
459  xrcond = 0.0;
460  if (info > 0)
461  chol_mat.resize (info - 1, info - 1);
462  else if (calc_cond)
463  {
464  octave_idx_type dpocon_info = 0;
465 
466  // Now calculate the condition number for non-singular matrix.
467  Array<double> z (dim_vector (3*n, 1));
468  double *pz = z.fortran_vec ();
469  Array<octave_idx_type> iz (dim_vector (n, 1));
470  octave_idx_type *piz = iz.fortran_vec ();
471  if (is_upper)
472  F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h,
473  n, anorm, xrcond, pz, piz, dpocon_info
474  F77_CHAR_ARG_LEN (1)));
475  else
476  F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 ("L", 1), n, h,
477  n, anorm, xrcond, pz, piz, dpocon_info
478  F77_CHAR_ARG_LEN (1)));
479 
480  if (dpocon_info != 0)
481  info = -1;
482  }
483 
484  return info;
485  }
486 
487 #if defined (HAVE_QRUPDATE)
488 
489  template <>
490  void
491  chol<Matrix>::update (const ColumnVector& u)
492  {
493  octave_idx_type n = chol_mat.rows ();
494 
495  if (u.numel () != n)
496  (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
497 
498  ColumnVector utmp = u;
499 
500  OCTAVE_LOCAL_BUFFER (double, w, n);
501 
502  F77_XFCN (dch1up, DCH1UP, (n, chol_mat.fortran_vec (), chol_mat.rows (),
503  utmp.fortran_vec (), w));
504  }
505 
506  template <>
507  octave_idx_type
508  chol<Matrix>::downdate (const ColumnVector& u)
509  {
510  octave_idx_type info = -1;
511 
512  octave_idx_type n = chol_mat.rows ();
513 
514  if (u.numel () != n)
515  (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
516 
517  ColumnVector utmp = u;
518 
519  OCTAVE_LOCAL_BUFFER (double, w, n);
520 
521  F77_XFCN (dch1dn, DCH1DN, (n, chol_mat.fortran_vec (), chol_mat.rows (),
522  utmp.fortran_vec (), w, info));
523 
524  return info;
525  }
526 
527  template <>
528  octave_idx_type
529  chol<Matrix>::insert_sym (const ColumnVector& u, octave_idx_type j)
530  {
531  octave_idx_type info = -1;
532 
533  octave_idx_type n = chol_mat.rows ();
534 
535  if (u.numel () != n + 1)
536  (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
537  if (j < 0 || j > n)
538  (*current_liboctave_error_handler) ("cholinsert: index out of range");
539 
540  ColumnVector utmp = u;
541 
542  OCTAVE_LOCAL_BUFFER (double, w, n);
543 
544  chol_mat.resize (n+1, n+1);
545 
546  F77_XFCN (dchinx, DCHINX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
547  j + 1, utmp.fortran_vec (), w, info));
548 
549  return info;
550  }
551 
552  template <>
553  void
554  chol<Matrix>::delete_sym (octave_idx_type j)
555  {
556  octave_idx_type n = chol_mat.rows ();
557 
558  if (j < 0 || j > n-1)
559  (*current_liboctave_error_handler) ("choldelete: index out of range");
560 
561  OCTAVE_LOCAL_BUFFER (double, w, n);
562 
563  F77_XFCN (dchdex, DCHDEX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
564  j + 1, w));
565 
566  chol_mat.resize (n-1, n-1);
567  }
568 
569  template <>
570  void
571  chol<Matrix>::shift_sym (octave_idx_type i, octave_idx_type j)
572  {
573  octave_idx_type n = chol_mat.rows ();
574 
575  if (i < 0 || i > n-1 || j < 0 || j > n-1)
576  (*current_liboctave_error_handler) ("cholshift: index out of range");
577 
578  OCTAVE_LOCAL_BUFFER (double, w, 2*n);
579 
580  F77_XFCN (dchshx, DCHSHX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
581  i + 1, j + 1, w));
582  }
583 
584 #endif
585 
586  template <>
587  octave_idx_type
588  chol<FloatMatrix>::init (const FloatMatrix& a, bool upper, bool calc_cond)
589  {
590  octave_idx_type a_nr = a.rows ();
591  octave_idx_type a_nc = a.cols ();
592 
593  if (a_nr != a_nc)
594  (*current_liboctave_error_handler) ("chol: requires square matrix");
595 
596  octave_idx_type n = a_nc;
597  octave_idx_type info;
598 
599  is_upper = upper;
600 
601  chol_mat.clear (n, n);
602  if (is_upper)
603  for (octave_idx_type j = 0; j < n; j++)
604  {
605  for (octave_idx_type i = 0; i <= j; i++)
606  chol_mat.xelem (i, j) = a(i, j);
607  for (octave_idx_type i = j+1; i < n; i++)
608  chol_mat.xelem (i, j) = 0.0f;
609  }
610  else
611  for (octave_idx_type j = 0; j < n; j++)
612  {
613  for (octave_idx_type i = 0; i < j; i++)
614  chol_mat.xelem (i, j) = 0.0f;
615  for (octave_idx_type i = j; i < n; i++)
616  chol_mat.xelem (i, j) = a(i, j);
617  }
618  float *h = chol_mat.fortran_vec ();
619 
620  // Calculate the norm of the matrix, for later use.
621  float anorm = 0;
622  if (calc_cond)
623  anorm = xnorm (a, 1);
624 
625  if (is_upper)
626  F77_XFCN (spotrf, SPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1), n, h, n, info
627  F77_CHAR_ARG_LEN (1)));
628  else
629  F77_XFCN (spotrf, SPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), n, h, n, info
630  F77_CHAR_ARG_LEN (1)));
631 
632  xrcond = 0.0;
633  if (info > 0)
634  chol_mat.resize (info - 1, info - 1);
635  else if (calc_cond)
636  {
637  octave_idx_type spocon_info = 0;
638 
639  // Now calculate the condition number for non-singular matrix.
640  Array<float> z (dim_vector (3*n, 1));
641  float *pz = z.fortran_vec ();
642  Array<octave_idx_type> iz (dim_vector (n, 1));
643  octave_idx_type *piz = iz.fortran_vec ();
644  if (is_upper)
645  F77_XFCN (spocon, SPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h,
646  n, anorm, xrcond, pz, piz, spocon_info
647  F77_CHAR_ARG_LEN (1)));
648  else
649  F77_XFCN (spocon, SPOCON, (F77_CONST_CHAR_ARG2 ("L", 1), n, h,
650  n, anorm, xrcond, pz, piz, spocon_info
651  F77_CHAR_ARG_LEN (1)));
652 
653  if (spocon_info != 0)
654  info = -1;
655  }
656 
657  return info;
658  }
659 
660 #if defined (HAVE_QRUPDATE)
661 
662  template <>
663  void
664  chol<FloatMatrix>::update (const FloatColumnVector& u)
665  {
666  octave_idx_type n = chol_mat.rows ();
667 
668  if (u.numel () != n)
669  (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
670 
671  FloatColumnVector utmp = u;
672 
673  OCTAVE_LOCAL_BUFFER (float, w, n);
674 
675  F77_XFCN (sch1up, SCH1UP, (n, chol_mat.fortran_vec (), chol_mat.rows (),
676  utmp.fortran_vec (), w));
677  }
678 
679  template <>
680  octave_idx_type
681  chol<FloatMatrix>::downdate (const FloatColumnVector& u)
682  {
683  octave_idx_type info = -1;
684 
685  octave_idx_type n = chol_mat.rows ();
686 
687  if (u.numel () != n)
688  (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
689 
690  FloatColumnVector utmp = u;
691 
692  OCTAVE_LOCAL_BUFFER (float, w, n);
693 
694  F77_XFCN (sch1dn, SCH1DN, (n, chol_mat.fortran_vec (), chol_mat.rows (),
695  utmp.fortran_vec (), w, info));
696 
697  return info;
698  }
699 
700  template <>
701  octave_idx_type
702  chol<FloatMatrix>::insert_sym (const FloatColumnVector& u, octave_idx_type j)
703  {
704  octave_idx_type info = -1;
705 
706  octave_idx_type n = chol_mat.rows ();
707 
708  if (u.numel () != n + 1)
709  (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
710  if (j < 0 || j > n)
711  (*current_liboctave_error_handler) ("cholinsert: index out of range");
712 
713  FloatColumnVector utmp = u;
714 
715  OCTAVE_LOCAL_BUFFER (float, w, n);
716 
717  chol_mat.resize (n+1, n+1);
718 
719  F77_XFCN (schinx, SCHINX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
720  j + 1, utmp.fortran_vec (), w, info));
721 
722  return info;
723  }
724 
725  template <>
726  void
727  chol<FloatMatrix>::delete_sym (octave_idx_type j)
728  {
729  octave_idx_type n = chol_mat.rows ();
730 
731  if (j < 0 || j > n-1)
732  (*current_liboctave_error_handler) ("choldelete: index out of range");
733 
734  OCTAVE_LOCAL_BUFFER (float, w, n);
735 
736  F77_XFCN (schdex, SCHDEX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
737  j + 1, w));
738 
739  chol_mat.resize (n-1, n-1);
740  }
741 
742  template <>
743  void
744  chol<FloatMatrix>::shift_sym (octave_idx_type i, octave_idx_type j)
745  {
746  octave_idx_type n = chol_mat.rows ();
747 
748  if (i < 0 || i > n-1 || j < 0 || j > n-1)
749  (*current_liboctave_error_handler) ("cholshift: index out of range");
750 
751  OCTAVE_LOCAL_BUFFER (float, w, 2*n);
752 
753  F77_XFCN (schshx, SCHSHX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
754  i + 1, j + 1, w));
755  }
756 
757 #endif
758 
759  template <>
760  octave_idx_type
761  chol<ComplexMatrix>::init (const ComplexMatrix& a, bool upper, bool calc_cond)
762  {
763  octave_idx_type a_nr = a.rows ();
764  octave_idx_type a_nc = a.cols ();
765 
766  if (a_nr != a_nc)
767  (*current_liboctave_error_handler) ("chol: requires square matrix");
768 
769  octave_idx_type n = a_nc;
770  octave_idx_type info;
771 
772  is_upper = upper;
773 
774  chol_mat.clear (n, n);
775  if (is_upper)
776  for (octave_idx_type j = 0; j < n; j++)
777  {
778  for (octave_idx_type i = 0; i <= j; i++)
779  chol_mat.xelem (i, j) = a(i, j);
780  for (octave_idx_type i = j+1; i < n; i++)
781  chol_mat.xelem (i, j) = 0.0;
782  }
783  else
784  for (octave_idx_type j = 0; j < n; j++)
785  {
786  for (octave_idx_type i = 0; i < j; i++)
787  chol_mat.xelem (i, j) = 0.0;
788  for (octave_idx_type i = j; i < n; i++)
789  chol_mat.xelem (i, j) = a(i, j);
790  }
791  Complex *h = chol_mat.fortran_vec ();
792 
793  // Calculate the norm of the matrix, for later use.
794  double anorm = 0;
795  if (calc_cond)
796  anorm = xnorm (a, 1);
797 
798  if (is_upper)
799  F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1), n,
800  F77_DBLE_CMPLX_ARG (h), n, info
801  F77_CHAR_ARG_LEN (1)));
802  else
803  F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), n,
804  F77_DBLE_CMPLX_ARG (h), n, info
805  F77_CHAR_ARG_LEN (1)));
806 
807  xrcond = 0.0;
808  if (info > 0)
809  chol_mat.resize (info - 1, info - 1);
810  else if (calc_cond)
811  {
812  octave_idx_type zpocon_info = 0;
813 
814  // Now calculate the condition number for non-singular matrix.
815  Array<Complex> z (dim_vector (2*n, 1));
816  Complex *pz = z.fortran_vec ();
817  Array<double> rz (dim_vector (n, 1));
818  double *prz = rz.fortran_vec ();
819  F77_XFCN (zpocon, ZPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n,
820  F77_DBLE_CMPLX_ARG (h),
821  n, anorm, xrcond, F77_DBLE_CMPLX_ARG (pz), prz, zpocon_info
822  F77_CHAR_ARG_LEN (1)));
823 
824  if (zpocon_info != 0)
825  info = -1;
826  }
827 
828  return info;
829  }
830 
831 #if defined (HAVE_QRUPDATE)
832 
833  template <>
834  void
835  chol<ComplexMatrix>::update (const ComplexColumnVector& u)
836  {
837  octave_idx_type n = chol_mat.rows ();
838 
839  if (u.numel () != n)
840  (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
841 
842  ComplexColumnVector utmp = u;
843 
844  OCTAVE_LOCAL_BUFFER (double, rw, n);
845 
846  F77_XFCN (zch1up, ZCH1UP, (n, F77_DBLE_CMPLX_ARG (chol_mat.fortran_vec ()),
847  chol_mat.rows (),
848  F77_DBLE_CMPLX_ARG (utmp.fortran_vec ()), rw));
849  }
850 
851  template <>
852  octave_idx_type
853  chol<ComplexMatrix>::downdate (const ComplexColumnVector& u)
854  {
855  octave_idx_type info = -1;
856 
857  octave_idx_type n = chol_mat.rows ();
858 
859  if (u.numel () != n)
860  (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
861 
862  ComplexColumnVector utmp = u;
863 
864  OCTAVE_LOCAL_BUFFER (double, rw, n);
865 
866  F77_XFCN (zch1dn, ZCH1DN, (n, F77_DBLE_CMPLX_ARG (chol_mat.fortran_vec ()),
867  chol_mat.rows (),
868  F77_DBLE_CMPLX_ARG (utmp.fortran_vec ()), rw, info));
869 
870  return info;
871  }
872 
873  template <>
874  octave_idx_type
875  chol<ComplexMatrix>::insert_sym (const ComplexColumnVector& u,
876  octave_idx_type j)
877  {
878  octave_idx_type info = -1;
879 
880  octave_idx_type n = chol_mat.rows ();
881 
882  if (u.numel () != n + 1)
883  (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
884  if (j < 0 || j > n)
885  (*current_liboctave_error_handler) ("cholinsert: index out of range");
886 
887  ComplexColumnVector utmp = u;
888 
889  OCTAVE_LOCAL_BUFFER (double, rw, n);
890 
891  chol_mat.resize (n+1, n+1);
892 
893  F77_XFCN (zchinx, ZCHINX, (n, F77_DBLE_CMPLX_ARG (chol_mat.fortran_vec ()),
894  chol_mat.rows (),
895  j + 1, F77_DBLE_CMPLX_ARG (utmp.fortran_vec ()), rw, info));
896 
897  return info;
898  }
899 
900  template <>
901  void
902  chol<ComplexMatrix>::delete_sym (octave_idx_type j)
903  {
904  octave_idx_type n = chol_mat.rows ();
905 
906  if (j < 0 || j > n-1)
907  (*current_liboctave_error_handler) ("choldelete: index out of range");
908 
909  OCTAVE_LOCAL_BUFFER (double, rw, n);
910 
911  F77_XFCN (zchdex, ZCHDEX, (n, F77_DBLE_CMPLX_ARG (chol_mat.fortran_vec ()),
912  chol_mat.rows (),
913  j + 1, rw));
914 
915  chol_mat.resize (n-1, n-1);
916  }
917 
918  template <>
919  void
920  chol<ComplexMatrix>::shift_sym (octave_idx_type i, octave_idx_type j)
921  {
922  octave_idx_type n = chol_mat.rows ();
923 
924  if (i < 0 || i > n-1 || j < 0 || j > n-1)
925  (*current_liboctave_error_handler) ("cholshift: index out of range");
926 
927  OCTAVE_LOCAL_BUFFER (Complex, w, n);
928  OCTAVE_LOCAL_BUFFER (double, rw, n);
929 
930  F77_XFCN (zchshx, ZCHSHX, (n, F77_DBLE_CMPLX_ARG (chol_mat.fortran_vec ()),
931  chol_mat.rows (),
932  i + 1, j + 1, F77_DBLE_CMPLX_ARG (w), rw));
933  }
934 
935 #endif
936 
937  template <>
938  octave_idx_type
939  chol<FloatComplexMatrix>::init (const FloatComplexMatrix& a, bool upper,
940  bool calc_cond)
941  {
942  octave_idx_type a_nr = a.rows ();
943  octave_idx_type a_nc = a.cols ();
944 
945  if (a_nr != a_nc)
946  (*current_liboctave_error_handler) ("chol: requires square matrix");
947 
948  octave_idx_type n = a_nc;
949  octave_idx_type info;
950 
951  is_upper = upper;
952 
953  chol_mat.clear (n, n);
954  if (is_upper)
955  for (octave_idx_type j = 0; j < n; j++)
956  {
957  for (octave_idx_type i = 0; i <= j; i++)
958  chol_mat.xelem (i, j) = a(i, j);
959  for (octave_idx_type i = j+1; i < n; i++)
960  chol_mat.xelem (i, j) = 0.0f;
961  }
962  else
963  for (octave_idx_type j = 0; j < n; j++)
964  {
965  for (octave_idx_type i = 0; i < j; i++)
966  chol_mat.xelem (i, j) = 0.0f;
967  for (octave_idx_type i = j; i < n; i++)
968  chol_mat.xelem (i, j) = a(i, j);
969  }
970  FloatComplex *h = chol_mat.fortran_vec ();
971 
972  // Calculate the norm of the matrix, for later use.
973  float anorm = 0;
974  if (calc_cond)
975  anorm = xnorm (a, 1);
976 
977  if (is_upper)
978  F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1), n, F77_CMPLX_ARG (h),
979  n, info
980  F77_CHAR_ARG_LEN (1)));
981  else
982  F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), n, F77_CMPLX_ARG (h),
983  n, info
984  F77_CHAR_ARG_LEN (1)));
985 
986  xrcond = 0.0;
987  if (info > 0)
988  chol_mat.resize (info - 1, info - 1);
989  else if (calc_cond)
990  {
991  octave_idx_type cpocon_info = 0;
992 
993  // Now calculate the condition number for non-singular matrix.
994  Array<FloatComplex> z (dim_vector (2*n, 1));
995  FloatComplex *pz = z.fortran_vec ();
996  Array<float> rz (dim_vector (n, 1));
997  float *prz = rz.fortran_vec ();
998  F77_XFCN (cpocon, CPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, F77_CMPLX_ARG (h),
999  n, anorm, xrcond, F77_CMPLX_ARG (pz), prz, cpocon_info
1000  F77_CHAR_ARG_LEN (1)));
1001 
1002  if (cpocon_info != 0)
1003  info = -1;
1004  }
1005 
1006  return info;
1007  }
1008 
1009 #if defined (HAVE_QRUPDATE)
1010 
1011  template <>
1012  void
1013  chol<FloatComplexMatrix>::update (const FloatComplexColumnVector& u)
1014  {
1015  octave_idx_type n = chol_mat.rows ();
1016 
1017  if (u.numel () != n)
1018  (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
1019 
1020  FloatComplexColumnVector utmp = u;
1021 
1022  OCTAVE_LOCAL_BUFFER (float, rw, n);
1023 
1024  F77_XFCN (cch1up, CCH1UP, (n, F77_CMPLX_ARG (chol_mat.fortran_vec ()),
1025  chol_mat.rows (),
1026  F77_CMPLX_ARG (utmp.fortran_vec ()), rw));
1027  }
1028 
1029  template <>
1030  octave_idx_type
1031  chol<FloatComplexMatrix>::downdate (const FloatComplexColumnVector& u)
1032  {
1033  octave_idx_type info = -1;
1034 
1035  octave_idx_type n = chol_mat.rows ();
1036 
1037  if (u.numel () != n)
1038  (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
1039 
1040  FloatComplexColumnVector utmp = u;
1041 
1042  OCTAVE_LOCAL_BUFFER (float, rw, n);
1043 
1044  F77_XFCN (cch1dn, CCH1DN, (n, F77_CMPLX_ARG (chol_mat.fortran_vec ()),
1045  chol_mat.rows (),
1046  F77_CMPLX_ARG (utmp.fortran_vec ()), rw, info));
1047 
1048  return info;
1049  }
1050 
1051  template <>
1052  octave_idx_type
1053  chol<FloatComplexMatrix>::insert_sym (const FloatComplexColumnVector& u,
1054  octave_idx_type j)
1055  {
1056  octave_idx_type info = -1;
1057 
1058  octave_idx_type n = chol_mat.rows ();
1059 
1060  if (u.numel () != n + 1)
1061  (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
1062  if (j < 0 || j > n)
1063  (*current_liboctave_error_handler) ("cholinsert: index out of range");
1064 
1065  FloatComplexColumnVector utmp = u;
1066 
1067  OCTAVE_LOCAL_BUFFER (float, rw, n);
1068 
1069  chol_mat.resize (n+1, n+1);
1070 
1071  F77_XFCN (cchinx, CCHINX, (n, F77_CMPLX_ARG (chol_mat.fortran_vec ()),
1072  chol_mat.rows (),
1073  j + 1, F77_CMPLX_ARG (utmp.fortran_vec ()), rw, info));
1074 
1075  return info;
1076  }
1077 
1078  template <>
1079  void
1080  chol<FloatComplexMatrix>::delete_sym (octave_idx_type j)
1081  {
1082  octave_idx_type n = chol_mat.rows ();
1083 
1084  if (j < 0 || j > n-1)
1085  (*current_liboctave_error_handler) ("choldelete: index out of range");
1086 
1087  OCTAVE_LOCAL_BUFFER (float, rw, n);
1088 
1089  F77_XFCN (cchdex, CCHDEX, (n, F77_CMPLX_ARG (chol_mat.fortran_vec ()),
1090  chol_mat.rows (),
1091  j + 1, rw));
1092 
1093  chol_mat.resize (n-1, n-1);
1094  }
1095 
1096  template <>
1097  void
1098  chol<FloatComplexMatrix>::shift_sym (octave_idx_type i, octave_idx_type j)
1099  {
1100  octave_idx_type n = chol_mat.rows ();
1101 
1102  if (i < 0 || i > n-1 || j < 0 || j > n-1)
1103  (*current_liboctave_error_handler) ("cholshift: index out of range");
1104 
1105  OCTAVE_LOCAL_BUFFER (FloatComplex, w, n);
1106  OCTAVE_LOCAL_BUFFER (float, rw, n);
1107 
1108  F77_XFCN (cchshx, CCHSHX, (n, F77_CMPLX_ARG (chol_mat.fortran_vec ()),
1109  chol_mat.rows (),
1110  i + 1, j + 1, F77_CMPLX_ARG (w), rw));
1111  }
1112 
1113 #endif
1114 
1115  // Instantiations we need.
1116 
1117  template class chol<Matrix>;
1118 
1119  template class chol<FloatMatrix>;
1120 
1121  template class chol<ComplexMatrix>;
1122 
1123  template class chol<FloatComplexMatrix>;
1124 
1125  template Matrix
1126  chol2inv<Matrix> (const Matrix& r);
1127 
1128  template ComplexMatrix
1129  chol2inv<ComplexMatrix> (const ComplexMatrix& r);
1130 
1131  template FloatMatrix
1132  chol2inv<FloatMatrix> (const FloatMatrix& r);
1133 
1134  template FloatComplexMatrix
1135  chol2inv<FloatComplexMatrix> (const FloatComplexMatrix& r);
1136  }
1137 }
octave::math::chol2inv< FloatComplexMatrix >
template FloatComplexMatrix chol2inv< FloatComplexMatrix >(const FloatComplexMatrix &r)
octave
Octave interface to the compression and uncompression libraries.
Definition: aepbalance.cc:47
octave::math::chol::init
octave_idx_type init(const T &a, bool upper, bool calc_cond)
Array::numel
octave_idx_type numel(void) const
Number of elements in the array.
Definition: Array.h:363
chol
Definition: mx-defs.h:65
current_liboctave_error_handler
OCTAVE_NORETURN liboctave_error_handler current_liboctave_error_handler
Definition: lo-error.c:38
F77_DBLE_CMPLX_ARG
#define F77_DBLE_CMPLX_ARG(x)
Definition: f77-fcn.h:345
octave::math::conj
double conj(double x)
Definition: lo-mappers.h:93
k
for large enough k
Definition: lu.cc:606
octave::math::chol2inv< Matrix >
template Matrix chol2inv< Matrix >(const Matrix &r)
zero
is greater than zero
Definition: load-path.cc:2339
octave::math::chol::insert_sym
octave_idx_type insert_sym(const VT &u, octave_idx_type j)
u
u
Definition: lu.cc:138
F77_XFCN
#define F77_XFCN(f, F, args)
Definition: f77-fcn.h:52
a_nc
octave_idx_type a_nc
Definition: sylvester.cc:72
octave::math::warn_qrupdate_once
void warn_qrupdate_once(void)
Array::rows
octave_idx_type rows(void) const
Definition: Array.h:401
a
calling an anonymous function involves an overhead quite comparable to the overhead of an m file function Passing a handle to a built in function is because the interpreter is not involved in the internal loop For a
Definition: cellfun.cc:398
octave::math::chol::inverse
T inverse(void) const
Definition: chol.cc:255
h
double h
Definition: graphics.cc:11205
conj
ComplexColumnVector conj(const ComplexColumnVector &a)
Definition: CColVector.cc:216
a_nr
octave_idx_type a_nr
Definition: sylvester.cc:71
octave::math::chol2inv
T chol2inv(const T &r)
Definition: chol.cc:247
Faddeeva::w
std::complex< double > w(std::complex< double > z, double relerr=0)
octave::math::chol::downdate
octave_idx_type downdate(const VT &u)
chol2inv_internal
static Matrix chol2inv_internal(const Matrix &r, bool is_upper=true)
Definition: chol.cc:54
tmp
double tmp
Definition: data.cc:6300
octave::math::chol::delete_sym
void delete_sym(octave_idx_type j)
retval
octave_value retval
Definition: data.cc:6294
xnorm
OCTAVE_API double xnorm(const ColumnVector &x, double p)
Definition: oct-norm.cc:549
octave::math::chol2inv< ComplexMatrix >
template ComplexMatrix chol2inv< ComplexMatrix >(const ComplexMatrix &r)
octave::math::imag
double imag(double)
Definition: lo-mappers.h:103
Matrix
Definition: dMatrix.h:37
octave::math::chol::set
void set(const T &R)
Definition: chol.cc:262
octave::math::chol::update
void update(const VT &u)
Array::xelem
T & xelem(octave_idx_type n)
Definition: Array.h:455
F77_CMPLX_ARG
#define F77_CMPLX_ARG(x)
Definition: f77-fcn.h:339
i
=val(i)}if ode{val(i)}occurs in table i
Definition: lookup.cc:239
p
p
Definition: lu.cc:138
octave::math::chol::shift_sym
void shift_sym(octave_idx_type i, octave_idx_type j)
octave::math::chol2inv< FloatMatrix >
template FloatMatrix chol2inv< FloatMatrix >(const FloatMatrix &r)
OCTAVE_LOCAL_BUFFER
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition: oct-locbuf.h:200
octave_idx_type
idx_vector
Definition: idx-vector.h:52
FloatComplex
std::complex< float > FloatComplex
Definition: oct-cmplx.h:32
Complex
std::complex< double > Complex
Definition: oct-cmplx.h:31
Array::fortran_vec
const T * fortran_vec(void) const
Definition: Array.h:584
Array::cols
octave_idx_type cols(void) const
Definition: Array.h:409
dim_vector
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:87
Generated on Wed May 10 2017 15:42:57 for GNU Octave by   doxygen 1.8.8