GNU Octave 10.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-2025 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
40OCTAVE_BEGIN_NAMESPACE(octave)
41
42OCTAVE_BEGIN_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
50template <typename T>
51void
52umfpack_defaults (double *Control);
53
54template <typename T>
55void
56umfpack_free_numeric (void **Numeric);
57
58template <typename T>
59void
60umfpack_free_symbolic (void **Symbolic);
61
62template <typename T>
63octave_idx_type
64umfpack_get_lunz (octave_idx_type *lnz, octave_idx_type *unz,
65 void *Numeric);
66
67template <typename T>
68octave_idx_type
69umfpack_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
77template <typename T>
78octave_idx_type
79umfpack_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
84template <typename T>
85octave_idx_type
86umfpack_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
92template <typename T>
93void
94umfpack_report_control (const double *Control);
95
96template <typename T>
97void
98umfpack_report_info (const double *Control, const double *Info);
99
100template <typename T>
101void
102umfpack_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
108template <typename T>
109void
110umfpack_report_numeric (void *Numeric, const double *Control);
111
112template <typename T>
113void
114umfpack_report_perm (octave_idx_type np, const octave_idx_type *Perm,
115 const double *Control);
116
117template <typename T>
118void
119umfpack_report_status (double *Control, octave_idx_type status);
120
121template <typename T>
122void
123umfpack_report_symbolic (void *Symbolic, const double *Control);
124
125#if defined (HAVE_UMFPACK)
126
127// SparseMatrix Specialization.
128
129template <>
130inline OCTAVE_API void
131umfpack_defaults<double> (double *Control)
132{
133 UMFPACK_DNAME (defaults) (Control);
134}
135
136template <>
137inline OCTAVE_API void
138umfpack_free_numeric<double> (void **Numeric)
139{
140 UMFPACK_DNAME (free_numeric) (Numeric);
141}
142
143template <>
144inline OCTAVE_API void
145umfpack_free_symbolic<double> (void **Symbolic)
146{
147 UMFPACK_DNAME (free_symbolic) (Symbolic);
148}
149
150template <>
151inline OCTAVE_API octave_idx_type
152umfpack_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
162template <>
163inline OCTAVE_API octave_idx_type
164umfpack_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
180template <>
181inline OCTAVE_API octave_idx_type
182umfpack_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
192template <>
193inline OCTAVE_API octave_idx_type
194umfpack_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
207template <>
208inline OCTAVE_API void
209umfpack_report_control<double> (const double *Control)
210{
211 UMFPACK_DNAME (report_control) (Control);
212}
213
214template <>
215inline OCTAVE_API void
216umfpack_report_info<double> (const double *Control, const double *Info)
217{
218 UMFPACK_DNAME (report_info) (Control, Info);
219}
220
221template <>
222inline OCTAVE_API void
223umfpack_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
234template <>
235inline OCTAVE_API void
236umfpack_report_numeric<double> (void *Numeric, const double *Control)
237{
238 UMFPACK_DNAME (report_numeric) (Numeric, Control);
239}
240
241template <>
242inline OCTAVE_API void
243umfpack_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
249template <>
250inline OCTAVE_API void
251umfpack_report_status<double> (double *Control, octave_idx_type status)
252{
253 UMFPACK_DNAME (report_status) (Control, status);
254}
255
256template <>
257inline OCTAVE_API void
258umfpack_report_symbolic<double> (void *Symbolic, const double *Control)
259{
260 UMFPACK_DNAME (report_symbolic) (Symbolic, Control);
261}
262
263// SparseComplexMatrix specialization.
264
265template <>
266inline OCTAVE_API void
267umfpack_defaults<Complex> (double *Control)
268{
269 UMFPACK_ZNAME (defaults) (Control);
270}
271
272template <>
273inline OCTAVE_API void
274umfpack_free_numeric<Complex> (void **Numeric)
275{
276 UMFPACK_ZNAME (free_numeric) (Numeric);
277}
278
279template <>
280inline OCTAVE_API void
281umfpack_free_symbolic<Complex> (void **Symbolic)
282{
283 UMFPACK_ZNAME (free_symbolic) (Symbolic);
284}
285
286template <>
287inline OCTAVE_API octave_idx_type
288umfpack_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
298template <>
299inline OCTAVE_API octave_idx_type
300umfpack_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
319template <>
320inline OCTAVE_API octave_idx_type
321umfpack_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
332template <>
333inline OCTAVE_API octave_idx_type
334umfpack_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
348template <>
349inline OCTAVE_API void
350umfpack_report_control<Complex> (const double *Control)
351{
352 UMFPACK_ZNAME (report_control) (Control);
353}
354
355template <>
356inline OCTAVE_API void
357umfpack_report_info<Complex> (const double *Control, const double *Info)
358{
359 UMFPACK_ZNAME (report_info) (Control, Info);
360}
361
362template <>
363inline OCTAVE_API void
364umfpack_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
376template <>
377inline OCTAVE_API void
378umfpack_report_numeric<Complex> (void *Numeric, const double *Control)
379{
380 UMFPACK_ZNAME (report_numeric) (Numeric, Control);
381}
382
383template <>
384inline OCTAVE_API void
385umfpack_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
391template <>
392inline OCTAVE_API void
393umfpack_report_status<Complex> (double *Control, octave_idx_type status)
394{
395 UMFPACK_ZNAME (report_status) (Control, status);
396}
397
398template <>
399inline OCTAVE_API void
400umfpack_report_symbolic <Complex> (void *Symbolic, const double *Control)
401{
402 UMFPACK_ZNAME (report_symbolic) (Symbolic, Control);
403}
404
405#endif
406
407template <typename lu_type>
408sparse_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.rwdata ();
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.rwdata ();
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.rwdata ();
556
557 m_Q.resize (dim_vector (nc, 1));
558 octave_idx_type *q = m_Q.rwdata ();
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
617template <typename lu_type>
618sparse_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.rwdata ();
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.rwdata ();
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.rwdata ();
794
795 m_Q.resize (dim_vector (nc, 1));
796 octave_idx_type *q = m_Q.rwdata ();
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
864template <typename lu_type>
865lu_type
866sparse_lu<lu_type>::Y () 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
901template <typename lu_type>
902SparseMatrix
903sparse_lu<lu_type>::Pr () 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
921template <typename lu_type>
922ColumnVector
923sparse_lu<lu_type>::Pr_vec () 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
935template <typename lu_type>
936PermMatrix
937sparse_lu<lu_type>::Pr_mat () const
938{
939 return PermMatrix (m_P, false);
940}
941
942template <typename lu_type>
943SparseMatrix
944sparse_lu<lu_type>::Pc () 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
962template <typename lu_type>
963ColumnVector
964sparse_lu<lu_type>::Pc_vec () 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
976template <typename lu_type>
977PermMatrix
978sparse_lu<lu_type>::Pc_mat () const
979{
980 return PermMatrix (m_Q, true);
981}
982
983// Instantiations we need.
984template class OCTAVE_API sparse_lu<SparseMatrix>;
985
986template class OCTAVE_API sparse_lu<SparseComplexMatrix>;
987
988OCTAVE_END_NAMESPACE(math)
989OCTAVE_END_NAMESPACE(octave)
Array::xelem
T & xelem(octave_idx_type n)
Size of the specified dimension.
Definition Array.h:525
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:1787
Array::numel
octave_idx_type numel() const
Number of elements in the array.
Definition Array.h:418
Sparse::cidx
octave_idx_type * cidx()
Definition Sparse.h:593
Sparse::data
T * data()
Definition Sparse.h:571
Sparse::ridx
octave_idx_type * ridx()
Definition Sparse.h:580
Sparse::xcidx
octave_idx_type * xcidx()
Definition Sparse.h:599
Sparse::xridx
octave_idx_type * xridx()
Definition Sparse.h:586
dim_vector
Vector representing the dimensions (size) of an Array.
Definition dim-vector.h:90
octave_idx_type
sparse_lu
Definition sparse-lu.h:48
sparse_lu::Y
lu_type Y() const
Definition sparse-lu.cc:866
sparse_lu::Pr
SparseMatrix Pr() const
Definition sparse-lu.cc:903
sparse_lu::sparse_lu
sparse_lu()
Definition sparse-lu.h:53
sparse_lu::m_L
lu_type m_L
Definition sparse-lu.h:115
sparse_lu::m_cond
double m_cond
Definition sparse-lu.h:119
sparse_lu::m_Q
MArray< octave_idx_type > m_Q
Definition sparse-lu.h:122
sparse_lu::m_R
SparseMatrix m_R
Definition sparse-lu.h:117
sparse_lu::m_U
lu_type m_U
Definition sparse-lu.h:116
sparse_lu::Pr_mat
PermMatrix Pr_mat() const
Definition sparse-lu.cc:937
sparse_lu::lu_elt_type
lu_type::element_type lu_elt_type
Definition sparse-lu.h:51
sparse_lu::Pc_mat
PermMatrix Pc_mat() const
Definition sparse-lu.cc:978
sparse_lu::m_P
MArray< octave_idx_type > m_P
Definition sparse-lu.h:121
sparse_lu::Pc
SparseMatrix Pc() const
Definition sparse-lu.cc:944
sparse_lu::Pc_vec
ColumnVector Pc_vec() const
Definition sparse-lu.cc:964
sparse_lu::Pr_vec
ColumnVector Pr_vec() const
Definition sparse-lu.cc:923
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:5731
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:216
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:153
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:145
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:301
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:183
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:335
umfpack_report_status< double >
void umfpack_report_status< double >(double *Control, octave_idx_type status)
Definition sparse-lu.cc:251
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:131
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:386
umfpack_report_control< Complex >
void umfpack_report_control< Complex >(const double *Control)
Definition sparse-lu.cc:350
umfpack_free_symbolic< Complex >
void umfpack_free_symbolic< Complex >(void **Symbolic)
Definition sparse-lu.cc:281
umfpack_defaults< Complex >
void umfpack_defaults< Complex >(double *Control)
Definition sparse-lu.cc:267
umfpack_free_numeric< double >
void umfpack_free_numeric< double >(void **Numeric)
Definition sparse-lu.cc:138
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:365
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:224
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:258
umfpack_report_control< double >
void umfpack_report_control< double >(const double *Control)
Definition sparse-lu.cc:209
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:289
umfpack_report_numeric< double >
void umfpack_report_numeric< double >(void *Numeric, const double *Control)
Definition sparse-lu.cc:236
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:165
umfpack_free_numeric< Complex >
void umfpack_free_numeric< Complex >(void **Numeric)
Definition sparse-lu.cc:274
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:393
umfpack_report_numeric< Complex >
void umfpack_report_numeric< Complex >(void *Numeric, const double *Control)
Definition sparse-lu.cc:378
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:195
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:357
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:322
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:244
umfpack_report_symbolic< Complex >
void umfpack_report_symbolic< Complex >(void *Symbolic, const double *Control)
Definition sparse-lu.cc:400