Mercurial > hg > octave-lyh

comparison src/DLD-FUNCTIONS/lu.cc @ 11586:12df7854fa7c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
strip trailing whitespace from source files
author John W. Eaton <jwe@octave.org>
date Thu, 20 Jan 2011 17:24:59 -0500
parents 7d6d8c1e471f
children 990762e784fe
comparison
equal deleted inserted replaced
11585:1473d0cf86d2 11586:12df7854fa7c
142 octave_value_list retval; 142 octave_value_list retval;
143 int nargin = args.length (); 143 int nargin = args.length ();
144 bool issparse = (nargin > 0 && args(0).is_sparse_type ()); 144 bool issparse = (nargin > 0 && args(0).is_sparse_type ());
145 bool scale = (nargout == 5); 145 bool scale = (nargout == 5);
146 146
147 if (nargin < 1 || (issparse && (nargin > 3 || nargout > 5)) 147 if (nargin < 1 || (issparse && (nargin > 3 || nargout > 5))
148 || (!issparse && (nargin > 2 || nargout > 3))) 148 || (!issparse && (nargin > 2 || nargout > 3)))
149 { 149 {
150 print_usage (); 150 print_usage ();
151 return retval; 151 return retval;
152 } 152 }
229 retval (0) = fact.Y (); 229 retval (0) = fact.Y ();
230 else 230 else
231 { 231 {
232 PermMatrix P = fact.Pr_mat (); 232 PermMatrix P = fact.Pr_mat ();
233 SparseMatrix L = P.transpose () * fact.L (); 233 SparseMatrix L = P.transpose () * fact.L ();
234 retval(1) = octave_value (fact.U (), 234 retval(1) = octave_value (fact.U (),
235 MatrixType (MatrixType::Upper)); 235 MatrixType (MatrixType::Upper));
236 236
237 retval(0) = octave_value (L, 237 retval(0) = octave_value (L,
238 MatrixType (MatrixType::Permuted_Lower, 238 MatrixType (MatrixType::Permuted_Lower,
239 nr, fact.row_perm ())); 239 nr, fact.row_perm ()));
240 } 240 }
241 } 241 }
242 break; 242 break;
243 243
248 if (vecout) 248 if (vecout)
249 retval (2) = fact.Pr_vec (); 249 retval (2) = fact.Pr_vec ();
250 else 250 else
251 retval(2) = fact.Pr_mat (); 251 retval(2) = fact.Pr_mat ();
252 252
253 retval(1) = octave_value (fact.U (), 253 retval(1) = octave_value (fact.U (),
254 MatrixType (MatrixType::Upper)); 254 MatrixType (MatrixType::Upper));
255 retval(0) = octave_value (fact.L (), 255 retval(0) = octave_value (fact.L (),
256 MatrixType (MatrixType::Lower)); 256 MatrixType (MatrixType::Lower));
257 } 257 }
258 break; 258 break;
259 259
260 case 4: 260 case 4:
273 else 273 else
274 { 274 {
275 retval(3) = fact.Pc_mat (); 275 retval(3) = fact.Pc_mat ();
276 retval(2) = fact.Pr_mat (); 276 retval(2) = fact.Pr_mat ();
277 } 277 }
278 retval(1) = octave_value (fact.U (), 278 retval(1) = octave_value (fact.U (),
279 MatrixType (MatrixType::Upper)); 279 MatrixType (MatrixType::Upper));
280 retval(0) = octave_value (fact.L (), 280 retval(0) = octave_value (fact.L (),
281 MatrixType (MatrixType::Lower)); 281 MatrixType (MatrixType::Lower));
282 } 282 }
283 break; 283 break;
284 } 284 }
285 } 285 }
299 retval (0) = fact.Y (); 299 retval (0) = fact.Y ();
300 else 300 else
301 { 301 {
302 PermMatrix P = fact.Pr_mat (); 302 PermMatrix P = fact.Pr_mat ();
303 SparseComplexMatrix L = P.transpose () * fact.L (); 303 SparseComplexMatrix L = P.transpose () * fact.L ();
304 retval(1) = octave_value (fact.U (), 304 retval(1) = octave_value (fact.U (),
305 MatrixType (MatrixType::Upper)); 305 MatrixType (MatrixType::Upper));
306 306
307 retval(0) = octave_value (L, 307 retval(0) = octave_value (L,
308 MatrixType (MatrixType::Permuted_Lower, 308 MatrixType (MatrixType::Permuted_Lower,
309 nr, fact.row_perm ())); 309 nr, fact.row_perm ()));
310 } 310 }
311 } 311 }
312 break; 312 break;
313 313
318 if (vecout) 318 if (vecout)
319 retval (2) = fact.Pr_vec (); 319 retval (2) = fact.Pr_vec ();
320 else 320 else
321 retval(2) = fact.Pr_mat (); 321 retval(2) = fact.Pr_mat ();
322 322
323 retval(1) = octave_value (fact.U (), 323 retval(1) = octave_value (fact.U (),
324 MatrixType (MatrixType::Upper)); 324 MatrixType (MatrixType::Upper));
325 retval(0) = octave_value (fact.L (), 325 retval(0) = octave_value (fact.L (),
326 MatrixType (MatrixType::Lower)); 326 MatrixType (MatrixType::Lower));
327 } 327 }
328 break; 328 break;
329 329
330 case 4: 330 case 4:
343 else 343 else
344 { 344 {
345 retval(3) = fact.Pc_mat (); 345 retval(3) = fact.Pc_mat ();
346 retval(2) = fact.Pr_mat (); 346 retval(2) = fact.Pr_mat ();
347 } 347 }
348 retval(1) = octave_value (fact.U (), 348 retval(1) = octave_value (fact.U (),
349 MatrixType (MatrixType::Upper)); 349 MatrixType (MatrixType::Upper));
350 retval(0) = octave_value (fact.L (), 350 retval(0) = octave_value (fact.L (),
351 MatrixType (MatrixType::Lower)); 351 MatrixType (MatrixType::Lower));
352 } 352 }
353 break; 353 break;
354 } 354 }
355 } 355 }
655 octave_value argu = args(1); 655 octave_value argu = args(1);
656 octave_value argp = pivoted ? args(2) : octave_value (); 656 octave_value argp = pivoted ? args(2) : octave_value ();
657 octave_value argx = args(2 + pivoted); 657 octave_value argx = args(2 + pivoted);
658 octave_value argy = args(3 + pivoted); 658 octave_value argy = args(3 + pivoted);
659 659
660 if (argl.is_numeric_type () && argu.is_numeric_type () 660 if (argl.is_numeric_type () && argu.is_numeric_type ()
661 && argx.is_numeric_type () && argy.is_numeric_type () 661 && argx.is_numeric_type () && argy.is_numeric_type ()
662 && (! pivoted || argp.is_perm_matrix ())) 662 && (! pivoted || argp.is_perm_matrix ()))
663 { 663 {
664 if (check_lu_dims (argl, argu, argp)) 664 if (check_lu_dims (argl, argu, argp))
665 { 665 {
666 PermMatrix P = (pivoted 666 PermMatrix P = (pivoted
667 ? argp.perm_matrix_value () 667 ? argp.perm_matrix_value ()
668 : PermMatrix::eye (argl.rows ())); 668 : PermMatrix::eye (argl.rows ()));
669 669
670 if (argl.is_real_type () 670 if (argl.is_real_type ()
671 && argu.is_real_type () 671 && argu.is_real_type ()
672 && argx.is_real_type () 672 && argx.is_real_type ()
673 && argy.is_real_type ()) 673 && argy.is_real_type ())
674 { 674 {
675 // all real case 675 // all real case
676 if (argl.is_single_type () 676 if (argl.is_single_type ()
677 || argu.is_single_type () 677 || argu.is_single_type ()
678 || argx.is_single_type () 678 || argx.is_single_type ()
679 || argy.is_single_type ()) 679 || argy.is_single_type ())
680 { 680 {
681 FloatMatrix L = argl.float_matrix_value (); 681 FloatMatrix L = argl.float_matrix_value ();
682 FloatMatrix U = argu.float_matrix_value (); 682 FloatMatrix U = argu.float_matrix_value ();
683 FloatMatrix x = argx.float_matrix_value (); 683 FloatMatrix x = argx.float_matrix_value ();
714 } 714 }
715 } 715 }
716 else 716 else
717 { 717 {
718 // complex case 718 // complex case
719 if (argl.is_single_type () 719 if (argl.is_single_type ()
720 || argu.is_single_type () 720 || argu.is_single_type ()
721 || argx.is_single_type () 721 || argx.is_single_type ()
722 || argy.is_single_type ()) 722 || argy.is_single_type ())
723 { 723 {
724 FloatComplexMatrix L = argl.float_complex_matrix_value (); 724 FloatComplexMatrix L = argl.float_complex_matrix_value ();
725 FloatComplexMatrix U = argu.float_complex_matrix_value (); 725 FloatComplexMatrix U = argu.float_complex_matrix_value ();
726 FloatComplexMatrix x = argx.float_complex_matrix_value (); 726 FloatComplexMatrix x = argx.float_complex_matrix_value ();
729 FloatComplexLU fact (L, U, P); 729 FloatComplexLU fact (L, U, P);
730 if (pivoted) 730 if (pivoted)
731 fact.update_piv (x, y); 731 fact.update_piv (x, y);
732 else 732 else
733 fact.update (x, y); 733 fact.update (x, y);
734 734
735 if (pivoted) 735 if (pivoted)
736 retval(2) = fact.P (); 736 retval(2) = fact.P ();
737 retval(1) = get_lu_u (fact); 737 retval(1) = get_lu_u (fact);
738 retval(0) = get_lu_l (fact); 738 retval(0) = get_lu_l (fact);
739 } 739 }
747 ComplexLU fact (L, U, P); 747 ComplexLU fact (L, U, P);
748 if (pivoted) 748 if (pivoted)
749 fact.update_piv (x, y); 749 fact.update_piv (x, y);
750 else 750 else
751 fact.update (x, y); 751 fact.update (x, y);
752 752
753 if (pivoted) 753 if (pivoted)
754 retval(2) = fact.P (); 754 retval(2) = fact.P ();
755 retval(1) = get_lu_u (fact); 755 retval(1) = get_lu_u (fact);
756 retval(0) = get_lu_l (fact); 756 retval(0) = get_lu_l (fact);
757 } 757 }
772 %! 0.594638 0.425302 0.603537; 772 %! 0.594638 0.425302 0.603537;
773 %! 0.383594 0.291238 0.085574; 773 %! 0.383594 0.291238 0.085574;
774 %! 0.265712 0.268003 0.238409; 774 %! 0.265712 0.268003 0.238409;
775 %! 0.669966 0.743851 0.445057 ]; 775 %! 0.669966 0.743851 0.445057 ];
776 %! 776 %!
777 %! u = [0.85082; 777 %! u = [0.85082;
778 %! 0.76426; 778 %! 0.76426;
779 %! 0.42883; 779 %! 0.42883;
780 %! 0.53010; 780 %! 0.53010;
781 %! 0.80683 ]; 781 %! 0.80683 ];
782 %! 782 %!
783 %! v = [0.98810; 783 %! v = [0.98810;
784 %! 0.24295; 784 %! 0.24295;
785 %! 0.43167 ]; 785 %! 0.43167 ];
805 %! [L,U,P] = lu(A); 805 %! [L,U,P] = lu(A);
806 %! [L,U] = luupdate(L,U,P*u,v); 806 %! [L,U] = luupdate(L,U,P*u,v);
807 %! assert(norm(vec(tril(L)-L),Inf) == 0) 807 %! assert(norm(vec(tril(L)-L),Inf) == 0)
808 %! assert(norm(vec(triu(U)-U),Inf) == 0) 808 %! assert(norm(vec(triu(U)-U),Inf) == 0)
809 %! assert(norm(vec(P'*L*U - A - u*v.'),Inf) < norm(A)*1e1*eps) 809 %! assert(norm(vec(P'*L*U - A - u*v.'),Inf) < norm(A)*1e1*eps)
810 %! 810 %!
811 %!testif HAVE_QRUPDATE_LUU 811 %!testif HAVE_QRUPDATE_LUU
812 %! [L,U,P] = lu(Ac); 812 %! [L,U,P] = lu(Ac);
813 %! [L,U] = luupdate(L,U,P*uc,vc); 813 %! [L,U] = luupdate(L,U,P*uc,vc);
814 %! assert(norm(vec(tril(L)-L),Inf) == 0) 814 %! assert(norm(vec(tril(L)-L),Inf) == 0)
815 %! assert(norm(vec(triu(U)-U),Inf) == 0) 815 %! assert(norm(vec(triu(U)-U),Inf) == 0)
819 %! [L,U,P] = lu(single(A)); 819 %! [L,U,P] = lu(single(A));
820 %! [L,U] = luupdate(L,U,P*single(u),single(v)); 820 %! [L,U] = luupdate(L,U,P*single(u),single(v));
821 %! assert(norm(vec(tril(L)-L),Inf) == 0) 821 %! assert(norm(vec(tril(L)-L),Inf) == 0)
822 %! assert(norm(vec(triu(U)-U),Inf) == 0) 822 %! assert(norm(vec(triu(U)-U),Inf) == 0)
823 %! assert(norm(vec(P'*L*U - single(A) - single(u)*single(v).'),Inf) < norm(single(A))*1e1*eps('single')) 823 %! assert(norm(vec(P'*L*U - single(A) - single(u)*single(v).'),Inf) < norm(single(A))*1e1*eps('single'))
824 %! 824 %!
825 %!testif HAVE_QRUPDATE_LUU 825 %!testif HAVE_QRUPDATE_LUU
826 %! [L,U,P] = lu(single(Ac)); 826 %! [L,U,P] = lu(single(Ac));
827 %! [L,U] = luupdate(L,U,P*single(uc),single(vc)); 827 %! [L,U] = luupdate(L,U,P*single(uc),single(vc));
828 %! assert(norm(vec(tril(L)-L),Inf) == 0) 828 %! assert(norm(vec(tril(L)-L),Inf) == 0)
829 %! assert(norm(vec(triu(U)-U),Inf) == 0) 829 %! assert(norm(vec(triu(U)-U),Inf) == 0)
833 %! [L,U,P] = lu(A); 833 %! [L,U,P] = lu(A);
834 %! [L,U,P] = luupdate(L,U,P,u,v); 834 %! [L,U,P] = luupdate(L,U,P,u,v);
835 %! assert(norm(vec(tril(L)-L),Inf) == 0) 835 %! assert(norm(vec(tril(L)-L),Inf) == 0)
836 %! assert(norm(vec(triu(U)-U),Inf) == 0) 836 %! assert(norm(vec(triu(U)-U),Inf) == 0)
837 %! assert(norm(vec(P'*L*U - A - u*v.'),Inf) < norm(A)*1e1*eps) 837 %! assert(norm(vec(P'*L*U - A - u*v.'),Inf) < norm(A)*1e1*eps)
838 %! 838 %!
839 %!testif HAVE_QRUPDATE_LUU 839 %!testif HAVE_QRUPDATE_LUU
840 %! [L,U,P] = lu(Ac); 840 %! [L,U,P] = lu(Ac);
841 %! [L,U,P] = luupdate(L,U,P,uc,vc); 841 %! [L,U,P] = luupdate(L,U,P,uc,vc);
842 %! assert(norm(vec(tril(L)-L),Inf) == 0) 842 %! assert(norm(vec(tril(L)-L),Inf) == 0)
843 %! assert(norm(vec(triu(U)-U),Inf) == 0) 843 %! assert(norm(vec(triu(U)-U),Inf) == 0)
847 %! [L,U,P] = lu(single(A)); 847 %! [L,U,P] = lu(single(A));
848 %! [L,U,P] = luupdate(L,U,P,single(u),single(v)); 848 %! [L,U,P] = luupdate(L,U,P,single(u),single(v));
849 %! assert(norm(vec(tril(L)-L),Inf) == 0) 849 %! assert(norm(vec(tril(L)-L),Inf) == 0)
850 %! assert(norm(vec(triu(U)-U),Inf) == 0) 850 %! assert(norm(vec(triu(U)-U),Inf) == 0)
851 %! assert(norm(vec(P'*L*U - single(A) - single(u)*single(v).'),Inf) < norm(single(A))*1e1*eps('single')) 851 %! assert(norm(vec(P'*L*U - single(A) - single(u)*single(v).'),Inf) < norm(single(A))*1e1*eps('single'))
852 %! 852 %!
853 %!testif HAVE_QRUPDATE_LUU 853 %!testif HAVE_QRUPDATE_LUU
854 %! [L,U,P] = lu(single(Ac)); 854 %! [L,U,P] = lu(single(Ac));
855 %! [L,U,P] = luupdate(L,U,P,single(uc),single(vc)); 855 %! [L,U,P] = luupdate(L,U,P,single(uc),single(vc));
856 %! assert(norm(vec(tril(L)-L),Inf) == 0) 856 %! assert(norm(vec(tril(L)-L),Inf) == 0)
857 %! assert(norm(vec(triu(U)-U),Inf) == 0) 857 %! assert(norm(vec(triu(U)-U),Inf) == 0)