Mercurial > octave-nkf

comparison libinterp/dldfcn/chol.cc @ 15195:2fc554ffbc28

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
split libinterp from src * libinterp: New directory. Move all files from src directory here except Makefile.am, main.cc, main-cli.cc, mkoctfile.in.cc, mkoctfilr.in.sh, octave-config.in.cc, octave-config.in.sh. * libinterp/Makefile.am: New file, extracted from src/Makefile.am. * src/Makefile.am: Delete everything except targets and definitions needed to build and link main and utility programs. * Makefile.am (SUBDIRS): Include libinterp in the list. * autogen.sh: Run config-module.sh in libinterp/dldfcn directory, not src/dldfcn directory. * configure.ac (AC_CONFIG_SRCDIR): Use libinterp/octave.cc, not src/octave.cc. (DL_LDFLAGS, LIBOCTINTERP): Use libinterp, not src. (AC_CONFIG_FILES): Include libinterp/Makefile in the list. * find-docstring-files.sh: Look in libinterp, not src. * gui/src/Makefile.am (liboctgui_la_CPPFLAGS): Find header files in libinterp, not src.
author John W. Eaton <jwe@octave.org>
date Sat, 18 Aug 2012 16:23:39 -0400
parents src/dldfcn/chol.cc@82d51d6e470c
children 94cdf82d4a0c
comparison
equal deleted inserted replaced
15194:0f0b795044c3 15195:2fc554ffbc28
1 /*
2
3 Copyright (C) 1996-2012 John W. Eaton
4 Copyright (C) 2008-2009 Jaroslav Hajek
5 Copyright (C) 2008-2009 VZLU Prague
6
7 This file is part of Octave.
8
9 Octave is free software; you can redistribute it and/or modify it
10 under the terms of the GNU General Public License as published by the
11 Free Software Foundation; either version 3 of the License, or (at your
12 option) any later version.
13
14 Octave is distributed in the hope that it will be useful, but WITHOUT
15 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
16 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
17 for more details.
18
19 You should have received a copy of the GNU General Public License
20 along with Octave; see the file COPYING. If not, see
21 <http://www.gnu.org/licenses/>.
22
23 */
24
25
26 #ifdef HAVE_CONFIG_H
27 #include <config.h>
28 #endif
29
30 #include "CmplxCHOL.h"
31 #include "dbleCHOL.h"
32 #include "fCmplxCHOL.h"
33 #include "floatCHOL.h"
34 #include "SparseCmplxCHOL.h"
35 #include "SparsedbleCHOL.h"
36 #include "oct-spparms.h"
37 #include "sparse-util.h"
38
39 #include "ov-re-sparse.h"
40 #include "ov-cx-sparse.h"
41 #include "defun-dld.h"
42 #include "error.h"
43 #include "gripes.h"
44 #include "oct-obj.h"
45 #include "utils.h"
46
47 template <class CHOLT>
48 static octave_value
49 get_chol_r (const CHOLT& fact)
50 {
51 return octave_value (fact.chol_matrix (),
52 MatrixType (MatrixType::Upper));
53 }
54
55 template <class CHOLT>
56 static octave_value
57 get_chol_l (const CHOLT& fact)
58 {
59 return octave_value (fact.chol_matrix ().transpose (),
60 MatrixType (MatrixType::Lower));
61 }
62
63 DEFUN_DLD (chol, args, nargout,
64 "-*- texinfo -*-\n\
65 @deftypefn {Loadable Function} {@var{R} =} chol (@var{A})\n\
66 @deftypefnx {Loadable Function} {[@var{R}, @var{p}] =} chol (@var{A})\n\
67 @deftypefnx {Loadable Function} {[@var{R}, @var{p}, @var{Q}] =} chol (@var{S})\n\
68 @deftypefnx {Loadable Function} {[@var{R}, @var{p}, @var{Q}] =} chol (@var{S}, \"vector\")\n\
69 @deftypefnx {Loadable Function} {[@var{L}, @dots{}] =} chol (@dots{}, \"lower\")\n\
70 @deftypefnx {Loadable Function} {[@var{L}, @dots{}] =} chol (@dots{}, \"upper\")\n\
71 @cindex Cholesky factorization\n\
72 Compute the Cholesky@tie{}factor, @var{R}, of the symmetric positive definite\n\
73 matrix @var{A}, where\n\
74 @tex\n\
75 $ R^T R = A $.\n\
76 @end tex\n\
77 @ifnottex\n\
78 \n\
79 @example\n\
80 @var{R}' * @var{R} = @var{A}.\n\
81 @end example\n\
82 \n\
83 @end ifnottex\n\
84 \n\
85 Called with one output argument @code{chol} fails if @var{A} or @var{S} is\n\
86 not positive definite. With two or more output arguments @var{p} flags\n\
87 whether the matrix was positive definite and @code{chol} does not fail. A\n\
88 zero value indicated that the matrix was positive definite and the @var{R}\n\
89 gives the factorization, and @var{p} will have a positive value otherwise.\n\
90 \n\
91 If called with 3 outputs then a sparsity preserving row/column permutation\n\
92 is applied to @var{A} prior to the factorization. That is @var{R}\n\
93 is the factorization of @code{@var{A}(@var{Q},@var{Q})} such that\n\
94 @tex\n\
95 $ R^T R = Q^T A Q$.\n\
96 @end tex\n\
97 @ifnottex\n\
98 \n\
99 @example\n\
100 @var{R}' * @var{R} = @var{Q}' * @var{A} * @var{Q}.\n\
101 @end example\n\
102 \n\
103 @end ifnottex\n\
104 \n\
105 The sparsity preserving permutation is generally returned as a matrix.\n\
106 However, given the flag \"vector\", @var{Q} will be returned as a vector\n\
107 such that\n\
108 @tex\n\
109 $ R^T R = A (Q, Q)$.\n\
110 @end tex\n\
111 @ifnottex\n\
112 \n\
113 @example\n\
114 @var{R}' * @var{R} = @var{A}(@var{Q}, @var{Q}).\n\
115 @end example\n\
116 \n\
117 @end ifnottex\n\
118 \n\
119 Called with either a sparse or full matrix and using the \"lower\" flag,\n\
120 @code{chol} returns the lower triangular factorization such that\n\
121 @tex\n\
122 $ L L^T = A $.\n\
123 @end tex\n\
124 @ifnottex\n\
125 \n\
126 @example\n\
127 @var{L} * @var{L}' = @var{A}.\n\
128 @end example\n\
129 \n\
130 @end ifnottex\n\
131 \n\
132 For full matrices, if the \"lower\" flag is set only the lower triangular\n\
133 part of the matrix is used for the factorization, otherwise the upper\n\
134 triangular part is used.\n\
135 \n\
136 In general the lower triangular factorization is significantly faster for\n\
137 sparse matrices.\n\
138 @seealso{cholinv, chol2inv}\n\
139 @end deftypefn")
140 {
141 octave_value_list retval;
142 int nargin = args.length ();
143 bool LLt = false;
144 bool vecout = false;
145
146 if (nargin < 1 || nargin > 3 || nargout > 3
147 || (! args(0).is_sparse_type () && nargout > 2))
148 {
149 print_usage ();
150 return retval;
151 }
152
153 int n = 1;
154 while (n < nargin && ! error_state)
155 {
156 std::string tmp = args(n++).string_value ();
157
158 if (! error_state )
159 {
160 if (tmp.compare ("vector") == 0)
161 vecout = true;
162 else if (tmp.compare ("lower") == 0)
163 // FIXME currently the option "lower" is handled by transposing the
164 // matrix, factorizing it with the lapack function DPOTRF ('U', ...)
165 // and finally transposing the factor. It would be more efficient to use
166 // DPOTRF ('L', ...) in this case.
167 LLt = true;
168 else if (tmp.compare ("upper") == 0)
169 LLt = false;
170 else
171 error ("chol: unexpected second or third input");
172 }
173 else
174 error ("chol: expecting trailing string arguments");
175 }
176
177 if (! error_state)
178 {
179 octave_value arg = args(0);
180
181 octave_idx_type nr = arg.rows ();
182 octave_idx_type nc = arg.columns ();
183 bool natural = (nargout != 3);
184
185 int arg_is_empty = empty_arg ("chol", nr, nc);
186
187 if (arg_is_empty < 0)
188 return retval;
189 if (arg_is_empty > 0)
190 return octave_value (Matrix ());
191
192 if (arg.is_sparse_type ())
193 {
194 if (arg.is_real_type ())
195 {
196 SparseMatrix m = arg.sparse_matrix_value ();
197
198 if (! error_state)
199 {
200 octave_idx_type info = nargout;
201 SparseCHOL fact (m, info, natural);
202 if (nargout == 3)
203 {
204 if (vecout)
205 retval(2) = fact.perm ();
206 else
207 retval(2) = fact.Q ();
208 }
209
210 if (nargout > 1 || info == 0)
211 {
212 retval(1) = fact.P ();
213 if (LLt)
214 retval(0) = fact.L ();
215 else
216 retval(0) = fact.R ();
217 }
218 else
219 error ("chol: input matrix must be positive definite");
220 }
221 }
222 else if (arg.is_complex_type ())
223 {
224 SparseComplexMatrix m = arg.sparse_complex_matrix_value ();
225
226 if (! error_state)
227 {
228 octave_idx_type info = nargout;
229 SparseComplexCHOL fact (m, info, natural);
230
231 if (nargout == 3)
232 {
233 if (vecout)
234 retval(2) = fact.perm ();
235 else
236 retval(2) = fact.Q ();
237 }
238
239 if (nargout > 1 || info == 0)
240 {
241 retval(1) = fact.P ();
242 if (LLt)
243 retval(0) = fact.L ();
244 else
245 retval(0) = fact.R ();
246 }
247 else
248 error ("chol: input matrix must be positive definite");
249 }
250 }
251 else
252 gripe_wrong_type_arg ("chol", arg);
253 }
254 else if (arg.is_single_type ())
255 {
256 if (arg.is_real_type ())
257 {
258 FloatMatrix m = arg.float_matrix_value ();
259
260 if (! error_state)
261 {
262 octave_idx_type info;
263
264 FloatCHOL fact;
265 if (LLt)
266 fact = FloatCHOL (m.transpose (), info);
267 else
268 fact = FloatCHOL (m, info);
269
270 if (nargout == 2 || info == 0)
271 {
272 retval(1) = info;
273 if (LLt)
274 retval(0) = get_chol_l (fact);
275 else
276 retval(0) = get_chol_r (fact);
277 }
278 else
279 error ("chol: input matrix must be positive definite");
280 }
281 }
282 else if (arg.is_complex_type ())
283 {
284 FloatComplexMatrix m = arg.float_complex_matrix_value ();
285
286 if (! error_state)
287 {
288 octave_idx_type info;
289
290 FloatComplexCHOL fact;
291 if (LLt)
292 fact = FloatComplexCHOL (m.transpose (), info);
293 else
294 fact = FloatComplexCHOL (m, info);
295
296 if (nargout == 2 || info == 0)
297 {
298 retval(1) = info;
299 if (LLt)
300 retval(0) = get_chol_l (fact);
301 else
302 retval(0) = get_chol_r (fact);
303 }
304 else
305 error ("chol: input matrix must be positive definite");
306 }
307 }
308 else
309 gripe_wrong_type_arg ("chol", arg);
310 }
311 else
312 {
313 if (arg.is_real_type ())
314 {
315 Matrix m = arg.matrix_value ();
316
317 if (! error_state)
318 {
319 octave_idx_type info;
320
321 CHOL fact;
322 if (LLt)
323 fact = CHOL (m.transpose (), info);
324 else
325 fact = CHOL (m, info);
326
327 if (nargout == 2 || info == 0)
328 {
329 retval(1) = info;
330 if (LLt)
331 retval(0) = get_chol_l (fact);
332 else
333 retval(0) = get_chol_r (fact);
334 }
335 else
336 error ("chol: input matrix must be positive definite");
337 }
338 }
339 else if (arg.is_complex_type ())
340 {
341 ComplexMatrix m = arg.complex_matrix_value ();
342
343 if (! error_state)
344 {
345 octave_idx_type info;
346
347 ComplexCHOL fact;
348 if (LLt)
349 fact = ComplexCHOL (m.transpose (), info);
350 else
351 fact = ComplexCHOL (m, info);
352
353 if (nargout == 2 || info == 0)
354 {
355 retval(1) = info;
356 if (LLt)
357 retval(0) = get_chol_l (fact);
358 else
359 retval(0) = get_chol_r (fact);
360 }
361 else
362 error ("chol: input matrix must be positive definite");
363 }
364 }
365 else
366 gripe_wrong_type_arg ("chol", arg);
367 }
368 }
369
370 return retval;
371 }
372
373 /*
374 %!assert (chol ([2, 1; 1, 1]), [sqrt(2), 1/sqrt(2); 0, 1/sqrt(2)], sqrt (eps))
375 %!assert (chol (single ([2, 1; 1, 1])), single ([sqrt(2), 1/sqrt(2); 0, 1/sqrt(2)]), sqrt (eps ("single")))
376
377 %!error chol ()
378 %!error <matrix must be positive definite> chol ([1, 2; 3, 4])
379 %!error <requires square matrix> chol ([1, 2; 3, 4; 5, 6])
380 %!error <unexpected second or third input> chol (1, 2)
381 */
382
383 DEFUN_DLD (cholinv, args, ,
384 "-*- texinfo -*-\n\
385 @deftypefn {Loadable Function} {} cholinv (@var{A})\n\
386 Use the Cholesky@tie{}factorization to compute the inverse of the\n\
387 symmetric positive definite matrix @var{A}.\n\
388 @seealso{chol, chol2inv, inv}\n\
389 @end deftypefn")
390 {
391 octave_value retval;
392
393 int nargin = args.length ();
394
395 if (nargin == 1)
396 {
397 octave_value arg = args(0);
398
399 octave_idx_type nr = arg.rows ();
400 octave_idx_type nc = arg.columns ();
401
402 if (nr == 0 || nc == 0)
403 retval = Matrix ();
404 else
405 {
406 if (arg.is_sparse_type ())
407 {
408 if (arg.is_real_type ())
409 {
410 SparseMatrix m = arg.sparse_matrix_value ();
411
412 if (! error_state)
413 {
414 octave_idx_type info;
415 SparseCHOL chol (m, info);
416 if (info == 0)
417 retval = chol.inverse ();
418 else
419 error ("cholinv: A must be positive definite");
420 }
421 }
422 else if (arg.is_complex_type ())
423 {
424 SparseComplexMatrix m = arg.sparse_complex_matrix_value ();
425
426 if (! error_state)
427 {
428 octave_idx_type info;
429 SparseComplexCHOL chol (m, info);
430 if (info == 0)
431 retval = chol.inverse ();
432 else
433 error ("cholinv: A must be positive definite");
434 }
435 }
436 else
437 gripe_wrong_type_arg ("cholinv", arg);
438 }
439 else if (arg.is_single_type ())
440 {
441 if (arg.is_real_type ())
442 {
443 FloatMatrix m = arg.float_matrix_value ();
444
445 if (! error_state)
446 {
447 octave_idx_type info;
448 FloatCHOL chol (m, info);
449 if (info == 0)
450 retval = chol.inverse ();
451 else
452 error ("cholinv: A must be positive definite");
453 }
454 }
455 else if (arg.is_complex_type ())
456 {
457 FloatComplexMatrix m = arg.float_complex_matrix_value ();
458
459 if (! error_state)
460 {
461 octave_idx_type info;
462 FloatComplexCHOL chol (m, info);
463 if (info == 0)
464 retval = chol.inverse ();
465 else
466 error ("cholinv: A must be positive definite");
467 }
468 }
469 else
470 gripe_wrong_type_arg ("chol", arg);
471 }
472 else
473 {
474 if (arg.is_real_type ())
475 {
476 Matrix m = arg.matrix_value ();
477
478 if (! error_state)
479 {
480 octave_idx_type info;
481 CHOL chol (m, info);
482 if (info == 0)
483 retval = chol.inverse ();
484 else
485 error ("cholinv: A must be positive definite");
486 }
487 }
488 else if (arg.is_complex_type ())
489 {
490 ComplexMatrix m = arg.complex_matrix_value ();
491
492 if (! error_state)
493 {
494 octave_idx_type info;
495 ComplexCHOL chol (m, info);
496 if (info == 0)
497 retval = chol.inverse ();
498 else
499 error ("cholinv: A must be positive definite");
500 }
501 }
502 else
503 gripe_wrong_type_arg ("chol", arg);
504 }
505 }
506 }
507 else
508 print_usage ();
509
510 return retval;
511 }
512
513 /*
514 %!shared A, Ainv
515 %! A = [2,0.2;0.2,1];
516 %! Ainv = inv (A);
517 %!test
518 %! Ainv1 = cholinv (A);
519 %! assert (norm (Ainv-Ainv1), 0, 1e-10);
520 %!testif HAVE_CHOLMOD
521 %! Ainv2 = inv (sparse (A));
522 %! assert (norm (Ainv-Ainv2), 0, 1e-10);
523 %!testif HAVE_CHOLMOD
524 %! Ainv3 = cholinv (sparse (A));
525 %! assert (norm (Ainv-Ainv3), 0, 1e-10);
526 */
527
528 DEFUN_DLD (chol2inv, args, ,
529 "-*- texinfo -*-\n\
530 @deftypefn {Loadable Function} {} chol2inv (@var{U})\n\
531 Invert a symmetric, positive definite square matrix from its Cholesky\n\
532 decomposition, @var{U}. Note that @var{U} should be an upper-triangular\n\
533 matrix with positive diagonal elements. @code{chol2inv (@var{U})}\n\
534 provides @code{inv (@var{U}'*@var{U})} but it is much faster than\n\
535 using @code{inv}.\n\
536 @seealso{chol, cholinv, inv}\n\
537 @end deftypefn")
538 {
539 octave_value retval;
540
541 int nargin = args.length ();
542
543 if (nargin == 1)
544 {
545 octave_value arg = args(0);
546
547 octave_idx_type nr = arg.rows ();
548 octave_idx_type nc = arg.columns ();
549
550 if (nr == 0 || nc == 0)
551 retval = Matrix ();
552 else
553 {
554 if (arg.is_sparse_type ())
555 {
556 if (arg.is_real_type ())
557 {
558 SparseMatrix r = arg.sparse_matrix_value ();
559
560 if (! error_state)
561 retval = chol2inv (r);
562 }
563 else if (arg.is_complex_type ())
564 {
565 SparseComplexMatrix r = arg.sparse_complex_matrix_value ();
566
567 if (! error_state)
568 retval = chol2inv (r);
569 }
570 else
571 gripe_wrong_type_arg ("chol2inv", arg);
572 }
573 else if (arg.is_single_type ())
574 {
575 if (arg.is_real_type ())
576 {
577 FloatMatrix r = arg.float_matrix_value ();
578
579 if (! error_state)
580 retval = chol2inv (r);
581 }
582 else if (arg.is_complex_type ())
583 {
584 FloatComplexMatrix r = arg.float_complex_matrix_value ();
585
586 if (! error_state)
587 retval = chol2inv (r);
588 }
589 else
590 gripe_wrong_type_arg ("chol2inv", arg);
591
592 }
593 else
594 {
595 if (arg.is_real_type ())
596 {
597 Matrix r = arg.matrix_value ();
598
599 if (! error_state)
600 retval = chol2inv (r);
601 }
602 else if (arg.is_complex_type ())
603 {
604 ComplexMatrix r = arg.complex_matrix_value ();
605
606 if (! error_state)
607 retval = chol2inv (r);
608 }
609 else
610 gripe_wrong_type_arg ("chol2inv", arg);
611 }
612 }
613 }
614 else
615 print_usage ();
616
617 return retval;
618 }
619
620 DEFUN_DLD (cholupdate, args, nargout,
621 "-*- texinfo -*-\n\
622 @deftypefn {Loadable Function} {[@var{R1}, @var{info}] =} cholupdate (@var{R}, @var{u}, @var{op})\n\
623 Update or downdate a Cholesky@tie{}factorization. Given an upper triangular\n\
624 matrix @var{R} and a column vector @var{u}, attempt to determine another\n\
625 upper triangular matrix @var{R1} such that\n\
626 \n\
627 @itemize @bullet\n\
628 @item\n\
629 @var{R1}'*@var{R1} = @var{R}'*@var{R} + @var{u}*@var{u}'\n\
630 if @var{op} is \"+\"\n\
631 \n\
632 @item\n\
633 @var{R1}'*@var{R1} = @var{R}'*@var{R} - @var{u}*@var{u}'\n\
634 if @var{op} is \"-\"\n\
635 @end itemize\n\
636 \n\
637 If @var{op} is \"-\", @var{info} is set to\n\
638 \n\
639 @itemize\n\
640 @item 0 if the downdate was successful,\n\
641 \n\
642 @item 1 if @var{R}'*@var{R} - @var{u}*@var{u}' is not positive definite,\n\
643 \n\
644 @item 2 if @var{R} is singular.\n\
645 @end itemize\n\
646 \n\
647 If @var{info} is not present, an error message is printed in cases 1 and 2.\n\
648 @seealso{chol, qrupdate}\n\
649 @end deftypefn")
650 {
651 octave_idx_type nargin = args.length ();
652
653 octave_value_list retval;
654
655 if (nargin > 3 || nargin < 2)
656 {
657 print_usage ();
658 return retval;
659 }
660
661 octave_value argr = args(0);
662 octave_value argu = args(1);
663
664 if (argr.is_numeric_type () && argu.is_numeric_type ()
665 && (nargin < 3 || args(2).is_string ()))
666 {
667 octave_idx_type n = argr.rows ();
668
669 std::string op = (nargin < 3) ? "+" : args(2).string_value ();
670
671 bool down = op == "-";
672
673 if (down || op == "+")
674 if (argr.columns () == n && argu.rows () == n && argu.columns () == 1)
675 {
676 int err = 0;
677 if (argr.is_single_type () || argu.is_single_type ())
678 {
679 if (argr.is_real_type () && argu.is_real_type ())
680 {
681 // real case
682 FloatMatrix R = argr.float_matrix_value ();
683 FloatColumnVector u = argu.float_column_vector_value ();
684
685 FloatCHOL fact;
686 fact.set (R);
687
688 if (down)
689 err = fact.downdate (u);
690 else
691 fact.update (u);
692
693 retval(0) = get_chol_r (fact);
694 }
695 else
696 {
697 // complex case
698 FloatComplexMatrix R = argr.float_complex_matrix_value ();
699 FloatComplexColumnVector u = argu.float_complex_column_vector_value ();
700
701 FloatComplexCHOL fact;
702 fact.set (R);
703
704 if (down)
705 err = fact.downdate (u);
706 else
707 fact.update (u);
708
709 retval(0) = get_chol_r (fact);
710 }
711 }
712 else
713 {
714 if (argr.is_real_type () && argu.is_real_type ())
715 {
716 // real case
717 Matrix R = argr.matrix_value ();
718 ColumnVector u = argu.column_vector_value ();
719
720 CHOL fact;
721 fact.set (R);
722
723 if (down)
724 err = fact.downdate (u);
725 else
726 fact.update (u);
727
728 retval(0) = get_chol_r (fact);
729 }
730 else
731 {
732 // complex case
733 ComplexMatrix R = argr.complex_matrix_value ();
734 ComplexColumnVector u = argu.complex_column_vector_value ();
735
736 ComplexCHOL fact;
737 fact.set (R);
738
739 if (down)
740 err = fact.downdate (u);
741 else
742 fact.update (u);
743
744 retval(0) = get_chol_r (fact);
745 }
746 }
747
748 if (nargout > 1)
749 retval(1) = err;
750 else if (err == 1)
751 error ("cholupdate: downdate violates positiveness");
752 else if (err == 2)
753 error ("cholupdate: singular matrix");
754 }
755 else
756 error ("cholupdate: dimension mismatch between R and U");
757 else
758 error ("cholupdate: OP must be \"+\" or \"-\"");
759 }
760 else
761 print_usage ();
762
763 return retval;
764 }
765
766 /*
767 %!shared A, u, Ac, uc
768 %! A = [ 0.436997 -0.131721 0.124120 -0.061673 ;
769 %! -0.131721 0.738529 0.019851 -0.140295 ;
770 %! 0.124120 0.019851 0.354879 -0.059472 ;
771 %! -0.061673 -0.140295 -0.059472 0.600939 ];
772 %!
773 %! u = [ 0.98950 ;
774 %! 0.39844 ;
775 %! 0.63484 ;
776 %! 0.13351 ];
777 %! Ac = [ 0.5585528 + 0.0000000i -0.1662088 - 0.0315341i 0.0107873 + 0.0236411i -0.0276775 - 0.0186073i ;
778 %! -0.1662088 + 0.0315341i 0.6760061 + 0.0000000i 0.0011452 - 0.0475528i 0.0145967 + 0.0247641i ;
779 %! 0.0107873 - 0.0236411i 0.0011452 + 0.0475528i 0.6263149 - 0.0000000i -0.1585837 - 0.0719763i ;
780 %! -0.0276775 + 0.0186073i 0.0145967 - 0.0247641i -0.1585837 + 0.0719763i 0.6034234 - 0.0000000i ];
781 %!
782 %! uc = [ 0.54267 + 0.91519i ;
783 %! 0.99647 + 0.43141i ;
784 %! 0.83760 + 0.68977i ;
785 %! 0.39160 + 0.90378i ];
786
787 %!test
788 %! R = chol (A);
789 %! R1 = cholupdate (R, u);
790 %! assert (norm (triu (R1)-R1, Inf), 0);
791 %! assert (norm (R1'*R1 - R'*R - u*u', Inf) < 1e1*eps);
792 %!
793 %! R1 = cholupdate (R1, u, "-");
794 %! assert (norm (triu (R1)-R1, Inf), 0);
795 %! assert (norm (R1 - R, Inf) < 1e1*eps);
796
797 %!test
798 %! R = chol (Ac);
799 %! R1 = cholupdate (R, uc);
800 %! assert (norm (triu (R1)-R1, Inf), 0);
801 %! assert (norm (R1'*R1 - R'*R - uc*uc', Inf) < 1e1*eps);
802 %!
803 %! R1 = cholupdate (R1, uc, "-");
804 %! assert (norm (triu (R1)-R1, Inf), 0);
805 %! assert (norm (R1 - R, Inf) < 1e1*eps);
806
807 %!test
808 %! R = chol (single (A));
809 %! R1 = cholupdate (R, single (u));
810 %! assert (norm (triu (R1)-R1, Inf), single (0));
811 %! assert (norm (R1'*R1 - R'*R - single (u*u'), Inf) < 1e1*eps ("single"));
812 %!
813 %! R1 = cholupdate (R1, single (u), "-");
814 %! assert (norm (triu (R1)-R1, Inf), single (0));
815 %! assert (norm (R1 - R, Inf) < 2e1*eps ("single"));
816
817 %!test
818 %! R = chol (single (Ac));
819 %! R1 = cholupdate (R, single (uc));
820 %! assert (norm (triu (R1)-R1, Inf), single (0));
821 %! assert (norm (R1'*R1 - R'*R - single (uc*uc'), Inf) < 1e1*eps ("single"));
822 %!
823 %! R1 = cholupdate (R1, single (uc), "-");
824 %! assert (norm (triu (R1)-R1, Inf), single (0));
825 %! assert (norm (R1 - R, Inf) < 2e1*eps ("single"));
826 */
827
828 DEFUN_DLD (cholinsert, args, nargout,
829 "-*- texinfo -*-\n\
830 @deftypefn {Loadable Function} {@var{R1} =} cholinsert (@var{R}, @var{j}, @var{u})\n\
831 @deftypefnx {Loadable Function} {[@var{R1}, @var{info}] =} cholinsert (@var{R}, @var{j}, @var{u})\n\
832 Given a Cholesky@tie{}factorization of a real symmetric or complex Hermitian\n\
833 positive definite matrix @w{@var{A} = @var{R}'*@var{R}}, @var{R}@tie{}upper\n\
834 triangular, return the Cholesky@tie{}factorization of\n\
835 @var{A1}, where @w{A1(p,p) = A}, @w{A1(:,j) = A1(j,:)' = u} and\n\
836 @w{p = [1:j-1,j+1:n+1]}. @w{u(j)} should be positive.\n\
837 On return, @var{info} is set to\n\
838 \n\
839 @itemize\n\
840 @item 0 if the insertion was successful,\n\
841 \n\
842 @item 1 if @var{A1} is not positive definite,\n\
843 \n\
844 @item 2 if @var{R} is singular.\n\
845 @end itemize\n\
846 \n\
847 If @var{info} is not present, an error message is printed in cases 1 and 2.\n\
848 @seealso{chol, cholupdate, choldelete}\n\
849 @end deftypefn")
850 {
851 octave_idx_type nargin = args.length ();
852
853 octave_value_list retval;
854
855 if (nargin != 3)
856 {
857 print_usage ();
858 return retval;
859 }
860
861 octave_value argr = args(0);
862 octave_value argj = args(1);
863 octave_value argu = args(2);
864
865 if (argr.is_numeric_type () && argu.is_numeric_type ()
866 && argj.is_real_scalar ())
867 {
868 octave_idx_type n = argr.rows ();
869 octave_idx_type j = argj.scalar_value ();
870
871 if (argr.columns () == n && argu.rows () == n+1 && argu.columns () == 1)
872 {
873 if (j > 0 && j <= n+1)
874 {
875 int err = 0;
876 if (argr.is_single_type () || argu.is_single_type ())
877 {
878 if (argr.is_real_type () && argu.is_real_type ())
879 {
880 // real case
881 FloatMatrix R = argr.float_matrix_value ();
882 FloatColumnVector u = argu.float_column_vector_value ();
883
884 FloatCHOL fact;
885 fact.set (R);
886 err = fact.insert_sym (u, j-1);
887
888 retval(0) = get_chol_r (fact);
889 }
890 else
891 {
892 // complex case
893 FloatComplexMatrix R = argr.float_complex_matrix_value ();
894 FloatComplexColumnVector u = argu.float_complex_column_vector_value ();
895
896 FloatComplexCHOL fact;
897 fact.set (R);
898 err = fact.insert_sym (u, j-1);
899
900 retval(0) = get_chol_r (fact);
901 }
902 }
903 else
904 {
905 if (argr.is_real_type () && argu.is_real_type ())
906 {
907 // real case
908 Matrix R = argr.matrix_value ();
909 ColumnVector u = argu.column_vector_value ();
910
911 CHOL fact;
912 fact.set (R);
913 err = fact.insert_sym (u, j-1);
914
915 retval(0) = get_chol_r (fact);
916 }
917 else
918 {
919 // complex case
920 ComplexMatrix R = argr.complex_matrix_value ();
921 ComplexColumnVector u = argu.complex_column_vector_value ();
922
923 ComplexCHOL fact;
924 fact.set (R);
925 err = fact.insert_sym (u, j-1);
926
927 retval(0) = get_chol_r (fact);
928 }
929 }
930
931 if (nargout > 1)
932 retval(1) = err;
933 else if (err == 1)
934 error ("cholinsert: insertion violates positiveness");
935 else if (err == 2)
936 error ("cholinsert: singular matrix");
937 else if (err == 3)
938 error ("cholinsert: diagonal element must be real");
939 }
940 else
941 error ("cholinsert: index J out of range");
942 }
943 else
944 error ("cholinsert: dimension mismatch between R and U");
945 }
946 else
947 print_usage ();
948
949 return retval;
950 }
951
952 /*
953 %!test
954 %! u2 = [ 0.35080 ;
955 %! 0.63930 ;
956 %! 3.31057 ;
957 %! -0.13825 ;
958 %! 0.45266 ];
959 %!
960 %! R = chol (A);
961 %!
962 %! j = 3; p = [1:j-1, j+1:5];
963 %! R1 = cholinsert (R, j, u2);
964 %! A1 = R1'*R1;
965 %!
966 %! assert (norm (triu (R1)-R1, Inf), 0);
967 %! assert (norm (A1(p,p) - A, Inf) < 1e1*eps);
968
969 %!test
970 %! u2 = [ 0.35080 + 0.04298i;
971 %! 0.63930 + 0.23778i;
972 %! 3.31057 + 0.00000i;
973 %! -0.13825 + 0.19879i;
974 %! 0.45266 + 0.50020i];
975 %!
976 %! R = chol (Ac);
977 %!
978 %! j = 3; p = [1:j-1, j+1:5];
979 %! R1 = cholinsert (R, j, u2);
980 %! A1 = R1'*R1;
981 %!
982 %! assert (norm (triu (R1)-R1, Inf), 0);
983 %! assert (norm (A1(p,p) - Ac, Inf) < 1e1*eps);
984
985 %!test
986 %! u2 = single ([ 0.35080 ;
987 %! 0.63930 ;
988 %! 3.31057 ;
989 %! -0.13825 ;
990 %! 0.45266 ]);
991 %!
992 %! R = chol (single (A));
993 %!
994 %! j = 3; p = [1:j-1, j+1:5];
995 %! R1 = cholinsert (R, j, u2);
996 %! A1 = R1'*R1;
997 %!
998 %! assert (norm (triu (R1)-R1, Inf), single (0));
999 %! assert (norm (A1(p,p) - A, Inf) < 1e1*eps ("single"));
1000
1001 %!test
1002 %! u2 = single ([ 0.35080 + 0.04298i;
1003 %! 0.63930 + 0.23778i;
1004 %! 3.31057 + 0.00000i;
1005 %! -0.13825 + 0.19879i;
1006 %! 0.45266 + 0.50020i]);
1007 %!
1008 %! R = chol (single (Ac));
1009 %!
1010 %! j = 3; p = [1:j-1, j+1:5];
1011 %! R1 = cholinsert (R, j, u2);
1012 %! A1 = R1'*R1;
1013 %!
1014 %! assert (norm (triu (R1)-R1, Inf), single (0));
1015 %! assert (norm (A1(p,p) - single (Ac), Inf) < 2e1*eps ("single"));
1016
1017 %!test
1018 %! cu = chol (triu (A), "upper");
1019 %! cl = chol (tril (A), "lower");
1020 %! assert (cu, cl', eps);
1021
1022 %!test
1023 %! cca = chol (Ac);
1024 %!
1025 %! ccal = chol (Ac, "lower");
1026 %! ccal2 = chol (tril (Ac), "lower");
1027 %!
1028 %! ccau = chol (Ac, "upper");
1029 %! ccau2 = chol (triu (Ac), "upper");
1030 %!
1031 %! assert (cca'*cca, Ac, eps);
1032 %! assert (ccau'*ccau, Ac, eps);
1033 %! assert (ccau2'*ccau2, Ac, eps);
1034 %!
1035 %! assert (cca, ccal', eps);
1036 %! assert (cca, ccau, eps);
1037 %! assert (cca, ccal2', eps);
1038 %! assert (cca, ccau2, eps);
1039
1040 %!test
1041 %! cca = chol (single (Ac));
1042 %!
1043 %! ccal = chol (single (Ac), "lower");
1044 %! ccal2 = chol (tril (single (Ac)), "lower");
1045 %!
1046 %! ccau = chol (single (Ac), "upper");
1047 %! ccau2 = chol (triu (single (Ac)), "upper");
1048 %!
1049 %! assert (cca'*cca, single (Ac), eps ("single"));
1050 %! assert (ccau'*ccau, single (Ac), eps ("single"));
1051 %! assert (ccau2'*ccau2, single (Ac), eps ("single"));
1052 %!
1053 %! assert (cca, ccal', eps ("single"));
1054 %! assert (cca, ccau, eps ("single"));
1055 %! assert (cca, ccal2', eps ("single"));
1056 %! assert (cca, ccau2, eps ("single"));
1057
1058 %!test
1059 %! a = [12, 2, 3, 4;
1060 %! 2, 14, 5, 3;
1061 %! 3, 5, 16, 6;
1062 %! 4, 3, 6, 16];
1063 %!
1064 %! b = [0, 1, 2, 3;
1065 %! -1, 0, 1, 2;
1066 %! -2, -1, 0, 1;
1067 %! -3, -2, -1, 0];
1068 %!
1069 %! ca = a + i*b;
1070 %!
1071 %! cca = chol (ca);
1072 %!
1073 %! ccal = chol (ca, "lower");
1074 %! ccal2 = chol (tril (ca), "lower");
1075 %!
1076 %! ccau = chol (ca, "upper");
1077 %! ccau2 = chol (triu (ca), "upper");
1078 %!
1079 %! assert (cca'*cca, ca, 16*eps);
1080 %! assert (ccau'*ccau, ca, 16*eps);
1081 %! assert (ccau2'*ccau2, ca, 16*eps);
1082 %!
1083 %! assert (cca, ccal', 16*eps);
1084 %! assert (cca, ccau, 16*eps);
1085 %! assert (cca, ccal2', 16*eps);
1086 %! assert (cca, ccau2, 16*eps);
1087 */
1088
1089 DEFUN_DLD (choldelete, args, ,
1090 "-*- texinfo -*-\n\
1091 @deftypefn {Loadable Function} {@var{R1} =} choldelete (@var{R}, @var{j})\n\
1092 Given a Cholesky@tie{}factorization of a real symmetric or complex Hermitian\n\
1093 positive definite matrix @w{@var{A} = @var{R}'*@var{R}}, @var{R}@tie{}upper\n\
1094 triangular, return the Cholesky@tie{}factorization of @w{A(p,p)}, where\n\
1095 @w{p = [1:j-1,j+1:n+1]}.\n\
1096 @seealso{chol, cholupdate, cholinsert}\n\
1097 @end deftypefn")
1098 {
1099 octave_idx_type nargin = args.length ();
1100
1101 octave_value_list retval;
1102
1103 if (nargin != 2)
1104 {
1105 print_usage ();
1106 return retval;
1107 }
1108
1109 octave_value argr = args(0);
1110 octave_value argj = args(1);
1111
1112 if (argr.is_numeric_type () && argj.is_real_scalar ())
1113 {
1114 octave_idx_type n = argr.rows ();
1115 octave_idx_type j = argj.scalar_value ();
1116
1117 if (argr.columns () == n)
1118 {
1119 if (j > 0 && j <= n)
1120 {
1121 if (argr.is_single_type ())
1122 {
1123 if (argr.is_real_type ())
1124 {
1125 // real case
1126 FloatMatrix R = argr.float_matrix_value ();
1127
1128 FloatCHOL fact;
1129 fact.set (R);
1130 fact.delete_sym (j-1);
1131
1132 retval(0) = get_chol_r (fact);
1133 }
1134 else
1135 {
1136 // complex case
1137 FloatComplexMatrix R = argr.float_complex_matrix_value ();
1138
1139 FloatComplexCHOL fact;
1140 fact.set (R);
1141 fact.delete_sym (j-1);
1142
1143 retval(0) = get_chol_r (fact);
1144 }
1145 }
1146 else
1147 {
1148 if (argr.is_real_type ())
1149 {
1150 // real case
1151 Matrix R = argr.matrix_value ();
1152
1153 CHOL fact;
1154 fact.set (R);
1155 fact.delete_sym (j-1);
1156
1157 retval(0) = get_chol_r (fact);
1158 }
1159 else
1160 {
1161 // complex case
1162 ComplexMatrix R = argr.complex_matrix_value ();
1163
1164 ComplexCHOL fact;
1165 fact.set (R);
1166 fact.delete_sym (j-1);
1167
1168 retval(0) = get_chol_r (fact);
1169 }
1170 }
1171 }
1172 else
1173 error ("choldelete: index J out of range");
1174 }
1175 else
1176 error ("choldelete: matrix R must be square");
1177 }
1178 else
1179 print_usage ();
1180
1181 return retval;
1182 }
1183
1184 /*
1185 %!test
1186 %! R = chol (A);
1187 %!
1188 %! j = 3; p = [1:j-1,j+1:4];
1189 %! R1 = choldelete (R, j);
1190 %!
1191 %! assert (norm (triu (R1)-R1, Inf), 0);
1192 %! assert (norm (R1'*R1 - A(p,p), Inf) < 1e1*eps);
1193
1194 %!test
1195 %! R = chol (Ac);
1196 %!
1197 %! j = 3; p = [1:j-1,j+1:4];
1198 %! R1 = choldelete (R, j);
1199 %!
1200 %! assert (norm (triu (R1)-R1, Inf), 0);
1201 %! assert (norm (R1'*R1 - Ac(p,p), Inf) < 1e1*eps);
1202
1203 %!test
1204 %! R = chol (single (A));
1205 %!
1206 %! j = 3; p = [1:j-1,j+1:4];
1207 %! R1 = choldelete (R, j);
1208 %!
1209 %! assert (norm (triu (R1)-R1, Inf), single (0));
1210 %! assert (norm (R1'*R1 - single (A(p,p)), Inf) < 1e1*eps ("single"));
1211
1212 %!test
1213 %! R = chol (single (Ac));
1214 %!
1215 %! j = 3; p = [1:j-1,j+1:4];
1216 %! R1 = choldelete (R,j);
1217 %!
1218 %! assert (norm (triu (R1)-R1, Inf), single (0));
1219 %! assert (norm (R1'*R1 - single (Ac(p,p)), Inf) < 1e1*eps ("single"));
1220 */
1221
1222 DEFUN_DLD (cholshift, args, ,
1223 "-*- texinfo -*-\n\
1224 @deftypefn {Loadable Function} {@var{R1} =} cholshift (@var{R}, @var{i}, @var{j})\n\
1225 Given a Cholesky@tie{}factorization of a real symmetric or complex Hermitian\n\
1226 positive definite matrix @w{@var{A} = @var{R}'*@var{R}}, @var{R}@tie{}upper\n\
1227 triangular, return the Cholesky@tie{}factorization of\n\
1228 @w{@var{A}(p,p)}, where @w{p} is the permutation @*\n\
1229 @code{p = [1:i-1, shift(i:j, 1), j+1:n]} if @w{@var{i} < @var{j}} @*\n\
1230 or @*\n\
1231 @code{p = [1:j-1, shift(j:i,-1), i+1:n]} if @w{@var{j} < @var{i}}. @*\n\
1232 \n\
1233 @seealso{chol, cholinsert, choldelete}\n\
1234 @end deftypefn")
1235 {
1236 octave_idx_type nargin = args.length ();
1237
1238 octave_value_list retval;
1239
1240 if (nargin != 3)
1241 {
1242 print_usage ();
1243 return retval;
1244 }
1245
1246 octave_value argr = args(0);
1247 octave_value argi = args(1);
1248 octave_value argj = args(2);
1249
1250 if (argr.is_numeric_type () && argi.is_real_scalar () && argj.is_real_scalar ())
1251 {
1252 octave_idx_type n = argr.rows ();
1253 octave_idx_type i = argi.scalar_value ();
1254 octave_idx_type j = argj.scalar_value ();
1255
1256 if (argr.columns () == n)
1257 {
1258 if (j > 0 && j <= n+1 && i > 0 && i <= n+1)
1259 {
1260
1261 if (argr.is_single_type () && argi.is_single_type () &&
1262 argj.is_single_type ())
1263 {
1264 if (argr.is_real_type ())
1265 {
1266 // real case
1267 FloatMatrix R = argr.float_matrix_value ();
1268
1269 FloatCHOL fact;
1270 fact.set (R);
1271 fact.shift_sym (i-1, j-1);
1272
1273 retval(0) = get_chol_r (fact);
1274 }
1275 else
1276 {
1277 // complex case
1278 FloatComplexMatrix R = argr.float_complex_matrix_value ();
1279
1280 FloatComplexCHOL fact;
1281 fact.set (R);
1282 fact.shift_sym (i-1, j-1);
1283
1284 retval(0) = get_chol_r (fact);
1285 }
1286 }
1287 else
1288 {
1289 if (argr.is_real_type ())
1290 {
1291 // real case
1292 Matrix R = argr.matrix_value ();
1293
1294 CHOL fact;
1295 fact.set (R);
1296 fact.shift_sym (i-1, j-1);
1297
1298 retval(0) = get_chol_r (fact);
1299 }
1300 else
1301 {
1302 // complex case
1303 ComplexMatrix R = argr.complex_matrix_value ();
1304
1305 ComplexCHOL fact;
1306 fact.set (R);
1307 fact.shift_sym (i-1, j-1);
1308
1309 retval(0) = get_chol_r (fact);
1310 }
1311 }
1312 }
1313 else
1314 error ("cholshift: index I or J is out of range");
1315 }
1316 else
1317 error ("cholshift: R must be a square matrix");
1318 }
1319 else
1320 print_usage ();
1321
1322 return retval;
1323 }
1324
1325 /*
1326 %!test
1327 %! R = chol (A);
1328 %!
1329 %! i = 1; j = 3; p = [1:i-1, shift(i:j,-1), j+1:4];
1330 %! R1 = cholshift (R, i, j);
1331 %!
1332 %! assert (norm (triu (R1)-R1, Inf), 0);
1333 %! assert (norm (R1'*R1 - A(p,p), Inf) < 1e1*eps);
1334 %!
1335 %! j = 1; i = 3; p = [1:j-1, shift(j:i,+1), i+1:4];
1336 %! R1 = cholshift (R, i, j);
1337 %!
1338 %! assert (norm (triu (R1) - R1, Inf), 0);
1339 %! assert (norm (R1'*R1 - A(p,p), Inf) < 1e1*eps);
1340
1341 %!test
1342 %! R = chol (Ac);
1343 %!
1344 %! i = 1; j = 3; p = [1:i-1, shift(i:j,-1), j+1:4];
1345 %! R1 = cholshift (R, i, j);
1346 %!
1347 %! assert (norm (triu (R1)-R1, Inf), 0);
1348 %! assert (norm (R1'*R1 - Ac(p,p), Inf) < 1e1*eps);
1349 %!
1350 %! j = 1; i = 3; p = [1:j-1, shift(j:i,+1), i+1:4];
1351 %! R1 = cholshift (R, i, j);
1352 %!
1353 %! assert (norm (triu (R1)-R1, Inf), 0);
1354 %! assert (norm (R1'*R1 - Ac(p,p), Inf) < 1e1*eps);
1355
1356 %!test
1357 %! R = chol (single (A));
1358 %!
1359 %! i = 1; j = 3; p = [1:i-1, shift(i:j,-1), j+1:4];
1360 %! R1 = cholshift (R, i, j);
1361 %!
1362 %! assert (norm (triu (R1)-R1, Inf), 0);
1363 %! assert (norm (R1'*R1 - single (A(p,p)), Inf) < 1e1*eps ("single"));
1364 %!
1365 %! j = 1; i = 3; p = [1:j-1, shift(j:i,+1), i+1:4];
1366 %! R1 = cholshift (R, i, j);
1367 %!
1368 %! assert (norm (triu (R1)-R1, Inf), 0);
1369 %! assert (norm (R1'*R1 - single (A(p,p)), Inf) < 1e1*eps ("single"));
1370
1371 %!test
1372 %! R = chol (single (Ac));
1373 %!
1374 %! i = 1; j = 3; p = [1:i-1, shift(i:j,-1), j+1:4];
1375 %! R1 = cholshift (R, i, j);
1376 %!
1377 %! assert (norm (triu (R1)-R1, Inf), 0);
1378 %! assert (norm (R1'*R1 - single (Ac(p,p)), Inf) < 1e1*eps ("single"));
1379 %!
1380 %! j = 1; i = 3; p = [1:j-1, shift(j:i,+1), i+1:4];
1381 %! R1 = cholshift (R, i, j);
1382 %!
1383 %! assert (norm (triu (R1)-R1, Inf), 0);
1384 %! assert (norm (R1'*R1 - single (Ac(p,p)), Inf) < 1e1*eps ("single"));
1385 */