Mercurial > hg > octave-lyh
annotate src/DLD-FUNCTIONS/symbfact.cc @ 9064:7c02ec148a3c
Check grammar on all .cc files
Same check as previously done on .m files
Attempt to enforce some conformity in documentation text for rules
such as two spaces after a period, commas around latin abbreviations, etc.
| author | Rik <rdrider0-list@yahoo.com> |
|---|---|
| date | Sat, 28 Mar 2009 13:57:22 -0700 |
| parents | eb63fbe60fab |
| children | 40dfc0c99116 |
| rev | line source |
|---|---|
| 5506 | 1 /* |
| 2 | |
| 8920 | 3 Copyright (C) 2005, 2006, 2007, 2008 David Bateman |
| 7016 | 4 Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005 Andy Adler |
| 5 | |
| 6 This file is part of Octave. | |
| 5506 | 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 | |
| 7016 | 10 Free Software Foundation; either version 3 of the License, or (at your |
| 11 option) any later version. | |
| 5506 | 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 | |
| 7016 | 19 along with Octave; see the file COPYING. If not, see |
| 20 <http://www.gnu.org/licenses/>. | |
| 5506 | 21 |
| 22 */ | |
| 23 | |
| 24 #ifdef HAVE_CONFIG_H | |
| 25 #include <config.h> | |
| 26 #endif | |
| 27 | |
|
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7036
diff
changeset
|
28 #include "SparseCmplxCHOL.h" |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7036
diff
changeset
|
29 #include "SparsedbleCHOL.h" |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7036
diff
changeset
|
30 #include "oct-spparms.h" |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7036
diff
changeset
|
31 #include "sparse-util.h" |
|
8377
25bc2d31e1bf
improve OCTAVE_LOCAL_BUFFER
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
32 #include "oct-locbuf.h" |
|
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7036
diff
changeset
|
33 |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7036
diff
changeset
|
34 #include "ov-re-sparse.h" |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7036
diff
changeset
|
35 #include "ov-cx-sparse.h" |
| 5506 | 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\ | |
| 6547 | 44 @deftypefn {Loadable Function} {[@var{count}, @var{h}, @var{parent}, @var{post}, @var{r}] =} symbfact (@var{s}, @var{typ}, @var{mode})\n\ |
| 5506 | 45 \n\ |
| 46 Performs a symbolic factorization analysis on the sparse matrix @var{s}.\n\ | |
| 47 Where\n\ | |
| 48 \n\ | |
| 49 @table @asis\n\ | |
| 50 @item @var{s}\n\ | |
| 51 @var{s} is a complex or real sparse matrix.\n\ | |
| 52 \n\ | |
| 53 @item @var{typ}\n\ | |
| 54 Is the type of the factorization and can be one of\n\ | |
| 55 \n\ | |
| 56 @table @code\n\ | |
| 57 @item sym\n\ | |
|
9064
7c02ec148a3c
Check grammar on all .cc files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
58 Factorize @var{s}. This is the default.\n\ |
| 5506 | 59 \n\ |
| 60 @item col\n\ | |
| 61 Factorize @code{@var{s}' * @var{s}}.\n\ | |
| 62 @item row\n\ | |
| 63 Factorize @code{@var{s} * @var{s}'}.\n\ | |
| 64 @item lo\n\ | |
| 65 Factorize @code{@var{s}'}\n\ | |
| 66 @end table\n\ | |
| 67 \n\ | |
| 68 @item @var{mode}\n\ | |
| 69 The default is to return the Cholesky factorization for @var{r}, and if\n\ | |
| 7001 | 70 @var{mode} is 'L', the conjugate transpose of the Cholesky factorization\n\ |
|
9064
7c02ec148a3c
Check grammar on all .cc files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
71 is returned. The conjugate transpose version is faster and uses less\n\ |
| 5506 | 72 memory, but returns the same values for @var{count}, @var{h}, @var{parent}\n\ |
| 73 and @var{post} outputs.\n\ | |
| 74 @end table\n\ | |
| 75 \n\ | |
| 76 The output variables are\n\ | |
| 77 \n\ | |
| 78 @table @asis\n\ | |
| 79 @item @var{count}\n\ | |
| 80 The row counts of the Cholesky factorization as determined by @var{typ}.\n\ | |
| 81 \n\ | |
| 82 @item @var{h}\n\ | |
| 83 The height of the elimination tree.\n\ | |
| 84 \n\ | |
| 85 @item @var{parent}\n\ | |
| 86 The elimination tree itself.\n\ | |
| 87 \n\ | |
| 88 @item @var{post}\n\ | |
| 89 A sparse boolean matrix whose structure is that of the Cholesky\n\ | |
| 90 factorization as determined by @var{typ}.\n\ | |
| 91 @end table\n\ | |
| 92 @end deftypefn") | |
| 93 { | |
| 94 octave_value_list retval; | |
| 95 int nargin = args.length (); | |
| 96 | |
| 97 if (nargin < 1 || nargin > 3 || nargout > 5) | |
| 98 { | |
| 5823 | 99 print_usage (); |
| 5506 | 100 return retval; |
| 101 } | |
| 102 | |
| 5512 | 103 #ifdef HAVE_CHOLMOD |
| 104 | |
| 5506 | 105 cholmod_common Common; |
| 106 cholmod_common *cm = &Common; | |
| 107 CHOLMOD_NAME(start) (cm); | |
| 108 | |
| 5893 | 109 double spu = octave_sparse_params::get_key ("spumoni"); |
| 5506 | 110 if (spu == 0.) |
| 111 { | |
| 112 cm->print = -1; | |
| 7520 | 113 cm->print_function = 0; |
| 5506 | 114 } |
| 115 else | |
| 116 { | |
| 5760 | 117 cm->print = static_cast<int> (spu) + 2; |
| 5506 | 118 cm->print_function =&SparseCholPrint; |
| 119 } | |
| 120 | |
| 121 cm->error_handler = &SparseCholError; | |
| 122 cm->complex_divide = CHOLMOD_NAME(divcomplex); | |
| 123 cm->hypotenuse = CHOLMOD_NAME(hypot); | |
| 124 | |
| 125 double dummy; | |
| 126 cholmod_sparse Astore; | |
| 127 cholmod_sparse *A = &Astore; | |
| 5527 | 128 A->packed = true; |
| 129 A->sorted = true; | |
| 7520 | 130 A->nz = 0; |
| 5506 | 131 #ifdef IDX_TYPE_LONG |
| 132 A->itype = CHOLMOD_LONG; | |
| 133 #else | |
| 134 A->itype = CHOLMOD_INT; | |
| 135 #endif | |
| 136 A->dtype = CHOLMOD_DOUBLE; | |
| 137 A->stype = 1; | |
| 138 A->x = &dummy; | |
| 139 | |
| 140 if (args(0).is_real_type ()) | |
| 141 { | |
| 142 const SparseMatrix a = args(0).sparse_matrix_value(); | |
| 143 A->nrow = a.rows(); | |
| 144 A->ncol = a.cols(); | |
| 145 A->p = a.cidx(); | |
| 146 A->i = a.ridx(); | |
| 5604 | 147 A->nzmax = a.nnz(); |
| 5506 | 148 A->xtype = CHOLMOD_REAL; |
| 149 | |
| 150 if (a.rows() > 0 && a.cols() > 0) | |
| 151 A->x = a.data(); | |
| 152 } | |
| 153 else if (args(0).is_complex_type ()) | |
| 154 { | |
| 155 const SparseComplexMatrix a = args(0).sparse_complex_matrix_value(); | |
| 156 A->nrow = a.rows(); | |
| 157 A->ncol = a.cols(); | |
| 158 A->p = a.cidx(); | |
| 159 A->i = a.ridx(); | |
| 5604 | 160 A->nzmax = a.nnz(); |
| 5506 | 161 A->xtype = CHOLMOD_COMPLEX; |
| 162 | |
| 163 if (a.rows() > 0 && a.cols() > 0) | |
| 164 A->x = a.data(); | |
| 165 } | |
| 166 else | |
|
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7520
diff
changeset
|
167 gripe_wrong_type_arg ("symbfact", args(0)); |
| 5506 | 168 |
| 5527 | 169 octave_idx_type coletree = false; |
| 5506 | 170 octave_idx_type n = A->nrow; |
| 171 | |
| 172 if (nargin > 1) | |
| 173 { | |
| 174 char ch; | |
| 175 std::string str = args(1).string_value(); | |
| 176 ch = tolower (str.c_str()[0]); | |
| 177 if (ch == 'r') | |
| 178 A->stype = 0; | |
| 179 else if (ch == 'c') | |
| 180 { | |
| 181 n = A->ncol; | |
| 5527 | 182 coletree = true; |
| 5506 | 183 A->stype = 0; |
| 184 } | |
| 185 else if (ch == 's') | |
| 186 A->stype = 1; | |
| 187 else if (ch == 's') | |
| 188 A->stype = -1; | |
| 189 else | |
| 190 error ("Unrecognized typ in symbolic factorization"); | |
| 191 } | |
| 192 | |
| 193 if (A->stype && A->nrow != A->ncol) | |
| 194 error ("Matrix must be square"); | |
| 195 | |
| 196 if (!error_state) | |
| 197 { | |
| 198 OCTAVE_LOCAL_BUFFER (octave_idx_type, Parent, n); | |
| 199 OCTAVE_LOCAL_BUFFER (octave_idx_type, Post, n); | |
| 200 OCTAVE_LOCAL_BUFFER (octave_idx_type, ColCount, n); | |
| 201 OCTAVE_LOCAL_BUFFER (octave_idx_type, First, n); | |
| 202 OCTAVE_LOCAL_BUFFER (octave_idx_type, Level, n); | |
| 203 | |
| 204 cholmod_sparse *F = CHOLMOD_NAME(transpose) (A, 0, cm); | |
| 205 cholmod_sparse *Aup, *Alo; | |
| 206 | |
| 207 if (A->stype == 1 || coletree) | |
| 208 { | |
| 209 Aup = A ; | |
| 210 Alo = F ; | |
| 211 } | |
| 212 else | |
| 213 { | |
| 214 Aup = F ; | |
| 215 Alo = A ; | |
| 216 } | |
| 217 | |
| 218 CHOLMOD_NAME(etree) (Aup, Parent, cm); | |
| 219 | |
| 220 if (cm->status < CHOLMOD_OK) | |
| 221 { | |
| 222 error("matrix corrupted"); | |
| 223 goto symbfact_error; | |
| 224 } | |
| 225 | |
| 7520 | 226 if (CHOLMOD_NAME(postorder) (Parent, n, 0, Post, cm) != n) |
| 5506 | 227 { |
| 228 error("postorder failed"); | |
| 229 goto symbfact_error; | |
| 230 } | |
| 231 | |
| 7520 | 232 CHOLMOD_NAME(rowcolcounts) (Alo, 0, 0, Parent, Post, 0, |
| 5506 | 233 ColCount, First, Level, cm); |
| 234 | |
| 235 if (cm->status < CHOLMOD_OK) | |
| 236 { | |
| 237 error("matrix corrupted"); | |
| 238 goto symbfact_error; | |
| 239 } | |
| 240 | |
| 241 if (nargout > 4) | |
| 242 { | |
| 243 cholmod_sparse *A1, *A2; | |
| 244 | |
| 245 if (A->stype == 1) | |
| 246 { | |
| 247 A1 = A; | |
| 7520 | 248 A2 = 0; |
| 5506 | 249 } |
| 250 else if (A->stype == -1) | |
| 251 { | |
| 252 A1 = F; | |
| 7520 | 253 A2 = 0; |
| 5506 | 254 } |
| 255 else if (coletree) | |
| 256 { | |
| 257 A1 = F; | |
| 258 A2 = A; | |
| 259 } | |
| 260 else | |
| 261 { | |
| 262 A1 = A; | |
| 263 A2 = F; | |
| 264 } | |
| 265 | |
| 266 // count the total number of entries in L | |
| 267 octave_idx_type lnz = 0 ; | |
| 268 for (octave_idx_type j = 0 ; j < n ; j++) | |
| 269 lnz += ColCount [j] ; | |
| 270 | |
| 271 | |
| 272 // allocate the output matrix L (pattern-only) | |
| 273 SparseBoolMatrix L (n, n, lnz); | |
| 274 | |
| 275 // initialize column pointers | |
| 276 lnz = 0; | |
| 277 for (octave_idx_type j = 0 ; j < n ; j++) | |
| 278 { | |
| 279 L.xcidx(j) = lnz; | |
| 280 lnz += ColCount [j]; | |
| 281 } | |
| 282 L.xcidx(n) = lnz; | |
| 283 | |
| 284 | |
| 285 /* create a copy of the column pointers */ | |
| 286 octave_idx_type *W = First; | |
| 287 for (octave_idx_type j = 0 ; j < n ; j++) | |
| 288 W [j] = L.xcidx(j); | |
| 289 | |
| 290 // get workspace for computing one row of L | |
| 5527 | 291 cholmod_sparse *R = cholmod_allocate_sparse (n, 1, n, false, true, |
| 5506 | 292 0, CHOLMOD_PATTERN, cm); |
| 293 octave_idx_type *Rp = static_cast<octave_idx_type *>(R->p); | |
| 294 octave_idx_type *Ri = static_cast<octave_idx_type *>(R->i); | |
| 295 | |
| 296 // compute L one row at a time | |
| 297 for (octave_idx_type k = 0 ; k < n ; k++) | |
| 298 { | |
| 299 // get the kth row of L and store in the columns of L | |
| 5717 | 300 CHOLMOD_NAME (row_subtree) (A1, A2, k, Parent, R, cm) ; |
| 5506 | 301 for (octave_idx_type p = 0 ; p < Rp [1] ; p++) |
| 302 L.xridx (W [Ri [p]]++) = k ; | |
| 303 | |
| 304 // add the diagonal entry | |
| 305 L.xridx (W [k]++) = k ; | |
| 306 } | |
| 307 | |
| 308 // free workspace | |
| 309 cholmod_free_sparse (&R, cm) ; | |
| 310 | |
| 311 | |
| 312 // transpose L to get R, or leave as is | |
| 313 if (nargin < 3) | |
| 314 L = L.transpose (); | |
| 315 | |
| 316 // fill numerical values of L with one's | |
| 317 for (octave_idx_type p = 0 ; p < lnz ; p++) | |
| 318 L.xdata(p) = true; | |
| 319 | |
| 320 retval(4) = L; | |
| 321 } | |
| 322 | |
| 323 ColumnVector tmp (n); | |
| 324 if (nargout > 3) | |
| 325 { | |
| 326 for (octave_idx_type i = 0; i < n; i++) | |
| 327 tmp(i) = Post[i] + 1; | |
| 328 retval(3) = tmp; | |
| 329 } | |
| 330 | |
| 331 if (nargout > 2) | |
| 332 { | |
| 333 for (octave_idx_type i = 0; i < n; i++) | |
| 334 tmp(i) = Parent[i] + 1; | |
| 335 retval(2) = tmp; | |
| 336 } | |
| 337 | |
| 338 if (nargout > 1) | |
| 339 { | |
| 340 /* compute the elimination tree height */ | |
| 341 octave_idx_type height = 0 ; | |
| 342 for (int i = 0 ; i < n ; i++) | |
| 343 height = (height > Level[i] ? height : Level[i]); | |
| 344 height++ ; | |
| 5760 | 345 retval(1) = static_cast<double> (height); |
| 5506 | 346 } |
| 347 | |
| 348 for (octave_idx_type i = 0; i < n; i++) | |
| 349 tmp(i) = ColCount[i]; | |
| 350 retval(0) = tmp; | |
| 351 } | |
| 352 | |
| 5512 | 353 symbfact_error: |
| 354 #else | |
| 355 error ("symbfact: not available in this version of Octave"); | |
| 356 #endif | |
| 357 | |
| 5506 | 358 return retval; |
| 359 } | |
| 360 | |
| 361 /* | |
| 362 ;;; Local Variables: *** | |
| 363 ;;; mode: C++ *** | |
| 364 ;;; End: *** | |
| 365 */ | |
| 366 |