24 #if defined (HAVE_CONFIG_H)
37 #if defined (HAVE_CXX_COMPLEX_SETTERS)
39 #elif defined (HAVE_CXX_COMPLEX_REFERENCE_ACCESSORS)
40 b.imag () = -b.imag ();
54 if (pivot.imag () != 0)
55 error (
"ichol: non-real pivot encountered. The matrix must be Hermitian.");
56 else if (pivot.real () < 0)
57 error (
"ichol: negative pivot encountered");
65 error (
"ichol: negative pivot encountered");
70 template <
typename octave_matrix_t,
typename T, T (*ichol_mult) (T, T),
71 bool (*ichol_checkpivot) (T)>
75 octave_idx_type j1, jend, j2, jrow, jjrow, j, jw,
i,
k, jj, r;
97 for (i = 0; i < n; i++)
106 for (k = 0; k < n; k++)
110 for (j = j1; j < j2; j++)
117 jjrow = Lfirst[jrow];
119 for (jj = jjrow; jj < jend; jj++)
123 tl = ichol_mult (data[jj], data[jjrow]);
136 if ((jjrow + 1) < jend)
141 Llist[j] = Llist[ridx[Lfirst[j]]];
142 Llist[ridx[Lfirst[j]]] = j;
149 data[j1] += dropsums[
k];
152 if (j1 == j2 || ridx[j1] != k)
153 error (
"ichol: encountered a pivot equal to 0");
155 if (! ichol_checkpivot (data[j1]))
158 data[cidx[
k]] = std::sqrt (data[j1]);
165 for (i = j1 + 1; i < j2; i++)
171 if ((Lfirst[k] + 1) < j2)
174 jjrow = ridx[Lfirst[
k]];
175 Llist[
k] = Llist[jjrow];
191 michol =
args(1).string_value ();
197 if (!
args(0).is_complex_type ())
201 sm =
feval (
"tril", arg_list)(0).sparse_matrix_value ();
211 sm =
feval (
"tril", arg_list)(0).sparse_complex_matrix_value ();
219 template <
typename octave_matrix_t,
typename T, T (*ichol_mult) (T, T),
220 bool (*ichol_checkpivot) (T)>
221 void ichol_t (
const octave_matrix_t& sm, octave_matrix_t& L,
const T* cols_norm,
238 T* data = sm.data ();
248 max_len += (0.1 * max_len) > n ? 0.1 * max_len : n;
261 std::vector <octave_idx_type> vec;
266 for (i = 0; i < n; i++)
276 for (k = 0; k < n; k++)
279 for (j = cidx[k]; j < cidx[k+1]; j++)
281 w_data[ridx[j]] = data[j];
291 jjrow = Lfirst[jrow];
292 jend = cidx_l[jrow+1];
293 for (jj = jjrow; jj < jend; jj++)
299 if (w_data[j] == zero)
304 w_data[j] -= ichol_mult (data_l[jj], data_l[jjrow]);
308 if ((jjrow + 1) < jend)
313 Llist[j] = Llist[ridx_l[Lfirst[j]]];
314 Llist[ridx_l[Lfirst[j]]] = j;
321 if ((max_len - total_len) < n)
323 max_len += (0.1 * max_len) > n ? 0.1 * max_len : n;
335 std::sort (vec.begin (), vec.begin () +
ind);
337 data_l[total_len] = w_data[
k];
338 ridx_l[total_len] =
k;
343 for (i = 0; i <
ind ; i++)
346 if (w_data[jrow] != zero)
348 if (
std::abs (w_data[jrow]) < (droptol * cols_norm[
k]))
352 col_drops[
k] += w_data[jrow];
353 col_drops[jrow] += w_data[jrow];
358 data_l[total_len + w_len] = w_data[jrow];
359 ridx_l[total_len + w_len] = jrow;
369 data_l[total_len] += col_drops[
k];
371 if (data_l[total_len] == zero)
372 error (
"ichol: encountered a pivot equal to 0");
374 if (! ichol_checkpivot (data_l[total_len]))
379 data_l[total_len] = std::sqrt (data_l[total_len]);
380 for (jj = total_len + 1; jj < (total_len + w_len); jj++)
381 data_l[jj] /= data_l[total_len];
386 error (
"ichol: integer overflow. Too many fill-in elements in L");
388 cidx_l[k+1] = cidx_l[
k] - cidx_l[0] + w_len;
393 Lfirst[
k] = cidx_l[
k];
394 if ((Lfirst[k] + 1) < cidx_l[k+1])
397 jjrow = ridx_l[Lfirst[
k]];
398 Llist[
k] = Llist[jjrow];
405 L = octave_matrix_t (n, n, total_len);
407 for (i = 0; i <= n; i++)
408 L.cidx (i) = cidx_l[
i];
410 for (i = 0; i < total_len; i++)
412 L.ridx (i) = ridx_l[
i];
413 L.data (i) = data_l[
i];
432 droptol =
args(1).double_value ();
435 michol =
args(2).string_value ();
438 if (!
args(0).is_complex_type ())
442 arg_list.
append (
args(0).sparse_matrix_value ());
444 feval (
"tril", arg_list)(0).sparse_matrix_value ();
447 arg_list(2) =
"cols";
448 cols_norm =
feval (
"norm", arg_list)(0).vector_value ();
452 (sm_l, L, cols_norm.
fortran_vec (), droptol, michol);
460 arg_list.
append (
args(0).sparse_complex_matrix_value ());
462 feval (
"tril", arg_list)(0).sparse_complex_matrix_value ();
465 arg_list(2) =
"cols";
466 cols_norm =
feval (
"norm", arg_list)(0).complex_vector_value ();
Complex ichol_mult_complex(Complex a, Complex b)
OCTAVE_EXPORT octave_value_list isa nd deftypefn *return ovl(args(0).is_integer_type())
octave_value_list & append(const octave_value &val)
#define DEFUN(name, args_name, nargout_name, doc)
void error(const char *fmt,...)
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")
void resize(const dim_vector &dv, const T &rfv)
bool ichol_checkpivot_complex(Complex pivot)
=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))<
double ichol_mult_real(double a, double b)
issues an error eealso double
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
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)
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