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 <cstddef>
31 
32 #include "CSparse.h"
33 #include "MatrixType.h"
34 #include "dRowVector.h"
35 #include "dSparse.h"
36 #include "lo-error.h"
37 #include "oct-cmplx.h"
38 #include "oct-refcount.h"
39 #include "oct-sparse.h"
40 #include "oct-spparms.h"
41 #include "quit.h"
42 #include "sparse-chol.h"
43 #include "sparse-util.h"
44 
45 namespace octave
46 {
47   namespace math
48   {
49     template <typename chol_type>
50     class sparse_chol<chol_type>::sparse_chol_rep
51     {
52     public:
53 
sparse_chol_rep(void)54       sparse_chol_rep (void)
55         : count (1), is_pd (false), minor_p (0), perms (), cond (0)
56 #if defined (HAVE_CHOLMOD)
57         , Lsparse (nullptr), Common ()
58 #endif
59       { }
60 
sparse_chol_rep(const chol_type & a,bool natural,bool force)61       sparse_chol_rep (const chol_type& a, bool natural, bool force)
62         : count (1), is_pd (false), minor_p (0), perms (), cond (0)
63 #if defined (HAVE_CHOLMOD)
64         , Lsparse (nullptr), Common ()
65 #endif
66       {
67         init (a, natural, force);
68       }
69 
sparse_chol_rep(const chol_type & a,octave_idx_type & info,bool natural,bool force)70       sparse_chol_rep (const chol_type& a, octave_idx_type& info,
71                        bool natural, bool force)
72         : count (1), is_pd (false), minor_p (0), perms (), cond (0)
73 #if defined (HAVE_CHOLMOD)
74         , Lsparse (nullptr), Common ()
75 #endif
76       {
77         info = init (a, natural, force);
78       }
79 
80       // No copying!
81 
82       sparse_chol_rep (const sparse_chol_rep&) = delete;
83 
84       sparse_chol_rep& operator = (const sparse_chol_rep&) = delete;
85 
~sparse_chol_rep(void)86       ~sparse_chol_rep (void)
87       {
88 #if defined (HAVE_CHOLMOD)
89         if (Lsparse)
90           CHOLMOD_NAME (free_sparse) (&Lsparse, &Common);
91 
92         CHOLMOD_NAME(finish) (&Common);
93 #endif
94       }
95 
96 #if defined (HAVE_CHOLMOD)
L(void) const97       cholmod_sparse * L (void) const
98       {
99         return Lsparse;
100       }
101 #endif
102 
P(void) const103       octave_idx_type P (void) const
104       {
105 #if defined (HAVE_CHOLMOD)
106         return (minor_p == static_cast<octave_idx_type> (Lsparse->ncol) ?
107                 0 : minor_p + 1);
108 #else
109         return 0;
110 #endif
111       }
112 
perm(void) const113       RowVector perm (void) const { return perms + 1; }
114 
115       SparseMatrix Q (void) const;
116 
is_positive_definite(void) const117       bool is_positive_definite (void) const { return is_pd; }
118 
rcond(void) const119       double rcond (void) const { return cond; }
120 
121       refcount<octave_idx_type> count;
122 
123     private:
124 
125       bool is_pd;
126 
127       octave_idx_type minor_p;
128 
129       RowVector perms;
130 
131       double cond;
132 
133 #if defined (HAVE_CHOLMOD)
134       cholmod_sparse *Lsparse;
135 
136       cholmod_common Common;
137 
138       void drop_zeros (const cholmod_sparse *S);
139 #endif
140 
141       octave_idx_type init (const chol_type& a, bool natural, bool force);
142     };
143 
144 #if defined (HAVE_CHOLMOD)
145 
146     // Can't use CHOLMOD_NAME(drop)(0.0, S, cm) because it doesn't treat
147     // complex matrices.
148 
149     template <typename chol_type>
150     void
drop_zeros(const cholmod_sparse * S)151     sparse_chol<chol_type>::sparse_chol_rep::drop_zeros (const cholmod_sparse *S)
152     {
153       if (! S)
154         return;
155 
156       octave_idx_type *Sp = static_cast<octave_idx_type *>(S->p);
157       octave_idx_type *Si = static_cast<octave_idx_type *>(S->i);
158       chol_elt *Sx = static_cast<chol_elt *>(S->x);
159 
160       octave_idx_type pdest = 0;
161       octave_idx_type ncol = S->ncol;
162 
163       for (octave_idx_type k = 0; k < ncol; k++)
164         {
165           octave_idx_type p = Sp[k];
166           octave_idx_type pend = Sp[k+1];
167           Sp[k] = pdest;
168 
169           for (; p < pend; p++)
170             {
171               chol_elt sik = Sx[p];
172 
173               if (CHOLMOD_IS_NONZERO (sik))
174                 {
175                   if (p != pdest)
176                     {
177                       Si[pdest] = Si[p];
178                       Sx[pdest] = sik;
179                     }
180 
181                   pdest++;
182                 }
183             }
184         }
185 
186       Sp[ncol] = pdest;
187     }
188 
189     // Must provide a specialization for this function.
190     template <typename T>
191     int
192     get_xtype (void);
193 
194     template <>
195     inline int
get_xtype(void)196     get_xtype<double> (void)
197     {
198       return CHOLMOD_REAL;
199     }
200 
201     template <>
202     inline int
get_xtype(void)203     get_xtype<Complex> (void)
204     {
205       return CHOLMOD_COMPLEX;
206     }
207 
208 #endif
209 
210     template <typename chol_type>
211     octave_idx_type
init(const chol_type & a,bool natural,bool force)212     sparse_chol<chol_type>::sparse_chol_rep::init (const chol_type& a,
213                                                    bool natural, bool force)
214     {
215       volatile octave_idx_type info = 0;
216 
217 #if defined (HAVE_CHOLMOD)
218 
219       octave_idx_type a_nr = a.rows ();
220       octave_idx_type a_nc = a.cols ();
221 
222       if (a_nr != a_nc)
223         (*current_liboctave_error_handler)
224           ("sparse_chol requires square matrix");
225 
226       cholmod_common *cm = &Common;
227 
228       // Setup initial parameters
229 
230       CHOLMOD_NAME(start) (cm);
231       cm->prefer_zomplex = false;
232 
233       double spu = octave_sparse_params::get_key ("spumoni");
234 
235       if (spu == 0.)
236         {
237           cm->print = -1;
238           SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, nullptr);
239         }
240       else
241         {
242           cm->print = static_cast<int> (spu) + 2;
243           SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function,
244                                    &SparseCholPrint);
245         }
246 
247       cm->error_handler = &SparseCholError;
248 
249       SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide,
250                                 divcomplex);
251 
252       SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
253 
254       cm->final_asis = false;
255       cm->final_super = false;
256       cm->final_ll = true;
257       cm->final_pack = true;
258       cm->final_monotonic = true;
259       cm->final_resymbol = false;
260 
261       cholmod_sparse A;
262       cholmod_sparse *ac = &A;
263       double dummy;
264 
265       ac->nrow = a_nr;
266       ac->ncol = a_nc;
267 
268       ac->p = a.cidx ();
269       ac->i = a.ridx ();
270       ac->nzmax = a.nnz ();
271       ac->packed = true;
272       ac->sorted = true;
273       ac->nz = nullptr;
274 #if defined (OCTAVE_ENABLE_64)
275       ac->itype = CHOLMOD_LONG;
276 #else
277       ac->itype = CHOLMOD_INT;
278 #endif
279       ac->dtype = CHOLMOD_DOUBLE;
280       ac->stype = 1;
281       ac->xtype = get_xtype<chol_elt> ();
282 
283       if (a_nr < 1)
284         ac->x = &dummy;
285       else
286         ac->x = a.data ();
287 
288       // use natural ordering if no q output parameter
289       if (natural)
290         {
291           cm->nmethods = 1;
292           cm->method[0].ordering = CHOLMOD_NATURAL;
293           cm->postorder = false;
294         }
295 
296       cholmod_factor *Lfactor;
297       BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
298       Lfactor = CHOLMOD_NAME(analyze) (ac, cm);
299       CHOLMOD_NAME(factorize) (ac, Lfactor, cm);
300       END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
301 
302       is_pd = cm->status == CHOLMOD_OK;
303       info = (is_pd ? 0 : cm->status);
304 
305       if (is_pd || force)
306         {
307           BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
308           cond = CHOLMOD_NAME(rcond) (Lfactor, cm);
309           END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
310 
311           minor_p = Lfactor->minor;
312 
313           BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
314           Lsparse = CHOLMOD_NAME(factor_to_sparse) (Lfactor, cm);
315           END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
316 
317           if (minor_p > 0 && minor_p < a_nr)
318             {
319               std::size_t n1 = a_nr + 1;
320               Lsparse->p = CHOLMOD_NAME(realloc) (minor_p+1,
321                                                   sizeof(octave_idx_type),
322                                                   Lsparse->p, &n1, cm);
323               BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
324               CHOLMOD_NAME(reallocate_sparse)
325                 (static_cast<octave_idx_type *>(Lsparse->p)[minor_p],
326                  Lsparse, cm);
327               END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
328 
329               Lsparse->ncol = minor_p;
330             }
331 
332           drop_zeros (Lsparse);
333 
334           if (! natural)
335             {
336               perms.resize (a_nr);
337               for (octave_idx_type i = 0; i < a_nr; i++)
338                 perms(i) = static_cast<octave_idx_type *>(Lfactor->Perm)[i];
339             }
340         }
341 
342       // NAME used to prefix statistics report from print_common
343       static char blank_name[] = " ";
344 
345       BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
346       CHOLMOD_NAME(print_common) (blank_name, cm);
347       CHOLMOD_NAME(free_factor) (&Lfactor, cm);
348       END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
349 
350       return info;
351 
352 #else
353 
354       octave_unused_parameter (a);
355       octave_unused_parameter (natural);
356       octave_unused_parameter (force);
357 
358       (*current_liboctave_error_handler)
359         ("support for CHOLMOD was unavailable or disabled when liboctave was built");
360 
361       return info;
362 
363 #endif
364     }
365 
366     template <typename chol_type>
367     SparseMatrix
Q(void) const368     sparse_chol<chol_type>::sparse_chol_rep::Q (void) const
369     {
370 #if defined (HAVE_CHOLMOD)
371 
372       octave_idx_type n = Lsparse->nrow;
373       SparseMatrix p (n, n, n);
374 
375       for (octave_idx_type i = 0; i < n; i++)
376         {
377           p.xcidx (i) = i;
378           p.xridx (i) = static_cast<octave_idx_type> (perms (i));
379           p.xdata (i) = 1;
380         }
381 
382       p.xcidx (n) = n;
383 
384       return p;
385 
386 #else
387 
388       return SparseMatrix ();
389 
390 #endif
391     }
392 
393     template <typename chol_type>
sparse_chol(void)394     sparse_chol<chol_type>::sparse_chol (void)
395       : rep (new typename sparse_chol<chol_type>::sparse_chol_rep ())
396     { }
397 
398     template <typename chol_type>
sparse_chol(const chol_type & a,bool natural,bool force)399     sparse_chol<chol_type>::sparse_chol (const chol_type& a, bool natural,
400                                          bool force)
401       : rep (new typename
402              sparse_chol<chol_type>::sparse_chol_rep (a, natural, force))
403     { }
404 
405     template <typename chol_type>
sparse_chol(const chol_type & a,octave_idx_type & info,bool natural,bool force)406     sparse_chol<chol_type>::sparse_chol (const chol_type& a,
407                                          octave_idx_type& info,
408                                          bool natural, bool force)
409       : rep (new typename
410              sparse_chol<chol_type>::sparse_chol_rep (a, info, natural, force))
411     { }
412 
413     template <typename chol_type>
sparse_chol(const chol_type & a,octave_idx_type & info,bool natural)414     sparse_chol<chol_type>::sparse_chol (const chol_type& a,
415                                          octave_idx_type& info,
416                                          bool natural)
417       : rep (new typename
418              sparse_chol<chol_type>::sparse_chol_rep (a, info, natural, false))
419     { }
420 
421     template <typename chol_type>
sparse_chol(const chol_type & a,octave_idx_type & info)422     sparse_chol<chol_type>::sparse_chol (const chol_type& a,
423                                          octave_idx_type& info)
424       : rep (new typename
425              sparse_chol<chol_type>::sparse_chol_rep (a, info, false, false))
426     { }
427 
428     template <typename chol_type>
sparse_chol(const sparse_chol<chol_type> & a)429     sparse_chol<chol_type>::sparse_chol (const sparse_chol<chol_type>& a)
430       : rep (a.rep)
431     {
432       rep->count++;
433     }
434 
435     template <typename chol_type>
~sparse_chol(void)436     sparse_chol<chol_type>::~sparse_chol (void)
437     {
438       if (--rep->count == 0)
439         delete rep;
440     }
441 
442     template <typename chol_type>
443     sparse_chol<chol_type>&
operator =(const sparse_chol & a)444     sparse_chol<chol_type>::operator = (const sparse_chol& a)
445     {
446       if (this != &a)
447         {
448           if (--rep->count == 0)
449             delete rep;
450 
451           rep = a.rep;
452           rep->count++;
453         }
454 
455       return *this;
456     }
457 
458     template <typename chol_type>
459     chol_type
L(void) const460     sparse_chol<chol_type>::L (void) const
461     {
462 #if defined (HAVE_CHOLMOD)
463 
464       cholmod_sparse *m = rep->L ();
465 
466       octave_idx_type nc = m->ncol;
467       octave_idx_type nnz = m->nzmax;
468 
469       chol_type ret (m->nrow, nc, nnz);
470 
471       for (octave_idx_type j = 0; j < nc+1; j++)
472         ret.xcidx (j) = static_cast<octave_idx_type *>(m->p)[j];
473 
474       for (octave_idx_type i = 0; i < nnz; i++)
475         {
476           ret.xridx (i) = static_cast<octave_idx_type *>(m->i)[i];
477           ret.xdata (i) = static_cast<chol_elt *>(m->x)[i];
478         }
479 
480       return ret;
481 
482 #else
483 
484       return chol_type ();
485 
486 #endif
487     }
488 
489     template <typename chol_type>
490     octave_idx_type
P(void) const491     sparse_chol<chol_type>::P (void) const
492     {
493       return rep->P ();
494     }
495 
496     template <typename chol_type>
497     RowVector
perm(void) const498     sparse_chol<chol_type>::perm (void) const
499     {
500       return rep->perm ();
501     }
502 
503     template <typename chol_type>
504     SparseMatrix
Q(void) const505     sparse_chol<chol_type>::Q (void) const
506     {
507       return rep->Q ();
508     }
509 
510     template <typename chol_type>
511     bool
is_positive_definite(void) const512     sparse_chol<chol_type>::is_positive_definite (void) const
513     {
514       return rep->is_positive_definite ();
515     }
516 
517     template <typename chol_type>
518     double
rcond(void) const519     sparse_chol<chol_type>::rcond (void) const
520     {
521       return rep->rcond ();
522     }
523 
524     template <typename chol_type>
525     chol_type
inverse(void) const526     sparse_chol<chol_type>::inverse (void) const
527     {
528       chol_type retval;
529 
530 #if defined (HAVE_CHOLMOD)
531 
532       cholmod_sparse *m = rep->L ();
533       octave_idx_type n = m->ncol;
534       RowVector perms = rep->perm ();
535       double rcond2;
536       octave_idx_type info;
537       MatrixType mattype (MatrixType::Upper);
538       chol_type linv = L ().hermitian ().inverse (mattype, info, rcond2, 1, 0);
539 
540       if (perms.numel () == n)
541         {
542           SparseMatrix Qc = Q ();
543 
544           retval = Qc * linv * linv.hermitian () * Qc.transpose ();
545         }
546       else
547         retval = linv * linv.hermitian ();
548 
549 #endif
550 
551       return retval;
552     }
553 
554     template <typename chol_type>
555     chol_type
chol2inv(const chol_type & r)556     chol2inv (const chol_type& r)
557     {
558       octave_idx_type r_nr = r.rows ();
559       octave_idx_type r_nc = r.cols ();
560       chol_type retval;
561 
562       if (r_nr != r_nc)
563         (*current_liboctave_error_handler) ("U must be a square matrix");
564 
565       MatrixType mattype (r);
566       int typ = mattype.type (false);
567       double rcond;
568       octave_idx_type info;
569       chol_type rtra, multip;
570 
571       if (typ == MatrixType::Upper)
572         {
573           rtra = r.transpose ();
574           multip = (rtra*r);
575         }
576       else if (typ == MatrixType::Lower)
577         {
578           rtra = r.transpose ();
579           multip = (r*rtra);
580         }
581       else
582         (*current_liboctave_error_handler) ("U must be a triangular matrix");
583 
584       MatrixType mattypenew (multip);
585       retval = multip.inverse (mattypenew, info, rcond, true, false);
586       return retval;
587     }
588 
589     // SparseComplexMatrix specialization (the value for the NATURAL
590     // parameter in the sparse_chol<T>::sparse_chol_rep constructor is
591     // different from the default).
592 
593     template <>
sparse_chol(const SparseComplexMatrix & a,octave_idx_type & info)594     sparse_chol<SparseComplexMatrix>::sparse_chol (const SparseComplexMatrix& a,
595                                                    octave_idx_type& info)
596       : rep (new sparse_chol<SparseComplexMatrix>::sparse_chol_rep (a, info,
597                                                                     true,
598                                                                     false))
599     { }
600 
601     // Instantiations we need.
602 
603     template class sparse_chol<SparseMatrix>;
604 
605     template class sparse_chol<SparseComplexMatrix>;
606 
607     template SparseMatrix
608     chol2inv<SparseMatrix> (const SparseMatrix& r);
609 
610     template SparseComplexMatrix
611     chol2inv<SparseComplexMatrix> (const SparseComplexMatrix& r);
612   }
613 }
614