Mercurial > octave-nkf

comparison liboctave/SparseCmplxLU.cc @ 11586:12df7854fa7c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
strip trailing whitespace from source files
author John W. Eaton <jwe@octave.org>
date Thu, 20 Jan 2011 17:24:59 -0500
parents a83bad07f7e3
children 72c96de7a403
comparison
equal deleted inserted replaced
11585:1473d0cf86d2 11586:12df7854fa7c
40 40
41 template class sparse_base_lu <SparseComplexMatrix, Complex, SparseMatrix, double>; 41 template class sparse_base_lu <SparseComplexMatrix, Complex, SparseMatrix, double>;
42 42
43 #include "oct-sparse.h" 43 #include "oct-sparse.h"
44 44
45 SparseComplexLU::SparseComplexLU (const SparseComplexMatrix& a, 45 SparseComplexLU::SparseComplexLU (const SparseComplexMatrix& a,
46 const Matrix& piv_thres, bool scale) 46 const Matrix& piv_thres, bool scale)
47 { 47 {
48 #ifdef HAVE_UMFPACK 48 #ifdef HAVE_UMFPACK
49 octave_idx_type nr = a.rows (); 49 octave_idx_type nr = a.rows ();
50 octave_idx_type nc = a.cols (); 50 octave_idx_type nc = a.cols ();
80 // Set whether we are allowed to modify Q or not 80 // Set whether we are allowed to modify Q or not
81 tmp = octave_sparse_params::get_key ("autoamd"); 81 tmp = octave_sparse_params::get_key ("autoamd");
82 if (!xisnan (tmp)) 82 if (!xisnan (tmp))
83 Control (UMFPACK_FIXQ) = tmp; 83 Control (UMFPACK_FIXQ) = tmp;
84 84
85 // Turn-off UMFPACK scaling for LU 85 // Turn-off UMFPACK scaling for LU
86 if (scale) 86 if (scale)
87 Control (UMFPACK_SCALE) = UMFPACK_SCALE_SUM; 87 Control (UMFPACK_SCALE) = UMFPACK_SCALE_SUM;
88 else 88 else
89 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE; 89 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
90 90
99 0, 1, control); 99 0, 1, control);
100 100
101 void *Symbolic; 101 void *Symbolic;
102 Matrix Info (1, UMFPACK_INFO); 102 Matrix Info (1, UMFPACK_INFO);
103 double *info = Info.fortran_vec (); 103 double *info = Info.fortran_vec ();
104 int status = UMFPACK_ZNAME (qsymbolic) (nr, nc, Ap, Ai, 104 int status = UMFPACK_ZNAME (qsymbolic) (nr, nc, Ap, Ai,
105 reinterpret_cast<const double *> (Ax), 105 reinterpret_cast<const double *> (Ax),
106 0, 0, 106 0, 0,
107 &Symbolic, control, info); 107 &Symbolic, control, info);
108 108
109 if (status < 0) 109 if (status < 0)
110 { 110 {
111 (*current_liboctave_error_handler) 111 (*current_liboctave_error_handler)
112 ("SparseComplexLU::SparseComplexLU symbolic factorization failed"); 112 ("SparseComplexLU::SparseComplexLU symbolic factorization failed");
113 113
114 UMFPACK_ZNAME (report_status) (control, status); 114 UMFPACK_ZNAME (report_status) (control, status);
115 UMFPACK_ZNAME (report_info) (control, info); 115 UMFPACK_ZNAME (report_info) (control, info);
116 116
129 129
130 cond = Info (UMFPACK_RCOND); 130 cond = Info (UMFPACK_RCOND);
131 131
132 if (status < 0) 132 if (status < 0)
133 { 133 {
134 (*current_liboctave_error_handler) 134 (*current_liboctave_error_handler)
135 ("SparseComplexLU::SparseComplexLU numeric factorization failed"); 135 ("SparseComplexLU::SparseComplexLU numeric factorization failed");
136 136
137 UMFPACK_ZNAME (report_status) (control, status); 137 UMFPACK_ZNAME (report_status) (control, status);
138 UMFPACK_ZNAME (report_info) (control, info); 138 UMFPACK_ZNAME (report_info) (control, info);
139 139
144 UMFPACK_ZNAME (report_numeric) (Numeric, control); 144 UMFPACK_ZNAME (report_numeric) (Numeric, control);
145 145
146 octave_idx_type lnz, unz, ignore1, ignore2, ignore3; 146 octave_idx_type lnz, unz, ignore1, ignore2, ignore3;
147 status = UMFPACK_ZNAME (get_lunz) (&lnz, &unz, &ignore1, 147 status = UMFPACK_ZNAME (get_lunz) (&lnz, &unz, &ignore1,
148 &ignore2, &ignore3, Numeric) ; 148 &ignore2, &ignore3, Numeric) ;
149 149
150 if (status < 0) 150 if (status < 0)
151 { 151 {
152 (*current_liboctave_error_handler) 152 (*current_liboctave_error_handler)
153 ("SparseComplexLU::SparseComplexLU extracting LU factors failed"); 153 ("SparseComplexLU::SparseComplexLU extracting LU factors failed");
154 154
155 UMFPACK_ZNAME (report_status) (control, status); 155 UMFPACK_ZNAME (report_status) (control, status);
156 UMFPACK_ZNAME (report_info) (control, info); 156 UMFPACK_ZNAME (report_info) (control, info);
157 157
178 Ufact = SparseComplexMatrix (n_inner, nc, unz); 178 Ufact = SparseComplexMatrix (n_inner, nc, unz);
179 179
180 octave_idx_type *Up = Ufact.cidx (); 180 octave_idx_type *Up = Ufact.cidx ();
181 octave_idx_type *Uj = Ufact.ridx (); 181 octave_idx_type *Uj = Ufact.ridx ();
182 Complex *Ux = Ufact.data (); 182 Complex *Ux = Ufact.data ();
183 183
184 Rfact = SparseMatrix (nr, nr, nr); 184 Rfact = SparseMatrix (nr, nr, nr);
185 for (octave_idx_type i = 0; i < nr; i++) 185 for (octave_idx_type i = 0; i < nr; i++)
186 { 186 {
187 Rfact.xridx (i) = i; 187 Rfact.xridx (i) = i;
188 Rfact.xcidx (i) = i; 188 Rfact.xcidx (i) = i;
206 206
207 UMFPACK_ZNAME (free_numeric) (&Numeric) ; 207 UMFPACK_ZNAME (free_numeric) (&Numeric) ;
208 208
209 if (status < 0) 209 if (status < 0)
210 { 210 {
211 (*current_liboctave_error_handler) 211 (*current_liboctave_error_handler)
212 ("SparseComplexLU::SparseComplexLU extracting LU factors failed"); 212 ("SparseComplexLU::SparseComplexLU extracting LU factors failed");
213 213
214 UMFPACK_ZNAME (report_status) (control, status); 214 UMFPACK_ZNAME (report_status) (control, status);
215 } 215 }
216 else 216 else
220 if (do_recip) 220 if (do_recip)
221 for (octave_idx_type i = 0; i < nr; i++) 221 for (octave_idx_type i = 0; i < nr; i++)
222 Rx[i] = 1.0 / Rx[i]; 222 Rx[i] = 1.0 / Rx[i];
223 223
224 UMFPACK_ZNAME (report_matrix) (nr, n_inner, 224 UMFPACK_ZNAME (report_matrix) (nr, n_inner,
225 Lfact.cidx (), Lfact.ridx (), 225 Lfact.cidx (), Lfact.ridx (),
226 reinterpret_cast<double *> (Lfact.data()), 226 reinterpret_cast<double *> (Lfact.data()),
227 0, 1, control); 227 0, 1, control);
228 228
229 UMFPACK_ZNAME (report_matrix) (n_inner, nc, 229 UMFPACK_ZNAME (report_matrix) (n_inner, nc,
230 Ufact.cidx (), Ufact.ridx (), 230 Ufact.cidx (), Ufact.ridx (),
231 reinterpret_cast<double *> (Ufact.data()), 231 reinterpret_cast<double *> (Ufact.data()),
232 0, 1, control); 232 0, 1, control);
233 UMFPACK_ZNAME (report_perm) (nr, p, control); 233 UMFPACK_ZNAME (report_perm) (nr, p, control);
234 UMFPACK_ZNAME (report_perm) (nc, q, control); 234 UMFPACK_ZNAME (report_perm) (nc, q, control);
235 } 235 }
236 236
241 #else 241 #else
242 (*current_liboctave_error_handler) ("UMFPACK not installed"); 242 (*current_liboctave_error_handler) ("UMFPACK not installed");
243 #endif 243 #endif
244 } 244 }
245 245
246 SparseComplexLU::SparseComplexLU (const SparseComplexMatrix& a, 246 SparseComplexLU::SparseComplexLU (const SparseComplexMatrix& a,
247 const ColumnVector& Qinit, 247 const ColumnVector& Qinit,
248 const Matrix& piv_thres, bool scale, 248 const Matrix& piv_thres, bool scale,
249 bool FixedQ, double droptol, 249 bool FixedQ, double droptol,
250 bool milu, bool udiag) 250 bool milu, bool udiag)
251 { 251 {
252 #ifdef HAVE_UMFPACK 252 #ifdef HAVE_UMFPACK
253 if (milu) 253 if (milu)
254 (*current_liboctave_error_handler) 254 (*current_liboctave_error_handler)
255 ("Modified incomplete LU not implemented"); 255 ("Modified incomplete LU not implemented");
256 else 256 else
257 { 257 {
258 octave_idx_type nr = a.rows (); 258 octave_idx_type nr = a.rows ();
259 octave_idx_type nc = a.cols (); 259 octave_idx_type nc = a.cols ();
260 260
297 tmp = octave_sparse_params::get_key ("autoamd"); 297 tmp = octave_sparse_params::get_key ("autoamd");
298 if (!xisnan (tmp)) 298 if (!xisnan (tmp))
299 Control (UMFPACK_FIXQ) = tmp; 299 Control (UMFPACK_FIXQ) = tmp;
300 } 300 }
301 301
302 // Turn-off UMFPACK scaling for LU 302 // Turn-off UMFPACK scaling for LU
303 if (scale) 303 if (scale)
304 Control (UMFPACK_SCALE) = UMFPACK_SCALE_SUM; 304 Control (UMFPACK_SCALE) = UMFPACK_SCALE_SUM;
305 else 305 else
306 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE; 306 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
307 307
309 309
310 const octave_idx_type *Ap = a.cidx (); 310 const octave_idx_type *Ap = a.cidx ();
311 const octave_idx_type *Ai = a.ridx (); 311 const octave_idx_type *Ai = a.ridx ();
312 const Complex *Ax = a.data (); 312 const Complex *Ax = a.data ();
313 313
314 UMFPACK_ZNAME (report_matrix) (nr, nc, Ap, Ai, 314 UMFPACK_ZNAME (report_matrix) (nr, nc, Ap, Ai,
315 reinterpret_cast<const double *> (Ax), 0, 315 reinterpret_cast<const double *> (Ax), 0,
316 1, control); 316 1, control);
317 317
318 void *Symbolic; 318 void *Symbolic;
319 Matrix Info (1, UMFPACK_INFO); 319 Matrix Info (1, UMFPACK_INFO);
326 OCTAVE_LOCAL_BUFFER (octave_idx_type, qinit, nc); 326 OCTAVE_LOCAL_BUFFER (octave_idx_type, qinit, nc);
327 327
328 for (octave_idx_type i = 0; i < nc; i++) 328 for (octave_idx_type i = 0; i < nc; i++)
329 qinit [i] = static_cast<octave_idx_type> (Qinit (i)); 329 qinit [i] = static_cast<octave_idx_type> (Qinit (i));
330 330
331 status = UMFPACK_ZNAME (qsymbolic) (nr, nc, Ap, Ai, 331 status = UMFPACK_ZNAME (qsymbolic) (nr, nc, Ap, Ai,
332 reinterpret_cast<const double *> (Ax), 332 reinterpret_cast<const double *> (Ax),
333 0, qinit, &Symbolic, control, 333 0, qinit, &Symbolic, control,
334 info); 334 info);
335 } while (0); 335 } while (0);
336 336
337 if (status < 0) 337 if (status < 0)
338 { 338 {
339 (*current_liboctave_error_handler) 339 (*current_liboctave_error_handler)
340 ("SparseComplexLU::SparseComplexLU symbolic factorization failed"); 340 ("SparseComplexLU::SparseComplexLU symbolic factorization failed");
341 341
342 UMFPACK_ZNAME (report_status) (control, status); 342 UMFPACK_ZNAME (report_status) (control, status);
343 UMFPACK_ZNAME (report_info) (control, info); 343 UMFPACK_ZNAME (report_info) (control, info);
344 344
347 else 347 else
348 { 348 {
349 UMFPACK_ZNAME (report_symbolic) (Symbolic, control); 349 UMFPACK_ZNAME (report_symbolic) (Symbolic, control);
350 350
351 void *Numeric; 351 void *Numeric;
352 status = UMFPACK_ZNAME (numeric) (Ap, Ai, 352 status = UMFPACK_ZNAME (numeric) (Ap, Ai,
353 reinterpret_cast<const double *> (Ax), 0, 353 reinterpret_cast<const double *> (Ax), 0,
354 Symbolic, &Numeric, control, info) ; 354 Symbolic, &Numeric, control, info) ;
355 UMFPACK_ZNAME (free_symbolic) (&Symbolic) ; 355 UMFPACK_ZNAME (free_symbolic) (&Symbolic) ;
356 356
357 cond = Info (UMFPACK_RCOND); 357 cond = Info (UMFPACK_RCOND);
358 358
359 if (status < 0) 359 if (status < 0)
360 { 360 {
361 (*current_liboctave_error_handler) 361 (*current_liboctave_error_handler)
362 ("SparseComplexLU::SparseComplexLU numeric factorization failed"); 362 ("SparseComplexLU::SparseComplexLU numeric factorization failed");
363 363
364 UMFPACK_ZNAME (report_status) (control, status); 364 UMFPACK_ZNAME (report_status) (control, status);
365 UMFPACK_ZNAME (report_info) (control, info); 365 UMFPACK_ZNAME (report_info) (control, info);
366 366
371 UMFPACK_ZNAME (report_numeric) (Numeric, control); 371 UMFPACK_ZNAME (report_numeric) (Numeric, control);
372 372
373 octave_idx_type lnz, unz, ignore1, ignore2, ignore3; 373 octave_idx_type lnz, unz, ignore1, ignore2, ignore3;
374 status = UMFPACK_ZNAME (get_lunz) (&lnz, &unz, 374 status = UMFPACK_ZNAME (get_lunz) (&lnz, &unz,
375 &ignore1, &ignore2, &ignore3, Numeric); 375 &ignore1, &ignore2, &ignore3, Numeric);
376 376
377 if (status < 0) 377 if (status < 0)
378 { 378 {
379 (*current_liboctave_error_handler) 379 (*current_liboctave_error_handler)
380 ("SparseComplexLU::SparseComplexLU extracting LU factors failed"); 380 ("SparseComplexLU::SparseComplexLU extracting LU factors failed");
381 381
382 UMFPACK_ZNAME (report_status) (control, status); 382 UMFPACK_ZNAME (report_status) (control, status);
383 UMFPACK_ZNAME (report_info) (control, info); 383 UMFPACK_ZNAME (report_info) (control, info);
384 384
405 Ufact = SparseComplexMatrix (n_inner, nc, unz); 405 Ufact = SparseComplexMatrix (n_inner, nc, unz);
406 406
407 octave_idx_type *Up = Ufact.cidx (); 407 octave_idx_type *Up = Ufact.cidx ();
408 octave_idx_type *Uj = Ufact.ridx (); 408 octave_idx_type *Uj = Ufact.ridx ();
409 Complex *Ux = Ufact.data (); 409 Complex *Ux = Ufact.data ();
410 410
411 Rfact = SparseMatrix (nr, nr, nr); 411 Rfact = SparseMatrix (nr, nr, nr);
412 for (octave_idx_type i = 0; i < nr; i++) 412 for (octave_idx_type i = 0; i < nr; i++)
413 { 413 {
414 Rfact.xridx (i) = i; 414 Rfact.xridx (i) = i;
415 Rfact.xcidx (i) = i; 415 Rfact.xcidx (i) = i;
422 422
423 Q.resize (dim_vector (nc, 1)); 423 Q.resize (dim_vector (nc, 1));
424 octave_idx_type *q = Q.fortran_vec (); 424 octave_idx_type *q = Q.fortran_vec ();
425 425
426 octave_idx_type do_recip; 426 octave_idx_type do_recip;
427 status = 427 status =
428 UMFPACK_ZNAME (get_numeric) (Ltp, Ltj, 428 UMFPACK_ZNAME (get_numeric) (Ltp, Ltj,
429 reinterpret_cast<double *> (Ltx), 429 reinterpret_cast<double *> (Ltx),
430 0, Up, Uj, 430 0, Up, Uj,
431 reinterpret_cast<double *> (Ux), 431 reinterpret_cast<double *> (Ux),
432 0, p, q, 0, 0, 432 0, p, q, 0, 0,
433 &do_recip, Rx, Numeric) ; 433 &do_recip, Rx, Numeric) ;
434 434
435 UMFPACK_ZNAME (free_numeric) (&Numeric) ; 435 UMFPACK_ZNAME (free_numeric) (&Numeric) ;
436 436
437 if (status < 0) 437 if (status < 0)
438 { 438 {
439 (*current_liboctave_error_handler) 439 (*current_liboctave_error_handler)
440 ("SparseComplexLU::SparseComplexLU extracting LU factors failed"); 440 ("SparseComplexLU::SparseComplexLU extracting LU factors failed");
441 441
442 UMFPACK_ZNAME (report_status) (control, status); 442 UMFPACK_ZNAME (report_status) (control, status);
443 } 443 }
444 else 444 else
447 447
448 if (do_recip) 448 if (do_recip)
449 for (octave_idx_type i = 0; i < nr; i++) 449 for (octave_idx_type i = 0; i < nr; i++)
450 Rx[i] = 1.0 / Rx[i]; 450 Rx[i] = 1.0 / Rx[i];
451 451
452 UMFPACK_ZNAME (report_matrix) (nr, n_inner, 452 UMFPACK_ZNAME (report_matrix) (nr, n_inner,
453 Lfact.cidx (), 453 Lfact.cidx (),
454 Lfact.ridx (), 454 Lfact.ridx (),
455 reinterpret_cast<double *> (Lfact.data()), 455 reinterpret_cast<double *> (Lfact.data()),
456 0, 1, control); 456 0, 1, control);
457 457
458 UMFPACK_ZNAME (report_matrix) (n_inner, nc, 458 UMFPACK_ZNAME (report_matrix) (n_inner, nc,
459 Ufact.cidx (), 459 Ufact.cidx (),
460 Ufact.ridx (), 460 Ufact.ridx (),
461 reinterpret_cast<double *> (Ufact.data()), 461 reinterpret_cast<double *> (Ufact.data()),
462 0, 1, control); 462 0, 1, control);
463 UMFPACK_ZNAME (report_perm) (nr, p, control); 463 UMFPACK_ZNAME (report_perm) (nr, p, control);
464 UMFPACK_ZNAME (report_perm) (nc, q, control); 464 UMFPACK_ZNAME (report_perm) (nc, q, control);
465 } 465 }
466 466
468 } 468 }
469 } 469 }
470 } 470 }
471 471
472 if (udiag) 472 if (udiag)
473 (*current_liboctave_error_handler) 473 (*current_liboctave_error_handler)
474 ("Option udiag of incomplete LU not implemented"); 474 ("Option udiag of incomplete LU not implemented");
475 } 475 }
476 #else 476 #else
477 (*current_liboctave_error_handler) ("UMFPACK not installed"); 477 (*current_liboctave_error_handler) ("UMFPACK not installed");
478 #endif 478 #endif
479 } 479 }