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