Mercurial > hg > octave-nkf
annotate liboctave/SparsedbleLU.cc @ 8924:3c3cbe8f18e0 ss-3-1-54
bump version info for snapshot
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Sat, 07 Mar 2009 13:28:03 -0500 |
| parents | eb63fbe60fab |
| children | 97de6c916498 |
| rev | line source |
|---|---|
| 5164 | 1 /* |
| 2 | |
| 8920 | 3 Copyright (C) 2004, 2005, 2006, 2007, 2008 David Bateman |
| 7016 | 4 Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004 Andy Adler |
| 5 | |
| 6 This file is part of Octave. | |
| 5164 | 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. | |
| 5164 | 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/>. | |
| 5164 | 21 |
| 22 */ | |
| 23 | |
| 24 #ifdef HAVE_CONFIG_H | |
| 25 #include <config.h> | |
| 26 #endif | |
| 27 | |
| 28 #include <vector> | |
| 29 | |
| 30 #include "lo-error.h" | |
|
8377
25bc2d31e1bf
improve OCTAVE_LOCAL_BUFFER
Jaroslav Hajek <highegg@gmail.com>
parents:
7520
diff
changeset
|
31 #include "oct-locbuf.h" |
| 5164 | 32 |
| 33 #include "SparsedbleLU.h" | |
| 34 #include "oct-spparms.h" | |
| 35 | |
| 36 // Instantiate the base LU class for the types we need. | |
| 37 | |
| 38 #include "sparse-base-lu.h" | |
| 39 #include "sparse-base-lu.cc" | |
| 40 | |
| 41 template class sparse_base_lu <SparseMatrix, double, SparseMatrix, double>; | |
| 42 | |
| 5451 | 43 #include "oct-sparse.h" |
| 5164 | 44 |
|
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
45 SparseLU::SparseLU (const SparseMatrix& a, const Matrix& piv_thres, bool scale) |
| 5164 | 46 { |
| 5203 | 47 #ifdef HAVE_UMFPACK |
| 5275 | 48 octave_idx_type nr = a.rows (); |
| 49 octave_idx_type nc = a.cols (); | |
| 5164 | 50 |
| 51 // Setup the control parameters | |
| 52 Matrix Control (UMFPACK_CONTROL, 1); | |
| 53 double *control = Control.fortran_vec (); | |
| 5322 | 54 UMFPACK_DNAME (defaults) (control); |
| 5164 | 55 |
| 5893 | 56 double tmp = octave_sparse_params::get_key ("spumoni"); |
| 5164 | 57 if (!xisnan (tmp)) |
| 58 Control (UMFPACK_PRL) = tmp; | |
| 59 | |
|
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
60 if (piv_thres.nelem() != 2) |
| 5164 | 61 { |
|
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
62 tmp = (piv_thres (0) > 1. ? 1. : piv_thres (0)); |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
63 if (!xisnan (tmp)) |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
64 Control (UMFPACK_PIVOT_TOLERANCE) = tmp; |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
65 tmp = (piv_thres (1) > 1. ? 1. : piv_thres (1)); |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
66 if (!xisnan (tmp)) |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
67 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp; |
| 5164 | 68 } |
| 69 else | |
| 70 { | |
| 5893 | 71 tmp = octave_sparse_params::get_key ("piv_tol"); |
| 5164 | 72 if (!xisnan (tmp)) |
|
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
73 Control (UMFPACK_PIVOT_TOLERANCE) = tmp; |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
74 |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
75 tmp = octave_sparse_params::get_key ("sym_tol"); |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
76 if (!xisnan (tmp)) |
| 5164 | 77 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp; |
| 78 } | |
| 79 | |
| 80 // Set whether we are allowed to modify Q or not | |
| 5893 | 81 tmp = octave_sparse_params::get_key ("autoamd"); |
| 5164 | 82 if (!xisnan (tmp)) |
| 83 Control (UMFPACK_FIXQ) = tmp; | |
| 84 | |
|
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
85 if (scale) |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
86 Control (UMFPACK_SCALE) = UMFPACK_SCALE_SUM; |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
87 else |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
88 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE; |
| 5164 | 89 |
| 5322 | 90 UMFPACK_DNAME (report_control) (control); |
| 5164 | 91 |
| 5275 | 92 const octave_idx_type *Ap = a.cidx (); |
| 93 const octave_idx_type *Ai = a.ridx (); | |
| 5164 | 94 const double *Ax = a.data (); |
| 95 | |
| 5322 | 96 UMFPACK_DNAME (report_matrix) (nr, nc, Ap, Ai, Ax, 1, control); |
| 5164 | 97 |
| 98 void *Symbolic; | |
| 99 Matrix Info (1, UMFPACK_INFO); | |
| 100 double *info = Info.fortran_vec (); | |
| 7520 | 101 int status = UMFPACK_DNAME (qsymbolic) (nr, nc, Ap, Ai, Ax, 0, |
| 5164 | 102 &Symbolic, control, info); |
| 103 | |
| 104 if (status < 0) | |
| 105 { | |
| 106 (*current_liboctave_error_handler) | |
| 107 ("SparseLU::SparseLU symbolic factorization failed"); | |
| 108 | |
| 5322 | 109 UMFPACK_DNAME (report_status) (control, status); |
| 110 UMFPACK_DNAME (report_info) (control, info); | |
| 5164 | 111 |
| 5322 | 112 UMFPACK_DNAME (free_symbolic) (&Symbolic) ; |
| 5164 | 113 } |
| 114 else | |
| 115 { | |
| 5322 | 116 UMFPACK_DNAME (report_symbolic) (Symbolic, control); |
| 5164 | 117 |
| 118 void *Numeric; | |
| 5322 | 119 status = UMFPACK_DNAME (numeric) (Ap, Ai, Ax, Symbolic, |
| 120 &Numeric, control, info) ; | |
| 121 UMFPACK_DNAME (free_symbolic) (&Symbolic) ; | |
| 5164 | 122 |
| 123 cond = Info (UMFPACK_RCOND); | |
| 124 | |
| 125 if (status < 0) | |
| 126 { | |
| 127 (*current_liboctave_error_handler) | |
| 128 ("SparseLU::SparseLU numeric factorization failed"); | |
| 129 | |
| 5322 | 130 UMFPACK_DNAME (report_status) (control, status); |
| 131 UMFPACK_DNAME (report_info) (control, info); | |
| 5164 | 132 |
| 5322 | 133 UMFPACK_DNAME (free_numeric) (&Numeric); |
| 5164 | 134 } |
| 135 else | |
| 136 { | |
| 5322 | 137 UMFPACK_DNAME (report_numeric) (Numeric, control); |
| 5164 | 138 |
| 5322 | 139 octave_idx_type lnz, unz, ignore1, ignore2, ignore3; |
| 140 status = UMFPACK_DNAME (get_lunz) (&lnz, &unz, &ignore1, | |
| 141 &ignore2, &ignore3, Numeric) ; | |
| 5164 | 142 |
| 143 if (status < 0) | |
| 144 { | |
| 145 (*current_liboctave_error_handler) | |
| 146 ("SparseLU::SparseLU extracting LU factors failed"); | |
| 147 | |
| 5322 | 148 UMFPACK_DNAME (report_status) (control, status); |
| 149 UMFPACK_DNAME (report_info) (control, info); | |
| 5164 | 150 |
| 5322 | 151 UMFPACK_DNAME (free_numeric) (&Numeric); |
| 5164 | 152 } |
| 153 else | |
| 154 { | |
| 5322 | 155 octave_idx_type n_inner = (nr < nc ? nr : nc); |
| 5164 | 156 |
| 157 if (lnz < 1) | |
| 5322 | 158 Lfact = SparseMatrix (n_inner, nr, |
| 5275 | 159 static_cast<octave_idx_type> (1)); |
| 5164 | 160 else |
| 5322 | 161 Lfact = SparseMatrix (n_inner, nr, lnz); |
| 5164 | 162 |
| 5275 | 163 octave_idx_type *Ltp = Lfact.cidx (); |
| 164 octave_idx_type *Ltj = Lfact.ridx (); | |
| 5164 | 165 double *Ltx = Lfact.data (); |
| 166 | |
| 167 if (unz < 1) | |
| 5322 | 168 Ufact = SparseMatrix (n_inner, nc, |
| 5275 | 169 static_cast<octave_idx_type> (1)); |
| 5164 | 170 else |
| 5322 | 171 Ufact = SparseMatrix (n_inner, nc, unz); |
| 5164 | 172 |
| 5275 | 173 octave_idx_type *Up = Ufact.cidx (); |
| 174 octave_idx_type *Uj = Ufact.ridx (); | |
| 5164 | 175 double *Ux = Ufact.data (); |
| 176 | |
|
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
177 Rfact = SparseMatrix (nr, nr, nr); |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
178 for (octave_idx_type i = 0; i < nr; i++) |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
179 { |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
180 Rfact.xridx (i) = i; |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
181 Rfact.xcidx (i) = i; |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
182 } |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
183 Rfact.xcidx (nr) = nr; |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
184 double *Rx = Rfact.data (); |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
185 |
| 5164 | 186 P.resize (nr); |
| 5322 | 187 octave_idx_type *p = P.fortran_vec (); |
| 5164 | 188 |
| 189 Q.resize (nc); | |
| 5322 | 190 octave_idx_type *q = Q.fortran_vec (); |
| 5164 | 191 |
| 5322 | 192 octave_idx_type do_recip; |
| 193 status = UMFPACK_DNAME (get_numeric) (Ltp, Ltj, Ltx, | |
| 7520 | 194 Up, Uj, Ux, p, q, 0, |
|
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
195 &do_recip, Rx, |
| 5164 | 196 Numeric) ; |
| 197 | |
| 5322 | 198 UMFPACK_DNAME (free_numeric) (&Numeric) ; |
| 5164 | 199 |
|
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
200 if (status < 0) |
| 5164 | 201 { |
| 202 (*current_liboctave_error_handler) | |
| 203 ("SparseLU::SparseLU extracting LU factors failed"); | |
| 204 | |
| 5322 | 205 UMFPACK_DNAME (report_status) (control, status); |
| 5164 | 206 } |
| 207 else | |
| 208 { | |
| 209 Lfact = Lfact.transpose (); | |
| 210 | |
|
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
211 if (do_recip) |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
212 for (octave_idx_type i = 0; i < nr; i++) |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
213 Rx[i] = 1.0 / Rx[i]; |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
214 |
| 5322 | 215 UMFPACK_DNAME (report_matrix) (nr, n_inner, |
| 216 Lfact.cidx (), Lfact.ridx (), | |
| 217 Lfact.data (), 1, control); | |
| 218 UMFPACK_DNAME (report_matrix) (n_inner, nc, | |
| 219 Ufact.cidx (), Ufact.ridx (), | |
| 220 Ufact.data (), 1, control); | |
| 221 UMFPACK_DNAME (report_perm) (nr, p, control); | |
| 222 UMFPACK_DNAME (report_perm) (nc, q, control); | |
| 5164 | 223 } |
| 224 | |
| 5322 | 225 UMFPACK_DNAME (report_info) (control, info); |
| 5164 | 226 } |
| 227 } | |
| 228 } | |
| 5203 | 229 #else |
| 230 (*current_liboctave_error_handler) ("UMFPACK not installed"); | |
| 231 #endif | |
| 5164 | 232 } |
| 233 | |
| 234 SparseLU::SparseLU (const SparseMatrix& a, const ColumnVector& Qinit, | |
|
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
235 const Matrix& piv_thres, bool scale, bool FixedQ, |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
236 double droptol, bool milu, bool udiag) |
| 5164 | 237 { |
| 5203 | 238 #ifdef HAVE_UMFPACK |
| 5282 | 239 if (milu) |
| 240 (*current_liboctave_error_handler) | |
| 241 ("Modified incomplete LU not implemented"); | |
| 5164 | 242 else |
| 243 { | |
| 5282 | 244 octave_idx_type nr = a.rows (); |
| 245 octave_idx_type nc = a.cols (); | |
| 5164 | 246 |
| 5282 | 247 // Setup the control parameters |
| 248 Matrix Control (UMFPACK_CONTROL, 1); | |
| 249 double *control = Control.fortran_vec (); | |
| 5322 | 250 UMFPACK_DNAME (defaults) (control); |
| 5164 | 251 |
| 5893 | 252 double tmp = octave_sparse_params::get_key ("spumoni"); |
| 5282 | 253 if (!xisnan (tmp)) |
| 254 Control (UMFPACK_PRL) = tmp; | |
|
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
255 |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
256 if (piv_thres.nelem() != 2) |
| 5282 | 257 { |
|
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
258 tmp = (piv_thres (0) > 1. ? 1. : piv_thres (0)); |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
259 if (!xisnan (tmp)) |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
260 Control (UMFPACK_PIVOT_TOLERANCE) = tmp; |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
261 tmp = (piv_thres (1) > 1. ? 1. : piv_thres (1)); |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
262 if (!xisnan (tmp)) |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
263 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp; |
| 5282 | 264 } |
| 265 else | |
| 266 { | |
| 5893 | 267 tmp = octave_sparse_params::get_key ("piv_tol"); |
| 5282 | 268 if (!xisnan (tmp)) |
|
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
269 Control (UMFPACK_PIVOT_TOLERANCE) = tmp; |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
270 |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
271 tmp = octave_sparse_params::get_key ("sym_tol"); |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
272 if (!xisnan (tmp)) |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
273 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp; |
| 5282 | 274 } |
| 5164 | 275 |
| 5282 | 276 if (droptol >= 0.) |
| 277 Control (UMFPACK_DROPTOL) = droptol; | |
| 5164 | 278 |
| 279 | |
| 5282 | 280 // Set whether we are allowed to modify Q or not |
| 281 if (FixedQ) | |
| 282 Control (UMFPACK_FIXQ) = 1.0; | |
| 283 else | |
| 284 { | |
| 5893 | 285 tmp = octave_sparse_params::get_key ("autoamd"); |
| 5282 | 286 if (!xisnan (tmp)) |
| 287 Control (UMFPACK_FIXQ) = tmp; | |
| 288 } | |
| 5164 | 289 |
|
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
290 if (scale) |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
291 Control (UMFPACK_SCALE) = UMFPACK_SCALE_SUM; |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
292 else |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
293 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE; |
| 5164 | 294 |
| 5322 | 295 UMFPACK_DNAME (report_control) (control); |
| 5164 | 296 |
| 5282 | 297 const octave_idx_type *Ap = a.cidx (); |
| 298 const octave_idx_type *Ai = a.ridx (); | |
| 299 const double *Ax = a.data (); | |
| 300 | |
| 5322 | 301 UMFPACK_DNAME (report_matrix) (nr, nc, Ap, Ai, Ax, 1, |
| 302 control); | |
| 5282 | 303 |
| 304 void *Symbolic; | |
| 305 Matrix Info (1, UMFPACK_INFO); | |
| 306 double *info = Info.fortran_vec (); | |
| 307 int status; | |
| 5164 | 308 |
| 5282 | 309 // Null loop so that qinit is imediately deallocated when not needed |
| 310 do { | |
| 5322 | 311 OCTAVE_LOCAL_BUFFER (octave_idx_type, qinit, nc); |
| 5164 | 312 |
| 5322 | 313 for (octave_idx_type i = 0; i < nc; i++) |
| 314 qinit [i] = static_cast<octave_idx_type> (Qinit (i)); | |
| 5164 | 315 |
| 5322 | 316 status = UMFPACK_DNAME (qsymbolic) (nr, nc, Ap, Ai, Ax, |
| 317 qinit, &Symbolic, control, info); | |
| 5282 | 318 } while (0); |
| 5164 | 319 |
| 320 if (status < 0) | |
| 321 { | |
| 322 (*current_liboctave_error_handler) | |
| 5282 | 323 ("SparseLU::SparseLU symbolic factorization failed"); |
| 5164 | 324 |
| 5322 | 325 UMFPACK_DNAME (report_status) (control, status); |
| 326 UMFPACK_DNAME (report_info) (control, info); | |
| 5164 | 327 |
| 5322 | 328 UMFPACK_DNAME (free_symbolic) (&Symbolic) ; |
| 5164 | 329 } |
| 330 else | |
| 331 { | |
| 5322 | 332 UMFPACK_DNAME (report_symbolic) (Symbolic, control); |
| 5164 | 333 |
| 5282 | 334 void *Numeric; |
| 5322 | 335 status = UMFPACK_DNAME (numeric) (Ap, Ai, Ax, Symbolic, |
| 336 &Numeric, control, info) ; | |
| 337 UMFPACK_DNAME (free_symbolic) (&Symbolic) ; | |
| 5282 | 338 |
| 339 cond = Info (UMFPACK_RCOND); | |
| 340 | |
| 5164 | 341 if (status < 0) |
| 342 { | |
| 343 (*current_liboctave_error_handler) | |
| 5282 | 344 ("SparseLU::SparseLU numeric factorization failed"); |
| 5164 | 345 |
| 5322 | 346 UMFPACK_DNAME (report_status) (control, status); |
| 347 UMFPACK_DNAME (report_info) (control, info); | |
| 5164 | 348 |
| 5322 | 349 UMFPACK_DNAME (free_numeric) (&Numeric); |
| 5164 | 350 } |
| 351 else | |
| 352 { | |
| 5322 | 353 UMFPACK_DNAME (report_numeric) (Numeric, control); |
| 5164 | 354 |
| 5322 | 355 octave_idx_type lnz, unz, ignore1, ignore2, ignore3; |
| 356 status = UMFPACK_DNAME (get_lunz) (&lnz, &unz, &ignore1, &ignore2, | |
| 357 &ignore3, Numeric) ; | |
| 5282 | 358 |
| 359 if (status < 0) | |
| 5164 | 360 { |
| 361 (*current_liboctave_error_handler) | |
| 362 ("SparseLU::SparseLU extracting LU factors failed"); | |
| 363 | |
| 5322 | 364 UMFPACK_DNAME (report_status) (control, status); |
| 365 UMFPACK_DNAME (report_info) (control, info); | |
| 5282 | 366 |
| 5322 | 367 UMFPACK_DNAME (free_numeric) (&Numeric); |
| 5164 | 368 } |
| 369 else | |
| 370 { | |
| 5322 | 371 octave_idx_type n_inner = (nr < nc ? nr : nc); |
| 5282 | 372 |
| 373 if (lnz < 1) | |
| 5322 | 374 Lfact = SparseMatrix (n_inner, nr, |
| 375 static_cast<octave_idx_type> (1)); | |
| 5282 | 376 else |
| 5322 | 377 Lfact = SparseMatrix (n_inner, nr, lnz); |
| 5282 | 378 |
| 379 octave_idx_type *Ltp = Lfact.cidx (); | |
| 380 octave_idx_type *Ltj = Lfact.ridx (); | |
| 381 double *Ltx = Lfact.data (); | |
| 382 | |
| 383 if (unz < 1) | |
| 5322 | 384 Ufact = SparseMatrix (n_inner, nc, |
| 385 static_cast<octave_idx_type> (1)); | |
| 5282 | 386 else |
| 5322 | 387 Ufact = SparseMatrix (n_inner, nc, unz); |
| 5282 | 388 |
| 389 octave_idx_type *Up = Ufact.cidx (); | |
| 390 octave_idx_type *Uj = Ufact.ridx (); | |
| 391 double *Ux = Ufact.data (); | |
| 392 | |
|
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
393 Rfact = SparseMatrix (nr, nr, nr); |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
394 for (octave_idx_type i = 0; i < nr; i++) |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
395 { |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
396 Rfact.xridx (i) = i; |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
397 Rfact.xcidx (i) = i; |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
398 } |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
399 Rfact.xcidx (nr) = nr; |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
400 double *Rx = Rfact.data (); |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
401 |
| 5282 | 402 P.resize (nr); |
| 5322 | 403 octave_idx_type *p = P.fortran_vec (); |
| 5282 | 404 |
| 405 Q.resize (nc); | |
| 5322 | 406 octave_idx_type *q = Q.fortran_vec (); |
| 5282 | 407 |
| 5322 | 408 octave_idx_type do_recip; |
| 409 status = UMFPACK_DNAME (get_numeric) (Ltp, Ltj, | |
| 410 Ltx, Up, Uj, Ux, p, q, | |
| 7520 | 411 0, &do_recip, |
|
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
412 Rx, Numeric) ; |
| 5282 | 413 |
| 5322 | 414 UMFPACK_DNAME (free_numeric) (&Numeric) ; |
| 5282 | 415 |
|
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
416 if (status < 0) |
| 5282 | 417 { |
| 418 (*current_liboctave_error_handler) | |
| 419 ("SparseLU::SparseLU extracting LU factors failed"); | |
| 420 | |
| 5322 | 421 UMFPACK_DNAME (report_status) (control, status); |
| 5282 | 422 } |
| 423 else | |
| 424 { | |
| 425 Lfact = Lfact.transpose (); | |
|
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
426 |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
427 if (do_recip) |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
428 for (octave_idx_type i = 0; i < nr; i++) |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
429 Rx[i] = 1.0 / Rx[i]; |
|
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
430 |
| 5322 | 431 UMFPACK_DNAME (report_matrix) (nr, n_inner, |
| 5282 | 432 Lfact.cidx (), |
| 433 Lfact.ridx (), | |
| 434 Lfact.data (), | |
| 435 1, control); | |
| 5322 | 436 UMFPACK_DNAME (report_matrix) (n_inner, nc, |
| 5282 | 437 Ufact.cidx (), |
| 438 Ufact.ridx (), | |
| 439 Ufact.data (), | |
| 440 1, control); | |
| 5322 | 441 UMFPACK_DNAME (report_perm) (nr, p, control); |
| 442 UMFPACK_DNAME (report_perm) (nc, q, control); | |
| 5282 | 443 } |
| 444 | |
| 5322 | 445 UMFPACK_DNAME (report_info) (control, info); |
| 5164 | 446 } |
| 447 } | |
| 448 } | |
| 5282 | 449 |
| 450 if (udiag) | |
| 451 (*current_liboctave_error_handler) | |
| 452 ("Option udiag of incomplete LU not implemented"); | |
| 5164 | 453 } |
| 5203 | 454 #else |
| 455 (*current_liboctave_error_handler) ("UMFPACK not installed"); | |
| 456 #endif | |
| 5164 | 457 } |
| 458 | |
| 459 /* | |
| 460 ;;; Local Variables: *** | |
| 461 ;;; mode: C++ *** | |
| 462 ;;; End: *** | |
| 463 */ | |
| 464 |