40 template <
typename octave_matrix_t,
typename T>
41 void ilu_0 (octave_matrix_t& sm,
const std::string milu =
"off")
57 else if (milu ==
"col")
65 for (i = 0; i < n; i++)
67 for (k = 0; k < n; k++)
72 for (j = j1; j <= j2; j++)
79 while ((jrow < k) && (j <= j2))
83 tl = data[j] / data[uptr[jrow]];
86 for (jj = uptr[jrow] + 1; jj < cidx[jrow+1]; jj++)
91 data[jw] -= tl * data[jj];
93 data[jw] -= data[j] * data[jj];
100 r += data[j] * data[jj];
109 for (jj = uptr[k] + 1; jj < cidx[k+1]; jj++)
110 data[jj] /= data[uptr[k]];
113 error (
"ilu: A has a zero on the diagonal");
119 error (
"ilu: encountered a pivot equal to 0");
122 for (i = j1; i <= j2; i++)
126 sm = sm.transpose ();
129 DEFUN (__ilu0__, args, nargout,
131 @deftypefn {Built-in Function} {[@var{L}, @var{U}] =} __ilu0__ (@var{A})\n\
132 @deftypefnx {Built-in Function} {[@var{L}, @var{U}] =} __ilu0__ (@var{A}, @var{milu})\n\
133 @deftypefnx {Built-in Function} {[@var{L}, @var{U}, @var{P}] =} __ilu0__ (@var{A}, @dots{})\n\
134 Undocumented internal function.\n\
139 int nargin = args.
length ();
142 if (nargout > 2 || nargin < 1 || nargin > 2)
152 if (! args(0).is_complex_type ())
155 ilu_0 <SparseMatrix, double> (sm, milu);
159 retval(1) =
feval (
"triu", param_list)(0).sparse_matrix_value ();
165 feval (
"tril", param_list)(0).sparse_matrix_value ();
171 ilu_0 <SparseComplexMatrix, Complex> (sm, milu);
176 feval (
"triu", param_list)(0).sparse_complex_matrix_value ();
179 octave_value (sm.cols ())))(0).sparse_complex_matrix_value ();
182 eye +
feval (
"tril", param_list)(0).sparse_complex_matrix_value ();
189 template <
typename octave_matrix_t,
typename T>
190 void ilu_crout (octave_matrix_t& sm_l, octave_matrix_t& sm_u,
191 octave_matrix_t& L, octave_matrix_t& U, T* cols_norm,
192 T* rows_norm,
const T droptol = 0,
193 const std::string milu =
"off")
198 enum {OFF, ROW, COL};
201 else if (milu ==
"col")
206 octave_idx_type jrow, i, j, k, jj, total_len_l, total_len_u, max_len_u,
207 max_len_l, w_len_u, w_len_l, cols_list_len, rows_list_len;
210 sm_u = sm_u.transpose ();
212 max_len_u = sm_u.nnz ();
213 max_len_u += (0.1 * max_len_u) > n ? 0.1 * max_len_u : n;
214 max_len_l = sm_l.nnz ();
215 max_len_l += (0.1 * max_len_l) > n ? 0.1 * max_len_l : n;
219 T* data_in_u = sm_u.data ();
222 T* data_in_l = sm_l.data ();
249 cidx_u[0] = cidx_in_u[0];
250 cidx_l[0] = cidx_in_l[0];
251 for (i = 0; i < n; i++)
263 for (k = 0; k < n; k++)
266 for (i = cidx_in_l[k]; i < cidx_in_l[k+1]; i++)
267 w_data_l[ridx_in_l[i]] = data_in_l[i];
269 for (i = cidx_in_u[k]; i < cidx_in_u[k+1]; i++)
270 w_data_u[ridx_in_u[i]] = data_in_u[i];
273 for (j = 0; j < rows_list_len; j++)
275 if ((Ufirst[rows_list[j]] != -1))
276 for (jj = Ufirst[rows_list[j]]; jj < cidx_u[rows_list[j]+1]; jj++)
279 w_data_u[jrow] -= data_u[jj] * data_l[Lfirst[rows_list[j]]];
283 for (j = 0; j < cols_list_len; j++)
285 if (Lfirst[cols_list[j]] != -1)
286 for (jj = Lfirst[cols_list[j]]; jj < cidx_l[cols_list[j]+1]; jj++)
289 w_data_l[jrow] -= data_l[jj] * data_u[Ufirst[cols_list[j]]];
293 if ((max_len_u - total_len_u) < n)
295 max_len_u += (0.1 * max_len_u) > n ? 0.1 * max_len_u : n;
302 if ((max_len_l - total_len_l) < n)
304 max_len_l += (0.1 * max_len_l) > n ? 0.1 * max_len_l : n;
313 data_u[total_len_u] = w_data_u[k];
314 ridx_u[total_len_u] = k;
316 for (i = k + 1; i < n; i++)
318 if (w_data_u[i] != zero)
320 if (
std::abs (w_data_u[i]) < (droptol * rows_norm[k]))
323 cr_sum[k] += w_data_u[i];
325 cr_sum[i] += w_data_u[i];
329 data_u[total_len_u + w_len_u] = w_data_u[i];
330 ridx_u[total_len_u + w_len_u] = i;
336 if (w_data_l[i] != zero)
338 if (
std::abs (w_data_l[i]) < (droptol * cols_norm[k]))
341 cr_sum[k] += w_data_l[i];
343 cr_sum[i] += w_data_l[i];
347 data_l[total_len_l + w_len_l] = w_data_l[i];
348 ridx_l[total_len_l + w_len_l] = i;
357 if (opt == COL || opt == ROW)
358 data_u[total_len_u] += cr_sum[k];
361 if (data_u[total_len_u] == zero)
363 error (
"ilu: encountered a pivot equal to 0");
368 for (i = total_len_l ; i < (total_len_l + w_len_l); i++)
369 data_l[i] /= data_u[total_len_u];
372 total_len_u += w_len_u;
373 total_len_l += w_len_l;
376 if (total_len_l < 0 || total_len_u < 0)
378 error (
"ilu: integer overflow. Too many fill-in elements in L or U");
381 cidx_u[k+1] = cidx_u[k] - cidx_u[0] + w_len_u;
382 cidx_l[k+1] = cidx_l[k] - cidx_l[0] + w_len_l;
395 Ufirst[k] = cidx_u[k];
399 Lfirst[k] = cidx_l[k];
404 for (i = 0; i <= k; i++)
408 jj = ridx_u[Ufirst[i]];
411 if (Ufirst[i] < (cidx_u[i+1]))
414 if (Ufirst[i] == cidx_u[i+1])
417 jj = ridx_u[Ufirst[i]];
422 cols_list[cols_list_len] = i;
429 jj = ridx_l[Lfirst[i]];
431 if (Lfirst[i] < (cidx_l[i+1]))
434 if (Lfirst[i] == cidx_l[i+1])
437 jj = ridx_l[Lfirst[i]];
441 rows_list[rows_list_len] = i;
452 L = octave_matrix_t (n, n, total_len_l);
453 U = octave_matrix_t (n, n, total_len_u);
454 for (i = 0; i <= n; i++)
455 L.cidx (i) = cidx_l[i];
456 for (i = 0; i < total_len_l; i++)
458 L.ridx (i) = ridx_l[i];
459 L.data (i) = data_l[i];
461 for (i = 0; i <= n; i++)
462 U.cidx (i) = cidx_u[i];
463 for (i = 0; i < total_len_u; i++)
465 U.ridx (i) = ridx_u[i];
466 U.data (i) = data_u[i];
472 DEFUN (__iluc__, args, nargout,
474 @deftypefn {Built-in Function} {[@var{L}, @var{U}] =} __iluc__ (@var{A})\n\
475 @deftypefnx {Built-in Function} {[@var{L}, @var{U}] =} __iluc__ (@var{A}, @var{droptol})\n\
476 @deftypefnx {Built-in Function} {[@var{L}, @var{U}] =} __iluc__ (@var{A}, @var{droptol}, @var{milu})\n\
477 @deftypefnx {Built-in Function} {[@var{L}, @var{U}, @var{P}] =} __iluc__ (@var{A}, @dots{})\n\
478 Undocumented internal function.\n\
482 int nargin = args.
length ();
483 std::string milu =
"off";
486 if (nargout != 2 || nargin < 1 || nargin > 3)
494 droptol = args(1).double_value ();
497 milu = args(2).string_value ();
500 if (! args(0).is_complex_type ())
503 param_list.
append (args(0).sparse_matrix_value ());
507 param_list(1) =
"rows";
508 rows_norm =
feval (
"norm", param_list)(0).vector_value ();
509 param_list(1) =
"cols";
510 cols_norm =
feval (
"norm", param_list)(0).vector_value ();
514 ilu_crout <SparseMatrix, double> (sm_l, sm_u, L, U,
522 feval (
"speye", param_list)(0).sparse_matrix_value ();
530 param_list.
append (args(0).sparse_complex_matrix_value ());
532 feval(
"triu", param_list)(0).sparse_complex_matrix_value ();
535 feval(
"tril", param_list)(0).sparse_complex_matrix_value ();
536 param_list(1) =
"rows";
537 rows_norm =
feval (
"norm", param_list)(0).complex_vector_value ();
538 param_list(1) =
"cols";
539 cols_norm =
feval (
"norm", param_list)(0).complex_vector_value ();
543 ilu_crout < SparseComplexMatrix, Complex >
550 feval (
"speye", param_list)(0).sparse_complex_matrix_value ();
567 template <
typename octave_matrix_t,
typename T>
568 void ilu_tp (octave_matrix_t& sm, octave_matrix_t& L, octave_matrix_t& U,
571 const T thresh = T(0),
const std::string milu =
"off",
572 const double udiag = 0)
575 enum {OFF, ROW, COL};
578 else if (milu ==
"col")
592 T* data_in = sm.data ();
593 octave_idx_type jrow, i, j, k, jj, c, total_len_l, total_len_u, p_perm,
594 max_ind, max_len_l, max_len_u;
597 T tl = zero, aux, maximum;
600 max_len_u += (0.1 * max_len_u) > n ? 0.1 * max_len_u : n;
602 max_len_l += (0.1 * max_len_l) > n ? 0.1 * max_len_l : n;
620 T total_sum, partial_col_sum = zero, partial_row_sum = zero;
621 std::set <octave_idx_type> iw_l;
622 std::set <octave_idx_type> iw_u;
623 std::set <octave_idx_type>::iterator it, it2;
630 cidx_l[0] = cidx_in[0];
631 cidx_u[0] = cidx_in[0];
632 for (i = 0; i < n; i++)
641 for (k = 0; k < n; k++)
644 for (j = cidx_in[k]; j < cidx_in[k+1]; j++)
646 p_perm = iperm[ridx_in[j]];
647 w_data[iperm[ridx_in[j]]] = data_in[j];
649 iw_l.insert (ridx_in[j]);
651 iw_u.insert (p_perm);
657 while ((jrow < k) && (it != iw_u.end ()))
660 partial_col_sum = w_data[jrow];
661 if (w_data[jrow] != zero)
665 partial_row_sum = w_data[jrow];
666 tl = w_data[jrow] / data_u[uptr[jrow]];
668 for (jj = cidx_l[jrow]; jj < cidx_l[jrow+1]; jj++)
670 p_perm = iperm[ridx_l[jj]];
671 aux = w_data[p_perm];
674 w_data[p_perm] -= tl * data_l[jj];
675 partial_row_sum += tl * data_l[jj];
679 tl = data_l[jj] * w_data[jrow];
680 w_data[p_perm] -= tl;
682 partial_col_sum += tl;
685 if ((aux == zero) && (w_data[p_perm] != zero))
688 iw_l.insert (ridx_l[jj]);
690 iw_u.insert (p_perm);
696 if ((
std::abs (w_data[jrow]) < (droptol * cols_norm[k]))
697 && (w_data[jrow] != zero))
700 total_sum += partial_col_sum;
702 total_sum += partial_row_sum;
721 if ((thresh > zero) && (k < (n - 1)))
723 maximum =
std::abs (w_data[k]) / thresh;
725 for (it = iw_l.begin (); it != iw_l.end (); ++it)
728 if (
std::abs (w_data[p_perm]) > maximum)
730 maximum =
std::abs (w_data[p_perm]);
736 p_perm = iperm[max_ind];
737 if (max_ind != perm[k])
740 if (w_data[k] != zero)
741 iw_l.insert (perm[k]);
746 iperm[perm[p_perm]] = k;
747 iperm[perm[k]] = p_perm;
749 perm[k] = perm[p_perm];
751 w_data[k] = w_data[p_perm];
752 w_data[p_perm] = aux;
760 while (it != iw_l.end ())
764 if (
std::abs (w_data[p_perm]) < (droptol * cols_norm[k]))
767 total_sum += w_data[p_perm];
768 w_data[p_perm] = zero;
782 if ((total_sum > zero) && (w_data[k] == zero))
784 w_data[k] += total_sum;
792 if (w_data[k] == zero)
801 error (
"ilu: encountered a pivot equal to 0");
808 for (it = iw_l.begin (); it != iw_l.end (); ++it)
811 w_data[p_perm] = w_data[p_perm] / w_data[k];
815 if ((max_len_u - total_len_u) < n)
817 max_len_u += (0.1 * max_len_u) > n ? 0.1 * max_len_u : n;
824 if ((max_len_l - total_len_l) < n)
826 max_len_l += (0.1 * max_len_l) > n ? 0.1 * max_len_l : n;
835 for (it = iw_u.begin (); it != iw_u.end (); ++it)
837 if (w_data[*it] != zero)
839 data_u[total_len_u + w_len_u] = w_data[*it];
840 ridx_u[total_len_u + w_len_u] = *it;
847 for (it = iw_l.begin (); it != iw_l.end (); ++it)
850 if (w_data[p_perm] != zero)
852 data_l[total_len_l + w_len_l] = w_data[p_perm];
853 ridx_l[total_len_l + w_len_l] = *it;
858 total_len_u += w_len_u;
859 total_len_l += w_len_l;
862 if (total_len_l < 0 || total_len_u < 0)
864 error (
"ilu: Integer overflow. Too many fill-in elements in L or U");
868 uptr[k] = total_len_u - 1;
869 cidx_u[k+1] = cidx_u[k] - cidx_u[0] + w_len_u;
870 cidx_l[k+1] = cidx_l[k] - cidx_l[0] + w_len_l;
878 octave_matrix_t *L_ptr;
879 octave_matrix_t *U_ptr;
880 octave_matrix_t diag (n, n, n);
889 L = octave_matrix_t (n, n, total_len_u - n);
890 U = octave_matrix_t (n, n, total_len_l);
896 L = octave_matrix_t (n, n, total_len_l);
897 U = octave_matrix_t (n, n, total_len_u);
900 for (i = 0; i <= n; i++)
902 L_ptr->cidx (i) = cidx_l[i];
903 U_ptr->cidx (i) = cidx_u[i];
905 U_ptr->cidx (i) -= i;
908 for (i = 0; i < n; i++)
911 diag.elem (i,i) = data_u[uptr[i]];
914 while (j < cidx_l[i+1])
916 L_ptr->ridx (j) = ridx_l[j];
917 L_ptr->data (j) = data_l[j];
922 while (j < cidx_u[i+1])
938 U_ptr->data (c) = data_u[j];
939 U_ptr->ridx (c) = ridx_u[j];
954 DEFUN (__ilutp__, args, nargout,
956 @deftypefn {Built-in Function} {[@var{L}, @var{U}] =} __ilutp__ (@var{A})\n\
957 @deftypefnx {Built-in Function} {[@var{L}, @var{U}] =} __ilutp__ (@var{A}, @var{droptol})\n\
958 @deftypefnx {Built-in Function} {[@var{L}, @var{U}] =} __ilutp__ (@var{A}, @var{droptol}, @var{thresh})\n\
959 @deftypefnx {Built-in Function} {[@var{L}, @var{U}] =} __ilutp__ (@var{A}, @var{droptol}, @var{thresh}, @var{milu})\n\
960 @deftypefnx {Built-in Function} {[@var{L}, @var{U}] =} __ilutp__ (@var{A}, @var{droptol}, @var{thresh}, @var{milu}, @var{udiag})\n\
961 @deftypefnx {Built-in Function} {[@var{L}, @var{U}, @var{P}] =} __ilutp__ (@var{A}, @dots{})\n\
962 Undocumented internal function.\n\
967 int nargin = args.
length ();
968 std::string milu =
"";
969 double droptol = 0, thresh = 1;
972 if (nargout < 2 || nargout > 3 || nargin < 1 || nargin > 5)
980 droptol = args(1).double_value ();
983 thresh = args(2).double_value ();
986 milu = args(3).string_value ();
989 udiag = args(4).double_value ();
993 if (! args(0).is_complex_type ())
998 nnz_u = (
feval (
"triu", param_list)(0).sparse_matrix_value ()).nnz ();
1000 nnz_l = (
feval (
"tril", param_list)(0).sparse_matrix_value ()).nnz ();
1002 param_list (1) =
"rows";
1004 param_list (1) =
"cols";
1005 rc_norm =
feval (
"norm", param_list)(0).vector_value ();
1006 param_list.
clear ();
1010 ilu_tp <SparseMatrix, double> (sm, L, U, nnz_u, nnz_l,
1012 perm, droptol, thresh, milu, udiag);
1017 feval (
"speye", param_list)(0).sparse_matrix_value ();
1025 else if (nargout == 2)
1027 retval(0) = L + eye;
1051 feval (
"triu", param_list)(0).sparse_complex_matrix_value ().nnz ();
1054 feval (
"tril", param_list)(0).sparse_complex_matrix_value ().nnz ();
1056 param_list(1) =
"rows";
1058 param_list(1) =
"cols";
1059 rc_norm =
feval (
"norm", param_list)(0).complex_vector_value ();
1061 param_list.
clear ();
1064 ilu_tp < SparseComplexMatrix, Complex>
1065 (sm, L, U, nnz_u, nnz_l, rc_norm.
fortran_vec (), perm,
1072 feval (
"speye", param_list)(0).sparse_complex_matrix_value ();
1080 else if (nargout == 2)
1082 retval(0) = L + eye;
octave_idx_type cols(void) const
static const idx_vector colon
OCTINTERP_API void print_usage(void)
octave_idx_type length(void) const
octave_value_list & append(const octave_value &val)
#define DEFUN(name, args_name, nargout_name, doc)
void error(const char *fmt,...)
Sparse< T > index(const idx_vector &i, bool resize_ok=false) const
octave_value_list feval(const std::string &name, const octave_value_list &args, int nargout)
void resize(const dim_vector &dv, const T &rfv)
void ilu_tp(octave_matrix_t &sm, octave_matrix_t &L, octave_matrix_t &U, octave_idx_type nnz_u, octave_idx_type nnz_l, T *cols_norm, Array< octave_idx_type > &perm_vec, const T droptol=T(0), const T thresh=T(0), const std::string milu="off", const double udiag=0)
void ilu_0(octave_matrix_t &sm, const std::string milu="off")
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
std::complex< double > Complex
const T * fortran_vec(void) const
void ilu_crout(octave_matrix_t &sm_l, octave_matrix_t &sm_u, octave_matrix_t &L, octave_matrix_t &U, T *cols_norm, T *rows_norm, const T droptol=0, const std::string milu="off")