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