Mercurial > octave-nkf

comparison src/DLD-FUNCTIONS/chol.cc @ 7554:40574114c514

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
implement Cholesky factorization updating
author Jaroslav Hajek <highegg@gmail.com>
date Tue, 04 Mar 2008 22:25:50 -0500
parents f3c00dc0912b
children 07522d7dcdf8
comparison
equal deleted inserted replaced
7553:56be6f31dd4e 7554:40574114c514
18 You should have received a copy of the GNU General Public License 18 You should have received a copy of the GNU General Public License
19 along with Octave; see the file COPYING. If not, see 19 along with Octave; see the file COPYING. If not, see
20 <http://www.gnu.org/licenses/>. 20 <http://www.gnu.org/licenses/>.
21 21
22 */ 22 */
23
24 // The cholupdate function was written by Jaroslav Hajek
25 // <highegg@gmail.com>, Copyright (C) 2008 VZLU Prague, a.s., Czech
26 // Republic.
23 27
24 #ifdef HAVE_CONFIG_H 28 #ifdef HAVE_CONFIG_H
25 #include <config.h> 29 #include <config.h>
26 #endif 30 #endif
27 31
439 print_usage (); 443 print_usage ();
440 444
441 return retval; 445 return retval;
442 } 446 }
443 447
448 DEFUN_DLD (cholupdate, args, nargout,
449 "-*- texinfo -*-\n\
450 @deftypefn {Loadable Function} {[@var{R1}, @var{info}]} = cholupdate (@var{R}, @var{u}, @var{op})\n\
451 Update or downdate a Cholesky factorization. Given an upper triangular\n\
452 matrix @var{R} and a column vector @var{u}, attempt to determine another\n\
453 upper triangular matrix @var{R1} such that\n\
454 @itemize @bullet\n\
455 @item\n\
456 @var{R1}'*@var{R1} = @var{R}'*@var{R} + @var{u}*@var{u}'\n\
457 if @var{op} is \"+\"\n\
458 @item\n\
459 @var{R1}'*@var{R1} = @var{R}'*@var{R} - @var{u}*@var{u}'\n\
460 if @var{op} is \"-\"\n\
461 @end itemize\n\
462 \n\
463 If @var{op} is \"-\", @var{info} is set to\n\
464 @itemize\n\
465 @item 0 if the downdate was successful,\n\
466 @item 1 if @var{R}'*@var{R} - @var{u}*@var{u}' is not positive definite,\n\
467 @item 2 if @var{R} is singular.\n\
468 @end itemize\n\
469 \n\
470 If @var{info} is not present, an error message is printed in cases 1 and 2.\n\
471 @seealso{chol, qrupdate}\n\
472 @end deftypefn")
473 {
474 int nargin = args.length ();
475
476 octave_value_list retval;
477
478 octave_value argR,argu,argop;
479
480 if ((nargin == 3 || nargin == 2)
481 && (argR = args(0), argR.is_matrix_type ())
482 && (argu = args(1), argu.is_matrix_type ())
483 && (nargin < 3 || (argop = args(2), argop.is_string ())))
484 {
485 octave_idx_type n = argR.rows ();
486
487 std::string op = (nargin < 3) ? "+" : argop.string_value();
488
489 bool down = false;
490
491 if (nargin < 3 || (op == "+") || (down = op == "-"))
492 if (argR.columns () == n
493 && argu.rows () == n && argu.columns () == 1)
494 {
495 if (argR.is_real_matrix () && argu.is_real_matrix ())
496 {
497 // real case
498 Matrix R = argR.matrix_value ();
499 Matrix u = argu.matrix_value ();
500
501 CHOL fact;
502 fact.set (R);
503 int err = 0;
504
505 if (down)
506 err = fact.downdate (u);
507 else
508 fact.update (u);
509
510 if (nargout > 1)
511 retval(1) = err;
512 else if (err)
513 error ("cholupdate: downdate violates positiveness");
514
515 retval(0) = fact.chol_matrix ();
516 }
517 else
518 {
519 // complex case
520 ComplexMatrix R = argR.complex_matrix_value ();
521 ComplexMatrix u = argu.complex_matrix_value ();
522
523 ComplexCHOL fact;
524 fact.set (R);
525 int err = 0;
526
527 if (down)
528 err = fact.downdate (u);
529 else
530 fact.update (u);
531
532 if (nargout > 1)
533 retval(1) = err;
534 else if (err)
535 error ("cholupdate: downdate violates positiveness");
536
537 retval(0) = fact.chol_matrix ();
538 }
539 }
540 else
541 error ("cholupdate: dimension mismatch");
542 else
543 error ("cholupdate: op must be \"+\" or \"-\"");
544 }
545 else
546 print_usage ();
547
548 return retval;
549 }
550
551 /*
552 %!test
553 %! A = [ 0.436997 -0.131721 0.124120 -0.061673 ;
554 %! -0.131721 0.738529 0.019851 -0.140295 ;
555 %! 0.124120 0.019851 0.354879 -0.059472 ;
556 %! -0.061673 -0.140295 -0.059472 0.600939 ];
557 %!
558 %! u = [ 0.98950 ;
559 %! 0.39844 ;
560 %! 0.63484 ;
561 %! 0.13351 ];
562 %!
563 %! R = chol(A);
564 %!
565 %! R1 = cholupdate(R,u);
566 %!
567 %! assert(norm(triu(R1)-R1,Inf) == 0)
568 %! assert(norm(R1'*R1 - R'*R - u*u',Inf) < 1e1*eps)
569 %!
570 %! R1 = cholupdate(R1,u,"-");
571 %!
572 %! assert(norm(triu(R1)-R1,Inf) == 0)
573 %! assert(norm(R1 - R,Inf) < 1e1*eps)
574 %!
575 %!test
576 %! A = [ 0.5585528 + 0.0000000i -0.1662088 - 0.0315341i 0.0107873 + 0.0236411i -0.0276775 - 0.0186073i ;
577 %! -0.1662088 + 0.0315341i 0.6760061 + 0.0000000i 0.0011452 - 0.0475528i 0.0145967 + 0.0247641i ;
578 %! 0.0107873 - 0.0236411i 0.0011452 + 0.0475528i 0.6263149 - 0.0000000i -0.1585837 - 0.0719763i ;
579 %! -0.0276775 + 0.0186073i 0.0145967 - 0.0247641i -0.1585837 + 0.0719763i 0.6034234 - 0.0000000i ];
580 %!
581 %! u = [ 0.54267 + 0.91519i ;
582 %! 0.99647 + 0.43141i ;
583 %! 0.83760 + 0.68977i ;
584 %! 0.39160 + 0.90378i ];
585 %!
586 %! R = chol(A);
587 %!
588 %! R1 = cholupdate(R,u);
589 %!
590 %! assert(norm(triu(R1)-R1,Inf) == 0)
591 %! assert(norm(R1'*R1 - R'*R - u*u',Inf) < 1e1*eps)
592 %!
593 %! R1 = cholupdate(R1,u,"-");
594 %!
595 %! assert(norm(triu(R1)-R1,Inf) == 0)
596 %! assert(norm(R1 - R,Inf) < 1e1*eps)
597 */
598
444 /* 599 /*
445 ;;; Local Variables: *** 600 ;;; Local Variables: ***
446 ;;; mode: C++ *** 601 ;;; mode: C++ ***
447 ;;; End: *** 602 ;;; End: ***
448 */ 603 */