Mercurial > octave-nkf

comparison src/DLD-FUNCTIONS/chol.cc @ 8547:d66c9b6e506a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
imported patch qrupdate.diff
author Jaroslav Hajek <highegg@gmail.com>
date Tue, 20 Jan 2009 21:16:42 +0100
parents 81d6ab3ac93c
children a6edd5c23cb5
comparison
equal deleted inserted replaced
8546:3d8a914c580e 8547:d66c9b6e506a
1 /* 1 /*
2 2
3 Copyright (C) 1996, 1997, 1999, 2000, 2002, 2005, 2006, 2007 3 Copyright (C) 1996, 1997, 1999, 2000, 2002, 2005, 2006, 2007
4 John W. Eaton 4 John W. Eaton
5 Copyright (C) 2008, 2009 Jaroslav Hajek
6 Copyright (C) 2008, 2009 VZLU Prague
5 7
6 This file is part of Octave. 8 This file is part of Octave.
7 9
8 Octave is free software; you can redistribute it and/or modify it 10 Octave is free software; you can redistribute it and/or modify it
9 under the terms of the GNU General Public License as published by the 11 under the terms of the GNU General Public License as published by the
19 along with Octave; see the file COPYING. If not, see 21 along with Octave; see the file COPYING. If not, see
20 <http://www.gnu.org/licenses/>. 22 <http://www.gnu.org/licenses/>.
21 23
22 */ 24 */
23 25
24 // The cholupdate, cholinsert, choldelete and cholshift functions were
25 // written by Jaroslav Hajek <highegg@gmail.com>, Copyright (C) 2008
26 // VZLU Prague, a.s., Czech Republic.
27 26
28 #ifdef HAVE_CONFIG_H 27 #ifdef HAVE_CONFIG_H
29 #include <config.h> 28 #include <config.h>
30 #endif 29 #endif
31 30
575 print_usage (); 574 print_usage ();
576 575
577 return retval; 576 return retval;
578 } 577 }
579 578
579 #ifdef HAVE_QRUPDATE
580
580 DEFUN_DLD (cholupdate, args, nargout, 581 DEFUN_DLD (cholupdate, args, nargout,
581 "-*- texinfo -*-\n\ 582 "-*- texinfo -*-\n\
582 @deftypefn {Loadable Function} {[@var{R1}, @var{info}] =} cholupdate (@var{R}, @var{u}, @var{op})\n\ 583 @deftypefn {Loadable Function} {[@var{R1}, @var{info}] =} cholupdate (@var{R}, @var{u}, @var{op})\n\
583 Update or downdate a Cholesky factorization. Given an upper triangular\n\ 584 Update or downdate a Cholesky factorization. Given an upper triangular\n\
584 matrix @var{R} and a column vector @var{u}, attempt to determine another\n\ 585 matrix @var{R} and a column vector @var{u}, attempt to determine another\n\
626 bool down = op == "-"; 627 bool down = op == "-";
627 628
628 if (down || op == "+") 629 if (down || op == "+")
629 if (argr.columns () == n && argu.rows () == n && argu.columns () == 1) 630 if (argr.columns () == n && argu.rows () == n && argu.columns () == 1)
630 { 631 {
632 int err = 0;
631 if (argr.is_single_type () || argu.is_single_type ()) 633 if (argr.is_single_type () || argu.is_single_type ())
632 { 634 {
633 if (argr.is_real_matrix () && argu.is_real_matrix ()) 635 if (argr.is_real_type () && argu.is_real_type ())
634 { 636 {
635 // real case 637 // real case
636 FloatMatrix R = argr.float_matrix_value (); 638 FloatMatrix R = argr.float_matrix_value ();
637 FloatMatrix u = argu.float_matrix_value (); 639 FloatColumnVector u = argu.float_column_vector_value ();
638 640
639 FloatCHOL fact; 641 FloatCHOL fact;
640 fact.set (R); 642 fact.set (R);
641 int err = 0;
642 643
643 if (down) 644 if (down)
644 err = fact.downdate (u); 645 err = fact.downdate (u);
645 else 646 else
646 fact.update (u); 647 fact.update (u);
647
648 if (nargout > 1)
649 retval(1) = err;
650 else if (err)
651 error ("cholupdate: downdate violates positiveness");
652 648
653 retval(0) = fact.chol_matrix (); 649 retval(0) = fact.chol_matrix ();
654 } 650 }
655 else 651 else
656 { 652 {
657 // complex case 653 // complex case
658 FloatComplexMatrix R = argr.float_complex_matrix_value (); 654 FloatComplexMatrix R = argr.float_complex_matrix_value ();
659 FloatComplexMatrix u = argu.float_complex_matrix_value (); 655 FloatComplexColumnVector u = argu.float_complex_column_vector_value ();
660 656
661 FloatComplexCHOL fact; 657 FloatComplexCHOL fact;
662 fact.set (R); 658 fact.set (R);
663 int err = 0;
664 659
665 if (down) 660 if (down)
666 err = fact.downdate (u); 661 err = fact.downdate (u);
667 else 662 else
668 fact.update (u); 663 fact.update (u);
669
670 if (nargout > 1)
671 retval(1) = err;
672 else if (err)
673 error ("cholupdate: downdate violates positiveness");
674 664
675 retval(0) = fact.chol_matrix (); 665 retval(0) = fact.chol_matrix ();
676 } 666 }
677 } 667 }
678 else 668 else
679 { 669 {
680 if (argr.is_real_matrix () && argu.is_real_matrix ()) 670 if (argr.is_real_type () && argu.is_real_type ())
681 { 671 {
682 // real case 672 // real case
683 Matrix R = argr.matrix_value (); 673 Matrix R = argr.matrix_value ();
684 Matrix u = argu.matrix_value (); 674 ColumnVector u = argu.column_vector_value ();
685 675
686 CHOL fact; 676 CHOL fact;
687 fact.set (R); 677 fact.set (R);
688 int err = 0;
689 678
690 if (down) 679 if (down)
691 err = fact.downdate (u); 680 err = fact.downdate (u);
692 else 681 else
693 fact.update (u); 682 fact.update (u);
694
695 if (nargout > 1)
696 retval(1) = err;
697 else if (err)
698 error ("cholupdate: downdate violates positiveness");
699 683
700 retval(0) = fact.chol_matrix (); 684 retval(0) = fact.chol_matrix ();
701 } 685 }
702 else 686 else
703 { 687 {
704 // complex case 688 // complex case
705 ComplexMatrix R = argr.complex_matrix_value (); 689 ComplexMatrix R = argr.complex_matrix_value ();
706 ComplexMatrix u = argu.complex_matrix_value (); 690 ComplexColumnVector u = argu.complex_column_vector_value ();
707 691
708 ComplexCHOL fact; 692 ComplexCHOL fact;
709 fact.set (R); 693 fact.set (R);
710 int err = 0;
711 694
712 if (down) 695 if (down)
713 err = fact.downdate (u); 696 err = fact.downdate (u);
714 else 697 else
715 fact.update (u); 698 fact.update (u);
716 699
717 if (nargout > 1)
718 retval(1) = err;
719 else if (err)
720 error ("cholupdate: downdate violates positiveness");
721
722 retval(0) = fact.chol_matrix (); 700 retval(0) = fact.chol_matrix ();
723 } 701 }
724 } 702 }
703
704 if (nargout > 1)
705 retval(1) = err;
706 else if (err == 1)
707 error ("cholupdate: downdate violates positiveness");
708 else if (err == 2)
709 error ("cholupdate: singular matrix");
725 } 710 }
726 else 711 else
727 error ("cholupdate: dimension mismatch"); 712 error ("cholupdate: dimension mismatch");
728 else 713 else
729 error ("cholupdate: op must be \"+\" or \"-\""); 714 error ("cholupdate: op must be \"+\" or \"-\"");
851 836
852 if (argr.columns () == n && argu.rows () == n+1 && argu.columns () == 1) 837 if (argr.columns () == n && argu.rows () == n+1 && argu.columns () == 1)
853 { 838 {
854 if (j > 0 && j <= n+1) 839 if (j > 0 && j <= n+1)
855 { 840 {
841 int err = 0;
856 if (argr.is_single_type () || argu.is_single_type ()) 842 if (argr.is_single_type () || argu.is_single_type ())
857 { 843 {
858 if (argr.is_real_matrix () && argu.is_real_matrix ()) 844 if (argr.is_real_type () && argu.is_real_type ())
859 { 845 {
860 // real case 846 // real case
861 FloatMatrix R = argr.float_matrix_value (); 847 FloatMatrix R = argr.float_matrix_value ();
862 FloatMatrix u = argu.float_matrix_value (); 848 FloatColumnVector u = argu.float_column_vector_value ();
863 849
864 FloatCHOL fact; 850 FloatCHOL fact;
865 fact.set (R); 851 fact.set (R);
866 int err = fact.insert_sym (u, j-1); 852 err = fact.insert_sym (u, j-1);
867
868 if (nargout > 1)
869 retval(1) = err;
870 else if (err)
871 error ("cholinsert: insertion violates positiveness");
872 853
873 retval(0) = fact.chol_matrix (); 854 retval(0) = fact.chol_matrix ();
874 } 855 }
875 else 856 else
876 { 857 {
877 // complex case 858 // complex case
878 FloatComplexMatrix R = argr.float_complex_matrix_value (); 859 FloatComplexMatrix R = argr.float_complex_matrix_value ();
879 FloatComplexMatrix u = argu.float_complex_matrix_value (); 860 FloatComplexColumnVector u = argu.float_complex_column_vector_value ();
880 861
881 FloatComplexCHOL fact; 862 FloatComplexCHOL fact;
882 fact.set (R); 863 fact.set (R);
883 int err = fact.insert_sym (u, j-1); 864 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 865
890 retval(0) = fact.chol_matrix (); 866 retval(0) = fact.chol_matrix ();
891 } 867 }
892 } 868 }
893 else 869 else
894 { 870 {
895 if (argr.is_real_matrix () && argu.is_real_matrix ()) 871 if (argr.is_real_type () && argu.is_real_type ())
896 { 872 {
897 // real case 873 // real case
898 Matrix R = argr.matrix_value (); 874 Matrix R = argr.matrix_value ();
899 Matrix u = argu.matrix_value (); 875 ColumnVector u = argu.column_vector_value ();
900 876
901 CHOL fact; 877 CHOL fact;
902 fact.set (R); 878 fact.set (R);
903 int err = fact.insert_sym (u, j-1); 879 err = fact.insert_sym (u, j-1);
904
905 if (nargout > 1)
906 retval(1) = err;
907 else if (err)
908 error ("cholinsert: insertion violates positiveness");
909 880
910 retval(0) = fact.chol_matrix (); 881 retval(0) = fact.chol_matrix ();
911 } 882 }
912 else 883 else
913 { 884 {
914 // complex case 885 // complex case
915 ComplexMatrix R = argr.complex_matrix_value (); 886 ComplexMatrix R = argr.complex_matrix_value ();
916 ComplexMatrix u = argu.complex_matrix_value (); 887 ComplexColumnVector u = argu.complex_column_vector_value ();
917 888
918 ComplexCHOL fact; 889 ComplexCHOL fact;
919 fact.set (R); 890 fact.set (R);
920 int err = fact.insert_sym (u, j-1); 891 err = fact.insert_sym (u, j-1);
921
922 if (nargout > 1)
923 retval(1) = err;
924 else if (err)
925 error ("cholinsert: insertion violates positiveness");
926 892
927 retval(0) = fact.chol_matrix (); 893 retval(0) = fact.chol_matrix ();
928 } 894 }
929 } 895 }
896
897 if (nargout > 1)
898 retval(1) = err;
899 else if (err == 1)
900 error ("cholinsert: insertion violates positiveness");
901 else if (err == 2)
902 error ("cholinsert: singular matrix");
903 else if (err == 3)
904 error ("cholinsert: diagonal element must be real");
930 } 905 }
931 else 906 else
932 error ("cholinsert: index out of range"); 907 error ("cholinsert: index out of range");
933 } 908 }
934 else 909 else
1035 { 1010 {
1036 if (j > 0 && j <= n) 1011 if (j > 0 && j <= n)
1037 { 1012 {
1038 if (argr.is_single_type ()) 1013 if (argr.is_single_type ())
1039 { 1014 {
1040 if (argr.is_real_matrix ()) 1015 if (argr.is_real_type ())
1041 { 1016 {
1042 // real case 1017 // real case
1043 FloatMatrix R = argr.float_matrix_value (); 1018 FloatMatrix R = argr.float_matrix_value ();
1044 1019
1045 FloatCHOL fact; 1020 FloatCHOL fact;
1060 retval(0) = fact.chol_matrix (); 1035 retval(0) = fact.chol_matrix ();
1061 } 1036 }
1062 } 1037 }
1063 else 1038 else
1064 { 1039 {
1065 if (argr.is_real_matrix ()) 1040 if (argr.is_real_type ())
1066 { 1041 {
1067 // real case 1042 // real case
1068 Matrix R = argr.matrix_value (); 1043 Matrix R = argr.matrix_value ();
1069 1044
1070 CHOL fact; 1045 CHOL fact;
1176 { 1151 {
1177 1152
1178 if (argr.is_single_type () && argi.is_single_type () && 1153 if (argr.is_single_type () && argi.is_single_type () &&
1179 argj.is_single_type ()) 1154 argj.is_single_type ())
1180 { 1155 {
1181 if (argr.is_real_matrix ()) 1156 if (argr.is_real_type ())
1182 { 1157 {
1183 // real case 1158 // real case
1184 FloatMatrix R = argr.float_matrix_value (); 1159 FloatMatrix R = argr.float_matrix_value ();
1185 1160
1186 FloatCHOL fact; 1161 FloatCHOL fact;
1201 retval(0) = fact.chol_matrix (); 1176 retval(0) = fact.chol_matrix ();
1202 } 1177 }
1203 } 1178 }
1204 else 1179 else
1205 { 1180 {
1206 if (argr.is_real_matrix ()) 1181 if (argr.is_real_type ())
1207 { 1182 {
1208 // real case 1183 // real case
1209 Matrix R = argr.matrix_value (); 1184 Matrix R = argr.matrix_value ();
1210 1185
1211 CHOL fact; 1186 CHOL fact;
1299 %! 1274 %!
1300 %! assert(norm(triu(R1)-R1,Inf) == 0) 1275 %! assert(norm(triu(R1)-R1,Inf) == 0)
1301 %! assert(norm(R1'*R1 - single(Ac(p,p)),Inf) < 1e1*eps('single')) 1276 %! assert(norm(R1'*R1 - single(Ac(p,p)),Inf) < 1e1*eps('single'))
1302 */ 1277 */
1303 1278
1279 #endif
1280
1304 /* 1281 /*
1305 ;;; Local Variables: *** 1282 ;;; Local Variables: ***
1306 ;;; mode: C++ *** 1283 ;;; mode: C++ ***
1307 ;;; End: *** 1284 ;;; End: ***
1308 */ 1285 */