GNU Octave  9.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-2024 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 
40 OCTAVE_BEGIN_NAMESPACE(octave)
41 
42 OCTAVE_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 
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 () 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 () 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 () 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 () const
938 {
939  return PermMatrix (m_P, false);
940 }
941 
942 template <typename lu_type>
943 SparseMatrix
944 sparse_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 
962 template <typename lu_type>
963 ColumnVector
964 sparse_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 
976 template <typename lu_type>
977 PermMatrix
978 sparse_lu<lu_type>::Pc_mat () 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 OCTAVE_END_NAMESPACE(math)
989 OCTAVE_END_NAMESPACE(octave)
Array::fortran_vec
T * fortran_vec()
Size of the specified dimension.
Definition: Array-base.cc:1764
Array::resize
void resize(const dim_vector &dv, const T &rfv)
Size of the specified dimension.
Definition: Array-base.cc:1023
Array::xelem
T & xelem(octave_idx_type n)
Size of the specified dimension.
Definition: Array.h:524
Array::numel
octave_idx_type numel() const
Number of elements in the array.
Definition: Array.h:414
Matrix
Definition: dMatrix.h:42
Sparse::xcidx
octave_idx_type * xcidx()
Definition: Sparse.h:602
Sparse::cidx
octave_idx_type * cidx()
Definition: Sparse.h:596
Sparse::data
T * data()
Definition: Sparse.h:574
Sparse::ridx
octave_idx_type * ridx()
Definition: Sparse.h:583
Sparse::xridx
octave_idx_type * xridx()
Definition: Sparse.h:589
dim_vector
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:94
octave_idx_type
sparse_lu
Definition: sparse-lu.h:50
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:55
sparse_lu::m_L
lu_type m_L
Definition: sparse-lu.h:117
sparse_lu::m_cond
double m_cond
Definition: sparse-lu.h:121
sparse_lu::m_Q
MArray< octave_idx_type > m_Q
Definition: sparse-lu.h:124
sparse_lu::m_R
SparseMatrix m_R
Definition: sparse-lu.h:119
sparse_lu::m_U
lu_type m_U
Definition: sparse-lu.h:118
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:53
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:123
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:5500
isnan
bool isnan(bool)
Definition: lo-mappers.h:178
OCTAVE_API
#define OCTAVE_API
Definition: main.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