24 #if defined (HAVE_CONFIG_H)
40 template <
typename octave_matrix_t,
typename T>
55 else if (milu ==
"col")
70 for (i = 0; i < n; i++)
74 for (k = 0; k < n; k++)
80 error (
"ilu: A has a zero on the diagonal");
82 for (j = j1; j < j2; j++)
88 while ((jrow < k) && (j < j2))
92 tl = data[j] / data[uptr[jrow]];
95 for (jj = uptr[jrow] + 1; jj < cidx[jrow+1]; jj++)
100 data[jw] -= tl * data[jj];
102 data[jw] -= data[j] * data[jj];
109 r += data[j] * data[jj];
119 for (jj = uptr[k] + 1; jj < cidx[k+1]; jj++)
120 data[jj] /= data[uptr[k]];
123 error (
"ilu: A has a zero on the diagonal");
126 error (
"ilu: encountered a pivot equal to 0");
128 for (i = j1; i < j2; i++)
133 sm = sm.transpose ();
146 if (nargin < 1 || nargin > 2)
157 if (!
args(0).is_complex_type ())
160 ilu_0 <SparseMatrix, double> (sm, milu);
163 retval(1) =
feval (
"triu", arg_list)(0).sparse_matrix_value ();
165 feval (
"speye",
ovl (sm.cols ()))(0).sparse_matrix_value ();
168 feval (
"tril", arg_list)(0).sparse_matrix_value ();
173 ilu_0 <SparseComplexMatrix, Complex> (sm, milu);
176 retval(1) =
feval (
"triu", arg_list)(0).sparse_complex_matrix_value ();
178 feval (
"speye",
ovl (sm.cols ()))(0).sparse_complex_matrix_value ();
181 feval (
"tril", arg_list)(0).sparse_complex_matrix_value ();
187 template <
typename octave_matrix_t,
typename T>
188 void ilu_crout (octave_matrix_t& sm_l, octave_matrix_t& sm_u,
189 octave_matrix_t& L, octave_matrix_t& U, T* cols_norm,
190 T* rows_norm,
const T droptol = 0,
195 enum {OFF, ROW, COL};
198 else if (milu ==
"col")
204 max_len_l, w_len_u, w_len_l, cols_list_len, rows_list_len;
207 sm_u = sm_u.transpose ();
209 max_len_u = sm_u.nnz ();
210 max_len_u += (0.1 * max_len_u) > n ? 0.1 * max_len_u : n;
211 max_len_l = sm_l.nnz ();
212 max_len_l += (0.1 * max_len_l) > n ? 0.1 * max_len_l : n;
217 T* data_in_u = sm_u.data ();
220 T* data_in_l = sm_l.data ();
226 T* data_l = data_out_l.fortran_vec ();
232 T* data_u = data_out_u.fortran_vec ();
248 cidx_u[0] = cidx_in_u[0];
249 cidx_l[0] = cidx_in_l[0];
250 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;
296 data_out_u.resize (
dim_vector (max_len_u, 1));
297 data_u = data_out_u.fortran_vec ();
298 ridx_out_u.resize (
dim_vector (max_len_u, 1));
299 ridx_u = ridx_out_u.fortran_vec ();
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;
305 data_out_l.resize (
dim_vector (max_len_l, 1));
306 data_l = data_out_l.fortran_vec ();
307 ridx_out_l.resize (
dim_vector (max_len_l, 1));
308 ridx_l = ridx_out_l.fortran_vec ();
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)
362 error (
"ilu: encountered a pivot equal to 0");
365 for (i = total_len_l ; i < (total_len_l + w_len_l); i++)
366 data_l[i] /= data_u[total_len_u];
368 total_len_u += w_len_u;
369 total_len_l += w_len_l;
372 if (total_len_l < 0 || total_len_u < 0)
373 error (
"ilu: integer overflow. Too many fill-in elements in L or U");
375 cidx_u[k+1] = cidx_u[
k] - cidx_u[0] + w_len_u;
376 cidx_l[k+1] = cidx_l[
k] - cidx_l[0] + w_len_l;
389 Ufirst[
k] = cidx_u[
k];
393 Lfirst[
k] = cidx_l[
k];
398 for (i = 0; i <=
k; i++)
402 jj = ridx_u[Ufirst[
i]];
405 if (Ufirst[i] < (cidx_u[i+1]))
408 if (Ufirst[i] == cidx_u[i+1])
411 jj = ridx_u[Ufirst[
i]];
416 cols_list[cols_list_len] =
i;
423 jj = ridx_l[Lfirst[
i]];
425 if (Lfirst[i] < (cidx_l[i+1]))
428 if (Lfirst[i] == cidx_l[i+1])
431 jj = ridx_l[Lfirst[
i]];
435 rows_list[rows_list_len] =
i;
444 L = octave_matrix_t (n, n, total_len_l);
445 U = octave_matrix_t (n, n, total_len_u);
448 for (i = 0; i <= n; i++)
449 L.cidx (i) = cidx_l[
i];
451 for (i = 0; i < total_len_l; i++)
453 L.ridx (i) = ridx_l[
i];
454 L.data (i) = data_l[
i];
457 for (i = 0; i <= n; i++)
458 U.cidx (i) = cidx_u[
i];
460 for (i = 0; i < total_len_u; i++)
462 U.ridx (i) = ridx_u[
i];
463 U.data (i) = data_u[
i];
479 if (nargin < 1 || nargin > 3)
487 droptol =
args(1).double_value ();
490 milu =
args(2).string_value ();
493 if (!
args(0).is_complex_type ())
496 arg_list.
append (
args(0).sparse_matrix_value ());
500 arg_list(1) =
"rows";
501 rows_norm =
feval (
"norm", arg_list)(0).vector_value ();
502 arg_list(1) =
"cols";
503 cols_norm =
feval (
"norm", arg_list)(0).vector_value ();
506 ilu_crout <SparseMatrix, double> (sm_l, sm_u, L, U,
512 return ovl (L + eye, U);
517 arg_list.
append (
args(0).sparse_complex_matrix_value ());
519 feval (
"triu", arg_list)(0).sparse_complex_matrix_value ();
522 feval (
"tril", arg_list)(0).sparse_complex_matrix_value ();
523 arg_list(1) =
"rows";
524 rows_norm =
feval (
"norm", arg_list)(0).complex_vector_value ();
525 arg_list(1) =
"cols";
526 cols_norm =
feval (
"norm", arg_list)(0).complex_vector_value ();
529 ilu_crout <SparseComplexMatrix, Complex>
535 feval (
"speye", arg_list)(0).sparse_complex_matrix_value ();
536 return ovl (L + eye, U);
548 template <
typename octave_matrix_t,
typename T>
549 void ilu_tp (octave_matrix_t& sm, octave_matrix_t& L, octave_matrix_t& U,
552 const T thresh = T(0),
const std::string milu =
"off",
553 const double udiag = 0)
556 enum {OFF, ROW, COL};
559 else if (milu ==
"col")
568 sm = sm.transpose ();
573 T* data_in = sm.data ();
575 max_ind, max_len_l, max_len_u;
578 T tl =
zero, aux, maximum;
581 max_len_u += (0.1 * max_len_u) > n ? 0.1 * max_len_u : n;
583 max_len_l += (0.1 * max_len_l) > n ? 0.1 * max_len_l : n;
591 T* data_l = data_out_l.fortran_vec ();
599 T* data_u = data_out_u.fortran_vec ();
603 T total_sum, partial_col_sum =
zero, partial_row_sum =
zero;
604 std::set <octave_idx_type> iw_l;
605 std::set <octave_idx_type> iw_u;
606 std::set <octave_idx_type>::iterator it, it2;
613 cidx_l[0] = cidx_in[0];
614 cidx_u[0] = cidx_in[0];
615 for (i = 0; i < n; i++)
625 for (k = 0; k < n; k++)
627 for (j = cidx_in[k]; j < cidx_in[k+1]; j++)
629 p_perm = iperm[ridx_in[j]];
630 w_data[iperm[ridx_in[j]]] = data_in[j];
632 iw_l.insert (ridx_in[j]);
634 iw_u.insert (p_perm);
640 while ((jrow < k) && (it != iw_u.end ()))
643 partial_col_sum = w_data[jrow];
644 if (w_data[jrow] != zero)
648 partial_row_sum = w_data[jrow];
649 tl = w_data[jrow] / data_u[uptr[jrow]];
651 for (jj = cidx_l[jrow]; jj < cidx_l[jrow+1]; jj++)
653 p_perm = iperm[ridx_l[jj]];
654 aux = w_data[p_perm];
657 w_data[p_perm] -= tl * data_l[jj];
658 partial_row_sum += tl * data_l[jj];
662 tl = data_l[jj] * w_data[jrow];
663 w_data[p_perm] -= tl;
665 partial_col_sum += tl;
668 if ((aux == zero) && (w_data[p_perm] !=
zero))
671 iw_l.insert (ridx_l[jj]);
673 iw_u.insert (p_perm);
679 if ((
std::abs (w_data[jrow]) < (droptol * cols_norm[k]))
680 && (w_data[jrow] != zero))
683 total_sum += partial_col_sum;
685 total_sum += partial_row_sum;
704 if ((thresh > zero) && (k < (n - 1)))
706 maximum =
std::abs (w_data[k]) / thresh;
708 for (it = iw_l.begin (); it != iw_l.end (); ++it)
711 if (
std::abs (w_data[p_perm]) > maximum)
713 maximum =
std::abs (w_data[p_perm]);
719 p_perm = iperm[max_ind];
720 if (max_ind != perm[k])
723 if (w_data[k] != zero)
724 iw_l.insert (perm[k]);
729 iperm[perm[p_perm]] =
k;
730 iperm[perm[
k]] = p_perm;
732 perm[
k] = perm[p_perm];
734 w_data[
k] = w_data[p_perm];
735 w_data[p_perm] = aux;
743 while (it != iw_l.end ())
747 if (
std::abs (w_data[p_perm]) < (droptol * cols_norm[
k]))
750 total_sum += w_data[p_perm];
751 w_data[p_perm] =
zero;
765 if ((total_sum > zero) && (w_data[k] == zero))
767 w_data[
k] += total_sum;
773 if (w_data[k] == zero)
776 error (
"ilu: encountered a pivot equal to 0");
784 for (it = iw_l.begin (); it != iw_l.end (); ++it)
787 w_data[p_perm] = w_data[p_perm] / w_data[
k];
790 if ((max_len_u - total_len_u) < n)
792 max_len_u += (0.1 * max_len_u) > n ? 0.1 * max_len_u : n;
793 data_out_u.resize (
dim_vector (max_len_u, 1));
794 data_u = data_out_u.fortran_vec ();
795 ridx_out_u.resize (
dim_vector (max_len_u, 1));
796 ridx_u = ridx_out_u.fortran_vec ();
799 if ((max_len_l - total_len_l) < n)
801 max_len_l += (0.1 * max_len_l) > n ? 0.1 * max_len_l : n;
802 data_out_l.resize (
dim_vector (max_len_l, 1));
803 data_l = data_out_l.fortran_vec ();
804 ridx_out_l.resize (
dim_vector (max_len_l, 1));
805 ridx_l = ridx_out_l.fortran_vec ();
810 for (it = iw_u.begin (); it != iw_u.end (); ++it)
812 if (w_data[*it] != zero)
814 data_u[total_len_u + w_len_u] = w_data[*it];
815 ridx_u[total_len_u + w_len_u] = *it;
823 for (it = iw_l.begin (); it != iw_l.end (); ++it)
826 if (w_data[p_perm] != zero)
828 data_l[total_len_l + w_len_l] = w_data[p_perm];
829 ridx_l[total_len_l + w_len_l] = *it;
834 total_len_u += w_len_u;
835 total_len_l += w_len_l;
839 if (total_len_l < 0 || total_len_u < 0)
840 error (
"ilu: Integer overflow. Too many fill-in elements in L or U");
843 uptr[
k] = total_len_u - 1;
845 cidx_u[k+1] = cidx_u[
k] - cidx_u[0] + w_len_u;
846 cidx_l[k+1] = cidx_l[
k] - cidx_l[0] + w_len_l;
852 octave_matrix_t *L_ptr;
853 octave_matrix_t *U_ptr;
854 octave_matrix_t diag (n, n, n);
863 L = octave_matrix_t (n, n, total_len_u - n);
864 U = octave_matrix_t (n, n, total_len_l);
870 L = octave_matrix_t (n, n, total_len_l);
871 U = octave_matrix_t (n, n, total_len_u);
874 for (i = 0; i <= n; i++)
876 L_ptr->cidx (i) = cidx_l[
i];
877 U_ptr->cidx (i) = cidx_u[
i];
879 U_ptr->cidx (i) -=
i;
882 for (i = 0; i < n; i++)
885 diag.elem (i,i) = data_u[uptr[
i]];
888 while (j < cidx_l[i+1])
890 L_ptr->ridx (j) = ridx_l[j];
891 L_ptr->data (j) = data_l[j];
896 while (j < cidx_u[i+1])
912 U_ptr->data (c) = data_u[j];
913 U_ptr->ridx (c) = ridx_u[j];
938 int nargin =
args.length ();
940 if (nargin < 1 || nargin > 5)
951 droptol =
args(1).double_value ();
954 thresh =
args(2).double_value ();
957 milu =
args(3).string_value ();
960 udiag =
args(4).double_value ();
964 if (!
args(0).is_complex_type ())
969 nnz_u = (
feval (
"triu", arg_list)(0).sparse_matrix_value ()).nnz ();
971 nnz_l = (
feval (
"tril", arg_list)(0).sparse_matrix_value ()).nnz ();
973 arg_list (1) =
"rows";
975 arg_list (1) =
"cols";
976 rc_norm =
feval (
"norm", arg_list)(0).vector_value ();
980 ilu_tp <SparseMatrix, double> (sm, L, U, nnz_u, nnz_l,
982 perm, droptol, thresh, milu, udiag);
985 feval (
"speye", arg_list)(0).sparse_matrix_value ();
1018 feval (
"triu", arg_list)(0).sparse_complex_matrix_value ().nnz ();
1021 feval (
"tril", arg_list)(0).sparse_complex_matrix_value ().nnz ();
1023 arg_list(1) =
"rows";
1025 arg_list(1) =
"cols";
1026 rc_norm =
feval (
"norm", arg_list)(0).complex_vector_value ();
1030 ilu_tp <SparseComplexMatrix, Complex>
1031 (sm, L, U, nnz_u, nnz_l, rc_norm.
fortran_vec (), perm,
1036 feval (
"speye", arg_list)(0).sparse_complex_matrix_value ();
octave_idx_type cols(void) const
static const idx_vector colon
OCTAVE_EXPORT octave_value_list isa nd deftypefn *return ovl(args(0).is_integer_type())
OCTINTERP_API void print_usage(void)
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_EXPORT octave_value_list return the number of command line arguments passed to Octave If called with the optional argument the function xample nargout(@histc)
the sparsity preserving column transformation such that that defines the pivoting threshold can be given in which case it defines the c
=val(i)}if ode{val(i)}occurs in table i
OCTAVE_EXPORT octave_value_list return the value of the option it must match the dimension of the state and the relative tolerance must also be a vector of the same length tem it must match the dimension of the state and the absolute tolerance must also be a vector of the same length The local error test applied at each integration step is xample roup abs(local error in x(i))<
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
Vector representing the dimensions (size) of an Array.
If this string is the system will ring the terminal sometimes it is useful to be able to print the original representation of the string