Mercurial > octave

comparison libinterp/dldfcn/symbfact.cc @ 21585:bad3ed83330d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
symbfact.cc: Overhaul dldfcn. * symbfact.cc: Redo docstring. Eliminate nargout input validation. Improve input validation for TYP and MODE variables. Add label cleanup: and use goto statements to guarantee that memory allocated in cholmod library is cleaned up. Use std::max rather than hand-coded max routine. Add BIST tests.
author Rik <rik@octave.org>
date Mon, 04 Apr 2016 09:00:22 -0700
parents ad0599a0acc6
children d7a268e68e69
comparison
equal deleted inserted replaced
21584:ee1a009dd60f 21585:bad3ed83330d
39 #include "ovl.h" 39 #include "ovl.h"
40 #include "utils.h" 40 #include "utils.h"
41 41
42 DEFUN_DLD (symbfact, args, nargout, 42 DEFUN_DLD (symbfact, args, nargout,
43 "-*- texinfo -*-\n\ 43 "-*- texinfo -*-\n\
44 @deftypefn {} {[@var{count}, @var{h}, @var{parent}, @var{post}, @var{r}] =} symbfact (@var{S})\n\ 44 @deftypefn {} {[@var{count}, @var{h}, @var{parent}, @var{post}, @var{R}] =} symbfact (@var{S})\n\
45 @deftypefnx {} {[@dots{}] =} symbfact (@var{S}, @var{typ})\n\ 45 @deftypefnx {} {[@dots{}] =} symbfact (@var{S}, @var{typ})\n\
46 @deftypefnx {} {[@dots{}] =} symbfact (@var{S}, @var{typ}, @var{mode})\n\ 46 @deftypefnx {} {[@dots{}] =} symbfact (@var{S}, @var{typ}, @var{mode})\n\
47 \n\ 47 \n\
48 Perform a symbolic factorization analysis on the sparse matrix @var{S}.\n\ 48 Perform a symbolic factorization analysis of the sparse matrix @var{S}.\n\
49 \n\ 49 \n\
50 The input variables are\n\ 50 The input variables are\n\
51 \n\ 51 \n\
52 @table @var\n\ 52 @table @var\n\
53 @item S\n\ 53 @item S\n\
54 @var{S} is a complex or real sparse matrix.\n\ 54 @var{S} is a real or complex sparse matrix.\n\
55 \n\ 55 \n\
56 @item typ\n\ 56 @item typ\n\
57 Is the type of the factorization and can be one of\n\ 57 Is the type of the factorization and can be one of\n\
58 \n\ 58 \n\
59 @table @samp\n\ 59 @table @asis\n\
60 @item sym\n\ 60 @item @qcode{\"sym\"} (default)\n\
61 Factorize @var{S}. This is the default.\n\ 61 Factorize @var{S}. Assumes @var{S} is symmetric and uses the upper\n\
62 \n\ 62 triangular portion of the matrix.\n\
63 @item col\n\ 63 \n\
64 Factorize @code{@var{S}' * @var{S}}.\n\ 64 @item @qcode{\"col\"}\n\
65 \n\ 65 Factorize @tcode{@var{S}' * @var{S}}.\n\
66 @item row\n\ 66 \n\
67 @item @qcode{\"row\"}\n\
67 Factorize @tcode{@var{S} * @var{S}'}.\n\ 68 Factorize @tcode{@var{S} * @var{S}'}.\n\
68 \n\ 69 \n\
69 @item lo\n\ 70 @item @qcode{\"lo\"}\n\
70 Factorize @tcode{@var{S}'}\n\ 71 Factorize @tcode{@var{S}'}. Assumes @var{S} is symmetric and uses the lower\n\
72 triangular portion of the matrix.\n\
71 @end table\n\ 73 @end table\n\
72 \n\ 74 \n\
73 @item mode\n\ 75 @item mode\n\
74 The default is to return the Cholesky@tie{}factorization for @var{r}, and if\n\ 76 When @var{mode} is unspecified return the Cholesky@tie{}factorization for\n\
75 @var{mode} is @qcode{'L'}, the conjugate transpose of the\n\ 77 @var{R}. If @var{mode} is @qcode{\"lower\"} or @qcode{\"L\"} then return\n\
76 Cholesky@tie{}factorization is returned. The conjugate transpose version is\n\ 78 the conjugate transpose @tcode{@var{R}'} which is a lower triangular factor.\n\
77 faster and uses less memory, but returns the same values for @var{count},\n\ 79 The conjugate transpose version is faster and uses less memory, but still\n\
78 @var{h}, @var{parent} and @var{post} outputs.\n\ 80 returns the same values for all other outputs: @var{count}, @var{h},\n\
81 @var{parent}, and @var{post}.\n\
79 @end table\n\ 82 @end table\n\
80 \n\ 83 \n\
81 The output variables are\n\ 84 The output variables are:\n\
82 \n\ 85 \n\
83 @table @var\n\ 86 @table @var\n\
84 @item count\n\ 87 @item count\n\
85 The row counts of the Cholesky@tie{}factorization as determined by\n\ 88 The row counts of the Cholesky@tie{}factorization as determined by\n\
86 @var{typ}.\n\ 89 @var{typ}. The computational difficulty of performing the true\n\
90 factorization using @code{chol} is @code{sum (@var{count} .^ 2)}.\n\
87 \n\ 91 \n\
88 @item h\n\ 92 @item h\n\
89 The height of the elimination tree.\n\ 93 The height of the elimination tree.\n\
90 \n\ 94 \n\
91 @item parent\n\ 95 @item parent\n\
92 The elimination tree itself.\n\ 96 The elimination tree itself.\n\
93 \n\ 97 \n\
94 @item post\n\ 98 @item post\n\
95 A sparse boolean matrix whose structure is that of the Cholesky\n\ 99 A sparse boolean matrix whose structure is that of the\n\
96 factorization as determined by @var{typ}.\n\ 100 Cholesky@tie{}factorization as determined by @var{typ}.\n\
97 @end table\n\ 101 @end table\n\
102 @seealso{chol, etree, treelayout}\n\
98 @end deftypefn") 103 @end deftypefn")
99 { 104 {
100 #ifdef HAVE_CHOLMOD 105 #ifdef HAVE_CHOLMOD
101 106
102 int nargin = args.length (); 107 int nargin = args.length ();
103 108
104 if (nargin < 1 || nargin > 3 || nargout > 5) 109 if (nargin < 1 || nargin > 3)
105 print_usage (); 110 print_usage ();
106 111
107 octave_value_list retval; 112 octave_value_list retval;
108 113
109 cholmod_common Common;
110 cholmod_common *cm = &Common;
111 CHOLMOD_NAME(start) (cm);
112
113 double spu = octave_sparse_params::get_key ("spumoni");
114 if (spu == 0.)
115 {
116 cm->print = -1;
117 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, 0);
118 }
119 else
120 {
121 cm->print = static_cast<int> (spu) + 2;
122 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &SparseCholPrint);
123 }
124
125 cm->error_handler = &SparseCholError;
126 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
127 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
128
129 double dummy; 114 double dummy;
130 cholmod_sparse Astore; 115 cholmod_sparse Astore;
131 cholmod_sparse *A = &Astore; 116 cholmod_sparse *A = &Astore;
132 A->packed = true; 117 A->packed = true;
133 A->sorted = true; 118 A->sorted = true;
173 octave_idx_type coletree = false; 158 octave_idx_type coletree = false;
174 octave_idx_type n = A->nrow; 159 octave_idx_type n = A->nrow;
175 160
176 if (nargin > 1) 161 if (nargin > 1)
177 { 162 {
163 std::string str = args(1).xstring_value ("TYP must be a string");
164 // FIXME: The input validation could be improved to use strncmp
178 char ch; 165 char ch;
179 std::string str = args(1).string_value (); 166 ch = tolower (str[0]);
180 ch = tolower (str.c_str ()[0]); 167 if (ch == 'r') // 'row'
181 if (ch == 'r')
182 A->stype = 0; 168 A->stype = 0;
183 else if (ch == 'c') 169 else if (ch == 'c') // 'col'
184 { 170 {
185 n = A->ncol; 171 n = A->ncol;
186 coletree = true; 172 coletree = true;
187 A->stype = 0; 173 A->stype = 0;
188 } 174 }
189 else if (ch == 's') 175 else if (ch == 's') // 'sym' (default)
190 A->stype = 1; 176 A->stype = 1;
191 else if (ch == 's') 177 else if (ch == 'l') // 'lo'
192 A->stype = -1; 178 A->stype = -1;
193 else 179 else
194 error ("symbfact: unrecognized TYP in symbolic factorization"); 180 error ("symbfact: unrecognized TYP \"%s\"", str.c_str ());
181 }
182
183 if (nargin == 3)
184 {
185 std::string str = args(2).xstring_value ("MODE must be a string");
186 // FIXME: The input validation could be improved to use strncmp
187 char ch;
188 ch = toupper (str[0]);
189 if (ch != 'L')
190 error ("symbfact: unrecognized MODE \"%s\"", str.c_str ());
195 } 191 }
196 192
197 if (A->stype && A->nrow != A->ncol) 193 if (A->stype && A->nrow != A->ncol)
198 err_square_matrix_required ("symbfact", "S"); 194 err_square_matrix_required ("symbfact", "S");
199 195
201 OCTAVE_LOCAL_BUFFER (octave_idx_type, Post, n); 197 OCTAVE_LOCAL_BUFFER (octave_idx_type, Post, n);
202 OCTAVE_LOCAL_BUFFER (octave_idx_type, ColCount, n); 198 OCTAVE_LOCAL_BUFFER (octave_idx_type, ColCount, n);
203 OCTAVE_LOCAL_BUFFER (octave_idx_type, First, n); 199 OCTAVE_LOCAL_BUFFER (octave_idx_type, First, n);
204 OCTAVE_LOCAL_BUFFER (octave_idx_type, Level, n); 200 OCTAVE_LOCAL_BUFFER (octave_idx_type, Level, n);
205 201
202 cholmod_common Common;
203 cholmod_common *cm = &Common;
204 CHOLMOD_NAME(start) (cm);
205
206 double spu = octave_sparse_params::get_key ("spumoni");
207 if (spu == 0.)
208 {
209 cm->print = -1;
210 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, 0);
211 }
212 else
213 {
214 cm->print = static_cast<int> (spu) + 2;
215 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &SparseCholPrint);
216 }
217
218 cm->error_handler = &SparseCholError;
219 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
220 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
221
206 cholmod_sparse *F = CHOLMOD_NAME(transpose) (A, 0, cm); 222 cholmod_sparse *F = CHOLMOD_NAME(transpose) (A, 0, cm);
207 cholmod_sparse *Aup, *Alo; 223 cholmod_sparse *Aup, *Alo;
208 224
209 if (A->stype == 1 || coletree) 225 if (A->stype == 1 || coletree)
210 { 226 {
211 Aup = A ; 227 Aup = A;
212 Alo = F ; 228 Alo = F;
213 } 229 }
214 else 230 else
215 { 231 {
216 Aup = F ; 232 Aup = F;
217 Alo = A ; 233 Alo = A;
218 } 234 }
219 235
220 CHOLMOD_NAME(etree) (Aup, Parent, cm); 236 CHOLMOD_NAME(etree) (Aup, Parent, cm);
221 237
238 ColumnVector tmp (n); // Declaration must precede any goto cleanup.
239 std::string err_msg;
240
222 if (cm->status < CHOLMOD_OK) 241 if (cm->status < CHOLMOD_OK)
223 error ("symbfact: matrix corrupted"); 242 {
243 err_msg = "symbfact: matrix corrupted";
244 goto cleanup;
245 }
224 246
225 if (CHOLMOD_NAME(postorder) (Parent, n, 0, Post, cm) != n) 247 if (CHOLMOD_NAME(postorder) (Parent, n, 0, Post, cm) != n)
226 error ("symbfact: postorder failed"); 248 {
249 err_msg = "symbfact: postorder failed";
250 goto cleanup;
251 }
227 252
228 CHOLMOD_NAME(rowcolcounts) (Alo, 0, 0, Parent, Post, 0, 253 CHOLMOD_NAME(rowcolcounts) (Alo, 0, 0, Parent, Post, 0,
229 ColCount, First, Level, cm); 254 ColCount, First, Level, cm);
230 255
231 if (cm->status < CHOLMOD_OK) 256 if (cm->status < CHOLMOD_OK)
232 error ("symbfact: matrix corrupted"); 257 {
258 err_msg = "symbfact: matrix corrupted";
259 goto cleanup;
260 }
233 261
234 if (nargout > 4) 262 if (nargout > 4)
235 { 263 {
236 cholmod_sparse *A1, *A2; 264 cholmod_sparse *A1, *A2;
237 265
255 A1 = A; 283 A1 = A;
256 A2 = F; 284 A2 = F;
257 } 285 }
258 286
259 // count the total number of entries in L 287 // count the total number of entries in L
260 octave_idx_type lnz = 0 ; 288 octave_idx_type lnz = 0;
261 for (octave_idx_type j = 0 ; j < n ; j++) 289 for (octave_idx_type j = 0 ; j < n ; j++)
262 lnz += ColCount[j]; 290 lnz += ColCount[j];
263
264 291
265 // allocate the output matrix L (pattern-only) 292 // allocate the output matrix L (pattern-only)
266 SparseBoolMatrix L (n, n, lnz); 293 SparseBoolMatrix L (n, n, lnz);
267 294
268 // initialize column pointers 295 // initialize column pointers
271 { 298 {
272 L.xcidx(j) = lnz; 299 L.xcidx(j) = lnz;
273 lnz += ColCount[j]; 300 lnz += ColCount[j];
274 } 301 }
275 L.xcidx(n) = lnz; 302 L.xcidx(n) = lnz;
276
277 303
278 // create a copy of the column pointers 304 // create a copy of the column pointers
279 octave_idx_type *W = First; 305 octave_idx_type *W = First;
280 for (octave_idx_type j = 0 ; j < n ; j++) 306 for (octave_idx_type j = 0 ; j < n ; j++)
281 W[j] = L.xcidx (j); 307 W[j] = L.xcidx (j);
282 308
283 // get workspace for computing one row of L 309 // get workspace for computing one row of L
284 cholmod_sparse *R 310 cholmod_sparse *R
285 = CHOLMOD_NAME (allocate_sparse) (n, 1, n, false, true, 311 = CHOLMOD_NAME(allocate_sparse) (n, 1, n, false, true,
286 0, CHOLMOD_PATTERN, cm); 312 0, CHOLMOD_PATTERN, cm);
287 octave_idx_type *Rp = static_cast<octave_idx_type *>(R->p); 313 octave_idx_type *Rp = static_cast<octave_idx_type *> (R->p);
288 octave_idx_type *Ri = static_cast<octave_idx_type *>(R->i); 314 octave_idx_type *Ri = static_cast<octave_idx_type *> (R->i);
289 315
290 // compute L one row at a time 316 // compute L one row at a time
291 for (octave_idx_type k = 0 ; k < n ; k++) 317 for (octave_idx_type k = 0 ; k < n ; k++)
292 { 318 {
293 // get the kth row of L and store in the columns of L 319 // get the kth row of L and store in the columns of L
294 CHOLMOD_NAME (row_subtree) (A1, A2, k, Parent, R, cm) ; 320 CHOLMOD_NAME(row_subtree) (A1, A2, k, Parent, R, cm);
295 for (octave_idx_type p = 0 ; p < Rp[1] ; p++) 321 for (octave_idx_type p = 0 ; p < Rp[1] ; p++)
296 L.xridx (W[Ri[p]]++) = k ; 322 L.xridx (W[Ri[p]]++) = k;
297 323
298 // add the diagonal entry 324 // add the diagonal entry
299 L.xridx (W[k]++) = k ; 325 L.xridx (W[k]++) = k;
300 } 326 }
301 327
302 // free workspace 328 // free workspace
303 CHOLMOD_NAME (free_sparse) (&R, cm) ; 329 CHOLMOD_NAME(free_sparse) (&R, cm);
304
305 330
306 // transpose L to get R, or leave as is 331 // transpose L to get R, or leave as is
307 if (nargin < 3) 332 if (nargin < 3)
308 L = L.transpose (); 333 L = L.transpose ();
309 334
312 L.xdata(p) = true; 337 L.xdata(p) = true;
313 338
314 retval(4) = L; 339 retval(4) = L;
315 } 340 }
316 341
317 ColumnVector tmp (n);
318 if (nargout > 3) 342 if (nargout > 3)
319 { 343 {
320 for (octave_idx_type i = 0; i < n; i++) 344 for (octave_idx_type i = 0; i < n; i++)
321 tmp(i) = Post[i] + 1; 345 tmp(i) = Post[i] + 1;
322 retval(3) = tmp; 346 retval(3) = tmp;
330 } 354 }
331 355
332 if (nargout > 1) 356 if (nargout > 1)
333 { 357 {
334 // compute the elimination tree height 358 // compute the elimination tree height
335 octave_idx_type height = 0 ; 359 octave_idx_type height = 0;
336 for (int i = 0 ; i < n ; i++) 360 for (int i = 0 ; i < n ; i++)
337 height = (height > Level[i] ? height : Level[i]); 361 height = std::max (height, Level[i]);
338 height++ ; 362 height++;
339 retval(1) = static_cast<double> (height); 363 retval(1) = static_cast<double> (height);
340 } 364 }
341 365
342 for (octave_idx_type i = 0; i < n; i++) 366 for (octave_idx_type i = 0; i < n; i++)
343 tmp(i) = ColCount[i]; 367 tmp(i) = ColCount[i];
344 retval(0) = tmp; 368 retval(0) = tmp;
369
370 cleanup:
371 CHOLMOD_NAME(free_sparse) (&F, cm);
372 CHOLMOD_NAME(finish) (cm);
373
374 if (! err_msg.empty ())
375 error (err_msg.c_str ());
345 376
346 return retval; 377 return retval;
347 378
348 #else 379 #else
349 err_disabled_feature ("symbfact", "CHOLMOD"); 380 err_disabled_feature ("symbfact", "CHOLMOD");
350 #endif 381 #endif
351 } 382 }
383
384 /*
385 %!testif HAVE_CHOLMOD
386 %! A = sparse (magic (3));
387 %! [count, h, parent, post, r] = symbfact (A);
388 %! assert (count, [3; 2; 1]);
389 %! assert (h, 3);
390 %! assert (parent, [2; 3; 0]);
391 %! assert (r, sparse (triu (true (3))));
392
393 %!testif HAVE_CHOLMOD
394 %! ## Test MODE "lower"
395 %! A = sparse (magic (3));
396 %! [~, ~, ~, ~, l] = symbfact (A, "sym", "lower");
397 %! assert (l, sparse (tril (true (3))));
398
399 %!testif HAVE_CHOLMOD
400 %! ## Bug #42587, singular matrix
401 %! A = sparse ([1 0 8;0 1 8;8 8 1]);
402 %! [count, h, parent, post, r] = symbfact (A);
403
404 ## Test input validation
405 %!testif HAVE_CHOLMOD
406 %! fail ("symbfact ()");
407 %! fail ("symbfact (1,2,3,4)");
408 %! fail ("symbfact ({1})", "wrong type argument 'cell'");
409 %! fail ("symbfact (sparse (1), {1})", "TYP must be a string");
410 %! fail ("symbfact (sparse (1), 'foobar')", 'unrecognized TYP "foobar"');
411 %! fail ("symbfact (sparse (1), 'sym', {'L'})", "MODE must be a string");
412 %! fail ('symbfact (sparse (1), "sym", "foobar")', 'unrecognized MODE "foobar"');
413 %! fail ("symbfact (sparse ([1, 2; 3, 4; 5, 6]))", "S must be a square matrix");
414
415 */