37 #if defined (HAVE_CXX_COMPLEX_SETTERS)
39 #elif defined (HAVE_CXX_COMPLEX_REFERENCE_ACCESSORS)
54 if (pivot.imag () != 0)
56 error (
"ichol: non-real pivot encountered. The matrix must be hermitian.");
59 else if (pivot.real () < 0)
61 error (
"ichol: negative pivot encountered");
71 error (
"ichol: negative pivot encountered");
77 template <
typename octave_matrix_t,
typename T, T (*ichol_mult) (T, T),
78 bool (*ichol_checkpivot) (T)>
79 void ichol_0 (octave_matrix_t& sm,
const std::string michol =
"off")
83 octave_idx_type j1, jend, j2, jrow, jjrow, j, jw, i, k, jj, r;
104 for (i = 0; i < n; i++)
113 for (k = 0; k < n; k++)
117 for (j = j1; j < j2; j++)
124 jjrow = Lfirst[jrow];
126 for (jj = jjrow; jj < jend; jj++)
130 tl = ichol_mult (data[jj], data[jjrow]);
143 if ((jjrow + 1) < jend)
148 Llist[j] = Llist[ridx[Lfirst[j]]];
149 Llist[ridx[Lfirst[j]]] = j;
156 data[j1] += dropsums[k];
160 error (
"ichol: encountered a pivot equal to 0");
164 if (! ichol_checkpivot (data[j1]))
167 data[cidx[k]] = std::sqrt (data[j1]);
174 for (i = j1 + 1; i < j2; i++)
180 if ((Lfirst[k] + 1) < j2)
183 jjrow = ridx[Lfirst[k]];
184 Llist[k] = Llist[jjrow];
191 DEFUN (__ichol0__, args, nargout,
193 @deftypefn {Built-in Function} {@var{L} =} __ichol0__ (@var{A})\n\
194 @deftypefnx {Built-in Function} {@var{L} =} __ichol0__ (@var{A}, @var{michol})\n\
195 Undocumented internal function.\n\
201 int nargin = args.
length ();
202 std::string michol =
"off";
204 if (nargout > 1 || nargin < 1 || nargin > 2)
211 michol = args(1).string_value ();
217 if (!args(0).is_complex_type ())
221 sm =
feval (
"tril", param_list)(0).sparse_matrix_value ();
231 sm =
feval (
"tril", param_list)(0).sparse_complex_matrix_value ();
241 template <
typename octave_matrix_t,
typename T, T (*ichol_mult) (T, T),
242 bool (*ichol_checkpivot) (T)>
243 void ichol_t (
const octave_matrix_t& sm, octave_matrix_t& L,
const T* cols_norm,
244 const T droptol,
const std::string michol =
"off")
261 T* data = sm.data ();
271 max_len += (0.1 * max_len) > n ? 0.1 * max_len : n;
284 std::vector <octave_idx_type> vec;
289 for (i = 0; i < n; i++)
299 for (k = 0; k < n; k++)
302 for (j = cidx[k]; j < cidx[k+1]; j++)
304 w_data[ridx[j]] = data[j];
314 jjrow = Lfirst[jrow];
315 jend = cidx_l[jrow+1];
316 for (jj = jjrow; jj < jend; jj++)
322 if (w_data[j] == zero)
327 w_data[j] -= ichol_mult (data_l[jj], data_l[jjrow]);
331 if ((jjrow + 1) < jend)
336 Llist[j] = Llist[ridx_l[Lfirst[j]]];
337 Llist[ridx_l[Lfirst[j]]] = j;
344 if ((max_len - total_len) < n)
346 max_len += (0.1 * max_len) > n ? 0.1 * max_len : n;
358 std::sort (vec.begin (), vec.begin () + ind);
360 data_l[total_len] = w_data[k];
361 ridx_l[total_len] = k;
366 for (i = 0; i < ind ; i++)
369 if (w_data[jrow] != zero)
371 if (
std::abs (w_data[jrow]) < (droptol * cols_norm[k]))
375 col_drops[k] += w_data[jrow];
376 col_drops[jrow] += w_data[jrow];
381 data_l[total_len + w_len] = w_data[jrow];
382 ridx_l[total_len + w_len] = jrow;
392 data_l[total_len] += col_drops[k];
394 if (data_l[total_len] == zero)
396 error (
"ichol: encountered a pivot equal to 0");
399 else if (! ichol_checkpivot (data_l[total_len]))
404 data_l[total_len] = std::sqrt (data_l[total_len]);
405 for (jj = total_len + 1; jj < (total_len + w_len); jj++)
406 data_l[jj] /= data_l[total_len];
412 error (
"ichol: integer overflow. Too many fill-in elements in L");
415 cidx_l[k+1] = cidx_l[k] - cidx_l[0] + w_len;
420 Lfirst[k] = cidx_l[k];
421 if ((Lfirst[k] + 1) < cidx_l[k+1])
424 jjrow = ridx_l[Lfirst[k]];
425 Llist[k] = Llist[jjrow];
434 L = octave_matrix_t (n, n, total_len);
435 for (i = 0; i <= n; i++)
436 L.cidx (i) = cidx_l[i];
437 for (i = 0; i < total_len; i++)
439 L.ridx (i) = ridx_l[i];
440 L.data (i) = data_l[i];
445 DEFUN (__icholt__, args, nargout,
447 @deftypefn {Built-in Function} {@var{L} =} __icholt__ (@var{A})\n\
448 @deftypefnx {Built-in Function} {@var{L} =} __icholt__ (@var{A}, @var{droptol})\n\
449 @deftypefnx {Built-in Function} {@var{L} =} __icholt__ (@var{A}, @var{droptol}, @var{michol})\n\
450 Undocumented internal function.\n\
454 int nargin = args.
length ();
456 std::string michol =
"off";
459 if (nargout > 1 || nargin < 1 || nargin > 3)
468 droptol = args(1).double_value ();
471 michol = args(2).string_value ();
474 if (! args(0).is_complex_type ())
478 param_list.
append (args(0).sparse_matrix_value ());
480 feval (
"tril", param_list)(0).sparse_matrix_value ();
481 param_list(0) = sm_l;
483 param_list(2) =
"cols";
484 cols_norm =
feval (
"norm", param_list)(0).vector_value ();
488 (sm_l, L, cols_norm.
fortran_vec (), droptol, michol);
496 param_list.
append (args(0).sparse_complex_matrix_value ());
498 feval (
"tril", param_list)(0).sparse_complex_matrix_value ();
499 param_list(0) = sm_l;
501 param_list(2) =
"cols";
502 cols_norm =
feval (
"norm", param_list)(0).complex_vector_value ();
Complex ichol_mult_complex(Complex a, Complex b)
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,...)
octave_value_list feval(const std::string &name, const octave_value_list &args, int nargout)
ComplexColumnVector conj(const ComplexColumnVector &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")
void resize(const dim_vector &dv, const T &rfv)
bool ichol_checkpivot_complex(Complex pivot)
double ichol_mult_real(double a, double b)
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
ColumnVector imag(const ComplexColumnVector &a)
void ichol_0(octave_matrix_t &sm, const std::string michol="off")
std::complex< double > Complex
const T * fortran_vec(void) const
bool ichol_checkpivot_real(double pivot)