Mercurial > octave

comparison libinterp/corefcn/__ilu__.cc @ 20892:c07bee629973

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
2015 Code Sprint: use ovl ().
author Rik <rik@octave.org>
date Mon, 14 Dec 2015 11:28:48 -0800
parents 1142cf6abc0d
children 03e4ddd49396
comparison
equal deleted inserted replaced
20891:95c0d4c07c56 20892:c07bee629973
132 @deftypefnx {} {[@var{L}, @var{U}] =} __ilu0__ (@var{A}, @var{milu})\n\ 132 @deftypefnx {} {[@var{L}, @var{U}] =} __ilu0__ (@var{A}, @var{milu})\n\
133 @deftypefnx {} {[@var{L}, @var{U}, @var{P}] =} __ilu0__ (@var{A}, @dots{})\n\ 133 @deftypefnx {} {[@var{L}, @var{U}, @var{P}] =} __ilu0__ (@var{A}, @dots{})\n\
134 Undocumented internal function.\n\ 134 Undocumented internal function.\n\
135 @end deftypefn") 135 @end deftypefn")
136 { 136 {
137 octave_value_list retval;
138
139 int nargin = args.length (); 137 int nargin = args.length ();
140 std::string milu;
141 138
142 if (nargout > 2 || nargin < 1 || nargin > 2) 139 if (nargout > 2 || nargin < 1 || nargin > 2)
143 print_usage (); 140 print_usage ();
141
142 octave_value_list retval (2);
143
144 std::string milu;
144 145
145 // In ILU0 algorithm the zero-pattern of the input matrix is preserved so 146 // In ILU0 algorithm the zero-pattern of the input matrix is preserved so
146 // it's structure does not change during the algorithm. The same input 147 // it's structure does not change during the algorithm. The same input
147 // matrix is used to build the output matrix due to that fact. 148 // matrix is used to build the output matrix due to that fact.
148 octave_value_list param_list; 149 octave_value_list param_list;
465 DEFUN (__iluc__, args, nargout, 466 DEFUN (__iluc__, args, nargout,
466 "-*- texinfo -*-\n\ 467 "-*- texinfo -*-\n\
467 @deftypefn {} {[@var{L}, @var{U}] =} __iluc__ (@var{A})\n\ 468 @deftypefn {} {[@var{L}, @var{U}] =} __iluc__ (@var{A})\n\
468 @deftypefnx {} {[@var{L}, @var{U}] =} __iluc__ (@var{A}, @var{droptol})\n\ 469 @deftypefnx {} {[@var{L}, @var{U}] =} __iluc__ (@var{A}, @var{droptol})\n\
469 @deftypefnx {} {[@var{L}, @var{U}] =} __iluc__ (@var{A}, @var{droptol}, @var{milu})\n\ 470 @deftypefnx {} {[@var{L}, @var{U}] =} __iluc__ (@var{A}, @var{droptol}, @var{milu})\n\
470 @deftypefnx {} {[@var{L}, @var{U}, @var{P}] =} __iluc__ (@var{A}, @dots{})\n\
471 Undocumented internal function.\n\ 471 Undocumented internal function.\n\
472 @end deftypefn") 472 @end deftypefn")
473 { 473 {
474 octave_value_list retval;
475 int nargin = args.length (); 474 int nargin = args.length ();
475
476 if (nargout != 2 || nargin < 1 || nargin > 3)
477 print_usage ();
478
476 std::string milu = "off"; 479 std::string milu = "off";
477 double droptol = 0; 480 double droptol = 0;
478
479 if (nargout != 2 || nargin < 1 || nargin > 3)
480 print_usage ();
481 481
482 // Don't repeat input validation of arguments done in ilu.m 482 // Don't repeat input validation of arguments done in ilu.m
483 if (nargin >= 2) 483 if (nargin >= 2)
484 droptol = args(1).double_value (); 484 droptol = args(1).double_value ();
485 485
497 param_list(1) = "rows"; 497 param_list(1) = "rows";
498 rows_norm = feval ("norm", param_list)(0).vector_value (); 498 rows_norm = feval ("norm", param_list)(0).vector_value ();
499 param_list(1) = "cols"; 499 param_list(1) = "cols";
500 cols_norm = feval ("norm", param_list)(0).vector_value (); 500 cols_norm = feval ("norm", param_list)(0).vector_value ();
501 param_list.clear (); 501 param_list.clear ();
502 SparseMatrix U; 502 SparseMatrix U, L;
503 SparseMatrix L;
504 ilu_crout <SparseMatrix, double> (sm_l, sm_u, L, U, 503 ilu_crout <SparseMatrix, double> (sm_l, sm_u, L, U,
505 cols_norm.fortran_vec (), 504 cols_norm.fortran_vec (),
506 rows_norm.fortran_vec (), 505 rows_norm.fortran_vec (),
507 droptol, milu); 506 droptol, milu);
508 param_list.append (octave_value (L.cols ())); 507 param_list.append (octave_value (L.cols ()));
509 SparseMatrix eye = 508 SparseMatrix eye =
510 feval ("speye", param_list)(0).sparse_matrix_value (); 509 feval ("speye", param_list)(0).sparse_matrix_value ();
511 retval(1) = U; 510 return ovl (L + eye, U);
512 retval(0) = L + eye;
513 } 511 }
514 else 512 else
515 { 513 {
516 Array<Complex> cols_norm, rows_norm; 514 Array<Complex> cols_norm, rows_norm;
517 param_list.append (args(0).sparse_complex_matrix_value ()); 515 param_list.append (args(0).sparse_complex_matrix_value ());
523 param_list(1) = "rows"; 521 param_list(1) = "rows";
524 rows_norm = feval ("norm", param_list)(0).complex_vector_value (); 522 rows_norm = feval ("norm", param_list)(0).complex_vector_value ();
525 param_list(1) = "cols"; 523 param_list(1) = "cols";
526 cols_norm = feval ("norm", param_list)(0).complex_vector_value (); 524 cols_norm = feval ("norm", param_list)(0).complex_vector_value ();
527 param_list.clear (); 525 param_list.clear ();
528 SparseComplexMatrix U; 526 SparseComplexMatrix U, L;
529 SparseComplexMatrix L;
530 ilu_crout < SparseComplexMatrix, Complex > 527 ilu_crout < SparseComplexMatrix, Complex >
531 (sm_l, sm_u, L, U, cols_norm.fortran_vec () , 528 (sm_l, sm_u, L, U, cols_norm.fortran_vec () ,
532 rows_norm.fortran_vec (), Complex (droptol), milu); 529 rows_norm.fortran_vec (), Complex (droptol), milu);
533 530
534 param_list.append (octave_value (L.cols ())); 531 param_list.append (octave_value (L.cols ()));
535 SparseComplexMatrix eye = 532 SparseComplexMatrix eye =
536 feval ("speye", param_list)(0).sparse_complex_matrix_value (); 533 feval ("speye", param_list)(0).sparse_complex_matrix_value ();
537 retval(1) = U; 534 return ovl (L + eye, U);
538 retval(0) = L + eye; 535 }
539 }
540
541 return retval;
542 } 536 }
543 537
544 // That function implements the IKJ and JKI variants of gaussian elimination 538 // That function implements the IKJ and JKI variants of gaussian elimination
545 // to perform the ILUTP decomposition. The behaviour is controlled by milu 539 // to perform the ILUTP decomposition. The behaviour is controlled by milu
546 // parameter. If milu = ['off'|'col'] the JKI version is performed taking 540 // parameter. If milu = ['off'|'col'] the JKI version is performed taking
942 @deftypefnx {} {[@var{L}, @var{U}] =} __ilutp__ (@var{A}, @var{droptol}, @var{thresh}, @var{milu}, @var{udiag})\n\ 936 @deftypefnx {} {[@var{L}, @var{U}] =} __ilutp__ (@var{A}, @var{droptol}, @var{thresh}, @var{milu}, @var{udiag})\n\
943 @deftypefnx {} {[@var{L}, @var{U}, @var{P}] =} __ilutp__ (@var{A}, @dots{})\n\ 937 @deftypefnx {} {[@var{L}, @var{U}, @var{P}] =} __ilutp__ (@var{A}, @dots{})\n\
944 Undocumented internal function.\n\ 938 Undocumented internal function.\n\
945 @end deftypefn") 939 @end deftypefn")
946 { 940 {
941 int nargin = args.length ();
942
943 if (nargout < 2 || nargout > 3 || nargin < 1 || nargin > 5)
944 print_usage ();
945
947 octave_value_list retval; 946 octave_value_list retval;
948
949 int nargin = args.length ();
950 std::string milu = ""; 947 std::string milu = "";
951 double droptol = 0; 948 double droptol = 0;
952 double thresh = 1; 949 double thresh = 1;
953 double udiag = 0; 950 double udiag = 0;
954
955 if (nargout < 2 || nargout > 3 || nargin < 1 || nargin > 5)
956 print_usage ();
957 951
958 // Don't repeat input validation of arguments done in ilu.m 952 // Don't repeat input validation of arguments done in ilu.m
959 if (nargin >= 2) 953 if (nargin >= 2)
960 droptol = args(1).double_value (); 954 droptol = args(1).double_value ();
961 955
983 else 977 else
984 param_list (1) = "cols"; 978 param_list (1) = "cols";
985 rc_norm = feval ("norm", param_list)(0).vector_value (); 979 rc_norm = feval ("norm", param_list)(0).vector_value ();
986 param_list.clear (); 980 param_list.clear ();
987 Array <octave_idx_type> perm (dim_vector (sm.cols (), 1)); 981 Array <octave_idx_type> perm (dim_vector (sm.cols (), 1));
988 SparseMatrix U; 982 SparseMatrix U, L;
989 SparseMatrix L;
990 ilu_tp <SparseMatrix, double> (sm, L, U, nnz_u, nnz_l, 983 ilu_tp <SparseMatrix, double> (sm, L, U, nnz_u, nnz_l,
991 rc_norm.fortran_vec (), 984 rc_norm.fortran_vec (),
992 perm, droptol, thresh, milu, udiag); 985 perm, droptol, thresh, milu, udiag);
993 param_list.append (octave_value (L.cols ())); 986 param_list.append (octave_value (L.cols ()));
994 SparseMatrix eye = 987 SparseMatrix eye =
1034 else 1027 else
1035 param_list(1) = "cols"; 1028 param_list(1) = "cols";
1036 rc_norm = feval ("norm", param_list)(0).complex_vector_value (); 1029 rc_norm = feval ("norm", param_list)(0).complex_vector_value ();
1037 Array <octave_idx_type> perm (dim_vector (sm.cols (), 1)); 1030 Array <octave_idx_type> perm (dim_vector (sm.cols (), 1));
1038 param_list.clear (); 1031 param_list.clear ();
1039 SparseComplexMatrix U; 1032 SparseComplexMatrix U, L;
1040 SparseComplexMatrix L;
1041 ilu_tp < SparseComplexMatrix, Complex> 1033 ilu_tp < SparseComplexMatrix, Complex>
1042 (sm, L, U, nnz_u, nnz_l, rc_norm.fortran_vec (), perm, 1034 (sm, L, U, nnz_u, nnz_l, rc_norm.fortran_vec (), perm,
1043 Complex (droptol), Complex (thresh), milu, udiag); 1035 Complex (droptol), Complex (thresh), milu, udiag);
1044 1036
1045 param_list.append (octave_value (L.cols ())); 1037 param_list.append (octave_value (L.cols ()));