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