GNU Octave  4.0.0
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-2015 Eduardo Ramos Fernández <eduradical951@gmail.com>
4 Copyright (C) 2013-2015 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 #ifdef 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 (-std::imag (b));
39 #elif defined (HAVE_CXX_COMPLEX_REFERENCE_ACCESSORS)
40  b.imag () = -std::imag (b);
41 #else
42  b = std::conj (b);
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  {
56  error ("ichol: non-real pivot encountered. The matrix must be hermitian.");
57  return false;
58  }
59  else if (pivot.real () < 0)
60  {
61  error ("ichol: negative pivot encountered");
62  return false;
63  }
64  return true;
65 }
66 
67 bool ichol_checkpivot_real (double pivot)
68 {
69  if (pivot < 0)
70  {
71  error ("ichol: negative pivot encountered");
72  return false;
73  }
74  return true;
75 }
76 
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")
80 {
81 
82  const octave_idx_type n = sm.cols ();
83  octave_idx_type j1, jend, j2, jrow, jjrow, j, jw, i, k, jj, r;
84  T tl;
85  char opt;
86  enum {OFF, ON};
87  if (michol == "on")
88  opt = ON;
89  else
90  opt = OFF;
91 
92  // Input matrix pointers
93  octave_idx_type* cidx = sm.cidx ();
94  octave_idx_type* ridx = sm.ridx ();
95  T* data = sm.data ();
96 
97  // Working arrays
101  OCTAVE_LOCAL_BUFFER (T, dropsums, n);
102 
103  // Initialize working arrays
104  for (i = 0; i < n; i++)
105  {
106  iw[i] = -1;
107  Llist[i] = -1;
108  Lfirst[i] = -1;
109  dropsums[i] = 0;
110  }
111 
112  // Main loop
113  for (k = 0; k < n; k++)
114  {
115  j1 = cidx[k];
116  j2 = cidx[k+1];
117  for (j = j1; j < j2; j++)
118  iw[ridx[j]] = j;
119 
120  jrow = Llist [k];
121  // Iterate over each non-zero element in the actual row.
122  while (jrow != -1)
123  {
124  jjrow = Lfirst[jrow];
125  jend = cidx[jrow+1];
126  for (jj = jjrow; jj < jend; jj++)
127  {
128  r = ridx[jj];
129  jw = iw[r];
130  tl = ichol_mult (data[jj], data[jjrow]);
131  if (jw != -1)
132  data[jw] -= tl;
133  else
134  // Because of the symmetry of the matrix, we know
135  // the drops in the column r are also in the column k.
136  if (opt == ON)
137  {
138  dropsums[r] -= tl;
139  dropsums[k] -= tl;
140  }
141  }
142  // Update the linked list and the first entry of the actual column.
143  if ((jjrow + 1) < jend)
144  {
145  Lfirst[jrow]++;
146  j = jrow;
147  jrow = Llist[jrow];
148  Llist[j] = Llist[ridx[Lfirst[j]]];
149  Llist[ridx[Lfirst[j]]] = j;
150  }
151  else
152  jrow = Llist[jrow];
153  }
154 
155  if (opt == ON)
156  data[j1] += dropsums[k];
157 
158  if (ridx[j1] != k)
159  {
160  error ("ichol: encountered a pivot equal to 0");
161  break;
162  }
163 
164  if (! ichol_checkpivot (data[j1]))
165  break;
166 
167  data[cidx[k]] = std::sqrt (data[j1]);
168 
169  // Update Llist and Lfirst with the k-column information. Also,
170  // scale the column elements by the pivot and reset the working array iw.
171  if (k < (n - 1))
172  {
173  iw[ridx[j1]] = -1;
174  for (i = j1 + 1; i < j2; i++)
175  {
176  iw[ridx[i]] = -1;
177  data[i] /= data[j1];
178  }
179  Lfirst[k] = j1;
180  if ((Lfirst[k] + 1) < j2)
181  {
182  Lfirst[k]++;
183  jjrow = ridx[Lfirst[k]];
184  Llist[k] = Llist[jjrow];
185  Llist[jjrow] = k;
186  }
187  }
188  }
189 }
190 
191 DEFUN (__ichol0__, args, nargout,
192  "-*- texinfo -*-\n\
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\
196 @end deftypefn")
197 
198 {
199  octave_value_list retval;
200 
201  int nargin = args.length ();
202  std::string michol = "off";
203 
204  if (nargout > 1 || nargin < 1 || nargin > 2)
205  {
206  print_usage ();
207  return retval;
208  }
209 
210  if (nargin == 2)
211  michol = args(1).string_value ();
212 
213  // In ICHOL0 algorithm the zero-pattern of the input matrix is preserved
214  // so it's structure does not change during the algorithm. The same input
215  // matrix is used to build the output matrix due to that fact.
216  octave_value_list param_list;
217  if (!args(0).is_complex_type ())
218  {
219  SparseMatrix sm = args(0).sparse_matrix_value ();
220  param_list.append (sm);
221  sm = feval ("tril", param_list)(0).sparse_matrix_value ();
223  ichol_checkpivot_real> (sm, michol);
224  if (! error_state)
225  retval(0) = sm;
226  }
227  else
228  {
229  SparseComplexMatrix sm = args(0).sparse_complex_matrix_value ();
230  param_list.append (sm);
231  sm = feval ("tril", param_list)(0).sparse_complex_matrix_value ();
233  ichol_checkpivot_complex> (sm, michol);
234  if (! error_state)
235  retval(0) = sm;
236  }
237 
238  return retval;
239 }
240 
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")
245 
246 {
247 
248  const octave_idx_type n = sm.cols ();
249  octave_idx_type j, jrow, jend, jjrow, i, k, jj, total_len,
250  w_len, max_len, ind;
251  char opt;
252  enum {OFF, ON};
253  if (michol == "on")
254  opt = ON;
255  else
256  opt = OFF;
257 
258  // Input matrix pointers
259  octave_idx_type* cidx = sm.cidx ();
260  octave_idx_type* ridx = sm.ridx ();
261  T* data = sm.data ();
262 
263  // Output matrix data structures. Because the final zero pattern pattern of
264  // the output matrix is not known due to fill-in elements, a heuristic
265  // approach has been adopted for memory allocation. The size of ridx_out_l
266  // and data_out_l is incremented 10% of their actual size (nnz (A) in the
267  // beginning). If that amount is less than n, their size is just incremented
268  // in n elements. This way the number of reallocations decreases throughout
269  // the process, obtaining a good performance.
270  max_len = sm.nnz ();
271  max_len += (0.1 * max_len) > n ? 0.1 * max_len : n;
272  Array <octave_idx_type> cidx_out_l (dim_vector (n + 1, 1));
273  octave_idx_type* cidx_l = cidx_out_l.fortran_vec ();
274  Array <octave_idx_type> ridx_out_l (dim_vector (max_len ,1));
275  octave_idx_type* ridx_l = ridx_out_l.fortran_vec ();
276  Array <T> data_out_l (dim_vector (max_len, 1));
277  T* data_l = data_out_l.fortran_vec ();
278 
279  // Working arrays
280  OCTAVE_LOCAL_BUFFER (T, w_data, n);
283  OCTAVE_LOCAL_BUFFER (T, col_drops, n);
284  std::vector <octave_idx_type> vec;
285  vec.resize (n);
286 
287  T zero = T (0);
288  cidx_l[0] = cidx[0];
289  for (i = 0; i < n; i++)
290  {
291  Llist[i] = -1;
292  Lfirst[i] = -1;
293  w_data[i] = 0;
294  col_drops[i] = zero;
295  vec[i] = 0;
296  }
297 
298  total_len = 0;
299  for (k = 0; k < n; k++)
300  {
301  ind = 0;
302  for (j = cidx[k]; j < cidx[k+1]; j++)
303  {
304  w_data[ridx[j]] = data[j];
305  if (ridx[j] != k)
306  {
307  vec[ind] = ridx[j];
308  ind++;
309  }
310  }
311  jrow = Llist[k];
312  while (jrow != -1)
313  {
314  jjrow = Lfirst[jrow];
315  jend = cidx_l[jrow+1];
316  for (jj = jjrow; jj < jend; jj++)
317  {
318  j = ridx_l[jj];
319  // If the element in the j position of the row is zero,
320  // then it will become non-zero, so we add it to the
321  // vector that tracks non-zero elements in the working row.
322  if (w_data[j] == zero)
323  {
324  vec[ind] = j;
325  ind++;
326  }
327  w_data[j] -= ichol_mult (data_l[jj], data_l[jjrow]);
328  }
329  // Update the actual column first element and
330  // update the linked list of the jrow row.
331  if ((jjrow + 1) < jend)
332  {
333  Lfirst[jrow]++;
334  j = jrow;
335  jrow = Llist[jrow];
336  Llist[j] = Llist[ridx_l[Lfirst[j]]];
337  Llist[ridx_l[Lfirst[j]]] = j;
338  }
339  else
340  jrow = Llist[jrow];
341  }
342 
343  // Resizing output arrays
344  if ((max_len - total_len) < n)
345  {
346  max_len += (0.1 * max_len) > n ? 0.1 * max_len : n;
347  data_out_l.resize (dim_vector (max_len, 1));
348  data_l = data_out_l.fortran_vec ();
349  ridx_out_l.resize (dim_vector (max_len, 1));
350  ridx_l = ridx_out_l.fortran_vec ();
351  }
352 
353  // The sorting of the non-zero elements of the working column can be
354  // handled in a couple of ways. The most efficient two I found, are
355  // keeping the elements in an ordered binary search tree dynamically or
356  // keep them unsorted in a vector and at the end of the outer iteration
357  // order them. The last approach exhibits lower execution times.
358  std::sort (vec.begin (), vec.begin () + ind);
359 
360  data_l[total_len] = w_data[k];
361  ridx_l[total_len] = k;
362  w_len = 1;
363 
364  // Extract the non-zero elements of working column and
365  // drop the elements that are lower than droptol * cols_norm[k].
366  for (i = 0; i < ind ; i++)
367  {
368  jrow = vec[i];
369  if (w_data[jrow] != zero)
370  {
371  if (std::abs (w_data[jrow]) < (droptol * cols_norm[k]))
372  {
373  if (opt == ON)
374  {
375  col_drops[k] += w_data[jrow];
376  col_drops[jrow] += w_data[jrow];
377  }
378  }
379  else
380  {
381  data_l[total_len + w_len] = w_data[jrow];
382  ridx_l[total_len + w_len] = jrow;
383  w_len++;
384  }
385  vec[i] = 0;
386  }
387  w_data[jrow] = zero;
388  }
389 
390  // Compensate column sums --> michol option
391  if (opt == ON)
392  data_l[total_len] += col_drops[k];
393 
394  if (data_l[total_len] == zero)
395  {
396  error ("ichol: encountered a pivot equal to 0");
397  break;
398  }
399  else if (! ichol_checkpivot (data_l[total_len]))
400  break;
401 
402  // Once elements are dropped and compensation of column sums are done,
403  // scale the elements by the pivot.
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];
407  total_len += w_len;
408  // Check if there are too many elements to be indexed with
409  // octave_idx_type type due to fill-in during the process.
410  if (total_len < 0)
411  {
412  error ("ichol: integer overflow. Too many fill-in elements in L");
413  break;
414  }
415  cidx_l[k+1] = cidx_l[k] - cidx_l[0] + w_len;
416 
417  // Update Llist and Lfirst with the k-column information.
418  if (k < (n - 1))
419  {
420  Lfirst[k] = cidx_l[k];
421  if ((Lfirst[k] + 1) < cidx_l[k+1])
422  {
423  Lfirst[k]++;
424  jjrow = ridx_l[Lfirst[k]];
425  Llist[k] = Llist[jjrow];
426  Llist[jjrow] = k;
427  }
428  }
429  }
430 
431  if (! error_state)
432  {
433  // Build the output matrices
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++)
438  {
439  L.ridx (i) = ridx_l[i];
440  L.data (i) = data_l[i];
441  }
442  }
443 }
444 
445 DEFUN (__icholt__, args, nargout,
446  "-*- texinfo -*-\n\
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\
451 @end deftypefn")
452 {
453  octave_value_list retval;
454  int nargin = args.length ();
455  // Default values of parameters
456  std::string michol = "off";
457  double droptol = 0;
458 
459  if (nargout > 1 || nargin < 1 || nargin > 3)
460  {
461  print_usage ();
462  return retval;
463  }
464 
465  // Don't repeat input validation of arguments done in ichol.m
466 
467  if (nargin >= 2)
468  droptol = args(1).double_value ();
469 
470  if (nargin == 3)
471  michol = args(2).string_value ();
472 
473  octave_value_list param_list;
474  if (! args(0).is_complex_type ())
475  {
476  Array <double> cols_norm;
477  SparseMatrix L;
478  param_list.append (args(0).sparse_matrix_value ());
479  SparseMatrix sm_l =
480  feval ("tril", param_list)(0).sparse_matrix_value ();
481  param_list(0) = sm_l;
482  param_list(1) = 1;
483  param_list(2) = "cols";
484  cols_norm = feval ("norm", param_list)(0).vector_value ();
485  param_list.clear ();
488  (sm_l, L, cols_norm.fortran_vec (), droptol, michol);
489  if (! error_state)
490  retval(0) = L;
491  }
492  else
493  {
494  Array <Complex> cols_norm;
496  param_list.append (args(0).sparse_complex_matrix_value ());
497  SparseComplexMatrix sm_l =
498  feval ("tril", param_list)(0).sparse_complex_matrix_value ();
499  param_list(0) = sm_l;
500  param_list(1) = 1;
501  param_list(2) = "cols";
502  cols_norm = feval ("norm", param_list)(0).complex_vector_value ();
503  param_list.clear ();
506  (sm_l, L, cols_norm.fortran_vec (),
507  Complex (droptol), michol);
508  if (! error_state)
509  retval(0) = L;
510  }
511 
512  return retval;
513 }
514 
515 /*
516 ## No test needed for internal helper function.
517 %!assert (1)
518 */
519 
void clear(void)
Definition: oct-obj.h:148
Complex ichol_mult_complex(Complex a, Complex b)
Definition: __ichol__.cc:35
OCTINTERP_API void print_usage(void)
Definition: defun.cc:51
octave_idx_type length(void) const
Definition: oct-obj.h:89
octave_value_list & append(const octave_value &val)
Definition: oct-obj.cc:85
#define DEFUN(name, args_name, nargout_name, doc)
Definition: defun.h:44
void error(const char *fmt,...)
Definition: error.cc:476
octave_value_list feval(const std::string &name, const octave_value_list &args, int nargout)
Definition: oct-parse.cc:8625
ComplexColumnVector conj(const ComplexColumnVector &a)
Definition: CColVector.cc:244
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:243
int error_state
Definition: error.cc:101
void resize(const dim_vector &dv, const T &rfv)
Definition: Array.cc:1033
bool ichol_checkpivot_complex(Complex pivot)
Definition: __ichol__.cc:52
double ichol_mult_real(double a, double b)
Definition: __ichol__.cc:47
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition: oct-locbuf.h:197
ColumnVector imag(const ComplexColumnVector &a)
Definition: dColVector.cc:162
void ichol_0(octave_matrix_t &sm, const std::string michol="off")
Definition: __ichol__.cc:79
std::complex< double > Complex
Definition: oct-cmplx.h:29
const T * fortran_vec(void) const
Definition: Array.h:481
bool ichol_checkpivot_real(double pivot)
Definition: __ichol__.cc:67
T abs(T x)
Definition: pr-output.cc:3062