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