GNU Octave 7.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
sparse-lu.cc
Go to the documentation of this file.
1////////////////////////////////////////////////////////////////////////
2//
3// Copyright (C) 1998-2022 The Octave Project Developers
4//
5// See the file COPYRIGHT.md in the top-level directory of this
6// distribution or <https://octave.org/copyright/>.
7//
8// This file is part of Octave.
9//
10// Octave is free software: you can redistribute it and/or modify it
11// under the terms of the GNU General Public License as published by
12// the Free Software Foundation, either version 3 of the License, or
13// (at your option) any later version.
14//
15// Octave is distributed in the hope that it will be useful, but
16// WITHOUT ANY WARRANTY; without even the implied warranty of
17// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18// GNU General Public License for more details.
19//
20// You should have received a copy of the GNU General Public License
21// along with Octave; see the file COPYING. If not, see
22// <https://www.gnu.org/licenses/>.
23//
24////////////////////////////////////////////////////////////////////////
25
26#if defined (HAVE_CONFIG_H)
27# include "config.h"
28#endif
29
30#include "CSparse.h"
31#include "PermMatrix.h"
32#include "dSparse.h"
33#include "lo-error.h"
34#include "lo-mappers.h"
35#include "oct-locbuf.h"
36#include "oct-sparse.h"
37#include "oct-spparms.h"
38#include "sparse-lu.h"
39
40namespace octave
41{
42 namespace math
43 {
44 // Wrappers for SuiteSparse (formerly UMFPACK) functions that have
45 // different names depending on the sparse matrix data type.
46 //
47 // All of these functions must be specialized to forward to the correct
48 // SuiteSparse functions.
49
50 template <typename T>
51 void
52 umfpack_defaults (double *Control);
53
54 template <typename T>
55 void
56 umfpack_free_numeric (void **Numeric);
57
58 template <typename T>
59 void
60 umfpack_free_symbolic (void **Symbolic);
61
62 template <typename T>
63 octave_idx_type
64 umfpack_get_lunz (octave_idx_type *lnz, octave_idx_type *unz,
65 void *Numeric);
66
67 template <typename T>
68 octave_idx_type
69 umfpack_get_numeric (octave_idx_type *Lp, octave_idx_type *Lj,
70 T *Lx, // Or Lz_packed
71 octave_idx_type *Up, octave_idx_type *Ui,
72 T *Ux, // Or Uz_packed
73 octave_idx_type *p, octave_idx_type *q,
74 double *Dz_packed, octave_idx_type *do_recip,
75 double *Rs, void *Numeric);
76
77 template <typename T>
78 octave_idx_type
79 umfpack_numeric (const octave_idx_type *Ap, const octave_idx_type *Ai,
80 const T *Ax, // Or Az_packed
81 void *Symbolic, void **Numeric,
82 const double *Control, double *Info);
83
84 template <typename T>
85 octave_idx_type
86 umfpack_qsymbolic (octave_idx_type n_row, octave_idx_type n_col,
87 const octave_idx_type *Ap, const octave_idx_type *Ai,
88 const T *Ax, // Or Az_packed
89 const octave_idx_type *Qinit, void **Symbolic,
90 const double *Control, double *Info);
91
92 template <typename T>
93 void
94 umfpack_report_control (const double *Control);
95
96 template <typename T>
97 void
98 umfpack_report_info (const double *Control, const double *Info);
99
100 template <typename T>
101 void
102 umfpack_report_matrix (octave_idx_type n_row, octave_idx_type n_col,
103 const octave_idx_type *Ap,
104 const octave_idx_type *Ai,
105 const T *Ax, // Or Az_packed
106 octave_idx_type col_form, const double *Control);
107
108 template <typename T>
109 void
110 umfpack_report_numeric (void *Numeric, const double *Control);
111
112 template <typename T>
113 void
114 umfpack_report_perm (octave_idx_type np, const octave_idx_type *Perm,
115 const double *Control);
116
117 template <typename T>
118 void
119 umfpack_report_status (double *Control, octave_idx_type status);
120
121 template <typename T>
122 void
123 umfpack_report_symbolic (void *Symbolic, const double *Control);
124
125#if defined (HAVE_UMFPACK)
126
127 // SparseMatrix Specialization.
128
129 template <>
130 inline OCTAVE_API void
131 umfpack_defaults<double> (double *Control)
132 {
133 UMFPACK_DNAME (defaults) (Control);
134 }
135
136 template <>
137 inline OCTAVE_API void
138 umfpack_free_numeric<double> (void **Numeric)
139 {
140 UMFPACK_DNAME (free_numeric) (Numeric);
141 }
142
143 template <>
144 inline OCTAVE_API void
145 umfpack_free_symbolic<double> (void **Symbolic)
146 {
147 UMFPACK_DNAME (free_symbolic) (Symbolic);
148 }
149
150 template <>
151 inline OCTAVE_API octave_idx_type
152 umfpack_get_lunz<double>
153 (octave_idx_type *lnz, octave_idx_type *unz, void *Numeric)
154 {
155 suitesparse_integer ignore1, ignore2, ignore3;
156
157 return UMFPACK_DNAME (get_lunz) (to_suitesparse_intptr (lnz),
158 to_suitesparse_intptr (unz),
159 &ignore1, &ignore2, &ignore3, Numeric);
160 }
161
162 template <>
163 inline OCTAVE_API octave_idx_type
164 umfpack_get_numeric<double>
165 (octave_idx_type *Lp, octave_idx_type *Lj, double *Lx,
166 octave_idx_type *Up, octave_idx_type *Ui, double *Ux,
167 octave_idx_type *p, octave_idx_type *q, double *Dx,
168 octave_idx_type *do_recip, double *Rs, void *Numeric)
169 {
170 return UMFPACK_DNAME (get_numeric) (to_suitesparse_intptr (Lp),
171 to_suitesparse_intptr (Lj),
172 Lx, to_suitesparse_intptr (Up),
173 to_suitesparse_intptr (Ui), Ux,
174 to_suitesparse_intptr (p),
175 to_suitesparse_intptr (q), Dx,
176 to_suitesparse_intptr (do_recip),
177 Rs, Numeric);
178 }
179
180 template <>
181 inline OCTAVE_API octave_idx_type
182 umfpack_numeric<double>
183 (const octave_idx_type *Ap, const octave_idx_type *Ai,
184 const double *Ax, void *Symbolic, void **Numeric,
185 const double *Control, double *Info)
186 {
187 return UMFPACK_DNAME (numeric) (to_suitesparse_intptr (Ap),
188 to_suitesparse_intptr (Ai),
189 Ax, Symbolic, Numeric, Control, Info);
190 }
191
192 template <>
193 inline OCTAVE_API octave_idx_type
194 umfpack_qsymbolic<double>
195 (octave_idx_type n_row, octave_idx_type n_col, const octave_idx_type *Ap,
196 const octave_idx_type *Ai, const double *Ax,
197 const octave_idx_type *Qinit, void **Symbolic,
198 const double *Control, double *Info)
199 {
200 return UMFPACK_DNAME (qsymbolic) (n_row, n_col,
201 to_suitesparse_intptr (Ap),
202 to_suitesparse_intptr (Ai), Ax,
203 to_suitesparse_intptr (Qinit),
204 Symbolic, Control, Info);
205 }
206
207 template <>
208 inline OCTAVE_API void
209 umfpack_report_control<double> (const double *Control)
210 {
211 UMFPACK_DNAME (report_control) (Control);
212 }
213
214 template <>
215 inline OCTAVE_API void
216 umfpack_report_info<double> (const double *Control, const double *Info)
217 {
218 UMFPACK_DNAME (report_info) (Control, Info);
219 }
220
221 template <>
222 inline OCTAVE_API void
223 umfpack_report_matrix<double>
224 (octave_idx_type n_row, octave_idx_type n_col, const octave_idx_type *Ap,
225 const octave_idx_type *Ai, const double *Ax, octave_idx_type col_form,
226 const double *Control)
227 {
228 UMFPACK_DNAME (report_matrix) (n_row, n_col,
229 to_suitesparse_intptr (Ap),
230 to_suitesparse_intptr (Ai), Ax,
231 col_form, Control);
232 }
233
234 template <>
235 inline OCTAVE_API void
236 umfpack_report_numeric<double> (void *Numeric, const double *Control)
237 {
238 UMFPACK_DNAME (report_numeric) (Numeric, Control);
239 }
240
241 template <>
242 inline OCTAVE_API void
243 umfpack_report_perm<double>
244 (octave_idx_type np, const octave_idx_type *Perm, const double *Control)
245 {
246 UMFPACK_DNAME (report_perm) (np, to_suitesparse_intptr (Perm), Control);
247 }
248
249 template <>
250 inline OCTAVE_API void
251 umfpack_report_status<double> (double *Control, octave_idx_type status)
252 {
253 UMFPACK_DNAME (report_status) (Control, status);
254 }
255
256 template <>
257 inline OCTAVE_API void
258 umfpack_report_symbolic<double> (void *Symbolic, const double *Control)
259 {
260 UMFPACK_DNAME (report_symbolic) (Symbolic, Control);
261 }
262
263 // SparseComplexMatrix specialization.
264
265 template <>
266 inline OCTAVE_API void
267 umfpack_defaults<Complex> (double *Control)
268 {
269 UMFPACK_ZNAME (defaults) (Control);
270 }
271
272 template <>
273 inline OCTAVE_API void
274 umfpack_free_numeric<Complex> (void **Numeric)
275 {
276 UMFPACK_ZNAME (free_numeric) (Numeric);
277 }
278
279 template <>
280 inline OCTAVE_API void
281 umfpack_free_symbolic<Complex> (void **Symbolic)
282 {
283 UMFPACK_ZNAME (free_symbolic) (Symbolic);
284 }
285
286 template <>
287 inline OCTAVE_API octave_idx_type
288 umfpack_get_lunz<Complex>
289 (octave_idx_type *lnz, octave_idx_type *unz, void *Numeric)
290 {
291 suitesparse_integer ignore1, ignore2, ignore3;
292
293 return UMFPACK_ZNAME (get_lunz) (to_suitesparse_intptr (lnz),
294 to_suitesparse_intptr (unz),
295 &ignore1, &ignore2, &ignore3, Numeric);
296 }
297
298 template <>
299 inline OCTAVE_API octave_idx_type
300 umfpack_get_numeric<Complex>
301 (octave_idx_type *Lp, octave_idx_type *Lj, Complex *Lz,
302 octave_idx_type *Up, octave_idx_type *Ui, Complex *Uz,
303 octave_idx_type *p, octave_idx_type *q, double *Dz,
304 octave_idx_type *do_recip, double *Rs, void *Numeric)
305 {
306 return UMFPACK_ZNAME (get_numeric) (to_suitesparse_intptr (Lp),
307 to_suitesparse_intptr (Lj),
308 reinterpret_cast<double *> (Lz),
309 nullptr, to_suitesparse_intptr (Up),
310 to_suitesparse_intptr (Ui),
311 reinterpret_cast<double *> (Uz),
312 nullptr, to_suitesparse_intptr (p),
313 to_suitesparse_intptr (q),
314 reinterpret_cast<double *> (Dz),
315 nullptr, to_suitesparse_intptr (do_recip),
316 Rs, Numeric);
317 }
318
319 template <>
320 inline OCTAVE_API octave_idx_type
321 umfpack_numeric<Complex>
322 (const octave_idx_type *Ap, const octave_idx_type *Ai,
323 const Complex *Az, void *Symbolic, void **Numeric,
324 const double *Control, double *Info)
325 {
326 return UMFPACK_ZNAME (numeric) (to_suitesparse_intptr (Ap),
327 to_suitesparse_intptr (Ai),
328 reinterpret_cast<const double *> (Az),
329 nullptr, Symbolic, Numeric, Control, Info);
330 }
331
332 template <>
333 inline OCTAVE_API octave_idx_type
334 umfpack_qsymbolic<Complex>
335 (octave_idx_type n_row, octave_idx_type n_col,
336 const octave_idx_type *Ap, const octave_idx_type *Ai,
337 const Complex *Az, const octave_idx_type *Qinit,
338 void **Symbolic, const double *Control, double *Info)
339 {
340 return UMFPACK_ZNAME (qsymbolic) (n_row, n_col,
341 to_suitesparse_intptr (Ap),
342 to_suitesparse_intptr (Ai),
343 reinterpret_cast<const double *> (Az),
344 nullptr, to_suitesparse_intptr (Qinit),
345 Symbolic, Control, Info);
346 }
347
348 template <>
349 inline OCTAVE_API void
350 umfpack_report_control<Complex> (const double *Control)
351 {
352 UMFPACK_ZNAME (report_control) (Control);
353 }
354
355 template <>
356 inline OCTAVE_API void
357 umfpack_report_info<Complex> (const double *Control, const double *Info)
358 {
359 UMFPACK_ZNAME (report_info) (Control, Info);
360 }
361
362 template <>
363 inline OCTAVE_API void
364 umfpack_report_matrix<Complex>
365 (octave_idx_type n_row, octave_idx_type n_col,
366 const octave_idx_type *Ap, const octave_idx_type *Ai,
367 const Complex *Az, octave_idx_type col_form, const double *Control)
368 {
369 UMFPACK_ZNAME (report_matrix) (n_row, n_col,
370 to_suitesparse_intptr (Ap),
371 to_suitesparse_intptr (Ai),
372 reinterpret_cast<const double *> (Az),
373 nullptr, col_form, Control);
374 }
375
376 template <>
377 inline OCTAVE_API void
378 umfpack_report_numeric<Complex> (void *Numeric, const double *Control)
379 {
380 UMFPACK_ZNAME (report_numeric) (Numeric, Control);
381 }
382
383 template <>
384 inline OCTAVE_API void
385 umfpack_report_perm<Complex>
386 (octave_idx_type np, const octave_idx_type *Perm, const double *Control)
387 {
388 UMFPACK_ZNAME (report_perm) (np, to_suitesparse_intptr (Perm), Control);
389 }
390
391 template <>
392 inline OCTAVE_API void
393 umfpack_report_status<Complex> (double *Control, octave_idx_type status)
394 {
395 UMFPACK_ZNAME (report_status) (Control, status);
396 }
397
398 template <>
399 inline OCTAVE_API void
400 umfpack_report_symbolic <Complex> (void *Symbolic, const double *Control)
401 {
402 UMFPACK_ZNAME (report_symbolic) (Symbolic, Control);
403 }
404
405#endif
406
407 template <typename lu_type>
408 sparse_lu<lu_type>::sparse_lu (const lu_type& a, const Matrix& piv_thres,
409 bool scale)
410 : m_L (), m_U (), m_R (), m_cond (0), m_P (), m_Q ()
411 {
412#if defined (HAVE_UMFPACK)
413 octave_idx_type nr = a.rows ();
414 octave_idx_type nc = a.cols ();
415
416 // Setup the control parameters
417 Matrix Control (UMFPACK_CONTROL, 1);
418 double *control = Control.fortran_vec ();
419 umfpack_defaults<lu_elt_type> (control);
420
421 double tmp = sparse_params::get_key ("spumoni");
422 if (! math::isnan (tmp))
423 Control (UMFPACK_PRL) = tmp;
424
425 if (piv_thres.numel () == 2)
426 {
427 tmp = (piv_thres (0) > 1. ? 1. : piv_thres (0));
428 if (! math::isnan (tmp))
429 Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
430
431 tmp = (piv_thres (1) > 1. ? 1. : piv_thres (1));
432 if (! math::isnan (tmp))
433 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
434 }
435 else
436 {
437 tmp = sparse_params::get_key ("piv_tol");
438 if (! math::isnan (tmp))
439 Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
440
441 tmp = sparse_params::get_key ("sym_tol");
442 if (! math::isnan (tmp))
443 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
444 }
445
446 // Set whether we are allowed to modify Q or not
447 tmp = sparse_params::get_key ("autoamd");
448 if (! math::isnan (tmp))
449 Control (UMFPACK_FIXQ) = tmp;
450
451 // Turn-off UMFPACK scaling for LU
452 if (scale)
453 Control (UMFPACK_SCALE) = UMFPACK_SCALE_SUM;
454 else
455 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
456
457 umfpack_report_control<lu_elt_type> (control);
458
459 const octave_idx_type *Ap = a.cidx ();
460 const octave_idx_type *Ai = a.ridx ();
461 const lu_elt_type *Ax = a.data ();
462
463 umfpack_report_matrix<lu_elt_type> (nr, nc, Ap, Ai, Ax,
464 static_cast<octave_idx_type> (1),
465 control);
466
467 void *Symbolic;
468 Matrix Info (1, UMFPACK_INFO);
469 double *info = Info.fortran_vec ();
470 int status = umfpack_qsymbolic<lu_elt_type> (nr, nc, Ap, Ai, Ax, nullptr,
471 &Symbolic, control, info);
472
473 if (status < 0)
474 {
475 umfpack_report_status<lu_elt_type> (control, status);
476 umfpack_report_info<lu_elt_type> (control, info);
477
478 umfpack_free_symbolic<lu_elt_type> (&Symbolic);
479
480 (*current_liboctave_error_handler)
481 ("sparse_lu: symbolic factorization failed");
482 }
483 else
484 {
485 umfpack_report_symbolic<lu_elt_type> (Symbolic, control);
486
487 void *Numeric;
488 status = umfpack_numeric<lu_elt_type> (Ap, Ai, Ax, Symbolic,
489 &Numeric, control, info);
490 umfpack_free_symbolic<lu_elt_type> (&Symbolic);
491
492 m_cond = Info (UMFPACK_RCOND);
493
494 if (status < 0)
495 {
496 umfpack_report_status<lu_elt_type> (control, status);
497 umfpack_report_info<lu_elt_type> (control, info);
498
499 umfpack_free_numeric<lu_elt_type> (&Numeric);
500
501 (*current_liboctave_error_handler)
502 ("sparse_lu: numeric factorization failed");
503 }
504 else
505 {
506 umfpack_report_numeric<lu_elt_type> (Numeric, control);
507
508 octave_idx_type lnz, unz;
509 status = umfpack_get_lunz<lu_elt_type> (&lnz, &unz, Numeric);
510
511 if (status < 0)
512 {
513 umfpack_report_status<lu_elt_type> (control, status);
514 umfpack_report_info<lu_elt_type> (control, info);
515
516 umfpack_free_numeric<lu_elt_type> (&Numeric);
517
518 (*current_liboctave_error_handler)
519 ("sparse_lu: extracting LU factors failed");
520 }
521 else
522 {
523 octave_idx_type n_inner = (nr < nc ? nr : nc);
524
525 if (lnz < 1)
526 m_L = lu_type (n_inner, nr,
527 static_cast<octave_idx_type> (1));
528 else
529 m_L = lu_type (n_inner, nr, lnz);
530
531 octave_idx_type *Ltp = m_L.cidx ();
532 octave_idx_type *Ltj = m_L.ridx ();
533 lu_elt_type *Ltx = m_L.data ();
534
535 if (unz < 1)
536 m_U = lu_type (n_inner, nc,
537 static_cast<octave_idx_type> (1));
538 else
539 m_U = lu_type (n_inner, nc, unz);
540
541 octave_idx_type *Up = m_U.cidx ();
542 octave_idx_type *Uj = m_U.ridx ();
543 lu_elt_type *Ux = m_U.data ();
544
545 m_R = SparseMatrix (nr, nr, nr);
546 for (octave_idx_type i = 0; i < nr; i++)
547 {
548 m_R.xridx (i) = i;
549 m_R.xcidx (i) = i;
550 }
551 m_R.xcidx (nr) = nr;
552 double *Rx = m_R.data ();
553
554 m_P.resize (dim_vector (nr, 1));
555 octave_idx_type *p = m_P.fortran_vec ();
556
557 m_Q.resize (dim_vector (nc, 1));
558 octave_idx_type *q = m_Q.fortran_vec ();
559
560 octave_idx_type do_recip;
561 status = umfpack_get_numeric<lu_elt_type> (Ltp, Ltj, Ltx,
562 Up, Uj, Ux,
563 p, q, nullptr,
564 &do_recip, Rx,
565 Numeric);
566
567 umfpack_free_numeric<lu_elt_type> (&Numeric);
568
569 if (status < 0)
570 {
571 umfpack_report_status<lu_elt_type> (control, status);
572
573 (*current_liboctave_error_handler)
574 ("sparse_lu: extracting LU factors failed");
575 }
576 else
577 {
578 m_L = m_L.transpose ();
579
580 if (do_recip)
581 for (octave_idx_type i = 0; i < nr; i++)
582 Rx[i] = 1.0 / Rx[i];
583
584 umfpack_report_matrix<lu_elt_type> (nr, n_inner,
585 m_L.cidx (),
586 m_L.ridx (),
587 m_L.data (),
588 static_cast<octave_idx_type> (1),
589 control);
590 umfpack_report_matrix<lu_elt_type> (n_inner, nc,
591 m_U.cidx (),
592 m_U.ridx (),
593 m_U.data (),
594 static_cast<octave_idx_type> (1),
595 control);
596 umfpack_report_perm<lu_elt_type> (nr, p, control);
597 umfpack_report_perm<lu_elt_type> (nc, q, control);
598 }
599
600 umfpack_report_info<lu_elt_type> (control, info);
601 }
602 }
603 }
604
605#else
606
607 octave_unused_parameter (a);
608 octave_unused_parameter (piv_thres);
609 octave_unused_parameter (scale);
610
611 (*current_liboctave_error_handler)
612 ("support for UMFPACK was unavailable or disabled when liboctave was built");
613
614#endif
615 }
616
617 template <typename lu_type>
618 sparse_lu<lu_type>::sparse_lu (const lu_type& a,
619 const ColumnVector& Qinit,
620 const Matrix& piv_thres, bool scale,
621 bool FixedQ, double droptol,
622 bool milu, bool udiag)
623 : m_L (), m_U (), m_R (), m_cond (0), m_P (), m_Q ()
624 {
625#if defined (HAVE_UMFPACK)
626
627 if (milu)
628 (*current_liboctave_error_handler)
629 ("Modified incomplete LU not implemented");
630
631 octave_idx_type nr = a.rows ();
632 octave_idx_type nc = a.cols ();
633
634 // Setup the control parameters
635 Matrix Control (UMFPACK_CONTROL, 1);
636 double *control = Control.fortran_vec ();
637 umfpack_defaults<lu_elt_type> (control);
638
639 double tmp = sparse_params::get_key ("spumoni");
640 if (! math::isnan (tmp))
641 Control (UMFPACK_PRL) = tmp;
642
643 if (piv_thres.numel () == 2)
644 {
645 tmp = (piv_thres (0) > 1. ? 1. : piv_thres (0));
646 if (! math::isnan (tmp))
647 Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
648 tmp = (piv_thres (1) > 1. ? 1. : piv_thres (1));
649 if (! math::isnan (tmp))
650 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
651 }
652 else
653 {
654 tmp = sparse_params::get_key ("piv_tol");
655 if (! math::isnan (tmp))
656 Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
657
658 tmp = sparse_params::get_key ("sym_tol");
659 if (! math::isnan (tmp))
660 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
661 }
662
663 if (droptol >= 0.)
664 Control (UMFPACK_DROPTOL) = droptol;
665
666 // Set whether we are allowed to modify Q or not
667 if (FixedQ)
668 Control (UMFPACK_FIXQ) = 1.0;
669 else
670 {
671 tmp = sparse_params::get_key ("autoamd");
672 if (! math::isnan (tmp))
673 Control (UMFPACK_FIXQ) = tmp;
674 }
675
676 // Turn-off UMFPACK scaling for LU
677 if (scale)
678 Control (UMFPACK_SCALE) = UMFPACK_SCALE_SUM;
679 else
680 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
681
682 umfpack_report_control<lu_elt_type> (control);
683
684 const octave_idx_type *Ap = a.cidx ();
685 const octave_idx_type *Ai = a.ridx ();
686 const lu_elt_type *Ax = a.data ();
687
688 umfpack_report_matrix<lu_elt_type> (nr, nc, Ap, Ai, Ax,
689 static_cast<octave_idx_type> (1),
690 control);
691
692 void *Symbolic;
693 Matrix Info (1, UMFPACK_INFO);
694 double *info = Info.fortran_vec ();
695 int status;
696
697 // Null loop so that qinit is immediately deallocated when not needed
698 do
699 {
700 OCTAVE_LOCAL_BUFFER (octave_idx_type, qinit, nc);
701
702 for (octave_idx_type i = 0; i < nc; i++)
703 qinit[i] = static_cast<octave_idx_type> (Qinit (i));
704
705 status = umfpack_qsymbolic<lu_elt_type> (nr, nc, Ap, Ai, Ax,
706 qinit, &Symbolic, control,
707 info);
708 }
709 while (0);
710
711 if (status < 0)
712 {
713 umfpack_report_status<lu_elt_type> (control, status);
714 umfpack_report_info<lu_elt_type> (control, info);
715
716 umfpack_free_symbolic<lu_elt_type> (&Symbolic);
717
718 (*current_liboctave_error_handler)
719 ("sparse_lu: symbolic factorization failed");
720 }
721 else
722 {
723 umfpack_report_symbolic<lu_elt_type> (Symbolic, control);
724
725 void *Numeric;
726 status = umfpack_numeric<lu_elt_type> (Ap, Ai, Ax, Symbolic,
727 &Numeric, control, info);
728 umfpack_free_symbolic<lu_elt_type> (&Symbolic);
729
730 m_cond = Info (UMFPACK_RCOND);
731
732 if (status < 0)
733 {
734 umfpack_report_status<lu_elt_type> (control, status);
735 umfpack_report_info<lu_elt_type> (control, info);
736
737 umfpack_free_numeric<lu_elt_type> (&Numeric);
738
739 (*current_liboctave_error_handler)
740 ("sparse_lu: numeric factorization failed");
741 }
742 else
743 {
744 umfpack_report_numeric<lu_elt_type> (Numeric, control);
745
746 octave_idx_type lnz, unz;
747 status = umfpack_get_lunz<lu_elt_type> (&lnz, &unz, Numeric);
748
749 if (status < 0)
750 {
751 umfpack_report_status<lu_elt_type> (control, status);
752 umfpack_report_info<lu_elt_type> (control, info);
753
754 umfpack_free_numeric<lu_elt_type> (&Numeric);
755
756 (*current_liboctave_error_handler)
757 ("sparse_lu: extracting LU factors failed");
758 }
759 else
760 {
761 octave_idx_type n_inner = (nr < nc ? nr : nc);
762
763 if (lnz < 1)
764 m_L = lu_type (n_inner, nr,
765 static_cast<octave_idx_type> (1));
766 else
767 m_L = lu_type (n_inner, nr, lnz);
768
769 octave_idx_type *Ltp = m_L.cidx ();
770 octave_idx_type *Ltj = m_L.ridx ();
771 lu_elt_type *Ltx = m_L.data ();
772
773 if (unz < 1)
774 m_U = lu_type (n_inner, nc,
775 static_cast<octave_idx_type> (1));
776 else
777 m_U = lu_type (n_inner, nc, unz);
778
779 octave_idx_type *Up = m_U.cidx ();
780 octave_idx_type *Uj = m_U.ridx ();
781 lu_elt_type *Ux = m_U.data ();
782
783 m_R = SparseMatrix (nr, nr, nr);
784 for (octave_idx_type i = 0; i < nr; i++)
785 {
786 m_R.xridx (i) = i;
787 m_R.xcidx (i) = i;
788 }
789 m_R.xcidx (nr) = nr;
790 double *Rx = m_R.data ();
791
792 m_P.resize (dim_vector (nr, 1));
793 octave_idx_type *p = m_P.fortran_vec ();
794
795 m_Q.resize (dim_vector (nc, 1));
796 octave_idx_type *q = m_Q.fortran_vec ();
797
798 octave_idx_type do_recip;
799 status = umfpack_get_numeric<lu_elt_type> (Ltp, Ltj, Ltx,
800 Up, Uj, Ux,
801 p, q, nullptr,
802 &do_recip,
803 Rx, Numeric);
804
805 umfpack_free_numeric<lu_elt_type> (&Numeric);
806
807 if (status < 0)
808 {
809 umfpack_report_status<lu_elt_type> (control, status);
810
811 (*current_liboctave_error_handler)
812 ("sparse_lu: extracting LU factors failed");
813 }
814 else
815 {
816 m_L = m_L.transpose ();
817
818 if (do_recip)
819 for (octave_idx_type i = 0; i < nr; i++)
820 Rx[i] = 1.0 / Rx[i];
821
822 umfpack_report_matrix<lu_elt_type> (nr, n_inner,
823 m_L.cidx (),
824 m_L.ridx (),
825 m_L.data (),
826 static_cast<octave_idx_type> (1),
827 control);
828 umfpack_report_matrix<lu_elt_type> (n_inner, nc,
829 m_U.cidx (),
830 m_U.ridx (),
831 m_U.data (),
832 static_cast<octave_idx_type> (1),
833 control);
834 umfpack_report_perm<lu_elt_type> (nr, p, control);
835 umfpack_report_perm<lu_elt_type> (nc, q, control);
836 }
837
838 umfpack_report_info<lu_elt_type> (control, info);
839 }
840 }
841 }
842
843 if (udiag)
844 (*current_liboctave_error_handler)
845 ("Option udiag of incomplete LU not implemented");
846
847#else
848
849 octave_unused_parameter (a);
850 octave_unused_parameter (Qinit);
851 octave_unused_parameter (piv_thres);
852 octave_unused_parameter (scale);
853 octave_unused_parameter (FixedQ);
854 octave_unused_parameter (droptol);
855 octave_unused_parameter (milu);
856 octave_unused_parameter (udiag);
857
858 (*current_liboctave_error_handler)
859 ("support for UMFPACK was unavailable or disabled when liboctave was built");
860
861#endif
862 }
863
864 template <typename lu_type>
865 lu_type
866 sparse_lu<lu_type>::Y (void) const
867 {
868 octave_idx_type nr = m_L.rows ();
869 octave_idx_type nz = m_L.cols ();
870 octave_idx_type nc = m_U.cols ();
871
872 lu_type Yout (nr, nc, m_L.nnz () + m_U.nnz () - (nr<nz?nr:nz));
873 octave_idx_type ii = 0;
874 Yout.xcidx (0) = 0;
875
876 for (octave_idx_type j = 0; j < nc; j++)
877 {
878 for (octave_idx_type i = m_U.cidx (j); i < m_U.cidx (j + 1); i++)
879 {
880 Yout.xridx (ii) = m_U.ridx (i);
881 Yout.xdata (ii++) = m_U.data (i);
882 }
883
884 if (j < nz)
885 {
886 // Note the +1 skips the 1.0 on the diagonal
887 for (octave_idx_type i = m_L.cidx (j) + 1;
888 i < m_L.cidx (j +1); i++)
889 {
890 Yout.xridx (ii) = m_L.ridx (i);
891 Yout.xdata (ii++) = m_L.data (i);
892 }
893 }
894
895 Yout.xcidx (j + 1) = ii;
896 }
897
898 return Yout;
899 }
900
901 template <typename lu_type>
902 SparseMatrix
903 sparse_lu<lu_type>::Pr (void) const
904 {
905 octave_idx_type nr = m_L.rows ();
906
907 SparseMatrix Pout (nr, nr, nr);
908
909 for (octave_idx_type i = 0; i < nr; i++)
910 {
911 Pout.cidx (i) = i;
912 Pout.ridx (m_P(i)) = i;
913 Pout.data (i) = 1;
914 }
915
916 Pout.cidx (nr) = nr;
917
918 return Pout;
919 }
920
921 template <typename lu_type>
922 ColumnVector
923 sparse_lu<lu_type>::Pr_vec (void) const
924 {
925 octave_idx_type nr = m_L.rows ();
926
927 ColumnVector Pout (nr);
928
929 for (octave_idx_type i = 0; i < nr; i++)
930 Pout.xelem (i) = static_cast<double> (m_P(i) + 1);
931
932 return Pout;
933 }
934
935 template <typename lu_type>
936 PermMatrix
937 sparse_lu<lu_type>::Pr_mat (void) const
938 {
939 return PermMatrix (m_P, false);
940 }
941
942 template <typename lu_type>
943 SparseMatrix
944 sparse_lu<lu_type>::Pc (void) const
945 {
946 octave_idx_type nc = m_U.cols ();
947
948 SparseMatrix Pout (nc, nc, nc);
949
950 for (octave_idx_type i = 0; i < nc; i++)
951 {
952 Pout.cidx (i) = i;
953 Pout.ridx (i) = m_Q(i);
954 Pout.data (i) = 1;
955 }
956
957 Pout.cidx (nc) = nc;
958
959 return Pout;
960 }
961
962 template <typename lu_type>
963 ColumnVector
964 sparse_lu<lu_type>::Pc_vec (void) const
965 {
966 octave_idx_type nc = m_U.cols ();
967
968 ColumnVector Pout (nc);
969
970 for (octave_idx_type i = 0; i < nc; i++)
971 Pout.xelem (i) = static_cast<double> (m_Q(i) + 1);
972
973 return Pout;
974 }
975
976 template <typename lu_type>
977 PermMatrix
978 sparse_lu<lu_type>::Pc_mat (void) const
979 {
980 return PermMatrix (m_Q, true);
981 }
982
983 // Instantiations we need.
984 template class OCTAVE_API sparse_lu<SparseMatrix>;
985
986 template class OCTAVE_API sparse_lu<SparseComplexMatrix>;
987 }
988}
Array::xelem
T & xelem(octave_idx_type n)
Size of the specified dimension.
Definition: Array.h:504
Array::numel
octave_idx_type numel(void) const
Number of elements in the array.
Definition: Array.h:411
Array::resize
OCTARRAY_API void resize(const dim_vector &dv, const T &rfv)
Size of the specified dimension.
Definition: Array.cc:1010
Array::fortran_vec
OCTARRAY_API T * fortran_vec(void)
Size of the specified dimension.
Definition: Array.cc:1744
Matrix
Definition: dMatrix.h:42
Sparse::data
T * data(void)
Definition: Sparse.h:574
Sparse::xridx
octave_idx_type * xridx(void)
Definition: Sparse.h:589
Sparse::ridx
octave_idx_type * ridx(void)
Definition: Sparse.h:583
Sparse::cidx
octave_idx_type * cidx(void)
Definition: Sparse.h:596
Sparse::xcidx
octave_idx_type * xcidx(void)
Definition: Sparse.h:602
dim_vector
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:94
octave::math::sparse_lu
Definition: sparse-lu.h:50
octave::math::sparse_lu::Pc_mat
OCTAVE_API PermMatrix Pc_mat(void) const
Definition: sparse-lu.cc:978
octave::math::sparse_lu::Pc_vec
OCTAVE_API ColumnVector Pc_vec(void) const
Definition: sparse-lu.cc:964
octave::math::sparse_lu::m_L
lu_type m_L
Definition: sparse-lu.h:117
octave::math::sparse_lu::lu_elt_type
lu_type::element_type lu_elt_type
Definition: sparse-lu.h:53
octave::math::sparse_lu::m_cond
double m_cond
Definition: sparse-lu.h:121
octave::math::sparse_lu::sparse_lu
sparse_lu(void)
Definition: sparse-lu.h:55
octave::math::sparse_lu::m_Q
MArray< octave_idx_type > m_Q
Definition: sparse-lu.h:124
octave::math::sparse_lu::m_R
SparseMatrix m_R
Definition: sparse-lu.h:119
octave::math::sparse_lu::Pr_mat
OCTAVE_API PermMatrix Pr_mat(void) const
Definition: sparse-lu.cc:937
octave::math::sparse_lu::Pr
OCTAVE_API SparseMatrix Pr(void) const
Definition: sparse-lu.cc:903
octave::math::sparse_lu::m_U
lu_type m_U
Definition: sparse-lu.h:118
octave::math::sparse_lu::Pr_vec
OCTAVE_API ColumnVector Pr_vec(void) const
Definition: sparse-lu.cc:923
octave::math::sparse_lu::Pc
OCTAVE_API SparseMatrix Pc(void) const
Definition: sparse-lu.cc:944
octave::math::sparse_lu::m_P
MArray< octave_idx_type > m_P
Definition: sparse-lu.h:123
octave::math::sparse_lu::Y
OCTAVE_API lu_type Y(void) const
Definition: sparse-lu.cc:866
octave::sparse_params::get_key
static double get_key(const std::string &key)
Definition: oct-spparms.cc:87
octave_idx_type
scale
void scale(Matrix &m, double x, double y, double z)
Definition: graphics.cc:5893
OCTAVE_API
#define OCTAVE_API
Definition: main.in.cc:55
SparseMatrix
class OCTAVE_API SparseMatrix
Definition: mx-fwd.h:55
PermMatrix
class OCTAVE_API PermMatrix
Definition: mx-fwd.h:64
octave::math::umfpack_report_symbolic< double >
OCTAVE_API void umfpack_report_symbolic< double >(void *Symbolic, const double *Control)
Definition: sparse-lu.cc:258
octave::math::umfpack_report_status< double >
OCTAVE_API void umfpack_report_status< double >(double *Control, octave_idx_type status)
Definition: sparse-lu.cc:251
octave::math::umfpack_defaults< Complex >
OCTAVE_API void umfpack_defaults< Complex >(double *Control)
Definition: sparse-lu.cc:267
octave::math::umfpack_report_info< double >
OCTAVE_API void umfpack_report_info< double >(const double *Control, const double *Info)
Definition: sparse-lu.cc:216
octave::math::umfpack_report_symbolic< Complex >
OCTAVE_API void umfpack_report_symbolic< Complex >(void *Symbolic, const double *Control)
Definition: sparse-lu.cc:400
octave::math::umfpack_free_numeric< Complex >
OCTAVE_API void umfpack_free_numeric< Complex >(void **Numeric)
Definition: sparse-lu.cc:274
octave::math::umfpack_qsymbolic< double >
OCTAVE_API octave_idx_type umfpack_qsymbolic< double >(octave_idx_type n_row, octave_idx_type n_col, const octave_idx_type *Ap, const octave_idx_type *Ai, const double *Ax, const octave_idx_type *Qinit, void **Symbolic, const double *Control, double *Info)
Definition: sparse-lu.cc:195
octave::math::umfpack_report_matrix
void umfpack_report_matrix(octave_idx_type n_row, octave_idx_type n_col, const octave_idx_type *Ap, const octave_idx_type *Ai, const T *Ax, octave_idx_type col_form, const double *Control)
octave::math::umfpack_report_status< Complex >
OCTAVE_API void umfpack_report_status< Complex >(double *Control, octave_idx_type status)
Definition: sparse-lu.cc:393
octave::math::umfpack_get_numeric< double >
OCTAVE_API octave_idx_type umfpack_get_numeric< double >(octave_idx_type *Lp, octave_idx_type *Lj, double *Lx, octave_idx_type *Up, octave_idx_type *Ui, double *Ux, octave_idx_type *p, octave_idx_type *q, double *Dx, octave_idx_type *do_recip, double *Rs, void *Numeric)
Definition: sparse-lu.cc:165
octave::math::umfpack_report_numeric
void umfpack_report_numeric(void *Numeric, const double *Control)
octave::math::umfpack_free_symbolic
void umfpack_free_symbolic(void **Symbolic)
octave::math::umfpack_report_matrix< Complex >
OCTAVE_API void umfpack_report_matrix< Complex >(octave_idx_type n_row, octave_idx_type n_col, const octave_idx_type *Ap, const octave_idx_type *Ai, const Complex *Az, octave_idx_type col_form, const double *Control)
Definition: sparse-lu.cc:365
octave::math::umfpack_free_symbolic< double >
OCTAVE_API void umfpack_free_symbolic< double >(void **Symbolic)
Definition: sparse-lu.cc:145
octave::math::isnan
bool isnan(bool)
Definition: lo-mappers.h:178
octave::math::umfpack_get_lunz< Complex >
OCTAVE_API octave_idx_type umfpack_get_lunz< Complex >(octave_idx_type *lnz, octave_idx_type *unz, void *Numeric)
Definition: sparse-lu.cc:289
octave::math::umfpack_report_perm< Complex >
OCTAVE_API void umfpack_report_perm< Complex >(octave_idx_type np, const octave_idx_type *Perm, const double *Control)
Definition: sparse-lu.cc:386
octave::math::umfpack_report_control
void umfpack_report_control(const double *Control)
octave::math::umfpack_qsymbolic< Complex >
OCTAVE_API octave_idx_type umfpack_qsymbolic< Complex >(octave_idx_type n_row, octave_idx_type n_col, const octave_idx_type *Ap, const octave_idx_type *Ai, const Complex *Az, const octave_idx_type *Qinit, void **Symbolic, const double *Control, double *Info)
Definition: sparse-lu.cc:335
octave::math::umfpack_get_lunz
octave_idx_type umfpack_get_lunz(octave_idx_type *lnz, octave_idx_type *unz, void *Numeric)
octave::math::umfpack_report_control< Complex >
OCTAVE_API void umfpack_report_control< Complex >(const double *Control)
Definition: sparse-lu.cc:350
octave::math::umfpack_report_numeric< Complex >
OCTAVE_API void umfpack_report_numeric< Complex >(void *Numeric, const double *Control)
Definition: sparse-lu.cc:378
octave::math::umfpack_get_numeric< Complex >
OCTAVE_API octave_idx_type umfpack_get_numeric< Complex >(octave_idx_type *Lp, octave_idx_type *Lj, Complex *Lz, octave_idx_type *Up, octave_idx_type *Ui, Complex *Uz, octave_idx_type *p, octave_idx_type *q, double *Dz, octave_idx_type *do_recip, double *Rs, void *Numeric)
Definition: sparse-lu.cc:301
octave::math::umfpack_report_numeric< double >
OCTAVE_API void umfpack_report_numeric< double >(void *Numeric, const double *Control)
Definition: sparse-lu.cc:236
octave::math::umfpack_report_control< double >
OCTAVE_API void umfpack_report_control< double >(const double *Control)
Definition: sparse-lu.cc:209
octave::math::umfpack_free_numeric< double >
OCTAVE_API void umfpack_free_numeric< double >(void **Numeric)
Definition: sparse-lu.cc:138
octave::math::umfpack_numeric< Complex >
OCTAVE_API octave_idx_type umfpack_numeric< Complex >(const octave_idx_type *Ap, const octave_idx_type *Ai, const Complex *Az, void *Symbolic, void **Numeric, const double *Control, double *Info)
Definition: sparse-lu.cc:322
octave::math::umfpack_defaults< double >
OCTAVE_API void umfpack_defaults< double >(double *Control)
Definition: sparse-lu.cc:131
octave::math::umfpack_free_numeric
void umfpack_free_numeric(void **Numeric)
octave::math::umfpack_numeric
octave_idx_type umfpack_numeric(const octave_idx_type *Ap, const octave_idx_type *Ai, const T *Ax, void *Symbolic, void **Numeric, const double *Control, double *Info)
octave::math::umfpack_get_lunz< double >
OCTAVE_API octave_idx_type umfpack_get_lunz< double >(octave_idx_type *lnz, octave_idx_type *unz, void *Numeric)
Definition: sparse-lu.cc:153
octave::math::umfpack_free_symbolic< Complex >
OCTAVE_API void umfpack_free_symbolic< Complex >(void **Symbolic)
Definition: sparse-lu.cc:281
octave::math::umfpack_report_matrix< double >
OCTAVE_API void umfpack_report_matrix< double >(octave_idx_type n_row, octave_idx_type n_col, const octave_idx_type *Ap, const octave_idx_type *Ai, const double *Ax, octave_idx_type col_form, const double *Control)
Definition: sparse-lu.cc:224
octave::math::umfpack_numeric< double >
OCTAVE_API octave_idx_type umfpack_numeric< double >(const octave_idx_type *Ap, const octave_idx_type *Ai, const double *Ax, void *Symbolic, void **Numeric, const double *Control, double *Info)
Definition: sparse-lu.cc:183
octave::math::umfpack_report_perm
void umfpack_report_perm(octave_idx_type np, const octave_idx_type *Perm, const double *Control)
octave::math::umfpack_report_perm< double >
OCTAVE_API void umfpack_report_perm< double >(octave_idx_type np, const octave_idx_type *Perm, const double *Control)
Definition: sparse-lu.cc:244
octave::math::umfpack_report_status
void umfpack_report_status(double *Control, octave_idx_type status)
octave::math::umfpack_defaults
void umfpack_defaults(double *Control)
octave::math::umfpack_report_info< Complex >
OCTAVE_API void umfpack_report_info< Complex >(const double *Control, const double *Info)
Definition: sparse-lu.cc:357
octave::math::umfpack_get_numeric
octave_idx_type umfpack_get_numeric(octave_idx_type *Lp, octave_idx_type *Lj, T *Lx, octave_idx_type *Up, octave_idx_type *Ui, T *Ux, octave_idx_type *p, octave_idx_type *q, double *Dz_packed, octave_idx_type *do_recip, double *Rs, void *Numeric)
octave::math::umfpack_qsymbolic
octave_idx_type umfpack_qsymbolic(octave_idx_type n_row, octave_idx_type n_col, const octave_idx_type *Ap, const octave_idx_type *Ai, const T *Ax, const octave_idx_type *Qinit, void **Symbolic, const double *Control, double *Info)
octave::math::umfpack_report_symbolic
void umfpack_report_symbolic(void *Symbolic, const double *Control)
octave::math::umfpack_report_info
void umfpack_report_info(const double *Control, const double *Info)
octave::to_suitesparse_intptr
suitesparse_integer * to_suitesparse_intptr(octave_idx_type *i)
Definition: oct-sparse.cc:51
Complex
std::complex< double > Complex
Definition: oct-cmplx.h:33
OCTAVE_LOCAL_BUFFER
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition: oct-locbuf.h:44
UMFPACK_ZNAME
#define UMFPACK_ZNAME(name)
Definition: oct-sparse.h:159
UMFPACK_DNAME
#define UMFPACK_DNAME(name)
Definition: oct-sparse.h:158