GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
chol.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (C) 1994-2024 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 OCTAVE_BEGIN_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 OCTAVE_BEGIN_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 () const
251 {
252  return chol2inv_internal (m_chol_mat, m_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  m_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 = m_chol_mat.rows ();
274 
275  if (u.numel () != n)
276  (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
277 
278  init (m_chol_mat.hermitian () * m_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 = m_chol_mat.rows ();
301 
302  if (u.numel () != n)
303  (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
304 
305  if (singular (m_chol_mat))
306  info = 2;
307  else
308  {
309  info = init (m_chol_mat.hermitian () * m_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 = m_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 (m_chol_mat))
335  info = 2;
336  else if (std::imag (u(j)) != zero)
337  info = 3;
338  else
339  {
340  T a = m_chol_mat.hermitian () * m_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 = m_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 = m_chol_mat.hermitian () * m_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 = m_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 = m_chol_mat.hermitian () * m_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_API 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  m_is_upper = upper;
422 
423  m_chol_mat.clear (n, n);
424  if (m_is_upper)
425  for (octave_idx_type j = 0; j < n; j++)
426  {
427  for (octave_idx_type i = 0; i <= j; i++)
428  m_chol_mat.xelem (i, j) = a(i, j);
429  for (octave_idx_type i = j+1; i < n; i++)
430  m_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  m_chol_mat.xelem (i, j) = 0.0;
437  for (octave_idx_type i = j; i < n; i++)
438  m_chol_mat.xelem (i, j) = a(i, j);
439  }
440  double *h = m_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 = octave::xnorm (a, 1);
446 
447  if (m_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  m_rcond = 0.0;
455  if (info > 0)
456  m_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 (m_is_upper)
466  F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h,
467  n, anorm, m_rcond, 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, m_rcond, 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 OCTAVE_API void
485 chol<Matrix>::update (const ColumnVector& u)
486 {
487  F77_INT n = to_f77_int (m_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, m_chol_mat.fortran_vec (), n,
497  utmp.fortran_vec (), w));
498 }
499 
500 template <>
501 OCTAVE_API octave_idx_type
502 chol<Matrix>::downdate (const ColumnVector& u)
503 {
504  F77_INT info = -1;
505 
506  F77_INT n = to_f77_int (m_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, m_chol_mat.fortran_vec (), n,
516  utmp.fortran_vec (), w, info));
517 
518  return info;
519 }
520 
521 template <>
522 OCTAVE_API 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 (m_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  m_chol_mat.resize (n+1, n+1);
540  F77_INT ldcm = to_f77_int (m_chol_mat.rows ());
541 
542  F77_XFCN (dchinx, DCHINX, (n, m_chol_mat.fortran_vec (), ldcm,
543  j + 1, utmp.fortran_vec (), w, info));
544 
545  return info;
546 }
547 
548 template <>
549 OCTAVE_API void
550 chol<Matrix>::delete_sym (octave_idx_type j_arg)
551 {
552  F77_INT n = to_f77_int (m_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, m_chol_mat.fortran_vec (), n, j + 1, w));
561 
562  m_chol_mat.resize (n-1, n-1);
563 }
564 
565 template <>
566 OCTAVE_API void
567 chol<Matrix>::shift_sym (octave_idx_type i_arg, octave_idx_type j_arg)
568 {
569  F77_INT n = to_f77_int (m_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, m_chol_mat.fortran_vec (), n,
579  i + 1, j + 1, w));
580 }
581 
582 #endif
583 
584 template <>
585 OCTAVE_API 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  m_is_upper = upper;
598 
599  m_chol_mat.clear (n, n);
600  if (m_is_upper)
601  for (octave_idx_type j = 0; j < n; j++)
602  {
603  for (octave_idx_type i = 0; i <= j; i++)
604  m_chol_mat.xelem (i, j) = a(i, j);
605  for (octave_idx_type i = j+1; i < n; i++)
606  m_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  m_chol_mat.xelem (i, j) = 0.0f;
613  for (octave_idx_type i = j; i < n; i++)
614  m_chol_mat.xelem (i, j) = a(i, j);
615  }
616  float *h = m_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 = octave::xnorm (a, 1);
622 
623  if (m_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  m_rcond = 0.0;
631  if (info > 0)
632  m_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 (m_is_upper)
642  F77_XFCN (spocon, SPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h,
643  n, anorm, m_rcond, 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, m_rcond, 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 OCTAVE_API void
661 chol<FloatMatrix>::update (const FloatColumnVector& u)
662 {
663  F77_INT n = to_f77_int (m_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, m_chol_mat.fortran_vec (), n,
673  utmp.fortran_vec (), w));
674 }
675 
676 template <>
677 OCTAVE_API octave_idx_type
678 chol<FloatMatrix>::downdate (const FloatColumnVector& u)
679 {
680  F77_INT info = -1;
681 
682  F77_INT n = to_f77_int (m_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, m_chol_mat.fortran_vec (), n,
692  utmp.fortran_vec (), w, info));
693 
694  return info;
695 }
696 
697 template <>
698 OCTAVE_API 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 (m_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  m_chol_mat.resize (n+1, n+1);
717  F77_INT ldcm = to_f77_int (m_chol_mat.rows ());
718 
719  F77_XFCN (schinx, SCHINX, (n, m_chol_mat.fortran_vec (), ldcm,
720  j + 1, utmp.fortran_vec (), w, info));
721 
722  return info;
723 }
724 
725 template <>
726 OCTAVE_API void
727 chol<FloatMatrix>::delete_sym (octave_idx_type j_arg)
728 {
729  F77_INT n = to_f77_int (m_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, m_chol_mat.fortran_vec (), n,
738  j + 1, w));
739 
740  m_chol_mat.resize (n-1, n-1);
741 }
742 
743 template <>
744 OCTAVE_API void
745 chol<FloatMatrix>::shift_sym (octave_idx_type i_arg, octave_idx_type j_arg)
746 {
747  F77_INT n = to_f77_int (m_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, m_chol_mat.fortran_vec (), n,
757  i + 1, j + 1, w));
758 }
759 
760 #endif
761 
762 template <>
763 OCTAVE_API 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  m_is_upper = upper;
776 
777  m_chol_mat.clear (n, n);
778  if (m_is_upper)
779  for (octave_idx_type j = 0; j < n; j++)
780  {
781  for (octave_idx_type i = 0; i <= j; i++)
782  m_chol_mat.xelem (i, j) = a(i, j);
783  for (octave_idx_type i = j+1; i < n; i++)
784  m_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  m_chol_mat.xelem (i, j) = 0.0;
791  for (octave_idx_type i = j; i < n; i++)
792  m_chol_mat.xelem (i, j) = a(i, j);
793  }
794  Complex *h = m_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 = octave::xnorm (a, 1);
800 
801  if (m_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  m_rcond = 0.0;
811  if (info > 0)
812  m_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, m_rcond,
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 OCTAVE_API void
838 chol<ComplexMatrix>::update (const ComplexColumnVector& u)
839 {
840  F77_INT n = to_f77_int (m_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 (m_chol_mat.fortran_vec ()),
851  n,
852  F77_DBLE_CMPLX_ARG (utmp.fortran_vec ()),
853  rw));
854 }
855 
856 template <>
857 OCTAVE_API octave_idx_type
858 chol<ComplexMatrix>::downdate (const ComplexColumnVector& u)
859 {
860  F77_INT info = -1;
861 
862  F77_INT n = to_f77_int (m_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 (m_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_API 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 (m_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  m_chol_mat.resize (n+1, n+1);
900  F77_INT ldcm = to_f77_int (m_chol_mat.rows ());
901 
902  F77_XFCN (zchinx, ZCHINX, (n,
903  F77_DBLE_CMPLX_ARG (m_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 OCTAVE_API void
913 chol<ComplexMatrix>::delete_sym (octave_idx_type j_arg)
914 {
915  F77_INT n = to_f77_int (m_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 (m_chol_mat.fortran_vec ()),
925  n, j + 1, rw));
926 
927  m_chol_mat.resize (n-1, n-1);
928 }
929 
930 template <>
931 OCTAVE_API 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 (m_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 (m_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_API 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  m_is_upper = upper;
968 
969  m_chol_mat.clear (n, n);
970  if (m_is_upper)
971  for (octave_idx_type j = 0; j < n; j++)
972  {
973  for (octave_idx_type i = 0; i <= j; i++)
974  m_chol_mat.xelem (i, j) = a(i, j);
975  for (octave_idx_type i = j+1; i < n; i++)
976  m_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  m_chol_mat.xelem (i, j) = 0.0f;
983  for (octave_idx_type i = j; i < n; i++)
984  m_chol_mat.xelem (i, j) = a(i, j);
985  }
986  FloatComplex *h = m_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 = octave::xnorm (a, 1);
992 
993  if (m_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  m_rcond = 0.0;
1003  if (info > 0)
1004  m_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, m_rcond,
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 OCTAVE_API void
1030 chol<FloatComplexMatrix>::update (const FloatComplexColumnVector& u)
1031 {
1032  F77_INT n = to_f77_int (m_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 (m_chol_mat.fortran_vec ()),
1042  n, F77_CMPLX_ARG (utmp.fortran_vec ()), rw));
1043 }
1044 
1045 template <>
1046 OCTAVE_API octave_idx_type
1047 chol<FloatComplexMatrix>::downdate (const FloatComplexColumnVector& u)
1048 {
1049  F77_INT info = -1;
1050 
1051  F77_INT n = to_f77_int (m_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 (m_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_API 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 (m_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  m_chol_mat.resize (n+1, n+1);
1087  F77_INT ldcm = to_f77_int (m_chol_mat.rows ());
1088 
1089  F77_XFCN (cchinx, CCHINX, (n, F77_CMPLX_ARG (m_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 OCTAVE_API void
1099 chol<FloatComplexMatrix>::delete_sym (octave_idx_type j_arg)
1100 {
1101  F77_INT n = to_f77_int (m_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 (m_chol_mat.fortran_vec ()),
1110  n, j + 1, rw));
1111 
1112  m_chol_mat.resize (n-1, n-1);
1113 }
1114 
1115 template <>
1116 OCTAVE_API 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 (m_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 (m_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 OCTAVE_API Matrix
1147 chol2inv<Matrix> (const Matrix& r);
1148 
1149 template OCTAVE_API ComplexMatrix
1150 chol2inv<ComplexMatrix> (const ComplexMatrix& r);
1151 
1152 template OCTAVE_API FloatMatrix
1153 chol2inv<FloatMatrix> (const FloatMatrix& r);
1154 
1155 template OCTAVE_API FloatComplexMatrix
1156 chol2inv<FloatComplexMatrix> (const FloatComplexMatrix& r);
1157 
1158 OCTAVE_END_NAMESPACE(math)
1159 OCTAVE_END_NAMESPACE(octave)
conj
ComplexColumnVector conj(const ComplexColumnVector &a)
Definition: CColVector.cc:217
Array::fortran_vec
T * fortran_vec()
Size of the specified dimension.
Definition: Array-base.cc:1764
Array::rows
octave_idx_type rows() const
Definition: Array.h:459
Array::cols
octave_idx_type cols() const
Definition: Array.h:469
Array::xelem
T & xelem(octave_idx_type n)
Size of the specified dimension.
Definition: Array.h:524
Array::numel
octave_idx_type numel() const
Number of elements in the array.
Definition: Array.h:414
Matrix
Definition: dMatrix.h:42
chol
Definition: chol.h:38
chol::set
void set(const T &R)
Definition: chol.cc:257
chol::shift_sym
void shift_sym(octave_idx_type i, octave_idx_type j)
chol::delete_sym
void delete_sym(octave_idx_type j)
chol::insert_sym
octave_idx_type insert_sym(const VT &u, octave_idx_type j)
chol::downdate
octave_idx_type downdate(const VT &u)
chol::update
void update(const VT &u)
chol::inverse
T inverse() const
Definition: chol.cc:250
dim_vector
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:94
idx_vector
Definition: idx-vector.h:60
octave_idx_type
imag
ColumnVector imag(const ComplexColumnVector &a)
Definition: dColVector.cc:143
OCTAVE_BEGIN_NAMESPACE
OCTAVE_BEGIN_NAMESPACE(octave) static octave_value daspk_fcn
F77_DBLE_CMPLX_ARG
#define F77_DBLE_CMPLX_ARG(x)
Definition: f77-fcn.h:316
F77_CMPLX_ARG
#define F77_CMPLX_ARG(x)
Definition: f77-fcn.h:310
F77_XFCN
#define F77_XFCN(f, F, args)
Definition: f77-fcn.h:45
F77_INT
octave_f77_int_type F77_INT
Definition: f77-fcn.h:306
chol2inv< Matrix >
template Matrix chol2inv< Matrix >(const Matrix &r)
chol2inv< FloatComplexMatrix >
template FloatComplexMatrix chol2inv< FloatComplexMatrix >(const FloatComplexMatrix &r)
chol2inv< ComplexMatrix >
template ComplexMatrix chol2inv< ComplexMatrix >(const ComplexMatrix &r)
chol2inv< FloatMatrix >
template FloatMatrix chol2inv< FloatMatrix >(const FloatMatrix &r)
chol2inv
T chol2inv(const T &r)
Definition: chol.cc:242
OCTAVE_API
#define OCTAVE_API
Definition: main.cc:55
n
octave_idx_type n
Definition: mx-inlines.cc:761
r
T * r
Definition: mx-inlines.cc:781
Faddeeva::w
std::complex< double > w(std::complex< double > z, double relerr=0)
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
double xnorm(const ColumnVector &x, double p)
Definition: oct-norm.cc:610
warn_qrupdate_once
void warn_qrupdate_once()