Mercurial > jwe > octave

comparison liboctave/numeric/sparse-lu.cc @ 21146:ea9c05014809

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
revamp sparse LU factorization classes * sparse-lu.h, sparse-lu.cc: Rename from sparse-base-lu.h and sparse-base-lu.cc, respectively. (class sparse_lu): Rename from sparse_base_lu. Incorporate code from SparseCmplxLU and SparsedbleLU classes into the sparse_lu template. * sparse-lu-inst.cc: New file. * SparseCmplxLU.cc, SparseCmplxLU.h, SparsedbleLU.cc, SparsedbleLU.h: Delete. * lu.cc, luinc.cc, CSparse.cc, dSparse.cc, eigs-base.cc: Change all uses of SparsedbleLU and SparseCmplxLU to use new sparse_lu template class. * liboctave/numeric/module.mk: Update.
author John W. Eaton <jwe@octave.org>
date Thu, 28 Jan 2016 00:15:33 -0500
parents liboctave/numeric/sparse-base-lu.cc@538b57866b90
children a10f60e13243
comparison
equal deleted inserted replaced
21145:307096fb67e1 21146:ea9c05014809
1 /*
2
3 Copyright (C) 2004-2015 David Bateman
4 Copyright (C) 1998-2004 Andy Adler
5
6 This file is part of Octave.
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
10 Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
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
19 along with Octave; see the file COPYING. If not, see
20 <http://www.gnu.org/licenses/>.
21
22 */
23
24 #ifdef HAVE_CONFIG_H
25 #include <config.h>
26 #endif
27
28 #include "CSparse.h"
29 #include "PermMatrix.h"
30 #include "dSparse.h"
31 #include "lo-error.h"
32 #include "oct-locbuf.h"
33 #include "oct-sparse.h"
34 #include "oct-spparms.h"
35 #include "sparse-lu.h"
36
37 // Wrappers for SuiteSparse (formerly UMFPACK) functions that have
38 // different names depending on the sparse matrix data type.
39 //
40 // All of these functions must be specialized to forward to the correct
41 // SuiteSparse functions.
42
43 template <typename T>
44 void
45 umfpack_defaults (double *Control);
46
47 template <typename T>
48 void
49 umfpack_free_numeric (void **Numeric);
50
51 template <typename T>
52 void
53 umfpack_free_symbolic (void **Symbolic);
54
55 template <typename T>
56 octave_idx_type
57 umfpack_get_lunz (octave_idx_type *lnz, octave_idx_type *unz, void *Numeric);
58
59 template <typename T>
60 octave_idx_type
61 umfpack_get_numeric (octave_idx_type *Lp, octave_idx_type *Lj,
62 T *Lx, // Or Lz_packed
63 octave_idx_type *Up, octave_idx_type *Ui,
64 T *Ux, // Or Uz_packed
65 octave_idx_type *p, octave_idx_type *q,
66 double *Dz_packed, octave_idx_type *do_recip,
67 double *Rs, void *Numeric);
68
69 template <typename T>
70 octave_idx_type
71 umfpack_numeric (const octave_idx_type *Ap, const octave_idx_type *Ai,
72 const T *Ax, // Or Az_packed
73 void *Symbolic, void **Numeric,
74 const double *Control, double *Info);
75
76 template <typename T>
77 octave_idx_type
78 umfpack_qsymbolic (octave_idx_type n_row, octave_idx_type n_col,
79 const octave_idx_type *Ap, const octave_idx_type *Ai,
80 const T *Ax, // Or Az_packed
81 const octave_idx_type *Qinit, void **Symbolic,
82 const double *Control, double *Info);
83
84 template <typename T>
85 void
86 umfpack_report_control (const double *Control);
87
88 template <typename T>
89 void
90 umfpack_report_info (const double *Control, const double *Info);
91
92 template <typename T>
93 void
94 umfpack_report_matrix (octave_idx_type n_row, octave_idx_type n_col,
95 const octave_idx_type *Ap, const octave_idx_type *Ai,
96 const T *Ax, // Or Az_packed
97 octave_idx_type col_form, const double *Control);
98
99 template <typename T>
100 void
101 umfpack_report_numeric (void *Numeric, const double *Control);
102
103 template <typename T>
104 void
105 umfpack_report_perm (octave_idx_type np, const octave_idx_type *Perm,
106 const double *Control);
107
108 template <typename T>
109 void
110 umfpack_report_status (double *Control, octave_idx_type status);
111
112 template <typename T>
113 void
114 umfpack_report_symbolic (void *Symbolic, const double *Control);
115
116 // SparseMatrix Specialization.
117
118 template <>
119 inline void
120 umfpack_defaults<double> (double *Control)
121 {
122 UMFPACK_DNAME (defaults) (Control);
123 }
124
125 template <>
126 inline void
127 umfpack_free_numeric<double> (void **Numeric)
128 {
129 UMFPACK_DNAME (free_numeric) (Numeric);
130 }
131
132 template <>
133 inline void
134 umfpack_free_symbolic<double> (void **Symbolic)
135 {
136 UMFPACK_DNAME (free_symbolic) (Symbolic);
137 }
138
139 template <>
140 inline octave_idx_type
141 umfpack_get_lunz<double>
142 (octave_idx_type *lnz, octave_idx_type *unz, void *Numeric)
143 {
144 octave_idx_type ignore1, ignore2, ignore3;
145
146 return UMFPACK_DNAME (get_lunz) (lnz, unz, &ignore1, &ignore2,
147 &ignore3, Numeric);
148 }
149
150 template <>
151 inline octave_idx_type
152 umfpack_get_numeric<double>
153 (octave_idx_type *Lp, octave_idx_type *Lj, double *Lx,
154 octave_idx_type *Up, octave_idx_type *Ui, double *Ux,
155 octave_idx_type *p, octave_idx_type *q, double *Dx,
156 octave_idx_type *do_recip, double *Rs, void *Numeric)
157 {
158 return UMFPACK_DNAME (get_numeric) (Lp, Lj, Lx, Up, Ui, Ux, p, q, Dx,
159 do_recip, Rs, Numeric);
160 }
161
162 template <>
163 inline octave_idx_type
164 umfpack_numeric<double>
165 (const octave_idx_type *Ap, const octave_idx_type *Ai,
166 const double *Ax, void *Symbolic, void **Numeric,
167 const double *Control, double *Info)
168 {
169 return UMFPACK_DNAME (numeric) (Ap, Ai, Ax, Symbolic, Numeric, Control,
170 Info);
171 }
172
173 template <>
174 inline octave_idx_type
175 umfpack_qsymbolic<double>
176 (octave_idx_type n_row, octave_idx_type n_col, const octave_idx_type *Ap,
177 const octave_idx_type *Ai, const double *Ax,
178 const octave_idx_type *Qinit, void **Symbolic,
179 const double *Control, double *Info)
180 {
181 return UMFPACK_DNAME (qsymbolic) (n_row, n_col, Ap, Ai, Ax, Qinit,
182 Symbolic, Control, Info);
183 }
184
185 template <>
186 inline void
187 umfpack_report_control<double> (const double *Control)
188 {
189 UMFPACK_DNAME (report_control) (Control);
190 }
191
192 template <>
193 inline void
194 umfpack_report_info<double> (const double *Control, const double *Info)
195 {
196 UMFPACK_DNAME (report_info) (Control, Info);
197 }
198
199 template <>
200 inline void
201 umfpack_report_matrix<double>
202 (octave_idx_type n_row, octave_idx_type n_col, const octave_idx_type *Ap,
203 const octave_idx_type *Ai, const double *Ax, octave_idx_type col_form,
204 const double *Control)
205 {
206 UMFPACK_DNAME (report_matrix) (n_row, n_col, Ap, Ai, Ax, col_form, Control);
207 }
208
209 template <>
210 inline void
211 umfpack_report_numeric<double> (void *Numeric, const double *Control)
212 {
213 UMFPACK_DNAME (report_numeric) (Numeric, Control);
214 }
215
216 template <>
217 inline void
218 umfpack_report_perm<double>
219 (octave_idx_type np, const octave_idx_type *Perm, const double *Control)
220 {
221 UMFPACK_DNAME (report_perm) (np, Perm, Control);
222 }
223
224 template <>
225 inline void
226 umfpack_report_status<double> (double *Control, octave_idx_type status)
227 {
228 UMFPACK_DNAME (report_status) (Control, status);
229 }
230
231 template <>
232 inline void
233 umfpack_report_symbolic<double> (void *Symbolic, const double *Control)
234 {
235 UMFPACK_DNAME (report_symbolic) (Symbolic, Control);
236 }
237
238 // SparseComplexMatrix specialization.
239
240 template <>
241 inline void
242 umfpack_defaults<Complex> (double *Control)
243 {
244 UMFPACK_ZNAME (defaults) (Control);
245 }
246
247 template <>
248 inline void
249 umfpack_free_numeric<Complex> (void **Numeric)
250 {
251 UMFPACK_ZNAME (free_numeric) (Numeric);
252 }
253
254 template <>
255 inline void
256 umfpack_free_symbolic<Complex> (void **Symbolic)
257 {
258 UMFPACK_ZNAME (free_symbolic) (Symbolic);
259 }
260
261 template <>
262 inline octave_idx_type
263 umfpack_get_lunz<Complex>
264 (octave_idx_type *lnz, octave_idx_type *unz, void *Numeric)
265 {
266 octave_idx_type ignore1, ignore2, ignore3;
267
268 return UMFPACK_ZNAME (get_lunz) (lnz, unz, &ignore1, &ignore2,
269 &ignore3, Numeric);
270 }
271
272 template <>
273 inline octave_idx_type
274 umfpack_get_numeric<Complex>
275 (octave_idx_type *Lp, octave_idx_type *Lj, Complex *Lz,
276 octave_idx_type *Up, octave_idx_type *Ui, Complex *Uz,
277 octave_idx_type *p, octave_idx_type *q, double *Dz,
278 octave_idx_type *do_recip, double *Rs, void *Numeric)
279 {
280 return UMFPACK_ZNAME (get_numeric) (Lp, Lj,
281 reinterpret_cast<double *> (Lz),
282 0, Up, Ui,
283 reinterpret_cast<double *> (Uz),
284 0, p, q,
285 reinterpret_cast<double *> (Dz),
286 0, do_recip, Rs, Numeric);
287 }
288
289 template <>
290 inline octave_idx_type
291 umfpack_numeric<Complex>
292 (const octave_idx_type *Ap, const octave_idx_type *Ai,
293 const Complex *Az, void *Symbolic, void **Numeric,
294 const double *Control, double *Info)
295 {
296 return UMFPACK_ZNAME (numeric) (Ap, Ai,
297 reinterpret_cast<const double *> (Az),
298 0, Symbolic, Numeric, Control, Info);
299 }
300
301
302 template <>
303 inline octave_idx_type
304 umfpack_qsymbolic<Complex>
305 (octave_idx_type n_row, octave_idx_type n_col,
306 const octave_idx_type *Ap, const octave_idx_type *Ai,
307 const Complex *Az, const octave_idx_type *Qinit,
308 void **Symbolic, const double *Control, double *Info)
309 {
310 return UMFPACK_ZNAME (qsymbolic) (n_row, n_col, Ap, Ai,
311 reinterpret_cast<const double *> (Az),
312 0, Qinit, Symbolic, Control, Info);
313 }
314
315
316 template <>
317 inline void
318 umfpack_report_control<Complex> (const double *Control)
319 {
320 UMFPACK_ZNAME (report_control) (Control);
321 }
322
323 template <>
324 inline void
325 umfpack_report_info<Complex> (const double *Control, const double *Info)
326 {
327 UMFPACK_ZNAME (report_info) (Control, Info);
328 }
329
330 template <>
331 inline void
332 umfpack_report_matrix<Complex>
333 (octave_idx_type n_row, octave_idx_type n_col,
334 const octave_idx_type *Ap, const octave_idx_type *Ai,
335 const Complex *Az, octave_idx_type col_form, const double *Control)
336 {
337 UMFPACK_ZNAME (report_matrix) (n_row, n_col, Ap, Ai,
338 reinterpret_cast<const double *> (Az),
339 0, col_form, Control);
340 }
341
342 template <>
343 inline void
344 umfpack_report_numeric<Complex> (void *Numeric, const double *Control)
345 {
346 UMFPACK_ZNAME (report_numeric) (Numeric, Control);
347 }
348
349 template <>
350 inline void
351 umfpack_report_perm<Complex>
352 (octave_idx_type np, const octave_idx_type *Perm, const double *Control)
353 {
354 UMFPACK_ZNAME (report_perm) (np, Perm, Control);
355 }
356
357 template <>
358 inline void
359 umfpack_report_status<Complex> (double *Control, octave_idx_type status)
360 {
361 UMFPACK_ZNAME (report_status) (Control, status);
362 }
363
364 template <>
365 inline void
366 umfpack_report_symbolic <Complex> (void *Symbolic, const double *Control)
367 {
368 UMFPACK_ZNAME (report_symbolic) (Symbolic, Control);
369 }
370
371 template <typename lu_type>
372 sparse_lu<lu_type>::sparse_lu (const lu_type& a, const Matrix& piv_thres,
373 bool scale)
374 : Lfact (), Ufact (), Rfact (), cond (0), P (), Q ()
375 {
376 #ifdef HAVE_UMFPACK
377 octave_idx_type nr = a.rows ();
378 octave_idx_type nc = a.cols ();
379
380 // Setup the control parameters
381 Matrix Control (UMFPACK_CONTROL, 1);
382 double *control = Control.fortran_vec ();
383 umfpack_defaults<lu_elt_type> (control);
384
385 double tmp = octave_sparse_params::get_key ("spumoni");
386 if (! xisnan (tmp))
387 Control (UMFPACK_PRL) = tmp;
388
389 if (piv_thres.numel () == 2)
390 {
391 tmp = (piv_thres (0) > 1. ? 1. : piv_thres (0));
392 if (! xisnan (tmp))
393 Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
394 tmp = (piv_thres (1) > 1. ? 1. : piv_thres (1));
395 if (! xisnan (tmp))
396 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
397 }
398 else
399 {
400 tmp = octave_sparse_params::get_key ("piv_tol");
401 if (! xisnan (tmp))
402 Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
403
404 tmp = octave_sparse_params::get_key ("sym_tol");
405 if (! xisnan (tmp))
406 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
407 }
408
409 // Set whether we are allowed to modify Q or not
410 tmp = octave_sparse_params::get_key ("autoamd");
411 if (! xisnan (tmp))
412 Control (UMFPACK_FIXQ) = tmp;
413
414 // Turn-off UMFPACK scaling for LU
415 if (scale)
416 Control (UMFPACK_SCALE) = UMFPACK_SCALE_SUM;
417 else
418 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
419
420 umfpack_report_control<lu_elt_type> (control);
421
422 const octave_idx_type *Ap = a.cidx ();
423 const octave_idx_type *Ai = a.ridx ();
424 const lu_elt_type *Ax = a.data ();
425
426 umfpack_report_matrix<lu_elt_type> (nr, nc, Ap, Ai, Ax, static_cast<octave_idx_type> (1), control);
427
428 void *Symbolic;
429 Matrix Info (1, UMFPACK_INFO);
430 double *info = Info.fortran_vec ();
431 int status = umfpack_qsymbolic<lu_elt_type> (nr, nc, Ap, Ai, Ax, 0, &Symbolic, control, info);
432
433 if (status < 0)
434 {
435 umfpack_report_status<lu_elt_type> (control, status);
436 umfpack_report_info<lu_elt_type> (control, info);
437
438 umfpack_free_symbolic<lu_elt_type> (&Symbolic);
439
440 (*current_liboctave_error_handler)
441 ("sparse_lu: symbolic factorization failed");
442 }
443 else
444 {
445 umfpack_report_symbolic<lu_elt_type> (Symbolic, control);
446
447 void *Numeric;
448 status = umfpack_numeric<lu_elt_type> (Ap, Ai, Ax, Symbolic, &Numeric, control, info);
449 umfpack_free_symbolic<lu_elt_type> (&Symbolic);
450
451 cond = Info (UMFPACK_RCOND);
452
453 if (status < 0)
454 {
455 umfpack_report_status<lu_elt_type> (control, status);
456 umfpack_report_info<lu_elt_type> (control, info);
457
458 umfpack_free_numeric<lu_elt_type> (&Numeric);
459
460 (*current_liboctave_error_handler)
461 ("sparse_lu: numeric factorization failed");
462 }
463 else
464 {
465 umfpack_report_numeric<lu_elt_type> (Numeric, control);
466
467 octave_idx_type lnz, unz;
468 status = umfpack_get_lunz<lu_elt_type> (&lnz, &unz, Numeric);
469
470 if (status < 0)
471 {
472 umfpack_report_status<lu_elt_type> (control, status);
473 umfpack_report_info<lu_elt_type> (control, info);
474
475 umfpack_free_numeric<lu_elt_type> (&Numeric);
476
477 (*current_liboctave_error_handler)
478 ("sparse_lu: extracting LU factors failed");
479 }
480 else
481 {
482 octave_idx_type n_inner = (nr < nc ? nr : nc);
483
484 if (lnz < 1)
485 Lfact = lu_type (n_inner, nr, static_cast<octave_idx_type> (1));
486 else
487 Lfact = lu_type (n_inner, nr, lnz);
488
489 octave_idx_type *Ltp = Lfact.cidx ();
490 octave_idx_type *Ltj = Lfact.ridx ();
491 lu_elt_type *Ltx = Lfact.data ();
492
493 if (unz < 1)
494 Ufact = lu_type (n_inner, nc, static_cast<octave_idx_type> (1));
495 else
496 Ufact = lu_type (n_inner, nc, unz);
497
498 octave_idx_type *Up = Ufact.cidx ();
499 octave_idx_type *Uj = Ufact.ridx ();
500 lu_elt_type *Ux = Ufact.data ();
501
502 Rfact = SparseMatrix (nr, nr, nr);
503 for (octave_idx_type i = 0; i < nr; i++)
504 {
505 Rfact.xridx (i) = i;
506 Rfact.xcidx (i) = i;
507 }
508 Rfact.xcidx (nr) = nr;
509 double *Rx = Rfact.data ();
510
511 P.resize (dim_vector (nr, 1));
512 octave_idx_type *p = P.fortran_vec ();
513
514 Q.resize (dim_vector (nc, 1));
515 octave_idx_type *q = Q.fortran_vec ();
516
517 octave_idx_type do_recip;
518 status = umfpack_get_numeric<lu_elt_type> (Ltp, Ltj, Ltx, Up, Uj, Ux, p, q, 0, &do_recip, Rx, Numeric);
519
520 umfpack_free_numeric<lu_elt_type> (&Numeric);
521
522 if (status < 0)
523 {
524 umfpack_report_status<lu_elt_type> (control, status);
525
526 (*current_liboctave_error_handler)
527 ("sparse_lu: extracting LU factors failed");
528 }
529 else
530 {
531 Lfact = Lfact.transpose ();
532
533 if (do_recip)
534 for (octave_idx_type i = 0; i < nr; i++)
535 Rx[i] = 1.0 / Rx[i];
536
537 umfpack_report_matrix<lu_elt_type> (nr, n_inner, Lfact.cidx (), Lfact.ridx (), Lfact.data (), static_cast<octave_idx_type> (1), control);
538 umfpack_report_matrix<lu_elt_type> (n_inner, nc, Ufact.cidx (), Ufact.ridx (), Ufact.data (), static_cast<octave_idx_type> (1), control);
539 umfpack_report_perm<lu_elt_type> (nr, p, control);
540 umfpack_report_perm<lu_elt_type> (nc, q, control);
541 }
542
543 umfpack_report_info<lu_elt_type> (control, info);
544 }
545 }
546 }
547
548 #else
549 (*current_liboctave_error_handler)
550 ("support for UMFPACK was unavailable or disabled when liboctave was built");
551 #endif
552 }
553
554 template <typename lu_type>
555 sparse_lu<lu_type>::sparse_lu (const lu_type& a,
556 const ColumnVector& Qinit,
557 const Matrix& piv_thres, bool scale,
558 bool FixedQ, double droptol,
559 bool milu, bool udiag)
560 : Lfact (), Ufact (), Rfact (), cond (0), P (), Q ()
561 {
562 #ifdef HAVE_UMFPACK
563 if (milu)
564 (*current_liboctave_error_handler)
565 ("Modified incomplete LU not implemented");
566
567 octave_idx_type nr = a.rows ();
568 octave_idx_type nc = a.cols ();
569
570 // Setup the control parameters
571 Matrix Control (UMFPACK_CONTROL, 1);
572 double *control = Control.fortran_vec ();
573 umfpack_defaults<lu_elt_type> (control);
574
575 double tmp = octave_sparse_params::get_key ("spumoni");
576 if (! xisnan (tmp))
577 Control (UMFPACK_PRL) = tmp;
578
579 if (piv_thres.numel () == 2)
580 {
581 tmp = (piv_thres (0) > 1. ? 1. : piv_thres (0));
582 if (! xisnan (tmp))
583 Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
584 tmp = (piv_thres (1) > 1. ? 1. : piv_thres (1));
585 if (! xisnan (tmp))
586 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
587 }
588 else
589 {
590 tmp = octave_sparse_params::get_key ("piv_tol");
591 if (! xisnan (tmp))
592 Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
593
594 tmp = octave_sparse_params::get_key ("sym_tol");
595 if (! xisnan (tmp))
596 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
597 }
598
599 if (droptol >= 0.)
600 Control (UMFPACK_DROPTOL) = droptol;
601
602 // Set whether we are allowed to modify Q or not
603 if (FixedQ)
604 Control (UMFPACK_FIXQ) = 1.0;
605 else
606 {
607 tmp = octave_sparse_params::get_key ("autoamd");
608 if (! xisnan (tmp))
609 Control (UMFPACK_FIXQ) = tmp;
610 }
611
612 // Turn-off UMFPACK scaling for LU
613 if (scale)
614 Control (UMFPACK_SCALE) = UMFPACK_SCALE_SUM;
615 else
616 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
617
618 umfpack_report_control<lu_elt_type> (control);
619
620 const octave_idx_type *Ap = a.cidx ();
621 const octave_idx_type *Ai = a.ridx ();
622 const lu_elt_type *Ax = a.data ();
623
624 umfpack_report_matrix<lu_elt_type> (nr, nc, Ap, Ai, Ax, static_cast<octave_idx_type> (1), control);
625
626 void *Symbolic;
627 Matrix Info (1, UMFPACK_INFO);
628 double *info = Info.fortran_vec ();
629 int status;
630
631 // Null loop so that qinit is imediately deallocated when not needed
632 do
633 {
634 OCTAVE_LOCAL_BUFFER (octave_idx_type, qinit, nc);
635
636 for (octave_idx_type i = 0; i < nc; i++)
637 qinit[i] = static_cast<octave_idx_type> (Qinit (i));
638
639 status = umfpack_qsymbolic<lu_elt_type> (nr, nc, Ap, Ai, Ax, qinit, &Symbolic, control, info);
640 }
641 while (0);
642
643 if (status < 0)
644 {
645 umfpack_report_status<lu_elt_type> (control, status);
646 umfpack_report_info<lu_elt_type> (control, info);
647
648 umfpack_free_symbolic<lu_elt_type> (&Symbolic);
649
650 (*current_liboctave_error_handler)
651 ("sparse_lu: symbolic factorization failed");
652 }
653 else
654 {
655 umfpack_report_symbolic<lu_elt_type> (Symbolic, control);
656
657 void *Numeric;
658 status = umfpack_numeric<lu_elt_type> (Ap, Ai, Ax, Symbolic, &Numeric, control, info);
659 umfpack_free_symbolic<lu_elt_type> (&Symbolic);
660
661 cond = Info (UMFPACK_RCOND);
662
663 if (status < 0)
664 {
665 umfpack_report_status<lu_elt_type> (control, status);
666 umfpack_report_info<lu_elt_type> (control, info);
667
668 umfpack_free_numeric<lu_elt_type> (&Numeric);
669
670 (*current_liboctave_error_handler)
671 ("sparse_lu: numeric factorization failed");
672 }
673 else
674 {
675 umfpack_report_numeric<lu_elt_type> (Numeric, control);
676
677 octave_idx_type lnz, unz;
678 status = umfpack_get_lunz<lu_elt_type> (&lnz, &unz, Numeric);
679
680 if (status < 0)
681 {
682 umfpack_report_status<lu_elt_type> (control, status);
683 umfpack_report_info<lu_elt_type> (control, info);
684
685 umfpack_free_numeric<lu_elt_type> (&Numeric);
686
687 (*current_liboctave_error_handler)
688 ("sparse_lu: extracting LU factors failed");
689 }
690 else
691 {
692 octave_idx_type n_inner = (nr < nc ? nr : nc);
693
694 if (lnz < 1)
695 Lfact = lu_type (n_inner, nr, static_cast<octave_idx_type> (1));
696 else
697 Lfact = lu_type (n_inner, nr, lnz);
698
699 octave_idx_type *Ltp = Lfact.cidx ();
700 octave_idx_type *Ltj = Lfact.ridx ();
701 lu_elt_type *Ltx = Lfact.data ();
702
703 if (unz < 1)
704 Ufact = lu_type (n_inner, nc, static_cast<octave_idx_type> (1));
705 else
706 Ufact = lu_type (n_inner, nc, unz);
707
708 octave_idx_type *Up = Ufact.cidx ();
709 octave_idx_type *Uj = Ufact.ridx ();
710 lu_elt_type *Ux = Ufact.data ();
711
712 Rfact = SparseMatrix (nr, nr, nr);
713 for (octave_idx_type i = 0; i < nr; i++)
714 {
715 Rfact.xridx (i) = i;
716 Rfact.xcidx (i) = i;
717 }
718 Rfact.xcidx (nr) = nr;
719 double *Rx = Rfact.data ();
720
721 P.resize (dim_vector (nr, 1));
722 octave_idx_type *p = P.fortran_vec ();
723
724 Q.resize (dim_vector (nc, 1));
725 octave_idx_type *q = Q.fortran_vec ();
726
727 octave_idx_type do_recip;
728 status = umfpack_get_numeric<lu_elt_type> (Ltp, Ltj, Ltx, Up, Uj, Ux, p, q, 0, &do_recip, Rx, Numeric);
729
730 umfpack_free_numeric<lu_elt_type> (&Numeric);
731
732 if (status < 0)
733 {
734 umfpack_report_status<lu_elt_type> (control, status);
735
736 (*current_liboctave_error_handler)
737 ("sparse_lu: extracting LU factors failed");
738 }
739 else
740 {
741 Lfact = Lfact.transpose ();
742
743 if (do_recip)
744 for (octave_idx_type i = 0; i < nr; i++)
745 Rx[i] = 1.0 / Rx[i];
746
747 umfpack_report_matrix<lu_elt_type> (nr, n_inner, Lfact.cidx (), Lfact.ridx (), Lfact.data (), static_cast<octave_idx_type> (1), control);
748 umfpack_report_matrix<lu_elt_type> (n_inner, nc, Ufact.cidx (), Ufact.ridx (), Ufact.data (), static_cast<octave_idx_type> (1), control);
749 umfpack_report_perm<lu_elt_type> (nr, p, control);
750 umfpack_report_perm<lu_elt_type> (nc, q, control);
751 }
752
753 umfpack_report_info<lu_elt_type> (control, info);
754 }
755 }
756 }
757
758 if (udiag)
759 (*current_liboctave_error_handler)
760 ("Option udiag of incomplete LU not implemented");
761
762 #else
763 (*current_liboctave_error_handler)
764 ("support for UMFPACK was unavailable or disabled when liboctave was built");
765 #endif
766 }
767
768 template <typename lu_type>
769 lu_type
770 sparse_lu<lu_type>::Y (void) const
771 {
772 octave_idx_type nr = Lfact.rows ();
773 octave_idx_type nz = Lfact.cols ();
774 octave_idx_type nc = Ufact.cols ();
775
776 lu_type Yout (nr, nc, Lfact.nnz () + Ufact.nnz () - (nr<nz?nr:nz));
777 octave_idx_type ii = 0;
778 Yout.xcidx (0) = 0;
779
780 for (octave_idx_type j = 0; j < nc; j++)
781 {
782 for (octave_idx_type i = Ufact.cidx (j); i < Ufact.cidx (j + 1); i++)
783 {
784 Yout.xridx (ii) = Ufact.ridx (i);
785 Yout.xdata (ii++) = Ufact.data (i);
786 }
787 if (j < nz)
788 {
789 // Note the +1 skips the 1.0 on the diagonal
790 for (octave_idx_type i = Lfact.cidx (j) + 1;
791 i < Lfact.cidx (j +1); i++)
792 {
793 Yout.xridx (ii) = Lfact.ridx (i);
794 Yout.xdata (ii++) = Lfact.data (i);
795 }
796 }
797 Yout.xcidx (j + 1) = ii;
798 }
799
800 return Yout;
801 }
802
803 template <typename lu_type>
804 SparseMatrix
805 sparse_lu<lu_type>::Pr (void) const
806 {
807
808 octave_idx_type nr = Lfact.rows ();
809
810 SparseMatrix Pout (nr, nr, nr);
811
812 for (octave_idx_type i = 0; i < nr; i++)
813 {
814 Pout.cidx (i) = i;
815 Pout.ridx (P (i)) = i;
816 Pout.data (i) = 1;
817 }
818 Pout.cidx (nr) = nr;
819
820 return Pout;
821 }
822
823 template <typename lu_type>
824 ColumnVector
825 sparse_lu<lu_type>::Pr_vec (void) const
826 {
827
828 octave_idx_type nr = Lfact.rows ();
829
830 ColumnVector Pout (nr);
831
832 for (octave_idx_type i = 0; i < nr; i++)
833 Pout.xelem (i) = static_cast<double> (P(i) + 1);
834
835 return Pout;
836 }
837
838 template <typename lu_type>
839 PermMatrix
840 sparse_lu<lu_type>::Pr_mat (void) const
841 {
842 return PermMatrix (P, false);
843 }
844
845 template <typename lu_type>
846 SparseMatrix
847 sparse_lu<lu_type>::Pc (void) const
848 {
849 octave_idx_type nc = Ufact.cols ();
850
851 SparseMatrix Pout (nc, nc, nc);
852
853 for (octave_idx_type i = 0; i < nc; i++)
854 {
855 Pout.cidx (i) = i;
856 Pout.ridx (i) = Q (i);
857 Pout.data (i) = 1;
858 }
859 Pout.cidx (nc) = nc;
860
861 return Pout;
862 }
863
864 template <typename lu_type>
865 ColumnVector
866 sparse_lu<lu_type>::Pc_vec (void) const
867 {
868
869 octave_idx_type nc = Ufact.cols ();
870
871 ColumnVector Pout (nc);
872
873 for (octave_idx_type i = 0; i < nc; i++)
874 Pout.xelem (i) = static_cast<double> (Q(i) + 1);
875
876 return Pout;
877 }
878
879 template <typename lu_type>
880 PermMatrix
881 sparse_lu<lu_type>::Pc_mat (void) const
882 {
883 return PermMatrix (Q, true);
884 }