GNU Octave  6.2.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-2021 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 namespace octave
41 {
42  namespace math
43  {
44  // Wrappers for SuiteSparse (formerly UMFPACK) functions that have
45  // different names depending on the sparse matrix data type.
46  //
47  // All of these functions must be specialized to forward to the correct
48  // SuiteSparse functions.
49 
50  template <typename T>
51  void
52  umfpack_defaults (double *Control);
53 
54  template <typename T>
55  void
56  umfpack_free_numeric (void **Numeric);
57 
58  template <typename T>
59  void
60  umfpack_free_symbolic (void **Symbolic);
61 
62  template <typename T>
63  octave_idx_type
64  umfpack_get_lunz (octave_idx_type *lnz, octave_idx_type *unz,
65  void *Numeric);
66 
67  template <typename T>
68  octave_idx_type
69  umfpack_get_numeric (octave_idx_type *Lp, octave_idx_type *Lj,
70  T *Lx, // Or Lz_packed
71  octave_idx_type *Up, octave_idx_type *Ui,
72  T *Ux, // Or Uz_packed
73  octave_idx_type *p, octave_idx_type *q,
74  double *Dz_packed, octave_idx_type *do_recip,
75  double *Rs, void *Numeric);
76 
77  template <typename T>
78  octave_idx_type
79  umfpack_numeric (const octave_idx_type *Ap, const octave_idx_type *Ai,
80  const T *Ax, // Or Az_packed
81  void *Symbolic, void **Numeric,
82  const double *Control, double *Info);
83 
84  template <typename T>
85  octave_idx_type
86  umfpack_qsymbolic (octave_idx_type n_row, octave_idx_type n_col,
87  const octave_idx_type *Ap, const octave_idx_type *Ai,
88  const T *Ax, // Or Az_packed
89  const octave_idx_type *Qinit, void **Symbolic,
90  const double *Control, double *Info);
91 
92  template <typename T>
93  void
94  umfpack_report_control (const double *Control);
95 
96  template <typename T>
97  void
98  umfpack_report_info (const double *Control, const double *Info);
99 
100  template <typename T>
101  void
102  umfpack_report_matrix (octave_idx_type n_row, octave_idx_type n_col,
103  const octave_idx_type *Ap,
104  const octave_idx_type *Ai,
105  const T *Ax, // Or Az_packed
106  octave_idx_type col_form, const double *Control);
107 
108  template <typename T>
109  void
110  umfpack_report_numeric (void *Numeric, const double *Control);
111 
112  template <typename T>
113  void
114  umfpack_report_perm (octave_idx_type np, const octave_idx_type *Perm,
115  const double *Control);
116 
117  template <typename T>
118  void
119  umfpack_report_status (double *Control, octave_idx_type status);
120 
121  template <typename T>
122  void
123  umfpack_report_symbolic (void *Symbolic, const double *Control);
124 
125 #if defined (HAVE_UMFPACK)
126 
127  // SparseMatrix Specialization.
128 
129  template <>
130  inline void
131  umfpack_defaults<double> (double *Control)
132  {
133  UMFPACK_DNAME (defaults) (Control);
134  }
135 
136  template <>
137  inline void
138  umfpack_free_numeric<double> (void **Numeric)
139  {
140  UMFPACK_DNAME (free_numeric) (Numeric);
141  }
142 
143  template <>
144  inline void
145  umfpack_free_symbolic<double> (void **Symbolic)
146  {
147  UMFPACK_DNAME (free_symbolic) (Symbolic);
148  }
149 
150  template <>
151  inline 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_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_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_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 void
209  umfpack_report_control<double> (const double *Control)
210  {
211  UMFPACK_DNAME (report_control) (Control);
212  }
213 
214  template <>
215  inline 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 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 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 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 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 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 void
267  umfpack_defaults<Complex> (double *Control)
268  {
269  UMFPACK_ZNAME (defaults) (Control);
270  }
271 
272  template <>
273  inline void
274  umfpack_free_numeric<Complex> (void **Numeric)
275  {
276  UMFPACK_ZNAME (free_numeric) (Numeric);
277  }
278 
279  template <>
280  inline void
281  umfpack_free_symbolic<Complex> (void **Symbolic)
282  {
283  UMFPACK_ZNAME (free_symbolic) (Symbolic);
284  }
285 
286  template <>
287  inline 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_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_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_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 void
350  umfpack_report_control<Complex> (const double *Control)
351  {
352  UMFPACK_ZNAME (report_control) (Control);
353  }
354 
355  template <>
356  inline 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 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 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 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 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 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  : Lfact (), Ufact (), Rfact (), cond (0), P (), 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 = octave_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 = octave_sparse_params::get_key ("piv_tol");
438  if (! math::isnan (tmp))
439  Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
440 
441  tmp = octave_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 = octave_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  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  Lfact = lu_type (n_inner, nr,
527  static_cast<octave_idx_type> (1));
528  else
529  Lfact = lu_type (n_inner, nr, lnz);
530 
531  octave_idx_type *Ltp = Lfact.cidx ();
532  octave_idx_type *Ltj = Lfact.ridx ();
533  lu_elt_type *Ltx = Lfact.data ();
534 
535  if (unz < 1)
536  Ufact = lu_type (n_inner, nc,
537  static_cast<octave_idx_type> (1));
538  else
539  Ufact = lu_type (n_inner, nc, unz);
540 
541  octave_idx_type *Up = Ufact.cidx ();
542  octave_idx_type *Uj = Ufact.ridx ();
543  lu_elt_type *Ux = Ufact.data ();
544 
545  Rfact = SparseMatrix (nr, nr, nr);
546  for (octave_idx_type i = 0; i < nr; i++)
547  {
548  Rfact.xridx (i) = i;
549  Rfact.xcidx (i) = i;
550  }
551  Rfact.xcidx (nr) = nr;
552  double *Rx = Rfact.data ();
553 
554  P.resize (dim_vector (nr, 1));
555  octave_idx_type *p = P.fortran_vec ();
556 
557  Q.resize (dim_vector (nc, 1));
558  octave_idx_type *q = 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  Lfact = Lfact.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  Lfact.cidx (),
586  Lfact.ridx (),
587  Lfact.data (),
588  static_cast<octave_idx_type> (1),
589  control);
590  umfpack_report_matrix<lu_elt_type> (n_inner, nc,
591  Ufact.cidx (),
592  Ufact.ridx (),
593  Ufact.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  : Lfact (), Ufact (), Rfact (), cond (0), P (), 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 = octave_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 = octave_sparse_params::get_key ("piv_tol");
655  if (! math::isnan (tmp))
656  Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
657 
658  tmp = octave_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 = octave_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  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  Lfact = lu_type (n_inner, nr,
765  static_cast<octave_idx_type> (1));
766  else
767  Lfact = lu_type (n_inner, nr, lnz);
768 
769  octave_idx_type *Ltp = Lfact.cidx ();
770  octave_idx_type *Ltj = Lfact.ridx ();
771  lu_elt_type *Ltx = Lfact.data ();
772 
773  if (unz < 1)
774  Ufact = lu_type (n_inner, nc,
775  static_cast<octave_idx_type> (1));
776  else
777  Ufact = lu_type (n_inner, nc, unz);
778 
779  octave_idx_type *Up = Ufact.cidx ();
780  octave_idx_type *Uj = Ufact.ridx ();
781  lu_elt_type *Ux = Ufact.data ();
782 
783  Rfact = SparseMatrix (nr, nr, nr);
784  for (octave_idx_type i = 0; i < nr; i++)
785  {
786  Rfact.xridx (i) = i;
787  Rfact.xcidx (i) = i;
788  }
789  Rfact.xcidx (nr) = nr;
790  double *Rx = Rfact.data ();
791 
792  P.resize (dim_vector (nr, 1));
793  octave_idx_type *p = P.fortran_vec ();
794 
795  Q.resize (dim_vector (nc, 1));
796  octave_idx_type *q = 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  Lfact = Lfact.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  Lfact.cidx (),
824  Lfact.ridx (),
825  Lfact.data (),
826  static_cast<octave_idx_type> (1),
827  control);
828  umfpack_report_matrix<lu_elt_type> (n_inner, nc,
829  Ufact.cidx (),
830  Ufact.ridx (),
831  Ufact.data (),
832  static_cast<octave_idx_type> (1),
833  control);
834  umfpack_report_perm<lu_elt_type> (nr, p, control);
835  umfpack_report_perm<lu_elt_type> (nc, q, control);
836  }
837 
838  umfpack_report_info<lu_elt_type> (control, info);
839  }
840  }
841  }
842 
843  if (udiag)
844  (*current_liboctave_error_handler)
845  ("Option udiag of incomplete LU not implemented");
846 
847 #else
848 
849  octave_unused_parameter (a);
850  octave_unused_parameter (Qinit);
851  octave_unused_parameter (piv_thres);
852  octave_unused_parameter (scale);
853  octave_unused_parameter (FixedQ);
854  octave_unused_parameter (droptol);
855  octave_unused_parameter (milu);
856  octave_unused_parameter (udiag);
857 
858  (*current_liboctave_error_handler)
859  ("support for UMFPACK was unavailable or disabled when liboctave was built");
860 
861 #endif
862  }
863 
864  template <typename lu_type>
865  lu_type
866  sparse_lu<lu_type>::Y (void) const
867  {
868  octave_idx_type nr = Lfact.rows ();
869  octave_idx_type nz = Lfact.cols ();
870  octave_idx_type nc = Ufact.cols ();
871 
872  lu_type Yout (nr, nc, Lfact.nnz () + Ufact.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 = Ufact.cidx (j); i < Ufact.cidx (j + 1); i++)
879  {
880  Yout.xridx (ii) = Ufact.ridx (i);
881  Yout.xdata (ii++) = Ufact.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 = Lfact.cidx (j) + 1;
888  i < Lfact.cidx (j +1); i++)
889  {
890  Yout.xridx (ii) = Lfact.ridx (i);
891  Yout.xdata (ii++) = Lfact.data (i);
892  }
893  }
894 
895  Yout.xcidx (j + 1) = ii;
896  }
897 
898  return Yout;
899  }
900 
901  template <typename lu_type>
902  SparseMatrix
903  sparse_lu<lu_type>::Pr (void) const
904  {
905  octave_idx_type nr = Lfact.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 (P (i)) = i;
913  Pout.data (i) = 1;
914  }
915 
916  Pout.cidx (nr) = nr;
917 
918  return Pout;
919  }
920 
921  template <typename lu_type>
922  ColumnVector
923  sparse_lu<lu_type>::Pr_vec (void) const
924  {
925  octave_idx_type nr = Lfact.rows ();
926 
927  ColumnVector Pout (nr);
928 
929  for (octave_idx_type i = 0; i < nr; i++)
930  Pout.xelem (i) = static_cast<double> (P(i) + 1);
931 
932  return Pout;
933  }
934 
935  template <typename lu_type>
936  PermMatrix
937  sparse_lu<lu_type>::Pr_mat (void) const
938  {
939  return PermMatrix (P, false);
940  }
941 
942  template <typename lu_type>
943  SparseMatrix
944  sparse_lu<lu_type>::Pc (void) const
945  {
946  octave_idx_type nc = Ufact.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) = Q (i);
954  Pout.data (i) = 1;
955  }
956 
957  Pout.cidx (nc) = nc;
958 
959  return Pout;
960  }
961 
962  template <typename lu_type>
963  ColumnVector
964  sparse_lu<lu_type>::Pc_vec (void) const
965  {
966  octave_idx_type nc = Ufact.cols ();
967 
968  ColumnVector Pout (nc);
969 
970  for (octave_idx_type i = 0; i < nc; i++)
971  Pout.xelem (i) = static_cast<double> (Q(i) + 1);
972 
973  return Pout;
974  }
975 
976  template <typename lu_type>
977  PermMatrix
978  sparse_lu<lu_type>::Pc_mat (void) const
979  {
980  return PermMatrix (Q, true);
981  }
982 
983  // Instantiations we need.
984 
985  template class sparse_lu<SparseMatrix>;
986 
987  template class sparse_lu<SparseComplexMatrix>;
988  }
989 }
Array::resize
void resize(const dim_vector &dv, const T &rfv)
Size of the specified dimension.
Definition: Array.cc:1011
Array::xelem
T & xelem(octave_idx_type n)
Size of the specified dimension.
Definition: Array.h:469
Array::numel
octave_idx_type numel(void) const
Number of elements in the array.
Definition: Array.h:377
Array::fortran_vec
const T * fortran_vec(void) const
Size of the specified dimension.
Definition: Array.h:583
Matrix
Definition: dMatrix.h:42
Sparse::xridx
octave_idx_type * xridx(void)
Definition: Sparse.h:485
Sparse::data
T * data(void)
Definition: Sparse.h:470
Sparse::cidx
octave_idx_type * cidx(void)
Definition: Sparse.h:492
Sparse::xcidx
octave_idx_type * xcidx(void)
Definition: Sparse.h:498
Sparse::ridx
octave_idx_type * ridx(void)
Definition: Sparse.h:479
dim_vector
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:95
octave::math::sparse_lu
Definition: sparse-lu.h:50
octave::math::sparse_lu::Pc_mat
PermMatrix Pc_mat(void) const
Definition: sparse-lu.cc:978
octave::math::sparse_lu::Pc_vec
ColumnVector Pc_vec(void) const
Definition: sparse-lu.cc:964
octave::math::sparse_lu::lu_elt_type
lu_type::element_type lu_elt_type
Definition: sparse-lu.h:53
octave::math::sparse_lu::sparse_lu
sparse_lu(void)
Definition: sparse-lu.h:55
octave::math::sparse_lu::cond
double cond
Definition: sparse-lu.h:119
octave::math::sparse_lu::Q
MArray< octave_idx_type > Q
Definition: sparse-lu.h:122
octave::math::sparse_lu::Pr_mat
PermMatrix Pr_mat(void) const
Definition: sparse-lu.cc:937
octave::math::sparse_lu::Lfact
lu_type Lfact
Definition: sparse-lu.h:115
octave::math::sparse_lu::Pr
SparseMatrix Pr(void) const
Definition: sparse-lu.cc:903
octave::math::sparse_lu::Rfact
SparseMatrix Rfact
Definition: sparse-lu.h:117
octave::math::sparse_lu::Pr_vec
ColumnVector Pr_vec(void) const
Definition: sparse-lu.cc:923
octave::math::sparse_lu::Pc
SparseMatrix Pc(void) const
Definition: sparse-lu.cc:944
octave::math::sparse_lu::P
MArray< octave_idx_type > P
Definition: sparse-lu.h:121
octave::math::sparse_lu::Ufact
lu_type Ufact
Definition: sparse-lu.h:116
octave::math::sparse_lu::Y
lu_type Y(void) const
Definition: sparse-lu.cc:866
octave_idx_type
octave_sparse_params::get_key
static double get_key(const std::string &key)
Definition: oct-spparms.cc:93
scale
void scale(Matrix &m, double x, double y, double z)
Definition: graphics.cc:5846
Q
F77_RET_T const F77_INT const F77_INT const F77_INT F77_DBLE const F77_INT F77_DBLE const F77_INT F77_DBLE * Q
Definition: lo-lapack-proto.h:793
octave::math::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
octave::math::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
octave::math::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
octave::math::umfpack_report_control< Complex >
void umfpack_report_control< Complex >(const double *Control)
Definition: sparse-lu.cc:350
octave::math::umfpack_defaults< double >
void umfpack_defaults< double >(double *Control)
Definition: sparse-lu.cc:131
octave::math::umfpack_free_numeric< Complex >
void umfpack_free_numeric< Complex >(void **Numeric)
Definition: sparse-lu.cc:274
octave::math::umfpack_report_matrix
void umfpack_report_matrix(octave_idx_type n_row, octave_idx_type n_col, const octave_idx_type *Ap, const octave_idx_type *Ai, const T *Ax, octave_idx_type col_form, const double *Control)
octave::math::umfpack_report_symbolic< double >
void umfpack_report_symbolic< double >(void *Symbolic, const double *Control)
Definition: sparse-lu.cc:258
octave::math::umfpack_free_symbolic< Complex >
void umfpack_free_symbolic< Complex >(void **Symbolic)
Definition: sparse-lu.cc:281
octave::math::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
octave::math::umfpack_report_info< Complex >
void umfpack_report_info< Complex >(const double *Control, const double *Info)
Definition: sparse-lu.cc:357
octave::math::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
octave::math::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
octave::math::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
octave::math::umfpack_report_numeric
void umfpack_report_numeric(void *Numeric, const double *Control)
octave::math::umfpack_free_symbolic
void umfpack_free_symbolic(void **Symbolic)
octave::math::isnan
bool isnan(bool)
Definition: lo-mappers.h:178
octave::math::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
octave::math::umfpack_report_control< double >
void umfpack_report_control< double >(const double *Control)
Definition: sparse-lu.cc:209
octave::math::umfpack_report_numeric< double >
void umfpack_report_numeric< double >(void *Numeric, const double *Control)
Definition: sparse-lu.cc:236
octave::math::umfpack_report_control
void umfpack_report_control(const double *Control)
octave::math::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
octave::math::umfpack_defaults< Complex >
void umfpack_defaults< Complex >(double *Control)
Definition: sparse-lu.cc:267
octave::math::umfpack_get_lunz
octave_idx_type umfpack_get_lunz(octave_idx_type *lnz, octave_idx_type *unz, void *Numeric)
octave::math::umfpack_report_status< Complex >
void umfpack_report_status< Complex >(double *Control, octave_idx_type status)
Definition: sparse-lu.cc:393
octave::math::umfpack_report_numeric< Complex >
void umfpack_report_numeric< Complex >(void *Numeric, const double *Control)
Definition: sparse-lu.cc:378
octave::math::umfpack_free_symbolic< double >
void umfpack_free_symbolic< double >(void **Symbolic)
Definition: sparse-lu.cc:145
octave::math::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
octave::math::umfpack_free_numeric
void umfpack_free_numeric(void **Numeric)
octave::math::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
octave::math::umfpack_numeric
octave_idx_type umfpack_numeric(const octave_idx_type *Ap, const octave_idx_type *Ai, const T *Ax, void *Symbolic, void **Numeric, const double *Control, double *Info)
octave::math::umfpack_report_info< double >
void umfpack_report_info< double >(const double *Control, const double *Info)
Definition: sparse-lu.cc:216
octave::math::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
octave::math::umfpack_report_status< double >
void umfpack_report_status< double >(double *Control, octave_idx_type status)
Definition: sparse-lu.cc:251
octave::math::umfpack_report_symbolic< Complex >
void umfpack_report_symbolic< Complex >(void *Symbolic, const double *Control)
Definition: sparse-lu.cc:400
octave::math::umfpack_report_perm
void umfpack_report_perm(octave_idx_type np, const octave_idx_type *Perm, const double *Control)
octave::math::umfpack_report_status
void umfpack_report_status(double *Control, octave_idx_type status)
octave::math::umfpack_free_numeric< double >
void umfpack_free_numeric< double >(void **Numeric)
Definition: sparse-lu.cc:138
octave::math::umfpack_defaults
void umfpack_defaults(double *Control)
octave::math::umfpack_get_numeric
octave_idx_type umfpack_get_numeric(octave_idx_type *Lp, octave_idx_type *Lj, T *Lx, octave_idx_type *Up, octave_idx_type *Ui, T *Ux, octave_idx_type *p, octave_idx_type *q, double *Dz_packed, octave_idx_type *do_recip, double *Rs, void *Numeric)
octave::math::umfpack_qsymbolic
octave_idx_type umfpack_qsymbolic(octave_idx_type n_row, octave_idx_type n_col, const octave_idx_type *Ap, const octave_idx_type *Ai, const T *Ax, const octave_idx_type *Qinit, void **Symbolic, const double *Control, double *Info)
octave::math::umfpack_report_symbolic
void umfpack_report_symbolic(void *Symbolic, const double *Control)
octave::math::umfpack_report_info
void umfpack_report_info(const double *Control, const double *Info)
octave::suitesparse_integer
int suitesparse_integer
Definition: oct-sparse.h:171
octave::to_suitesparse_intptr
suitesparse_integer * to_suitesparse_intptr(octave_idx_type *i)
Definition: oct-sparse.cc:51
Complex
std::complex< double > Complex
Definition: oct-cmplx.h:33
OCTAVE_LOCAL_BUFFER
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition: oct-locbuf.h:44
UMFPACK_ZNAME
#define UMFPACK_ZNAME(name)
Definition: oct-sparse.h:158
UMFPACK_DNAME
#define UMFPACK_DNAME(name)
Definition: oct-sparse.h:157