24 #if defined (HAVE_CONFIG_H) 34 #include "builtin-defun-decls.h" 42 template <
typename octave_matrix_t,
typename T>
57 else if (milu ==
"col")
72 for (
i = 0;
i < n;
i++)
76 for (
k = 0;
k < n;
k++)
82 error (
"ilu: A has a zero on the diagonal");
84 for (j = j1; j < j2; j++)
90 while ((jrow <
k) && (j < j2))
94 tl = data[j] / data[uptr[jrow]];
97 for (jj = uptr[jrow] + 1; jj < cidx[jrow+1]; jj++)
102 data[jw] -= tl * data[jj];
104 data[jw] -= data[j] * data[jj];
111 r += data[j] * data[jj];
121 for (jj = uptr[
k] + 1; jj < cidx[
k+1]; jj++)
122 data[jj] /= data[uptr[
k]];
125 error (
"ilu: A has a zero on the diagonal");
128 error (
"ilu: encountered a pivot equal to 0");
130 for (
i = j1;
i < j2;
i++)
135 sm = sm.transpose ();
138 DEFUN (__ilu0__, args, ,
145 if (args.length () != 2)
156 if (! args(0).iscomplex ())
161 ilu_0 <SparseMatrix, double> (sm, milu);
163 retval(0) = speye +
Ftril (
ovl (sm, -1))(0).sparse_matrix_value ();
164 retval(1) = Ftriu (
ovl (sm))(0).sparse_matrix_value ();
171 ilu_0 <SparseComplexMatrix, Complex> (sm, milu);
174 Ftril (
ovl (sm, -1))(0).sparse_complex_matrix_value ();
175 retval(1) = Ftriu (
ovl (sm))(0).sparse_complex_matrix_value ();
181 template <
typename octave_matrix_t,
typename T>
182 void ilu_crout (octave_matrix_t& sm_l, octave_matrix_t& sm_u,
183 octave_matrix_t& L, octave_matrix_t& U, T *cols_norm,
184 T *rows_norm,
const T droptol = 0,
189 enum {OFF, ROW, COL};
192 else if (milu ==
"col")
198 max_len_l, w_len_u, w_len_l, cols_list_len, rows_list_len;
201 sm_u = sm_u.transpose ();
203 max_len_u = sm_u.nnz ();
204 max_len_u += (0.1 * max_len_u) > n ? 0.1 * max_len_u : n;
205 max_len_l = sm_l.nnz ();
206 max_len_l += (0.1 * max_len_l) > n ? 0.1 * max_len_l : n;
211 T *data_in_u = sm_u.data ();
214 T *data_in_l = sm_l.data ();
220 T *data_l = data_out_l.fortran_vec ();
226 T *data_u = data_out_u.fortran_vec ();
242 cidx_u[0] = cidx_in_u[0];
243 cidx_l[0] = cidx_in_l[0];
244 for (
i = 0;
i < n;
i++)
257 for (
k = 0;
k < n;
k++)
260 for (
i = cidx_in_l[
k];
i < cidx_in_l[
k+1];
i++)
261 w_data_l[ridx_in_l[
i]] = data_in_l[
i];
263 for (
i = cidx_in_u[
k];
i < cidx_in_u[
k+1];
i++)
264 w_data_u[ridx_in_u[
i]] = data_in_u[
i];
267 for (j = 0; j < rows_list_len; j++)
269 if ((Ufirst[rows_list[j]] != -1))
270 for (jj = Ufirst[rows_list[j]]; jj < cidx_u[rows_list[j]+1]; jj++)
273 w_data_u[jrow] -= data_u[jj] * data_l[Lfirst[rows_list[j]]];
277 for (j = 0; j < cols_list_len; j++)
279 if (Lfirst[cols_list[j]] != -1)
280 for (jj = Lfirst[cols_list[j]]; jj < cidx_l[cols_list[j]+1]; jj++)
283 w_data_l[jrow] -= data_l[jj] * data_u[Ufirst[cols_list[j]]];
287 if ((max_len_u - total_len_u) < n)
289 max_len_u += (0.1 * max_len_u) > n ? 0.1 * max_len_u : n;
290 data_out_u.resize (
dim_vector (max_len_u, 1));
291 data_u = data_out_u.fortran_vec ();
292 ridx_out_u.resize (
dim_vector (max_len_u, 1));
293 ridx_u = ridx_out_u.fortran_vec ();
296 if ((max_len_l - total_len_l) < n)
298 max_len_l += (0.1 * max_len_l) > n ? 0.1 * max_len_l : n;
299 data_out_l.resize (
dim_vector (max_len_l, 1));
300 data_l = data_out_l.fortran_vec ();
301 ridx_out_l.resize (
dim_vector (max_len_l, 1));
302 ridx_l = ridx_out_l.fortran_vec ();
307 data_u[total_len_u] = w_data_u[
k];
308 ridx_u[total_len_u] =
k;
310 for (
i =
k + 1;
i < n;
i++)
312 if (w_data_u[
i] !=
zero)
314 if (
std::abs (w_data_u[
i]) < (droptol * rows_norm[
k]))
317 cr_sum[
k] += w_data_u[
i];
319 cr_sum[
i] += w_data_u[
i];
323 data_u[total_len_u + w_len_u] = w_data_u[
i];
324 ridx_u[total_len_u + w_len_u] =
i;
330 if (w_data_l[
i] !=
zero)
332 if (
std::abs (w_data_l[
i]) < (droptol * cols_norm[
k]))
335 cr_sum[
k] += w_data_l[
i];
337 cr_sum[
i] += w_data_l[
i];
341 data_l[total_len_l + w_len_l] = w_data_l[
i];
342 ridx_l[total_len_l + w_len_l] =
i;
351 if (opt == COL || opt == ROW)
352 data_u[total_len_u] += cr_sum[
k];
355 if (data_u[total_len_u] ==
zero)
356 error (
"ilu: encountered a pivot equal to 0");
359 for (
i = total_len_l ;
i < (total_len_l + w_len_l);
i++)
360 data_l[
i] /= data_u[total_len_u];
362 total_len_u += w_len_u;
363 total_len_l += w_len_l;
366 if (total_len_l < 0 || total_len_u < 0)
367 error (
"ilu: integer overflow. Too many fill-in elements in L or U");
369 cidx_u[
k+1] = cidx_u[
k] - cidx_u[0] + w_len_u;
370 cidx_l[
k+1] = cidx_l[
k] - cidx_l[0] + w_len_l;
383 Ufirst[
k] = cidx_u[
k];
387 Lfirst[
k] = cidx_l[
k];
392 for (
i = 0;
i <=
k;
i++)
396 jj = ridx_u[Ufirst[
i]];
399 if (Ufirst[
i] < (cidx_u[
i+1]))
402 if (Ufirst[
i] == cidx_u[
i+1])
405 jj = ridx_u[Ufirst[
i]];
410 cols_list[cols_list_len] =
i;
417 jj = ridx_l[Lfirst[
i]];
419 if (Lfirst[
i] < (cidx_l[
i+1]))
422 if (Lfirst[
i] == cidx_l[
i+1])
425 jj = ridx_l[Lfirst[
i]];
429 rows_list[rows_list_len] =
i;
438 L = octave_matrix_t (n, n, total_len_l);
439 U = octave_matrix_t (n, n, total_len_u);
442 for (
i = 0;
i <= n;
i++)
443 L.cidx (
i) = cidx_l[
i];
445 for (
i = 0;
i < total_len_l;
i++)
447 L.ridx (
i) = ridx_l[
i];
448 L.data (
i) = data_l[
i];
451 for (
i = 0;
i <= n;
i++)
452 U.cidx (
i) = cidx_u[
i];
454 for (
i = 0;
i < total_len_u;
i++)
456 U.ridx (
i) = ridx_u[
i];
457 U.data (
i) = data_u[
i];
463 DEFUN (__iluc__, args, ,
469 if (args.length () != 3)
472 double droptol = args(1).double_value ();
475 if (! args(0).iscomplex ())
482 ilu_crout <SparseMatrix, double> (sm_l, sm_u, L, U,
489 return ovl (L + speye, U);
500 ilu_crout <SparseComplexMatrix, Complex> (sm_l, sm_u, L, U,
507 return ovl (L + speye, U);
519 template <
typename octave_matrix_t,
typename T>
520 void ilu_tp (octave_matrix_t& sm, octave_matrix_t& L, octave_matrix_t& U,
523 const T thresh = T(0),
const std::string milu =
"off",
524 const double udiag = 0)
527 enum {OFF, ROW, COL};
530 else if (milu ==
"col")
539 sm = sm.transpose ();
544 T *data_in = sm.data ();
546 max_ind, max_len_l, max_len_u;
549 T tl =
zero, aux, maximum;
552 max_len_u += (0.1 * max_len_u) > n ? 0.1 * max_len_u : n;
554 max_len_l += (0.1 * max_len_l) > n ? 0.1 * max_len_l : n;
562 T *data_l = data_out_l.fortran_vec ();
570 T *data_u = data_out_u.fortran_vec ();
574 T total_sum, partial_col_sum =
zero, partial_row_sum =
zero;
575 std::set <octave_idx_type> iw_l;
576 std::set <octave_idx_type> iw_u;
577 std::set <octave_idx_type>::iterator it, it2;
584 cidx_l[0] = cidx_in[0];
585 cidx_u[0] = cidx_in[0];
586 for (
i = 0;
i < n;
i++)
596 for (
k = 0;
k < n;
k++)
598 for (j = cidx_in[
k]; j < cidx_in[
k+1]; j++)
600 p_perm = iperm[ridx_in[j]];
601 w_data[iperm[ridx_in[j]]] = data_in[j];
603 iw_l.insert (ridx_in[j]);
605 iw_u.insert (p_perm);
611 while ((jrow <
k) && (it != iw_u.end ()))
614 partial_col_sum = w_data[jrow];
615 if (w_data[jrow] !=
zero)
619 partial_row_sum = w_data[jrow];
620 tl = w_data[jrow] / data_u[uptr[jrow]];
622 for (jj = cidx_l[jrow]; jj < cidx_l[jrow+1]; jj++)
624 p_perm = iperm[ridx_l[jj]];
625 aux = w_data[p_perm];
628 w_data[p_perm] -= tl * data_l[jj];
629 partial_row_sum += tl * data_l[jj];
633 tl = data_l[jj] * w_data[jrow];
634 w_data[p_perm] -= tl;
636 partial_col_sum += tl;
639 if ((aux ==
zero) && (w_data[p_perm] !=
zero))
642 iw_l.insert (ridx_l[jj]);
644 iw_u.insert (p_perm);
650 if ((
std::abs (w_data[jrow]) < (droptol * cols_norm[
k]))
651 && (w_data[jrow] !=
zero))
654 total_sum += partial_col_sum;
656 total_sum += partial_row_sum;
675 if ((thresh >
zero) && (
k < (n - 1)))
679 for (it = iw_l.begin (); it != iw_l.end (); ++it)
682 if (
std::abs (w_data[p_perm]) > maximum)
684 maximum =
std::abs (w_data[p_perm]);
690 p_perm = iperm[max_ind];
691 if (max_ind != perm[
k])
694 if (w_data[
k] !=
zero)
695 iw_l.insert (perm[
k]);
700 iperm[perm[p_perm]] =
k;
701 iperm[perm[
k]] = p_perm;
703 perm[
k] = perm[p_perm];
705 w_data[
k] = w_data[p_perm];
706 w_data[p_perm] = aux;
714 while (it != iw_l.end ())
718 if (
std::abs (w_data[p_perm]) < (droptol * cols_norm[
k]))
721 total_sum += w_data[p_perm];
722 w_data[p_perm] =
zero;
736 if ((total_sum >
zero) && (w_data[
k] ==
zero))
738 w_data[
k] += total_sum;
744 if (w_data[
k] ==
zero)
747 error (
"ilu: encountered a pivot equal to 0");
755 for (it = iw_l.begin (); it != iw_l.end (); ++it)
758 w_data[p_perm] = w_data[p_perm] / w_data[
k];
761 if ((max_len_u - total_len_u) < n)
763 max_len_u += (0.1 * max_len_u) > n ? 0.1 * max_len_u : n;
764 data_out_u.resize (
dim_vector (max_len_u, 1));
765 data_u = data_out_u.fortran_vec ();
766 ridx_out_u.resize (
dim_vector (max_len_u, 1));
767 ridx_u = ridx_out_u.fortran_vec ();
770 if ((max_len_l - total_len_l) < n)
772 max_len_l += (0.1 * max_len_l) > n ? 0.1 * max_len_l : n;
773 data_out_l.resize (
dim_vector (max_len_l, 1));
774 data_l = data_out_l.fortran_vec ();
775 ridx_out_l.resize (
dim_vector (max_len_l, 1));
776 ridx_l = ridx_out_l.fortran_vec ();
781 for (it = iw_u.begin (); it != iw_u.end (); ++it)
783 if (w_data[*it] !=
zero)
785 data_u[total_len_u + w_len_u] = w_data[*it];
786 ridx_u[total_len_u + w_len_u] = *it;
794 for (it = iw_l.begin (); it != iw_l.end (); ++it)
797 if (w_data[p_perm] !=
zero)
799 data_l[total_len_l + w_len_l] = w_data[p_perm];
800 ridx_l[total_len_l + w_len_l] = *it;
805 total_len_u += w_len_u;
806 total_len_l += w_len_l;
810 if (total_len_l < 0 || total_len_u < 0)
811 error (
"ilu: Integer overflow. Too many fill-in elements in L or U");
814 uptr[
k] = total_len_u - 1;
816 cidx_u[
k+1] = cidx_u[
k] - cidx_u[0] + w_len_u;
817 cidx_l[
k+1] = cidx_l[
k] - cidx_l[0] + w_len_l;
823 octave_matrix_t *L_ptr;
824 octave_matrix_t *U_ptr;
825 octave_matrix_t diag (n, n, n);
834 L = octave_matrix_t (n, n, total_len_u - n);
835 U = octave_matrix_t (n, n, total_len_l);
841 L = octave_matrix_t (n, n, total_len_l);
842 U = octave_matrix_t (n, n, total_len_u);
845 for (
i = 0;
i <= n;
i++)
847 L_ptr->cidx (
i) = cidx_l[
i];
848 U_ptr->cidx (
i) = cidx_u[
i];
850 U_ptr->cidx (
i) -=
i;
853 for (
i = 0;
i < n;
i++)
856 diag.elem (
i,
i) = data_u[uptr[
i]];
859 while (j < cidx_l[
i+1])
861 L_ptr->ridx (j) = ridx_l[j];
862 L_ptr->data (j) = data_l[j];
867 while (j < cidx_u[
i+1])
883 U_ptr->data (
c) = data_u[j];
884 U_ptr->ridx (
c) = ridx_u[j];
905 if (args.length () != 5)
909 double droptol = args(1).double_value ();
910 double thresh = args(2).double_value ();
912 double udiag = args(4).double_value ();
916 if (! args(0).iscomplex ())
920 nnz_u = (Ftriu (
ovl (sm))(0).sparse_matrix_value ()).nnz ();
921 nnz_l = (
Ftril (
ovl (sm, -1))(0).sparse_matrix_value ()).nnz ();
929 ilu_tp <SparseMatrix, double> (sm, L, U, nnz_u, nnz_l,
931 perm, droptol, thresh, milu, udiag);
961 nnz_u = (Ftriu (
ovl (sm))(0).sparse_complex_matrix_value ()).nnz ();
962 nnz_l = (
Ftril (
ovl (sm, -1))(0).sparse_complex_matrix_value ()).nnz ();
970 ilu_tp <SparseComplexMatrix, Complex>
971 (sm, L, U, nnz_u, nnz_l, rc_norm.
fortran_vec (), perm,
OCTAVE_API ColumnVector xrownorms(const Matrix &m, double p)
static const idx_vector colon
OCTINTERP_API void print_usage(void)
OCTAVE_EXPORT octave_value_list or N dimensional array whose elements are all equal to the IEEE symbol zero divided by zero($0/0$)
const T * fortran_vec(void) const
#define DEFUN(name, args_name, nargout_name, doc)
Macro to define a builtin function.
void error(const char *fmt,...)
OCTAVE_EXPORT octave_value_list Ftril(const octave_value_list &args, int) and setting all other elements to zero. The optional second argument specifies how many diagonals above or below the main diagonal should also be set to zero. The default value of ar
nd example oindent opens the file binary numeric values will be read assuming they are stored in IEEE format with the least significant bit and then converted to the native representation Opening a file that is already open simply opens it again and returns a separate file id It is not an error to open a file several though writing to the same file through several different file ids may produce unexpected results The possible values of text mode reading and writing automatically converts linefeeds to the appropriate line end character for the you may append a you must also open the file in binary mode The parameter conversions are currently only supported for and permissions will be set to and then everything is written in a single operation This is very efficient and improves performance c
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)
OCTAVE_API RowVector xcolnorms(const Matrix &m, double p)
octave_idx_type cols(void) const
OCTAVE_EXPORT octave_value_list isa nd deftypefn *return ovl(args(0).isinteger())
void ilu_0(octave_matrix_t &sm, const std::string milu="off")
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
std::complex< double > Complex
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