24 #if defined (HAVE_CONFIG_H) 34 #include "builtin-defun-decls.h" 39 #if defined (HAVE_CXX_COMPLEX_SETTERS) 41 #elif defined (HAVE_CXX_COMPLEX_REFERENCE_ACCESSORS) 42 b.imag () = -
b.imag ();
56 if (pivot.imag () != 0)
57 error (
"ichol: non-real pivot encountered. The matrix must be Hermitian.");
58 else if (pivot.real () < 0)
59 error (
"ichol: negative pivot encountered");
67 error (
"ichol: negative pivot encountered");
72 template <
typename octave_matrix_t,
typename T, T (*ichol_mult) (T, T),
73 bool (*ichol_checkpivot) (T)>
77 octave_idx_type j1, jend, j2, jrow, jjrow, j, jw,
i,
k, jj, r;
99 for (
i = 0;
i < n;
i++)
108 for (
k = 0;
k < n;
k++)
112 for (j = j1; j < j2; j++)
119 jjrow = Lfirst[jrow];
121 for (jj = jjrow; jj < jend; jj++)
125 tl = ichol_mult (data[jj], data[jjrow]);
138 if ((jjrow + 1) < jend)
143 Llist[j] = Llist[ridx[Lfirst[j]]];
144 Llist[ridx[Lfirst[j]]] = j;
151 data[j1] += dropsums[
k];
154 if (j1 == j2 || ridx[j1] !=
k)
155 error (
"ichol: encountered a pivot equal to 0");
157 if (! ichol_checkpivot (data[j1]))
160 data[cidx[
k]] = std::sqrt (data[j1]);
167 for (
i = j1 + 1;
i < j2;
i++)
173 if ((Lfirst[
k] + 1) < j2)
176 jjrow = ridx[Lfirst[
k]];
177 Llist[
k] = Llist[jjrow];
184 DEFUN (__ichol0__, args, ,
190 if (args.length () != 2)
198 if (! args(0).iscomplex ())
208 Ftril (args(0))(0).sparse_complex_matrix_value ();
215 template <
typename octave_matrix_t,
typename T, T (*ichol_mult) (T, T),
216 bool (*ichol_checkpivot) (T)>
217 void ichol_t (
const octave_matrix_t& sm, octave_matrix_t& L,
const T* cols_norm,
234 T* data = sm.data ();
244 max_len += (0.1 * max_len) > n ? 0.1 * max_len : n;
257 std::vector<octave_idx_type> vec (n, 0);
258 std::vector<bool> mark (n,
false);
262 for (
i = 0;
i < n;
i++)
271 for (
k = 0;
k < n;
k++)
274 for (j = cidx[
k]; j < cidx[
k+1]; j++)
276 w_data[ridx[j]] = data[j];
279 mark[ridx[j]] =
true;
289 jjrow = Lfirst[jrow];
290 jend = cidx_l[jrow+1];
291 for (jj = jjrow; jj < jend; jj++)
303 w_data[j] -= ichol_mult (data_l[jj], data_l[jjrow]);
307 if ((jjrow + 1) < jend)
312 Llist[j] = Llist[ridx_l[Lfirst[j]]];
313 Llist[ridx_l[Lfirst[j]]] = j;
320 if ((max_len - total_len) < n)
322 max_len += (0.1 * max_len) > n ? 0.1 * max_len : n;
334 std::sort (vec.begin (), vec.begin () +
ind);
336 data_l[total_len] = w_data[
k];
337 ridx_l[total_len] =
k;
342 for (
i = 0;
i <
ind ;
i++)
345 if (w_data[jrow] !=
zero)
347 if (
std::abs (w_data[jrow]) < (droptol * cols_norm[
k]))
351 col_drops[
k] += w_data[jrow];
352 col_drops[jrow] += w_data[jrow];
357 data_l[total_len + w_len] = w_data[jrow];
358 ridx_l[total_len + w_len] = jrow;
370 data_l[total_len] += col_drops[
k];
372 if (data_l[total_len] ==
zero)
373 error (
"ichol: encountered a pivot equal to 0");
375 if (! ichol_checkpivot (data_l[total_len]))
380 data_l[total_len] = std::sqrt (data_l[total_len]);
381 for (jj = total_len + 1; jj < (total_len + w_len); jj++)
382 data_l[jj] /= data_l[total_len];
387 error (
"ichol: integer overflow. Too many fill-in elements in L");
389 cidx_l[
k+1] = cidx_l[
k] - cidx_l[0] + w_len;
394 Lfirst[
k] = cidx_l[
k];
395 if ((Lfirst[
k] + 1) < cidx_l[
k+1])
398 jjrow = ridx_l[Lfirst[
k]];
399 Llist[
k] = Llist[jjrow];
406 L = octave_matrix_t (n, n, total_len);
408 for (
i = 0;
i <= n;
i++)
409 L.cidx (
i) = cidx_l[
i];
411 for (
i = 0;
i < total_len;
i++)
413 L.ridx (
i) = ridx_l[
i];
414 L.data (
i) = data_l[
i];
418 DEFUN (__icholt__, args, ,
424 if (args.length () != 3)
427 double droptol = args(1).double_value ();
430 if (! args(0).iscomplex ())
444 Ftril (args(0))(0).sparse_complex_matrix_value ();
Complex ichol_mult_complex(Complex a, Complex b)
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
calling an anonymous function involves an overhead quite comparable to the overhead of an m file function Passing a handle to a built in function is because the interpreter is not involved in the internal loop For a
void ichol_t(const octave_matrix_t &sm, octave_matrix_t &L, const T *cols_norm, const T droptol, const std::string michol="off")
OCTAVE_API RowVector xcolnorms(const Matrix &m, double p)
void resize(const dim_vector &dv, const T &rfv)
Resizing (with fill).
OCTAVE_EXPORT octave_value_list isa nd deftypefn *return ovl(args(0).isinteger())
bool ichol_checkpivot_complex(Complex pivot)
double ichol_mult_real(double a, double b)
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
void ichol_0(octave_matrix_t &sm, const std::string michol="off")
std::complex< double > Complex
bool ichol_checkpivot_real(double pivot)
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