GNU Octave 10.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 
Loading...
Searching...
No Matches
chol.cc
Go to the documentation of this file.
1////////////////////////////////////////////////////////////////////////
2//
3// Copyright (C) 1994-2025 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
50OCTAVE_BEGIN_NAMESPACE(octave)
51
52static Matrix
53chol2inv_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.rwdata ();
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
100static FloatMatrix
101chol2inv_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.rwdata ();
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
148static ComplexMatrix
149chol2inv_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.rwdata ()), 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.rwdata ()), 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
193static FloatComplexMatrix
194chol2inv_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.rwdata ()), 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.rwdata ()), 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
238OCTAVE_BEGIN_NAMESPACE(math)
239
240template <typename T>
241T
242chol2inv (const T& r)
243{
244 return chol2inv_internal (r);
245}
246
247// Compute the inverse of a matrix using the Cholesky factorization.
248template <typename T>
249T
250chol<T>::inverse () const
251{
252 return chol2inv_internal (m_chol_mat, m_is_upper);
253}
254
255template <typename T>
256void
257chol<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
267template <typename T>
268void
269chol<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
282template <typename T>
283bool
284singular (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
292template <typename T>
293octave_idx_type
294chol<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
317template <typename T>
318octave_idx_type
319chol<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
359template <typename T>
360void
361chol<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
376template <typename T>
377void
378chol<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
408template <>
409OCTAVE_API octave_idx_type
410chol<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.rwdata ();
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.rwdata ();
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
483template <>
484OCTAVE_API void
485chol<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.rwdata (), n,
497 utmp.rwdata (), w));
498}
499
500template <>
501OCTAVE_API octave_idx_type
502chol<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.rwdata (), n,
516 utmp.rwdata (), w, info));
517
518 return info;
519}
520
521template <>
522OCTAVE_API octave_idx_type
523chol<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.rwdata (), ldcm,
543 j + 1, utmp.rwdata (), w, info));
544
545 return info;
546}
547
548template <>
549OCTAVE_API void
550chol<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.rwdata (), n, j + 1, w));
561
562 m_chol_mat.resize (n-1, n-1);
563}
564
565template <>
566OCTAVE_API void
567chol<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.rwdata (), n,
579 i + 1, j + 1, w));
580}
581
582#endif
583
584template <>
585OCTAVE_API octave_idx_type
586chol<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.rwdata ();
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.rwdata ();
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
659template <>
660OCTAVE_API void
661chol<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.rwdata (), n,
673 utmp.rwdata (), w));
674}
675
676template <>
677OCTAVE_API octave_idx_type
678chol<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.rwdata (), n,
692 utmp.rwdata (), w, info));
693
694 return info;
695}
696
697template <>
698OCTAVE_API octave_idx_type
699chol<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.rwdata (), ldcm,
720 j + 1, utmp.rwdata (), w, info));
721
722 return info;
723}
724
725template <>
726OCTAVE_API void
727chol<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.rwdata (), n,
738 j + 1, w));
739
740 m_chol_mat.resize (n-1, n-1);
741}
742
743template <>
744OCTAVE_API void
745chol<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.rwdata (), n,
757 i + 1, j + 1, w));
758}
759
760#endif
761
762template <>
763OCTAVE_API octave_idx_type
764chol<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.rwdata ();
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.rwdata ();
820 Array<double> rz (dim_vector (n, 1));
821 double *prz = rz.rwdata ();
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
836template <>
837OCTAVE_API void
838chol<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.rwdata ()),
851 n,
852 F77_DBLE_CMPLX_ARG (utmp.rwdata ()),
853 rw));
854}
855
856template <>
857OCTAVE_API octave_idx_type
858chol<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.rwdata ()),
873 n,
874 F77_DBLE_CMPLX_ARG (utmp.rwdata ()),
875 rw, info));
876
877 return info;
878}
879
880template <>
881OCTAVE_API octave_idx_type
882chol<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.rwdata ()),
904 ldcm, j + 1,
905 F77_DBLE_CMPLX_ARG (utmp.rwdata ()),
906 rw, info));
907
908 return info;
909}
910
911template <>
912OCTAVE_API void
913chol<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.rwdata ()),
925 n, j + 1, rw));
926
927 m_chol_mat.resize (n-1, n-1);
928}
929
930template <>
931OCTAVE_API void
932chol<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.rwdata ()),
947 n, i + 1, j + 1,
948 F77_DBLE_CMPLX_ARG (w), rw));
949}
950
951#endif
952
953template <>
954OCTAVE_API octave_idx_type
955chol<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.rwdata ();
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.rwdata ();
1012 Array<float> rz (dim_vector (n, 1));
1013 float *prz = rz.rwdata ();
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
1028template <>
1029OCTAVE_API void
1030chol<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.rwdata ()),
1042 n, F77_CMPLX_ARG (utmp.rwdata ()), rw));
1043}
1044
1045template <>
1046OCTAVE_API octave_idx_type
1047chol<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.rwdata ()),
1061 n, F77_CMPLX_ARG (utmp.rwdata ()),
1062 rw, info));
1063
1064 return info;
1065}
1066
1067template <>
1068OCTAVE_API octave_idx_type
1069chol<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.rwdata ()),
1090 ldcm, j + 1,
1091 F77_CMPLX_ARG (utmp.rwdata ()),
1092 rw, info));
1093
1094 return info;
1095}
1096
1097template <>
1098OCTAVE_API void
1099chol<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.rwdata ()),
1110 n, j + 1, rw));
1111
1112 m_chol_mat.resize (n-1, n-1);
1113}
1114
1115template <>
1116OCTAVE_API void
1117chol<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.rwdata ()),
1131 n, i + 1, j + 1, F77_CMPLX_ARG (w), rw));
1132}
1133
1134#endif
1135
1136// Instantiations we need.
1137
1138template class chol<Matrix>;
1139
1140template class chol<FloatMatrix>;
1141
1142template class chol<ComplexMatrix>;
1143
1144template class chol<FloatComplexMatrix>;
1145
1146template OCTAVE_API Matrix
1147chol2inv<Matrix> (const Matrix& r);
1148
1149template OCTAVE_API ComplexMatrix
1150chol2inv<ComplexMatrix> (const ComplexMatrix& r);
1151
1152template OCTAVE_API FloatMatrix
1153chol2inv<FloatMatrix> (const FloatMatrix& r);
1154
1155template OCTAVE_API FloatComplexMatrix
1156chol2inv<FloatComplexMatrix> (const FloatComplexMatrix& r);
1157
1158OCTAVE_END_NAMESPACE(math)
1159OCTAVE_END_NAMESPACE(octave)
Array
N Dimensional Array with copy-on-write semantics.
Definition Array.h:130
Array::xelem
T & xelem(octave_idx_type n)
Size of the specified dimension.
Definition Array.h:525
Array::rows
octave_idx_type rows() const
Definition Array.h:463
Array::cols
octave_idx_type cols() const
Definition Array.h:473
Array::rwdata
T * rwdata()
Size of the specified dimension.
Definition Array-base.cc:1787
Array::numel
octave_idx_type numel() const
Number of elements in the array.
Definition Array.h:418
chol
Definition chol.h:37
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:90
idx_vector
Definition idx-vector.h:58
octave_idx_type
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.in.cc:55
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
warn_qrupdate_once
void warn_qrupdate_once()