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