Mercurial > octave-nkf

comparison src/DLD-FUNCTIONS/symbfact.cc @ 7515:f3c00dc0912b

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