GNU Octave  4.2.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
__ichol__.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2014-2017 Eduardo Ramos Fern├índez <eduradical951@gmail.com>
4 Copyright (C) 2013-2016 Kai T. Ohlhus <k.ohlhus@gmail.com>
5 
6 This file is part of Octave.
7 
8 Octave is free software; you can redistribute it and/or modify it
9 under the terms of the GNU General Public License as published by the
10 Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12 
13 Octave is distributed in the hope that it will be useful, but WITHOUT
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
16 for more details.
17 
18 You should have received a copy of the GNU General Public License
19 along with Octave; see the file COPYING. If not, see
20 <http://www.gnu.org/licenses/>.
21 
22 */
23 
24 #if defined (HAVE_CONFIG_H)
25 # include "config.h"
26 #endif
27 
28 #include "oct-locbuf.h"
29 
30 #include "defun.h"
31 #include "error.h"
32 #include "parse.h"
33 
34 // Secondary functions for complex and real case used in ichol algorithms.
36 {
37 #if defined (HAVE_CXX_COMPLEX_SETTERS)
38  b.imag (-b.imag ());
39 #elif defined (HAVE_CXX_COMPLEX_REFERENCE_ACCESSORS)
40  b.imag () = -b.imag ();
41 #else
42  b = b.conj ();
43 #endif
44  return a * b;
45 }
46 
47 double ichol_mult_real (double a, double b)
48 {
49  return a * b;
50 }
51 
53 {
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");
58 
59  return true;
60 }
61 
62 bool ichol_checkpivot_real (double pivot)
63 {
64  if (pivot < 0)
65  error ("ichol: negative pivot encountered");
66 
67  return true;
68 }
69 
70 template <typename octave_matrix_t, typename T, T (*ichol_mult) (T, T),
71  bool (*ichol_checkpivot) (T)>
72 void ichol_0 (octave_matrix_t& sm, const std::string michol = "off")
73 {
74  const octave_idx_type n = sm.cols ();
75  octave_idx_type j1, jend, j2, jrow, jjrow, j, jw, i, k, jj, r;
76  T tl;
77 
78  char opt;
79  enum {OFF, ON};
80  if (michol == "on")
81  opt = ON;
82  else
83  opt = OFF;
84 
85  // Input matrix pointers
86  octave_idx_type* cidx = sm.cidx ();
87  octave_idx_type* ridx = sm.ridx ();
88  T* data = sm.data ();
89 
90  // Working arrays
94  OCTAVE_LOCAL_BUFFER (T, dropsums, n);
95 
96  // Initialize working arrays
97  for (i = 0; i < n; i++)
98  {
99  iw[i] = -1;
100  Llist[i] = -1;
101  Lfirst[i] = -1;
102  dropsums[i] = 0;
103  }
104 
105  // Loop over all columns
106  for (k = 0; k < n; k++)
107  {
108  j1 = cidx[k];
109  j2 = cidx[k+1];
110  for (j = j1; j < j2; j++)
111  iw[ridx[j]] = j;
112 
113  jrow = Llist[k];
114  // Iterate over each non-zero element in the actual row.
115  while (jrow != -1)
116  {
117  jjrow = Lfirst[jrow];
118  jend = cidx[jrow+1];
119  for (jj = jjrow; jj < jend; jj++)
120  {
121  r = ridx[jj];
122  jw = iw[r];
123  tl = ichol_mult (data[jj], data[jjrow]);
124  if (jw != -1)
125  data[jw] -= tl;
126  else
127  // Because of the symmetry of the matrix, we know
128  // the drops in the column r are also in the column k.
129  if (opt == ON)
130  {
131  dropsums[r] -= tl;
132  dropsums[k] -= tl;
133  }
134  }
135  // Update the linked list and the first entry of the actual column.
136  if ((jjrow + 1) < jend)
137  {
138  Lfirst[jrow]++;
139  j = jrow;
140  jrow = Llist[jrow];
141  Llist[j] = Llist[ridx[Lfirst[j]]];
142  Llist[ridx[Lfirst[j]]] = j;
143  }
144  else
145  jrow = Llist[jrow];
146  }
147 
148  if (opt == ON)
149  data[j1] += dropsums[k];
150 
151  // Test for j1 == j2 must be first to avoid invalid ridx[j1] access
152  if (j1 == j2 || ridx[j1] != k)
153  error ("ichol: encountered a pivot equal to 0");
154 
155  if (! ichol_checkpivot (data[j1]))
156  break;
157 
158  data[cidx[k]] = std::sqrt (data[j1]);
159 
160  // Update Llist and Lfirst with the k-column information. Also,
161  // scale the column elements by the pivot and reset the working array iw.
162  if (k < (n - 1))
163  {
164  iw[ridx[j1]] = -1;
165  for (i = j1 + 1; i < j2; i++)
166  {
167  iw[ridx[i]] = -1;
168  data[i] /= data[j1];
169  }
170  Lfirst[k] = j1;
171  if ((Lfirst[k] + 1) < j2)
172  {
173  Lfirst[k]++;
174  jjrow = ridx[Lfirst[k]];
175  Llist[k] = Llist[jjrow];
176  Llist[jjrow] = k;
177  }
178  }
179  }
180 }
181 
182 DEFUN (__ichol0__, args, ,
183  doc: /* -*- texinfo -*-
184 @deftypefn {} {@var{L} =} __ichol0__ (@var{A})
185 @deftypefnx {} {@var{L} =} __ichol0__ (@var{A}, @var{michol})
186 Undocumented internal function.
187 @end deftypefn */)
188 {
189  std::string michol = "off";
190  if (args.length ())
191  michol = args(1).string_value ();
192 
193  // In ICHOL0 algorithm the zero-pattern of the input matrix is preserved
194  // so its structure does not change during the algorithm. The same input
195  // matrix is used to build the output matrix due to that fact.
196  octave_value_list arg_list;
197  if (! args(0).is_complex_type ())
198  {
199  SparseMatrix sm = args(0).sparse_matrix_value ();
200  arg_list.append (sm);
201  sm = feval ("tril", arg_list)(0).sparse_matrix_value ();
203  ichol_checkpivot_real> (sm, michol);
204 
205  return ovl (sm);
206  }
207  else
208  {
209  SparseComplexMatrix sm = args(0).sparse_complex_matrix_value ();
210  arg_list.append (sm);
211  sm = feval ("tril", arg_list)(0).sparse_complex_matrix_value ();
213  ichol_checkpivot_complex> (sm, michol);
214 
215  return ovl (sm);
216  }
217 }
218 
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,
222  const T droptol, const std::string michol = "off")
223 {
224 
225  const octave_idx_type n = sm.cols ();
226  octave_idx_type j, jrow, jend, jjrow, i, k, jj, total_len,
227  w_len, max_len, ind;
228  char opt;
229  enum {OFF, ON};
230  if (michol == "on")
231  opt = ON;
232  else
233  opt = OFF;
234 
235  // Input matrix pointers
236  octave_idx_type* cidx = sm.cidx ();
237  octave_idx_type* ridx = sm.ridx ();
238  T* data = sm.data ();
239 
240  // Output matrix data structures. Because the final zero pattern pattern of
241  // the output matrix is not known due to fill-in elements, a heuristic
242  // approach has been adopted for memory allocation. The size of ridx_out_l
243  // and data_out_l is incremented 10% of their actual size (nnz (A) in the
244  // beginning). If that amount is less than n, their size is just incremented
245  // in n elements. This way the number of reallocations decreases throughout
246  // the process, obtaining a good performance.
247  max_len = sm.nnz ();
248  max_len += (0.1 * max_len) > n ? 0.1 * max_len : n;
249  Array <octave_idx_type> cidx_out_l (dim_vector (n + 1, 1));
250  octave_idx_type* cidx_l = cidx_out_l.fortran_vec ();
251  Array <octave_idx_type> ridx_out_l (dim_vector (max_len ,1));
252  octave_idx_type* ridx_l = ridx_out_l.fortran_vec ();
253  Array <T> data_out_l (dim_vector (max_len, 1));
254  T* data_l = data_out_l.fortran_vec ();
255 
256  // Working arrays
257  OCTAVE_LOCAL_BUFFER (T, w_data, n);
260  OCTAVE_LOCAL_BUFFER (T, col_drops, n);
261  std::vector <octave_idx_type> vec;
262  vec.resize (n);
263 
264  T zero = T (0);
265  cidx_l[0] = cidx[0];
266  for (i = 0; i < n; i++)
267  {
268  Llist[i] = -1;
269  Lfirst[i] = -1;
270  w_data[i] = 0;
271  col_drops[i] = zero;
272  vec[i] = 0;
273  }
274 
275  total_len = 0;
276  for (k = 0; k < n; k++)
277  {
278  ind = 0;
279  for (j = cidx[k]; j < cidx[k+1]; j++)
280  {
281  w_data[ridx[j]] = data[j];
282  if (ridx[j] != k)
283  {
284  vec[ind] = ridx[j];
285  ind++;
286  }
287  }
288  jrow = Llist[k];
289  while (jrow != -1)
290  {
291  jjrow = Lfirst[jrow];
292  jend = cidx_l[jrow+1];
293  for (jj = jjrow; jj < jend; jj++)
294  {
295  j = ridx_l[jj];
296  // If the element in the j position of the row is zero,
297  // then it will become non-zero, so we add it to the
298  // vector that tracks non-zero elements in the working row.
299  if (w_data[j] == zero)
300  {
301  vec[ind] = j;
302  ind++;
303  }
304  w_data[j] -= ichol_mult (data_l[jj], data_l[jjrow]);
305  }
306  // Update the actual column first element and
307  // update the linked list of the jrow row.
308  if ((jjrow + 1) < jend)
309  {
310  Lfirst[jrow]++;
311  j = jrow;
312  jrow = Llist[jrow];
313  Llist[j] = Llist[ridx_l[Lfirst[j]]];
314  Llist[ridx_l[Lfirst[j]]] = j;
315  }
316  else
317  jrow = Llist[jrow];
318  }
319 
320  // Resizing output arrays
321  if ((max_len - total_len) < n)
322  {
323  max_len += (0.1 * max_len) > n ? 0.1 * max_len : n;
324  data_out_l.resize (dim_vector (max_len, 1));
325  data_l = data_out_l.fortran_vec ();
326  ridx_out_l.resize (dim_vector (max_len, 1));
327  ridx_l = ridx_out_l.fortran_vec ();
328  }
329 
330  // The sorting of the non-zero elements of the working column can be
331  // handled in a couple of ways. The most efficient two I found, are
332  // keeping the elements in an ordered binary search tree dynamically or
333  // keep them unsorted in a vector and at the end of the outer iteration
334  // order them. The last approach exhibits lower execution times.
335  std::sort (vec.begin (), vec.begin () + ind);
336 
337  data_l[total_len] = w_data[k];
338  ridx_l[total_len] = k;
339  w_len = 1;
340 
341  // Extract the non-zero elements of working column and
342  // drop the elements that are lower than droptol * cols_norm[k].
343  for (i = 0; i < ind ; i++)
344  {
345  jrow = vec[i];
346  if (w_data[jrow] != zero)
347  {
348  if (std::abs (w_data[jrow]) < (droptol * cols_norm[k]))
349  {
350  if (opt == ON)
351  {
352  col_drops[k] += w_data[jrow];
353  col_drops[jrow] += w_data[jrow];
354  }
355  }
356  else
357  {
358  data_l[total_len + w_len] = w_data[jrow];
359  ridx_l[total_len + w_len] = jrow;
360  w_len++;
361  }
362  vec[i] = 0;
363  }
364  w_data[jrow] = zero;
365  }
366 
367  // Compensate column sums --> michol option
368  if (opt == ON)
369  data_l[total_len] += col_drops[k];
370 
371  if (data_l[total_len] == zero)
372  error ("ichol: encountered a pivot equal to 0");
373 
374  if (! ichol_checkpivot (data_l[total_len]))
375  break;
376 
377  // Once elements are dropped and compensation of column sums are done,
378  // scale the elements by the pivot.
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];
382  total_len += w_len;
383  // Check if there are too many elements to be indexed with
384  // octave_idx_type type due to fill-in during the process.
385  if (total_len < 0)
386  error ("ichol: integer overflow. Too many fill-in elements in L");
387 
388  cidx_l[k+1] = cidx_l[k] - cidx_l[0] + w_len;
389 
390  // Update Llist and Lfirst with the k-column information.
391  if (k < (n - 1))
392  {
393  Lfirst[k] = cidx_l[k];
394  if ((Lfirst[k] + 1) < cidx_l[k+1])
395  {
396  Lfirst[k]++;
397  jjrow = ridx_l[Lfirst[k]];
398  Llist[k] = Llist[jjrow];
399  Llist[jjrow] = k;
400  }
401  }
402  }
403 
404  // Build the output matrices
405  L = octave_matrix_t (n, n, total_len);
406 
407  for (i = 0; i <= n; i++)
408  L.cidx (i) = cidx_l[i];
409 
410  for (i = 0; i < total_len; i++)
411  {
412  L.ridx (i) = ridx_l[i];
413  L.data (i) = data_l[i];
414  }
415 }
416 
417 DEFUN (__icholt__, args, ,
418  doc: /* -*- texinfo -*-
419 @deftypefn {} {@var{L} =} __icholt__ (@var{A})
420 @deftypefnx {} {@var{L} =} __icholt__ (@var{A}, @var{droptol})
421 @deftypefnx {} {@var{L} =} __icholt__ (@var{A}, @var{droptol}, @var{michol})
422 Undocumented internal function.
423 @end deftypefn */)
424 {
425  int nargin = args.length ();
426  // Default values of parameters
427  std::string michol = "off";
428  double droptol = 0;
429 
430  // Don't repeat input validation of arguments done in ichol.m
431  if (nargin >= 2)
432  droptol = args(1).double_value ();
433 
434  if (nargin == 3)
435  michol = args(2).string_value ();
436 
437  octave_value_list arg_list;
438  if (! args(0).is_complex_type ())
439  {
440  Array <double> cols_norm;
441  SparseMatrix L;
442  arg_list.append (args(0).sparse_matrix_value ());
443  SparseMatrix sm_l =
444  feval ("tril", arg_list)(0).sparse_matrix_value ();
445  arg_list(0) = sm_l;
446  arg_list(1) = 1;
447  arg_list(2) = "cols";
448  cols_norm = feval ("norm", arg_list)(0).vector_value ();
449  arg_list.clear ();
452  (sm_l, L, cols_norm.fortran_vec (), droptol, michol);
453 
454  return ovl (L);
455  }
456  else
457  {
458  Array <Complex> cols_norm;
460  arg_list.append (args(0).sparse_complex_matrix_value ());
461  SparseComplexMatrix sm_l =
462  feval ("tril", arg_list)(0).sparse_complex_matrix_value ();
463  arg_list(0) = sm_l;
464  arg_list(1) = 1;
465  arg_list(2) = "cols";
466  cols_norm = feval ("norm", arg_list)(0).complex_vector_value ();
467  arg_list.clear ();
470  (sm_l, L, cols_norm.fortran_vec (),
471  Complex (droptol), michol);
472 
473  return ovl (L);
474  }
475 }
476 
477 /*
478 ## No test needed for internal helper function.
479 %!assert (1)
480 */
void clear(void)
Definition: ovl.h:152
Complex ichol_mult_complex(Complex a, Complex b)
Definition: __ichol__.cc:35
ind
Definition: sub2ind.cc:107
OCTAVE_EXPORT octave_value_list isa nd deftypefn *return ovl(args(0).is_integer_type())
octave_value_list & append(const octave_value &val)
Definition: ovl.cc:85
for large enough k
Definition: lu.cc:606
#define DEFUN(name, args_name, nargout_name, doc)
Definition: defun.h:46
void error(const char *fmt,...)
Definition: error.cc:570
is greater than zero
Definition: load-path.cc:2339
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
Definition: cellfun.cc:398
JNIEnv void * args
Definition: ov-java.cc:67
void ichol_t(const octave_matrix_t &sm, octave_matrix_t &L, const T *cols_norm, const T droptol, const std::string michol="off")
Definition: __ichol__.cc:221
int nargin
Definition: graphics.cc:10115
void resize(const dim_vector &dv, const T &rfv)
Definition: Array.cc:1028
feval(ar{f}, 1) esult
Definition: oct-parse.cc:8829
bool ichol_checkpivot_complex(Complex pivot)
Definition: __ichol__.cc:52
=val(i)}if ode{val(i)}occurs in table i
Definition: lookup.cc:239
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)
Definition: __ichol__.cc:47
issues an error eealso double
Definition: ov-bool-mat.cc:594
b
Definition: cellfun.cc:398
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition: oct-locbuf.h:200
void ichol_0(octave_matrix_t &sm, const std::string michol="off")
Definition: __ichol__.cc:72
std::complex< double > Complex
Definition: oct-cmplx.h:31
const T * fortran_vec(void) const
Definition: Array.h:584
bool ichol_checkpivot_real(double pivot)
Definition: __ichol__.cc:62
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:87
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
Definition: utils.cc:854