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
symbfact.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2005-2015 David Bateman
4 Copyright (C) 1998-2005 Andy Adler
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 "SparseCmplxCHOL.h"
29 #include "SparsedbleCHOL.h"
30 #include "oct-spparms.h"
31 #include "sparse-util.h"
32 #include "oct-locbuf.h"
33 
34 #include "ov-re-sparse.h"
35 #include "ov-cx-sparse.h"
36 #include "defun-dld.h"
37 #include "error.h"
38 #include "gripes.h"
39 #include "oct-obj.h"
40 #include "utils.h"
41 
42 DEFUN_DLD (symbfact, args, nargout,
43  "-*- texinfo -*-\n\
44 @deftypefn {Loadable Function} {[@var{count}, @var{h}, @var{parent}, @var{post}, @var{r}] =} symbfact (@var{S})\n\
45 @deftypefnx {Loadable Function} {[@dots{}] =} symbfact (@var{S}, @var{typ})\n\
46 @deftypefnx {Loadable Function} {[@dots{}] =} symbfact (@var{S}, @var{typ}, @var{mode})\n\
47 \n\
48 Perform a symbolic factorization analysis on the sparse matrix @var{S}.\n\
49 \n\
50 The input variables are\n\
51 \n\
52 @table @var\n\
53 @item S\n\
54 @var{S} is a complex or real sparse matrix.\n\
55 \n\
56 @item typ\n\
57 Is the type of the factorization and can be one of\n\
58 \n\
59 @table @samp\n\
60 @item sym\n\
61 Factorize @var{S}. This is the default.\n\
62 \n\
63 @item col\n\
64 Factorize @code{@var{S}' * @var{S}}.\n\
65 \n\
66 @item row\n\
67 Factorize @tcode{@var{S} * @var{S}'}.\n\
68 \n\
69 @item lo\n\
70 Factorize @tcode{@var{S}'}\n\
71 @end table\n\
72 \n\
73 @item mode\n\
74 The default is to return the Cholesky@tie{}factorization for @var{r}, and if\n\
75 @var{mode} is @qcode{'L'}, the conjugate transpose of the\n\
76 Cholesky@tie{}factorization is returned. The conjugate transpose version is\n\
77 faster and uses less memory, but returns the same values for @var{count},\n\
78 @var{h}, @var{parent} and @var{post} outputs.\n\
79 @end table\n\
80 \n\
81 The output variables are\n\
82 \n\
83 @table @var\n\
84 @item count\n\
85 The row counts of the Cholesky@tie{}factorization as determined by @var{typ}.\n\
86 \n\
87 @item h\n\
88 The height of the elimination tree.\n\
89 \n\
90 @item parent\n\
91 The elimination tree itself.\n\
92 \n\
93 @item post\n\
94 A sparse boolean matrix whose structure is that of the Cholesky\n\
95 factorization as determined by @var{typ}.\n\
96 @end table\n\
97 @end deftypefn")
98 {
99  octave_value_list retval;
100  int nargin = args.length ();
101 
102  if (nargin < 1 || nargin > 3 || nargout > 5)
103  {
104  print_usage ();
105  return retval;
106  }
107 
108 #ifdef HAVE_CHOLMOD
109 
110  cholmod_common Common;
111  cholmod_common *cm = &Common;
112  CHOLMOD_NAME(start) (cm);
113 
114  double spu = octave_sparse_params::get_key ("spumoni");
115  if (spu == 0.)
116  {
117  cm->print = -1;
118  SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, 0);
119  }
120  else
121  {
122  cm->print = static_cast<int> (spu) + 2;
123  SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &SparseCholPrint);
124  }
125 
126  cm->error_handler = &SparseCholError;
127  SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
128  SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
129 
130  double dummy;
131  cholmod_sparse Astore;
132  cholmod_sparse *A = &Astore;
133  A->packed = true;
134  A->sorted = true;
135  A->nz = 0;
136 #ifdef USE_64_BIT_IDX_T
137  A->itype = CHOLMOD_LONG;
138 #else
139  A->itype = CHOLMOD_INT;
140 #endif
141  A->dtype = CHOLMOD_DOUBLE;
142  A->stype = 1;
143  A->x = &dummy;
144 
145  if (args(0).is_real_type ())
146  {
147  const SparseMatrix a = args(0).sparse_matrix_value ();
148  A->nrow = a.rows ();
149  A->ncol = a.cols ();
150  A->p = a.cidx ();
151  A->i = a.ridx ();
152  A->nzmax = a.nnz ();
153  A->xtype = CHOLMOD_REAL;
154 
155  if (a.rows () > 0 && a.cols () > 0)
156  A->x = a.data ();
157  }
158  else if (args(0).is_complex_type ())
159  {
160  const SparseComplexMatrix a = args(0).sparse_complex_matrix_value ();
161  A->nrow = a.rows ();
162  A->ncol = a.cols ();
163  A->p = a.cidx ();
164  A->i = a.ridx ();
165  A->nzmax = a.nnz ();
166  A->xtype = CHOLMOD_COMPLEX;
167 
168  if (a.rows () > 0 && a.cols () > 0)
169  A->x = a.data ();
170  }
171  else
172  gripe_wrong_type_arg ("symbfact", args(0));
173 
174  octave_idx_type coletree = false;
175  octave_idx_type n = A->nrow;
176 
177  if (nargin > 1)
178  {
179  char ch;
180  std::string str = args(1).string_value ();
181  ch = tolower (str.c_str ()[0]);
182  if (ch == 'r')
183  A->stype = 0;
184  else if (ch == 'c')
185  {
186  n = A->ncol;
187  coletree = true;
188  A->stype = 0;
189  }
190  else if (ch == 's')
191  A->stype = 1;
192  else if (ch == 's')
193  A->stype = -1;
194  else
195  error ("symbfact: unrecognized TYP in symbolic factorization");
196  }
197 
198  if (A->stype && A->nrow != A->ncol)
199  error ("symbfact: S must be a square matrix");
200 
201  if (!error_state)
202  {
205  OCTAVE_LOCAL_BUFFER (octave_idx_type, ColCount, n);
208 
209  cholmod_sparse *F = CHOLMOD_NAME(transpose) (A, 0, cm);
210  cholmod_sparse *Aup, *Alo;
211 
212  if (A->stype == 1 || coletree)
213  {
214  Aup = A ;
215  Alo = F ;
216  }
217  else
218  {
219  Aup = F ;
220  Alo = A ;
221  }
222 
223  CHOLMOD_NAME(etree) (Aup, Parent, cm);
224 
225  if (cm->status < CHOLMOD_OK)
226  {
227  error ("matrix corrupted");
228  goto symbfact_error;
229  }
230 
231  if (CHOLMOD_NAME(postorder) (Parent, n, 0, Post, cm) != n)
232  {
233  error ("postorder failed");
234  goto symbfact_error;
235  }
236 
237  CHOLMOD_NAME(rowcolcounts) (Alo, 0, 0, Parent, Post, 0,
238  ColCount, First, Level, cm);
239 
240  if (cm->status < CHOLMOD_OK)
241  {
242  error ("matrix corrupted");
243  goto symbfact_error;
244  }
245 
246  if (nargout > 4)
247  {
248  cholmod_sparse *A1, *A2;
249 
250  if (A->stype == 1)
251  {
252  A1 = A;
253  A2 = 0;
254  }
255  else if (A->stype == -1)
256  {
257  A1 = F;
258  A2 = 0;
259  }
260  else if (coletree)
261  {
262  A1 = F;
263  A2 = A;
264  }
265  else
266  {
267  A1 = A;
268  A2 = F;
269  }
270 
271  // count the total number of entries in L
272  octave_idx_type lnz = 0 ;
273  for (octave_idx_type j = 0 ; j < n ; j++)
274  lnz += ColCount[j];
275 
276 
277  // allocate the output matrix L (pattern-only)
278  SparseBoolMatrix L (n, n, lnz);
279 
280  // initialize column pointers
281  lnz = 0;
282  for (octave_idx_type j = 0 ; j < n ; j++)
283  {
284  L.xcidx(j) = lnz;
285  lnz += ColCount[j];
286  }
287  L.xcidx(n) = lnz;
288 
289 
290  /* create a copy of the column pointers */
291  octave_idx_type *W = First;
292  for (octave_idx_type j = 0 ; j < n ; j++)
293  W[j] = L.xcidx (j);
294 
295  // get workspace for computing one row of L
296  cholmod_sparse *R
297  = CHOLMOD_NAME (allocate_sparse) (n, 1, n, false, true,
298  0, CHOLMOD_PATTERN, cm);
299  octave_idx_type *Rp = static_cast<octave_idx_type *>(R->p);
300  octave_idx_type *Ri = static_cast<octave_idx_type *>(R->i);
301 
302  // compute L one row at a time
303  for (octave_idx_type k = 0 ; k < n ; k++)
304  {
305  // get the kth row of L and store in the columns of L
306  CHOLMOD_NAME (row_subtree) (A1, A2, k, Parent, R, cm) ;
307  for (octave_idx_type p = 0 ; p < Rp[1] ; p++)
308  L.xridx (W[Ri[p]]++) = k ;
309 
310  // add the diagonal entry
311  L.xridx (W[k]++) = k ;
312  }
313 
314  // free workspace
315  CHOLMOD_NAME (free_sparse) (&R, cm) ;
316 
317 
318  // transpose L to get R, or leave as is
319  if (nargin < 3)
320  L = L.transpose ();
321 
322  // fill numerical values of L with one's
323  for (octave_idx_type p = 0 ; p < lnz ; p++)
324  L.xdata(p) = true;
325 
326  retval(4) = L;
327  }
328 
329  ColumnVector tmp (n);
330  if (nargout > 3)
331  {
332  for (octave_idx_type i = 0; i < n; i++)
333  tmp(i) = Post[i] + 1;
334  retval(3) = tmp;
335  }
336 
337  if (nargout > 2)
338  {
339  for (octave_idx_type i = 0; i < n; i++)
340  tmp(i) = Parent[i] + 1;
341  retval(2) = tmp;
342  }
343 
344  if (nargout > 1)
345  {
346  /* compute the elimination tree height */
347  octave_idx_type height = 0 ;
348  for (int i = 0 ; i < n ; i++)
349  height = (height > Level[i] ? height : Level[i]);
350  height++ ;
351  retval(1) = static_cast<double> (height);
352  }
353 
354  for (octave_idx_type i = 0; i < n; i++)
355  tmp(i) = ColCount[i];
356  retval(0) = tmp;
357  }
358 
359 symbfact_error:
360 #else
361  error ("symbfact: not available in this version of Octave");
362 #endif
363 
364  return retval;
365 }
octave_idx_type * xridx(void)
Definition: Sparse.h:524
octave_idx_type cols(void) const
Definition: Sparse.h:264
T * data(void)
Definition: Sparse.h:509
void gripe_wrong_type_arg(const char *name, const char *s, bool is_error)
Definition: gripes.cc:135
octave_idx_type rows(void) const
Definition: Sparse.h:263
OCTINTERP_API void print_usage(void)
Definition: defun.cc:51
octave_idx_type length(void) const
Definition: oct-obj.h:89
void SparseCholError(int status, char *file, int line, char *message)
Definition: sparse-util.cc:39
F77_RET_T const octave_idx_type Complex * A
Definition: CmplxGEPBAL.cc:39
static void transpose(octave_idx_type N, const octave_idx_type *ridx, const octave_idx_type *cidx, octave_idx_type *ridx2, octave_idx_type *cidx2)
Definition: symrcm.cc:382
octave_idx_type * xcidx(void)
Definition: Sparse.h:537
void error(const char *fmt,...)
Definition: error.cc:476
SparseBoolMatrix transpose(void) const
Definition: boolSparse.h:86
octave_idx_type * cidx(void)
Definition: Sparse.h:531
static void coletree(const octave_idx_type *ridx, const octave_idx_type *colbeg, octave_idx_type *colend, octave_idx_type *parent, octave_idx_type nr, octave_idx_type nc)
Definition: colamd.cc:166
static double get_key(const std::string &key)
Definition: oct-spparms.cc:99
octave_idx_type nnz(void) const
Definition: Sparse.h:248
int error_state
Definition: error.cc:101
int SparseCholPrint(const char *fmt,...)
Definition: sparse-util.cc:57
octave_idx_type * ridx(void)
Definition: Sparse.h:518
T * xdata(void)
Definition: Sparse.h:511
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition: oct-locbuf.h:197
#define DEFUN_DLD(name, args_name, nargout_name, doc)
Definition: defun-dld.h:59
void F(const TSRC *v, TRES *r, octave_idx_type m, octave_idx_type n)
Definition: mx-inlines.cc:527