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