Mercurial > octave-nkf

comparison src/DLD-FUNCTIONS/chol.cc @ 7814:87865ed7405f

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Second set of single precision test code and fix of resulting bugs
author David Bateman <dbateman@free.fr>
date Mon, 02 Jun 2008 16:57:45 +0200
parents 82be108cc558
children 4976f66d469b
comparison
equal deleted inserted replaced
7813:12a68443191c 7814:87865ed7405f
321 } 321 }
322 322
323 return retval; 323 return retval;
324 } 324 }
325 325
326 /*
327
328 %!assert(chol ([2, 1; 1, 1]), [sqrt(2), 1/sqrt(2); 0, 1/sqrt(2)], sqrt (eps));
329 %!assert(chol (single([2, 1; 1, 1])), single([sqrt(2), 1/sqrt(2); 0, 1/sqrt(2)]), sqrt (eps('single')));
330
331 %!error chol ([1, 2; 3, 4]);
332 %!error chol ([1, 2; 3, 4; 5, 6]);
333 %!error <Invalid call to chol.*> chol ();
334 %!error <unexpected second or third input.*> chol (1, 2);
335
336 */
337
326 DEFUN_DLD (cholinv, args, , 338 DEFUN_DLD (cholinv, args, ,
327 "-*- texinfo -*-\n\ 339 "-*- texinfo -*-\n\
328 @deftypefn {Loadable Function} {} cholinv (@var{a})\n\ 340 @deftypefn {Loadable Function} {} cholinv (@var{a})\n\
329 Use the Cholesky factorization to compute the inverse of the\n\ 341 Use the Cholesky factorization to compute the inverse of the\n\
330 symmetric positive definite matrix @var{a}.\n\ 342 symmetric positive definite matrix @var{a}.\n\
376 error ("cholinv: matrix not positive definite"); 388 error ("cholinv: matrix not positive definite");
377 } 389 }
378 } 390 }
379 else 391 else
380 gripe_wrong_type_arg ("cholinv", arg); 392 gripe_wrong_type_arg ("cholinv", arg);
393 }
394 else if (arg.is_single_type ())
395 {
396 if (arg.is_real_type ())
397 {
398 FloatMatrix m = arg.float_matrix_value ();
399
400 if (! error_state)
401 {
402 octave_idx_type info;
403 FloatCHOL chol (m, info);
404 if (info == 0)
405 retval = chol.inverse ();
406 else
407 error ("cholinv: matrix not positive definite");
408 }
409 }
410 else if (arg.is_complex_type ())
411 {
412 FloatComplexMatrix m = arg.float_complex_matrix_value ();
413
414 if (! error_state)
415 {
416 octave_idx_type info;
417 FloatComplexCHOL chol (m, info);
418 if (info == 0)
419 retval = chol.inverse ();
420 else
421 error ("cholinv: matrix not positive definite");
422 }
423 }
424 else
425 gripe_wrong_type_arg ("chol", arg);
381 } 426 }
382 else 427 else
383 { 428 {
384 if (arg.is_real_type ()) 429 if (arg.is_real_type ())
385 { 430 {
462 if (! error_state) 507 if (! error_state)
463 retval = chol2inv (r); 508 retval = chol2inv (r);
464 } 509 }
465 else 510 else
466 gripe_wrong_type_arg ("chol2inv", arg); 511 gripe_wrong_type_arg ("chol2inv", arg);
512 }
513 else if (arg.is_single_type ())
514 {
515 if (arg.is_real_type ())
516 {
517 FloatMatrix r = arg.float_matrix_value ();
518
519 if (! error_state)
520 retval = chol2inv (r);
521 }
522 else if (arg.is_complex_type ())
523 {
524 FloatComplexMatrix r = arg.float_complex_matrix_value ();
525
526 if (! error_state)
527 retval = chol2inv (r);
528 }
529 else
530 gripe_wrong_type_arg ("chol2inv", arg);
531
467 } 532 }
468 else 533 else
469 { 534 {
470 if (arg.is_real_type ()) 535 if (arg.is_real_type ())
471 { 536 {
541 bool down = op == "-"; 606 bool down = op == "-";
542 607
543 if (down || op == "+") 608 if (down || op == "+")
544 if (argr.columns () == n && argu.rows () == n && argu.columns () == 1) 609 if (argr.columns () == n && argu.rows () == n && argu.columns () == 1)
545 { 610 {
546 if (argr.is_real_matrix () && argu.is_real_matrix ()) 611 if (argr.is_single_type () || argu.is_single_type ())
547 { 612 {
548 // real case 613 if (argr.is_real_matrix () && argu.is_real_matrix ())
549 Matrix R = argr.matrix_value (); 614 {
550 Matrix u = argu.matrix_value (); 615 // real case
551 616 FloatMatrix R = argr.float_matrix_value ();
552 CHOL fact; 617 FloatMatrix u = argu.float_matrix_value ();
553 fact.set (R); 618
554 int err = 0; 619 FloatCHOL fact;
555 620 fact.set (R);
556 if (down) 621 int err = 0;
557 err = fact.downdate (u); 622
558 else 623 if (down)
559 fact.update (u); 624 err = fact.downdate (u);
560 625 else
561 if (nargout > 1) 626 fact.update (u);
562 retval(1) = err; 627
563 else if (err) 628 if (nargout > 1)
564 error ("cholupdate: downdate violates positiveness"); 629 retval(1) = err;
565 630 else if (err)
566 retval(0) = fact.chol_matrix (); 631 error ("cholupdate: downdate violates positiveness");
567 } 632
568 else 633 retval(0) = fact.chol_matrix ();
569 { 634 }
570 // complex case 635 else
571 ComplexMatrix R = argr.complex_matrix_value (); 636 {
572 ComplexMatrix u = argu.complex_matrix_value (); 637 // complex case
573 638 FloatComplexMatrix R = argr.float_complex_matrix_value ();
574 ComplexCHOL fact; 639 FloatComplexMatrix u = argu.float_complex_matrix_value ();
575 fact.set (R); 640
576 int err = 0; 641 FloatComplexCHOL fact;
577 642 fact.set (R);
578 if (down) 643 int err = 0;
579 err = fact.downdate (u); 644
580 else 645 if (down)
581 fact.update (u); 646 err = fact.downdate (u);
582 647 else
583 if (nargout > 1) 648 fact.update (u);
584 retval(1) = err; 649
585 else if (err) 650 if (nargout > 1)
586 error ("cholupdate: downdate violates positiveness"); 651 retval(1) = err;
587 652 else if (err)
588 retval(0) = fact.chol_matrix (); 653 error ("cholupdate: downdate violates positiveness");
589 } 654
655 retval(0) = fact.chol_matrix ();
656 }
657 }
658 else
659 {
660 if (argr.is_real_matrix () && argu.is_real_matrix ())
661 {
662 // real case
663 Matrix R = argr.matrix_value ();
664 Matrix u = argu.matrix_value ();
665
666 CHOL fact;
667 fact.set (R);
668 int err = 0;
669
670 if (down)
671 err = fact.downdate (u);
672 else
673 fact.update (u);
674
675 if (nargout > 1)
676 retval(1) = err;
677 else if (err)
678 error ("cholupdate: downdate violates positiveness");
679
680 retval(0) = fact.chol_matrix ();
681 }
682 else
683 {
684 // complex case
685 ComplexMatrix R = argr.complex_matrix_value ();
686 ComplexMatrix u = argu.complex_matrix_value ();
687
688 ComplexCHOL fact;
689 fact.set (R);
690 int err = 0;
691
692 if (down)
693 err = fact.downdate (u);
694 else
695 fact.update (u);
696
697 if (nargout > 1)
698 retval(1) = err;
699 else if (err)
700 error ("cholupdate: downdate violates positiveness");
701
702 retval(0) = fact.chol_matrix ();
703 }
704 }
590 } 705 }
591 else 706 else
592 error ("cholupdate: dimension mismatch"); 707 error ("cholupdate: dimension mismatch");
593 else 708 else
594 error ("cholupdate: op must be \"+\" or \"-\""); 709 error ("cholupdate: op must be \"+\" or \"-\"");
598 713
599 return retval; 714 return retval;
600 } 715 }
601 716
602 /* 717 /*
603 %!test 718 %!shared A, u, Ac, uc
604 %! A = [ 0.436997 -0.131721 0.124120 -0.061673 ; 719 %! A = [ 0.436997 -0.131721 0.124120 -0.061673 ;
605 %! -0.131721 0.738529 0.019851 -0.140295 ; 720 %! -0.131721 0.738529 0.019851 -0.140295 ;
606 %! 0.124120 0.019851 0.354879 -0.059472 ; 721 %! 0.124120 0.019851 0.354879 -0.059472 ;
607 %! -0.061673 -0.140295 -0.059472 0.600939 ]; 722 %! -0.061673 -0.140295 -0.059472 0.600939 ];
608 %! 723 %!
609 %! u = [ 0.98950 ; 724 %! u = [ 0.98950 ;
610 %! 0.39844 ; 725 %! 0.39844 ;
611 %! 0.63484 ; 726 %! 0.63484 ;
612 %! 0.13351 ]; 727 %! 0.13351 ];
613 %! 728 %! Ac = [ 0.5585528 + 0.0000000i -0.1662088 - 0.0315341i 0.0107873 + 0.0236411i -0.0276775 - 0.0186073i ;
729 %! -0.1662088 + 0.0315341i 0.6760061 + 0.0000000i 0.0011452 - 0.0475528i 0.0145967 + 0.0247641i ;
730 %! 0.0107873 - 0.0236411i 0.0011452 + 0.0475528i 0.6263149 - 0.0000000i -0.1585837 - 0.0719763i ;
731 %! -0.0276775 + 0.0186073i 0.0145967 - 0.0247641i -0.1585837 + 0.0719763i 0.6034234 - 0.0000000i ];
732 %!
733 %! uc = [ 0.54267 + 0.91519i ;
734 %! 0.99647 + 0.43141i ;
735 %! 0.83760 + 0.68977i ;
736 %! 0.39160 + 0.90378i ];
737
738
739
740 %!test
614 %! R = chol(A); 741 %! R = chol(A);
615 %! 742 %!
616 %! R1 = cholupdate(R,u); 743 %! R1 = cholupdate(R,u);
617 %! 744 %!
618 %! assert(norm(triu(R1)-R1,Inf) == 0) 745 %! assert(norm(triu(R1)-R1,Inf) == 0)
622 %! 749 %!
623 %! assert(norm(triu(R1)-R1,Inf) == 0) 750 %! assert(norm(triu(R1)-R1,Inf) == 0)
624 %! assert(norm(R1 - R,Inf) < 1e1*eps) 751 %! assert(norm(R1 - R,Inf) < 1e1*eps)
625 %! 752 %!
626 %!test 753 %!test
627 %! A = [ 0.5585528 + 0.0000000i -0.1662088 - 0.0315341i 0.0107873 + 0.0236411i -0.0276775 - 0.0186073i ; 754 %! R = chol(Ac);
628 %! -0.1662088 + 0.0315341i 0.6760061 + 0.0000000i 0.0011452 - 0.0475528i 0.0145967 + 0.0247641i ; 755 %!
629 %! 0.0107873 - 0.0236411i 0.0011452 + 0.0475528i 0.6263149 - 0.0000000i -0.1585837 - 0.0719763i ; 756 %! R1 = cholupdate(R,uc);
630 %! -0.0276775 + 0.0186073i 0.0145967 - 0.0247641i -0.1585837 + 0.0719763i 0.6034234 - 0.0000000i ]; 757 %!
631 %! 758 %! assert(norm(triu(R1)-R1,Inf) == 0)
632 %! u = [ 0.54267 + 0.91519i ; 759 %! assert(norm(R1'*R1 - R'*R - uc*uc',Inf) < 1e1*eps)
633 %! 0.99647 + 0.43141i ; 760 %!
634 %! 0.83760 + 0.68977i ; 761 %! R1 = cholupdate(R1,uc,"-");
635 %! 0.39160 + 0.90378i ];
636 %!
637 %! R = chol(A);
638 %!
639 %! R1 = cholupdate(R,u);
640 %!
641 %! assert(norm(triu(R1)-R1,Inf) == 0)
642 %! assert(norm(R1'*R1 - R'*R - u*u',Inf) < 1e1*eps)
643 %!
644 %! R1 = cholupdate(R1,u,"-");
645 %! 762 %!
646 %! assert(norm(triu(R1)-R1,Inf) == 0) 763 %! assert(norm(triu(R1)-R1,Inf) == 0)
647 %! assert(norm(R1 - R,Inf) < 1e1*eps) 764 %! assert(norm(R1 - R,Inf) < 1e1*eps)
765
766 %!test
767 %! R = chol(single(A));
768 %!
769 %! R1 = cholupdate(R,single(u));
770 %!
771 %! assert(norm(triu(R1)-R1,Inf) == 0)
772 %! assert(norm(R1'*R1 - R'*R - single(u*u'),Inf) < 1e1*eps('single'))
773 %!
774 %! R1 = cholupdate(R1,single(u),"-");
775 %!
776 %! assert(norm(triu(R1)-R1,Inf) == 0)
777 %! assert(norm(R1 - R,Inf) < 1e1*eps('single'))
778 %!
779 %!test
780 %! R = chol(single(Ac));
781 %!
782 %! R1 = cholupdate(R,single(uc));
783 %!
784 %! assert(norm(triu(R1)-R1,Inf) == 0)
785 %! assert(norm(R1'*R1 - R'*R - single(uc*uc'),Inf) < 1e1*eps('single'))
786 %!
787 %! R1 = cholupdate(R1,single(uc),"-");
788 %!
789 %! assert(norm(triu(R1)-R1,Inf) == 0)
790 %! assert(norm(R1 - R,Inf) < 1e1*eps('single'))
648 */ 791 */
649 792
650 DEFUN_DLD (cholinsert, args, nargout, 793 DEFUN_DLD (cholinsert, args, nargout,
651 "-*- texinfo -*-\n\ 794 "-*- texinfo -*-\n\
652 @deftypefn {Loadable Function} {[@var{R1}, @var{info}] =} cholinsert (@var{R}, @var{j}, @var{u})\n\ 795 @deftypefn {Loadable Function} {[@var{R1}, @var{info}] =} cholinsert (@var{R}, @var{j}, @var{u})\n\
688 831
689 if (argr.columns () == n && argu.rows () == n+1 && argu.columns () == 1) 832 if (argr.columns () == n && argu.rows () == n+1 && argu.columns () == 1)
690 { 833 {
691 if (j > 0 && j <= n+1) 834 if (j > 0 && j <= n+1)
692 { 835 {
693 if (argr.is_real_matrix () && argu.is_real_matrix ()) 836 if (argr.is_single_type () || argu.is_single_type ())
694 { 837 {
695 // real case 838 if (argr.is_real_matrix () && argu.is_real_matrix ())
696 Matrix R = argr.matrix_value (); 839 {
697 Matrix u = argu.matrix_value (); 840 // real case
698 841 FloatMatrix R = argr.float_matrix_value ();
699 CHOL fact; 842 FloatMatrix u = argu.float_matrix_value ();
700 fact.set (R); 843
701 int err = fact.insert_sym (u, j-1); 844 FloatCHOL fact;
702 845 fact.set (R);
703 if (nargout > 1) 846 int err = fact.insert_sym (u, j-1);
704 retval(1) = err; 847
705 else if (err) 848 if (nargout > 1)
706 error ("cholinsert: insertion violates positiveness"); 849 retval(1) = err;
707 850 else if (err)
708 retval(0) = fact.chol_matrix (); 851 error ("cholinsert: insertion violates positiveness");
709 } 852
710 else 853 retval(0) = fact.chol_matrix ();
711 { 854 }
712 // complex case 855 else
713 ComplexMatrix R = argr.complex_matrix_value (); 856 {
714 ComplexMatrix u = argu.complex_matrix_value (); 857 // complex case
715 858 FloatComplexMatrix R = argr.float_complex_matrix_value ();
716 ComplexCHOL fact; 859 FloatComplexMatrix u = argu.float_complex_matrix_value ();
717 fact.set (R); 860
718 int err = fact.insert_sym (u, j-1); 861 FloatComplexCHOL fact;
719 862 fact.set (R);
720 if (nargout > 1) 863 int err = fact.insert_sym (u, j-1);
721 retval(1) = err; 864
722 else if (err) 865 if (nargout > 1)
723 error ("cholinsert: insertion violates positiveness"); 866 retval(1) = err;
724 867 else if (err)
725 retval(0) = fact.chol_matrix (); 868 error ("cholinsert: insertion violates positiveness");
726 } 869
870 retval(0) = fact.chol_matrix ();
871 }
872 }
873 else
874 {
875 if (argr.is_real_matrix () && argu.is_real_matrix ())
876 {
877 // real case
878 Matrix R = argr.matrix_value ();
879 Matrix u = argu.matrix_value ();
880
881 CHOL fact;
882 fact.set (R);
883 int err = fact.insert_sym (u, j-1);
884
885 if (nargout > 1)
886 retval(1) = err;
887 else if (err)
888 error ("cholinsert: insertion violates positiveness");
889
890 retval(0) = fact.chol_matrix ();
891 }
892 else
893 {
894 // complex case
895 ComplexMatrix R = argr.complex_matrix_value ();
896 ComplexMatrix u = argu.complex_matrix_value ();
897
898 ComplexCHOL fact;
899 fact.set (R);
900 int err = fact.insert_sym (u, j-1);
901
902 if (nargout > 1)
903 retval(1) = err;
904 else if (err)
905 error ("cholinsert: insertion violates positiveness");
906
907 retval(0) = fact.chol_matrix ();
908 }
909 }
727 } 910 }
728 else 911 else
729 error ("cholinsert: index out of range"); 912 error ("cholinsert: index out of range");
730 } 913 }
731 else 914 else
737 return retval; 920 return retval;
738 } 921 }
739 922
740 /* 923 /*
741 %!test 924 %!test
742 %! A = [ 0.436997 -0.131721 0.124120 -0.061673 ; 925 %! u2 = [ 0.35080 ;
743 %! -0.131721 0.738529 0.019851 -0.140295 ; 926 %! 0.63930 ;
744 %! 0.124120 0.019851 0.354879 -0.059472 ; 927 %! 3.31057 ;
745 %! -0.061673 -0.140295 -0.059472 0.600939 ]; 928 %! -0.13825 ;
746 %! 929 %! 0.45266 ];
747 %! u = [ 0.35080 ;
748 %! 0.63930 ;
749 %! 3.31057 ;
750 %! -0.13825 ;
751 %! 0.45266 ];
752 %! 930 %!
753 %! R = chol(A); 931 %! R = chol(A);
754 %! 932 %!
755 %! j = 3; p = [1:j-1, j+1:5]; 933 %! j = 3; p = [1:j-1, j+1:5];
756 %! R1 = cholinsert(R,j,u); A1 = R1'*R1; 934 %! R1 = cholinsert(R,j,u2); A1 = R1'*R1;
757 %! 935 %!
758 %! assert(norm(triu(R1)-R1,Inf) == 0) 936 %! assert(norm(triu(R1)-R1,Inf) == 0)
759 %! assert(norm(A1(p,p) - A,Inf) < 1e1*eps) 937 %! assert(norm(A1(p,p) - A,Inf) < 1e1*eps)
760 %! 938 %!
761 %!test 939 %!test
762 %! A = [ 0.5585528 + 0.0000000i -0.1662088 - 0.0315341i 0.0107873 + 0.0236411i -0.0276775 - 0.0186073i ; 940 %! u2 = [ 0.35080 + 0.04298i;
763 %! -0.1662088 + 0.0315341i 0.6760061 + 0.0000000i 0.0011452 - 0.0475528i 0.0145967 + 0.0247641i ; 941 %! 0.63930 + 0.23778i;
764 %! 0.0107873 - 0.0236411i 0.0011452 + 0.0475528i 0.6263149 - 0.0000000i -0.1585837 - 0.0719763i ; 942 %! 3.31057 + 0.00000i;
765 %! -0.0276775 + 0.0186073i 0.0145967 - 0.0247641i -0.1585837 + 0.0719763i 0.6034234 - 0.0000000i ]; 943 %! -0.13825 + 0.19879i;
766 %! 944 %! 0.45266 + 0.50020i];
767 %! u = [ 0.35080 + 0.04298i;
768 %! 0.63930 + 0.23778i;
769 %! 3.31057 + 0.00000i;
770 %! -0.13825 + 0.19879i;
771 %! 0.45266 + 0.50020i];
772 %! 945 %!
773 %! R = chol(A); 946 %! R = chol(Ac);
774 %! 947 %!
775 %! j = 3; p = [1:j-1, j+1:5]; 948 %! j = 3; p = [1:j-1, j+1:5];
776 %! R1 = cholinsert(R,j,u); A1 = R1'*R1; 949 %! R1 = cholinsert(R,j,u2); A1 = R1'*R1;
777 %! 950 %!
778 %! assert(norm(triu(R1)-R1,Inf) == 0) 951 %! assert(norm(triu(R1)-R1,Inf) == 0)
779 %! assert(norm(A1(p,p) - A,Inf) < 1e1*eps) 952 %! assert(norm(A1(p,p) - Ac,Inf) < 1e1*eps)
953 %!
954
955 %!test
956 %! u2 = single ([ 0.35080 ;
957 %! 0.63930 ;
958 %! 3.31057 ;
959 %! -0.13825 ;
960 %! 0.45266 ]);
961 %!
962 %! R = chol(single(A));
963 %!
964 %! j = 3; p = [1:j-1, j+1:5];
965 %! R1 = cholinsert(R,j,u2); A1 = R1'*R1;
966 %!
967 %! assert(norm(triu(R1)-R1,Inf) == 0)
968 %! assert(norm(A1(p,p) - A,Inf) < 1e1*eps('single'))
969 %!
970 %!test
971 %! u2 = single ([ 0.35080 + 0.04298i;
972 %! 0.63930 + 0.23778i;
973 %! 3.31057 + 0.00000i;
974 %! -0.13825 + 0.19879i;
975 %! 0.45266 + 0.50020i]);
976 %!
977 %! R = chol(single(Ac));
978 %!
979 %! j = 3; p = [1:j-1, j+1:5];
980 %! R1 = cholinsert(R,j,u2); A1 = R1'*R1;
981 %!
982 %! assert(norm(triu(R1)-R1,Inf) == 0)
983 %! assert(norm(A1(p,p) - single(Ac),Inf) < 1e1*eps('single'))
780 %! 984 %!
781 */ 985 */
782 986
783 DEFUN_DLD (choldelete, args, nargout, 987 DEFUN_DLD (choldelete, args, nargout,
784 "-*- texinfo -*-\n\ 988 "-*- texinfo -*-\n\
809 1013
810 if (argr.columns () == n) 1014 if (argr.columns () == n)
811 { 1015 {
812 if (j > 0 && j <= n) 1016 if (j > 0 && j <= n)
813 { 1017 {
814 if (argr.is_real_matrix ()) 1018 if (argr.is_single_type ())
815 { 1019 {
816 // real case 1020 if (argr.is_real_matrix ())
817 Matrix R = argr.matrix_value (); 1021 {
818 1022 // real case
819 CHOL fact; 1023 FloatMatrix R = argr.float_matrix_value ();
820 fact.set (R); 1024
821 fact.delete_sym (j-1); 1025 FloatCHOL fact;
822 1026 fact.set (R);
823 retval(0) = fact.chol_matrix (); 1027 fact.delete_sym (j-1);
824 } 1028
825 else 1029 retval(0) = fact.chol_matrix ();
826 { 1030 }
827 // complex case 1031 else
828 ComplexMatrix R = argr.complex_matrix_value (); 1032 {
829 1033 // complex case
830 ComplexCHOL fact; 1034 FloatComplexMatrix R = argr.float_complex_matrix_value ();
831 fact.set (R); 1035
832 fact.delete_sym (j-1); 1036 FloatComplexCHOL fact;
833 1037 fact.set (R);
834 retval(0) = fact.chol_matrix (); 1038 fact.delete_sym (j-1);
835 } 1039
1040 retval(0) = fact.chol_matrix ();
1041 }
1042 }
1043 else
1044 {
1045 if (argr.is_real_matrix ())
1046 {
1047 // real case
1048 Matrix R = argr.matrix_value ();
1049
1050 CHOL fact;
1051 fact.set (R);
1052 fact.delete_sym (j-1);
1053
1054 retval(0) = fact.chol_matrix ();
1055 }
1056 else
1057 {
1058 // complex case
1059 ComplexMatrix R = argr.complex_matrix_value ();
1060
1061 ComplexCHOL fact;
1062 fact.set (R);
1063 fact.delete_sym (j-1);
1064
1065 retval(0) = fact.chol_matrix ();
1066 }
1067 }
836 } 1068 }
837 else 1069 else
838 error ("choldelete: index out of range"); 1070 error ("choldelete: index out of range");
839 } 1071 }
840 else 1072 else
846 return retval; 1078 return retval;
847 } 1079 }
848 1080
849 /* 1081 /*
850 %!test 1082 %!test
851 %! A = [ 0.436997 -0.131721 0.124120 -0.061673 ;
852 %! -0.131721 0.738529 0.019851 -0.140295 ;
853 %! 0.124120 0.019851 0.354879 -0.059472 ;
854 %! -0.061673 -0.140295 -0.059472 0.600939 ];
855 %!
856 %! R = chol(A); 1083 %! R = chol(A);
857 %! 1084 %!
858 %! j = 3; p = [1:j-1,j+1:4]; 1085 %! j = 3; p = [1:j-1,j+1:4];
859 %! R1 = choldelete(R,j); 1086 %! R1 = choldelete(R,j);
860 %! 1087 %!
861 %! assert(norm(triu(R1)-R1,Inf) == 0) 1088 %! assert(norm(triu(R1)-R1,Inf) == 0)
862 %! assert(norm(R1'*R1 - A(p,p),Inf) < 1e1*eps) 1089 %! assert(norm(R1'*R1 - A(p,p),Inf) < 1e1*eps)
863 %! 1090 %!
864 %!test 1091 %!test
865 %! A = [ 0.5585528 + 0.0000000i -0.1662088 - 0.0315341i 0.0107873 + 0.0236411i -0.0276775 - 0.0186073i ; 1092 %! R = chol(Ac);
866 %! -0.1662088 + 0.0315341i 0.6760061 + 0.0000000i 0.0011452 - 0.0475528i 0.0145967 + 0.0247641i ;
867 %! 0.0107873 - 0.0236411i 0.0011452 + 0.0475528i 0.6263149 - 0.0000000i -0.1585837 - 0.0719763i ;
868 %! -0.0276775 + 0.0186073i 0.0145967 - 0.0247641i -0.1585837 + 0.0719763i 0.6034234 - 0.0000000i ];
869 %!
870 %! R = chol(A);
871 %! 1093 %!
872 %! j = 3; p = [1:j-1,j+1:4]; 1094 %! j = 3; p = [1:j-1,j+1:4];
873 %! R1 = choldelete(R,j); 1095 %! R1 = choldelete(R,j);
874 %! 1096 %!
875 %! assert(norm(triu(R1)-R1,Inf) == 0) 1097 %! assert(norm(triu(R1)-R1,Inf) == 0)
876 %! assert(norm(R1'*R1 - A(p,p),Inf) < 1e1*eps) 1098 %! assert(norm(R1'*R1 - Ac(p,p),Inf) < 1e1*eps)
1099
1100 %!test
1101 %! R = chol(single(A));
1102 %!
1103 %! j = 3; p = [1:j-1,j+1:4];
1104 %! R1 = choldelete(R,j);
1105 %!
1106 %! assert(norm(triu(R1)-R1,Inf) == 0)
1107 %! assert(norm(R1'*R1 - single(A(p,p)),Inf) < 1e1*eps('single'))
1108 %!
1109 %!test
1110 %! R = chol(single(Ac));
1111 %!
1112 %! j = 3; p = [1:j-1,j+1:4];
1113 %! R1 = choldelete(R,j);
1114 %!
1115 %! assert(norm(triu(R1)-R1,Inf) == 0)
1116 %! assert(norm(R1'*R1 - single(Ac(p,p)),Inf) < 1e1*eps('single'))
877 */ 1117 */
878 1118
879 DEFUN_DLD (cholshift, args, nargout, 1119 DEFUN_DLD (cholshift, args, nargout,
880 "-*- texinfo -*-\n\ 1120 "-*- texinfo -*-\n\
881 @deftypefn {Loadable Function} {@var{R1} =} cholshift (@var{R}, @var{i}, @var{j})\n\ 1121 @deftypefn {Loadable Function} {@var{R1} =} cholshift (@var{R}, @var{i}, @var{j})\n\
912 1152
913 if (argr.columns () == n) 1153 if (argr.columns () == n)
914 { 1154 {
915 if (j > 0 && j <= n+1 && i > 0 && i <= n+1) 1155 if (j > 0 && j <= n+1 && i > 0 && i <= n+1)
916 { 1156 {
917 if (argr.is_real_matrix ()) 1157
918 { 1158 if (argr.is_single_type () && argi.is_single_type () &&
919 // real case 1159 argj.is_single_type ())
920 Matrix R = argr.matrix_value (); 1160 {
921 1161 if (argr.is_real_matrix ())
922 CHOL fact; 1162 {
923 fact.set (R); 1163 // real case
924 fact.shift_sym (i-1, j-1); 1164 FloatMatrix R = argr.float_matrix_value ();
925 1165
926 retval(0) = fact.chol_matrix (); 1166 FloatCHOL fact;
927 } 1167 fact.set (R);
928 else 1168 fact.shift_sym (i-1, j-1);
929 { 1169
930 // complex case 1170 retval(0) = fact.chol_matrix ();
931 ComplexMatrix R = argr.complex_matrix_value (); 1171 }
932 1172 else
933 ComplexCHOL fact; 1173 {
934 fact.set (R); 1174 // complex case
935 fact.shift_sym (i-1, j-1); 1175 FloatComplexMatrix R = argr.float_complex_matrix_value ();
936 1176
937 retval(0) = fact.chol_matrix (); 1177 FloatComplexCHOL fact;
938 } 1178 fact.set (R);
1179 fact.shift_sym (i-1, j-1);
1180
1181 retval(0) = fact.chol_matrix ();
1182 }
1183 }
1184 else
1185 {
1186 if (argr.is_real_matrix ())
1187 {
1188 // real case
1189 Matrix R = argr.matrix_value ();
1190
1191 CHOL fact;
1192 fact.set (R);
1193 fact.shift_sym (i-1, j-1);
1194
1195 retval(0) = fact.chol_matrix ();
1196 }
1197 else
1198 {
1199 // complex case
1200 ComplexMatrix R = argr.complex_matrix_value ();
1201
1202 ComplexCHOL fact;
1203 fact.set (R);
1204 fact.shift_sym (i-1, j-1);
1205
1206 retval(0) = fact.chol_matrix ();
1207 }
1208 }
939 } 1209 }
940 else 1210 else
941 error ("cholshift: index out of range"); 1211 error ("cholshift: index out of range");
942 } 1212 }
943 else 1213 else
949 return retval; 1219 return retval;
950 } 1220 }
951 1221
952 /* 1222 /*
953 %!test 1223 %!test
954 %! A = [ 0.436997 -0.131721 0.124120 -0.061673 ;
955 %! -0.131721 0.738529 0.019851 -0.140295 ;
956 %! 0.124120 0.019851 0.354879 -0.059472 ;
957 %! -0.061673 -0.140295 -0.059472 0.600939 ];
958 %!
959 %! R = chol(A); 1224 %! R = chol(A);
960 %! 1225 %!
961 %! i = 1; j = 3; p = [1:i-1, shift(i:j,-1), j+1:4]; 1226 %! i = 1; j = 3; p = [1:i-1, shift(i:j,-1), j+1:4];
962 %! R1 = cholshift(R,i,j); 1227 %! R1 = cholshift(R,i,j);
963 %! 1228 %!
969 %! 1234 %!
970 %! assert(norm(triu(R1)-R1,Inf) == 0) 1235 %! assert(norm(triu(R1)-R1,Inf) == 0)
971 %! assert(norm(R1'*R1 - A(p,p),Inf) < 1e1*eps) 1236 %! assert(norm(R1'*R1 - A(p,p),Inf) < 1e1*eps)
972 %! 1237 %!
973 %!test 1238 %!test
974 %! A = [ 0.5585528 + 0.0000000i -0.1662088 - 0.0315341i 0.0107873 + 0.0236411i -0.0276775 - 0.0186073i ; 1239 %! R = chol(Ac);
975 %! -0.1662088 + 0.0315341i 0.6760061 + 0.0000000i 0.0011452 - 0.0475528i 0.0145967 + 0.0247641i ;
976 %! 0.0107873 - 0.0236411i 0.0011452 + 0.0475528i 0.6263149 - 0.0000000i -0.1585837 - 0.0719763i ;
977 %! -0.0276775 + 0.0186073i 0.0145967 - 0.0247641i -0.1585837 + 0.0719763i 0.6034234 - 0.0000000i ];
978 %!
979 %! R = chol(A);
980 %! 1240 %!
981 %! i = 1; j = 3; p = [1:i-1, shift(i:j,-1), j+1:4]; 1241 %! i = 1; j = 3; p = [1:i-1, shift(i:j,-1), j+1:4];
982 %! R1 = cholshift(R,i,j); 1242 %! R1 = cholshift(R,i,j);
983 %! 1243 %!
984 %! assert(norm(triu(R1)-R1,Inf) == 0) 1244 %! assert(norm(triu(R1)-R1,Inf) == 0)
985 %! assert(norm(R1'*R1 - A(p,p),Inf) < 1e1*eps) 1245 %! assert(norm(R1'*R1 - Ac(p,p),Inf) < 1e1*eps)
986 %! 1246 %!
987 %! j = 1; i = 3; p = [1:j-1, shift(j:i,+1), i+1:4]; 1247 %! j = 1; i = 3; p = [1:j-1, shift(j:i,+1), i+1:4];
988 %! R1 = cholshift(R,i,j); 1248 %! R1 = cholshift(R,i,j);
989 %! 1249 %!
990 %! assert(norm(triu(R1)-R1,Inf) == 0) 1250 %! assert(norm(triu(R1)-R1,Inf) == 0)
991 %! assert(norm(R1'*R1 - A(p,p),Inf) < 1e1*eps) 1251 %! assert(norm(R1'*R1 - Ac(p,p),Inf) < 1e1*eps)
1252
1253 %!test
1254 %! R = chol(single(A));
1255 %!
1256 %! i = 1; j = 3; p = [1:i-1, shift(i:j,-1), j+1:4];
1257 %! R1 = cholshift(R,i,j);
1258 %!
1259 %! assert(norm(triu(R1)-R1,Inf) == 0)
1260 %! assert(norm(R1'*R1 - single(A(p,p)),Inf) < 1e1*eps('single'))
1261 %!
1262 %! j = 1; i = 3; p = [1:j-1, shift(j:i,+1), i+1:4];
1263 %! R1 = cholshift(R,i,j);
1264 %!
1265 %! assert(norm(triu(R1)-R1,Inf) == 0)
1266 %! assert(norm(R1'*R1 - single(A(p,p)),Inf) < 1e1*eps('single'))
1267 %!
1268 %!test
1269 %! R = chol(single(Ac));
1270 %!
1271 %! i = 1; j = 3; p = [1:i-1, shift(i:j,-1), j+1:4];
1272 %! R1 = cholshift(R,i,j);
1273 %!
1274 %! assert(norm(triu(R1)-R1,Inf) == 0)
1275 %! assert(norm(R1'*R1 - single(Ac(p,p)),Inf) < 1e1*eps('single'))
1276 %!
1277 %! j = 1; i = 3; p = [1:j-1, shift(j:i,+1), i+1:4];
1278 %! R1 = cholshift(R,i,j);
1279 %!
1280 %! assert(norm(triu(R1)-R1,Inf) == 0)
1281 %! assert(norm(R1'*R1 - single(Ac(p,p)),Inf) < 1e1*eps('single'))
992 */ 1282 */
993 1283
994 /* 1284 /*
995 ;;; Local Variables: *** 1285 ;;; Local Variables: ***
996 ;;; mode: C++ *** 1286 ;;; mode: C++ ***