Mercurial > octave-nkf

comparison liboctave/SparseCmplxLU.cc @ 7515:f3c00dc0912b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Eliminate the rest of the dispatched sparse functions
author David Bateman <dbateman@free.fr>
date Fri, 22 Feb 2008 15:50:51 +0100
parents a1dbe9d80eee
children b166043585a8
comparison
equal deleted inserted replaced
7514:4f6a73fd8df9 7515:f3c00dc0912b
40 template class sparse_base_lu <SparseComplexMatrix, Complex, SparseMatrix, double>; 40 template class sparse_base_lu <SparseComplexMatrix, Complex, SparseMatrix, double>;
41 41
42 #include "oct-sparse.h" 42 #include "oct-sparse.h"
43 43
44 SparseComplexLU::SparseComplexLU (const SparseComplexMatrix& a, 44 SparseComplexLU::SparseComplexLU (const SparseComplexMatrix& a,
45 double piv_thres) 45 const Matrix& piv_thres, bool scale)
46 { 46 {
47 #ifdef HAVE_UMFPACK 47 #ifdef HAVE_UMFPACK
48 octave_idx_type nr = a.rows (); 48 octave_idx_type nr = a.rows ();
49 octave_idx_type nc = a.cols (); 49 octave_idx_type nc = a.cols ();
50 50
54 UMFPACK_ZNAME (defaults) (control); 54 UMFPACK_ZNAME (defaults) (control);
55 55
56 double tmp = octave_sparse_params::get_key ("spumoni"); 56 double tmp = octave_sparse_params::get_key ("spumoni");
57 if (!xisnan (tmp)) 57 if (!xisnan (tmp))
58 Control (UMFPACK_PRL) = tmp; 58 Control (UMFPACK_PRL) = tmp;
59 if (piv_thres >= 0.) 59 if (piv_thres.nelem() != 2)
60 { 60 {
61 piv_thres = (piv_thres > 1. ? 1. : piv_thres); 61 tmp = (piv_thres (0) > 1. ? 1. : piv_thres (0));
62 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = piv_thres; 62 if (!xisnan (tmp))
63 Control (UMFPACK_PIVOT_TOLERANCE) = piv_thres; 63 Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
64 tmp = (piv_thres (1) > 1. ? 1. : piv_thres (1));
65 if (!xisnan (tmp))
66 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
64 } 67 }
65 else 68 else
66 { 69 {
67 tmp = octave_sparse_params::get_key ("piv_tol"); 70 tmp = octave_sparse_params::get_key ("piv_tol");
68 if (!xisnan (tmp)) 71 if (!xisnan (tmp))
69 { 72 Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
73
74 tmp = octave_sparse_params::get_key ("sym_tol");
75 if (!xisnan (tmp))
70 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp; 76 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
71 Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
72 }
73 } 77 }
74 78
75 // Set whether we are allowed to modify Q or not 79 // Set whether we are allowed to modify Q or not
76 tmp = octave_sparse_params::get_key ("autoamd"); 80 tmp = octave_sparse_params::get_key ("autoamd");
77 if (!xisnan (tmp)) 81 if (!xisnan (tmp))
78 Control (UMFPACK_FIXQ) = tmp; 82 Control (UMFPACK_FIXQ) = tmp;
79 83
80 // Turn-off UMFPACK scaling for LU 84 // Turn-off UMFPACK scaling for LU
81 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE; 85 if (scale)
86 Control (UMFPACK_SCALE) = UMFPACK_SCALE_SUM;
87 else
88 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
82 89
83 UMFPACK_ZNAME (report_control) (control); 90 UMFPACK_ZNAME (report_control) (control);
84 91
85 const octave_idx_type *Ap = a.cidx (); 92 const octave_idx_type *Ap = a.cidx ();
86 const octave_idx_type *Ai = a.ridx (); 93 const octave_idx_type *Ai = a.ridx ();
171 178
172 octave_idx_type *Up = Ufact.cidx (); 179 octave_idx_type *Up = Ufact.cidx ();
173 octave_idx_type *Uj = Ufact.ridx (); 180 octave_idx_type *Uj = Ufact.ridx ();
174 Complex *Ux = Ufact.data (); 181 Complex *Ux = Ufact.data ();
175 182
183 Rfact = SparseMatrix (nr, nr, nr);
184 for (octave_idx_type i = 0; i < nr; i++)
185 {
186 Rfact.xridx (i) = i;
187 Rfact.xcidx (i) = i;
188 }
189 Rfact.xcidx (nr) = nr;
190 double *Rx = Rfact.data ();
191
176 P.resize (nr); 192 P.resize (nr);
177 octave_idx_type *p = P.fortran_vec (); 193 octave_idx_type *p = P.fortran_vec ();
178 194
179 Q.resize (nc); 195 Q.resize (nc);
180 octave_idx_type *q = Q.fortran_vec (); 196 octave_idx_type *q = Q.fortran_vec ();
183 status = UMFPACK_ZNAME (get_numeric) (Ltp, Ltj, 199 status = UMFPACK_ZNAME (get_numeric) (Ltp, Ltj,
184 reinterpret_cast<double *> (Ltx), 200 reinterpret_cast<double *> (Ltx),
185 NULL, Up, Uj, 201 NULL, Up, Uj,
186 reinterpret_cast <double *> (Ux), 202 reinterpret_cast <double *> (Ux),
187 NULL, p, q, NULL, NULL, 203 NULL, p, q, NULL, NULL,
188 &do_recip, NULL, Numeric); 204 &do_recip, Rx, Numeric);
189 205
190 UMFPACK_ZNAME (free_numeric) (&Numeric) ; 206 UMFPACK_ZNAME (free_numeric) (&Numeric) ;
191 207
192 if (status < 0 || do_recip) 208 if (status < 0)
193 { 209 {
194 (*current_liboctave_error_handler) 210 (*current_liboctave_error_handler)
195 ("SparseComplexLU::SparseComplexLU extracting LU factors failed"); 211 ("SparseComplexLU::SparseComplexLU extracting LU factors failed");
196 212
197 UMFPACK_ZNAME (report_status) (control, status); 213 UMFPACK_ZNAME (report_status) (control, status);
198 } 214 }
199 else 215 else
200 { 216 {
201 Lfact = Lfact.transpose (); 217 Lfact = Lfact.transpose ();
218
219 if (do_recip)
220 for (octave_idx_type i = 0; i < nr; i++)
221 Rx[i] = 1.0 / Rx[i];
202 222
203 UMFPACK_ZNAME (report_matrix) (nr, n_inner, 223 UMFPACK_ZNAME (report_matrix) (nr, n_inner,
204 Lfact.cidx (), Lfact.ridx (), 224 Lfact.cidx (), Lfact.ridx (),
205 reinterpret_cast<double *> (Lfact.data()), 225 reinterpret_cast<double *> (Lfact.data()),
206 NULL, 1, control); 226 NULL, 1, control);
222 #endif 242 #endif
223 } 243 }
224 244
225 SparseComplexLU::SparseComplexLU (const SparseComplexMatrix& a, 245 SparseComplexLU::SparseComplexLU (const SparseComplexMatrix& a,
226 const ColumnVector& Qinit, 246 const ColumnVector& Qinit,
227 double piv_thres, bool FixedQ, 247 const Matrix& piv_thres, bool scale,
228 double droptol, bool milu, bool udiag) 248 bool FixedQ, double droptol,
249 bool milu, bool udiag)
229 { 250 {
230 #ifdef HAVE_UMFPACK 251 #ifdef HAVE_UMFPACK
231 if (milu) 252 if (milu)
232 (*current_liboctave_error_handler) 253 (*current_liboctave_error_handler)
233 ("Modified incomplete LU not implemented"); 254 ("Modified incomplete LU not implemented");
242 UMFPACK_ZNAME (defaults) (control); 263 UMFPACK_ZNAME (defaults) (control);
243 264
244 double tmp = octave_sparse_params::get_key ("spumoni"); 265 double tmp = octave_sparse_params::get_key ("spumoni");
245 if (!xisnan (tmp)) 266 if (!xisnan (tmp))
246 Control (UMFPACK_PRL) = tmp; 267 Control (UMFPACK_PRL) = tmp;
247 if (piv_thres >= 0.) 268 if (piv_thres.nelem() != 2)
248 { 269 {
249 piv_thres = (piv_thres > 1. ? 1. : piv_thres); 270 tmp = (piv_thres (0) > 1. ? 1. : piv_thres (0));
250 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = piv_thres; 271 if (!xisnan (tmp))
251 Control (UMFPACK_PIVOT_TOLERANCE) = piv_thres; 272 Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
273 tmp = (piv_thres (1) > 1. ? 1. : piv_thres (1));
274 if (!xisnan (tmp))
275 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
252 } 276 }
253 else 277 else
254 { 278 {
255 tmp = octave_sparse_params::get_key ("piv_tol"); 279 tmp = octave_sparse_params::get_key ("piv_tol");
256 if (!xisnan (tmp)) 280 if (!xisnan (tmp))
257 { 281 Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
258 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp; 282
259 Control (UMFPACK_PIVOT_TOLERANCE) = tmp; 283 tmp = octave_sparse_params::get_key ("sym_tol");
260 } 284 if (!xisnan (tmp))
285 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
261 } 286 }
262 287
263 if (droptol >= 0.) 288 if (droptol >= 0.)
264 Control (UMFPACK_DROPTOL) = droptol; 289 Control (UMFPACK_DROPTOL) = droptol;
265 290
272 if (!xisnan (tmp)) 297 if (!xisnan (tmp))
273 Control (UMFPACK_FIXQ) = tmp; 298 Control (UMFPACK_FIXQ) = tmp;
274 } 299 }
275 300
276 // Turn-off UMFPACK scaling for LU 301 // Turn-off UMFPACK scaling for LU
277 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE; 302 if (scale)
303 Control (UMFPACK_SCALE) = UMFPACK_SCALE_SUM;
304 else
305 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
278 306
279 UMFPACK_ZNAME (report_control) (control); 307 UMFPACK_ZNAME (report_control) (control);
280 308
281 const octave_idx_type *Ap = a.cidx (); 309 const octave_idx_type *Ap = a.cidx ();
282 const octave_idx_type *Ai = a.ridx (); 310 const octave_idx_type *Ai = a.ridx ();
377 405
378 octave_idx_type *Up = Ufact.cidx (); 406 octave_idx_type *Up = Ufact.cidx ();
379 octave_idx_type *Uj = Ufact.ridx (); 407 octave_idx_type *Uj = Ufact.ridx ();
380 Complex *Ux = Ufact.data (); 408 Complex *Ux = Ufact.data ();
381 409
410 Rfact = SparseMatrix (nr, nr, nr);
411 for (octave_idx_type i = 0; i < nr; i++)
412 {
413 Rfact.xridx (i) = i;
414 Rfact.xcidx (i) = i;
415 }
416 Rfact.xcidx (nr) = nr;
417 double *Rx = Rfact.data ();
418
382 P.resize (nr); 419 P.resize (nr);
383 octave_idx_type *p = P.fortran_vec (); 420 octave_idx_type *p = P.fortran_vec ();
384 421
385 Q.resize (nc); 422 Q.resize (nc);
386 octave_idx_type *q = Q.fortran_vec (); 423 octave_idx_type *q = Q.fortran_vec ();
390 UMFPACK_ZNAME (get_numeric) (Ltp, Ltj, 427 UMFPACK_ZNAME (get_numeric) (Ltp, Ltj,
391 reinterpret_cast<double *> (Ltx), 428 reinterpret_cast<double *> (Ltx),
392 NULL, Up, Uj, 429 NULL, Up, Uj,
393 reinterpret_cast<double *> (Ux), 430 reinterpret_cast<double *> (Ux),
394 NULL, p, q, NULL, NULL, 431 NULL, p, q, NULL, NULL,
395 &do_recip, NULL, Numeric) ; 432 &do_recip, Rx, Numeric) ;
396 433
397 UMFPACK_ZNAME (free_numeric) (&Numeric) ; 434 UMFPACK_ZNAME (free_numeric) (&Numeric) ;
398 435
399 if (status < 0 || do_recip) 436 if (status < 0)
400 { 437 {
401 (*current_liboctave_error_handler) 438 (*current_liboctave_error_handler)
402 ("SparseComplexLU::SparseComplexLU extracting LU factors failed"); 439 ("SparseComplexLU::SparseComplexLU extracting LU factors failed");
403 440
404 UMFPACK_ZNAME (report_status) (control, 441 UMFPACK_ZNAME (report_status) (control, status);
405 status);
406 } 442 }
407 else 443 else
408 { 444 {
409 Lfact = Lfact.transpose (); 445 Lfact = Lfact.transpose ();
446
447 if (do_recip)
448 for (octave_idx_type i = 0; i < nr; i++)
449 Rx[i] = 1.0 / Rx[i];
410 450
411 UMFPACK_ZNAME (report_matrix) (nr, n_inner, 451 UMFPACK_ZNAME (report_matrix) (nr, n_inner,
412 Lfact.cidx (), 452 Lfact.cidx (),
413 Lfact.ridx (), 453 Lfact.ridx (),
414 reinterpret_cast<double *> (Lfact.data()), 454 reinterpret_cast<double *> (Lfact.data()),