GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
sparse-base-chol.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2005-2013 David Bateman
4 Copyright (C) 1998-2005 Andy Adler
5 
6 This file is part of Octave.
7 
8 Octave is free software; you can redistribute it and/or modify it
9 under the terms of the GNU General Public License as published by the
10 Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12 
13 Octave is distributed in the hope that it will be useful, but WITHOUT
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
16 for more details.
17 
18 You should have received a copy of the GNU General Public License
19 along with Octave; see the file COPYING. If not, see
20 <http://www.gnu.org/licenses/>.
21 
22 */
23 
24 #ifdef HAVE_CONFIG_H
25 #include <config.h>
26 #endif
27 
28 #include "sparse-base-chol.h"
29 #include "sparse-util.h"
30 #include "lo-error.h"
31 #include "oct-sparse.h"
32 #include "oct-spparms.h"
33 #include "quit.h"
34 #include "MatrixType.h"
35 
36 #ifdef HAVE_CHOLMOD
37 // Can't use CHOLMOD_NAME(drop)(0.0, S, cm). It doesn't treat complex matrices
38 template <class chol_type, class chol_elt, class p_type>
39 void
41  (const cholmod_sparse* S)
42 {
43  chol_elt sik;
44  octave_idx_type *Sp, *Si;
45  chol_elt *Sx;
46  octave_idx_type pdest, k, ncol, p, pend;
47 
48  if (! S)
49  return;
50 
51  Sp = static_cast<octave_idx_type *>(S->p);
52  Si = static_cast<octave_idx_type *>(S->i);
53  Sx = static_cast<chol_elt *>(S->x);
54  pdest = 0;
55  ncol = S->ncol;
56 
57  for (k = 0; k < ncol; k++)
58  {
59  p = Sp[k];
60  pend = Sp[k+1];
61  Sp[k] = pdest;
62  for (; p < pend; p++)
63  {
64  sik = Sx[p];
65  if (CHOLMOD_IS_NONZERO (sik))
66  {
67  if (p != pdest)
68  {
69  Si[pdest] = Si[p];
70  Sx[pdest] = sik;
71  }
72  pdest++;
73  }
74  }
75  }
76  Sp[ncol] = pdest;
77 }
78 #endif
79 
80 template <class chol_type, class chol_elt, class p_type>
83  (const chol_type& a, bool natural, bool force)
84 {
85  volatile octave_idx_type info = 0;
86 
87 #ifdef HAVE_CHOLMOD
88  octave_idx_type a_nr = a.rows ();
89  octave_idx_type a_nc = a.cols ();
90 
91  if (a_nr != a_nc)
92  {
93  (*current_liboctave_error_handler)
94  ("SparseCHOL requires square matrix");
95  return -1;
96  }
97 
98  cholmod_common *cm = &Common;
99 
100  // Setup initial parameters
101  CHOLMOD_NAME(start) (cm);
102  cm->prefer_zomplex = false;
103 
104  double spu = octave_sparse_params::get_key ("spumoni");
105  if (spu == 0.)
106  {
107  cm->print = -1;
108  cm->print_function = 0;
109  }
110  else
111  {
112  cm->print = static_cast<int> (spu) + 2;
113  cm->print_function =&SparseCholPrint;
114  }
115 
116  cm->error_handler = &SparseCholError;
117  cm->complex_divide = CHOLMOD_NAME(divcomplex);
118  cm->hypotenuse = CHOLMOD_NAME(hypot);
119 
120  cm->final_asis = false;
121  cm->final_super = false;
122  cm->final_ll = true;
123  cm->final_pack = true;
124  cm->final_monotonic = true;
125  cm->final_resymbol = false;
126 
127  cholmod_sparse A;
128  cholmod_sparse *ac = &A;
129  double dummy;
130  ac->nrow = a_nr;
131  ac->ncol = a_nc;
132 
133  ac->p = a.cidx ();
134  ac->i = a.ridx ();
135  ac->nzmax = a.nnz ();
136  ac->packed = true;
137  ac->sorted = true;
138  ac->nz = 0;
139 #ifdef USE_64_BIT_IDX_T
140  ac->itype = CHOLMOD_LONG;
141 #else
142  ac->itype = CHOLMOD_INT;
143 #endif
144  ac->dtype = CHOLMOD_DOUBLE;
145  ac->stype = 1;
146 #ifdef OCTAVE_CHOLMOD_TYPE
147  ac->xtype = OCTAVE_CHOLMOD_TYPE;
148 #else
149  ac->xtype = CHOLMOD_REAL;
150 #endif
151 
152  if (a_nr < 1)
153  ac->x = &dummy;
154  else
155  ac->x = a.data ();
156 
157  // use natural ordering if no q output parameter
158  if (natural)
159  {
160  cm->nmethods = 1 ;
161  cm->method[0].ordering = CHOLMOD_NATURAL ;
162  cm->postorder = false ;
163  }
164 
165  cholmod_factor *Lfactor;
167  Lfactor = CHOLMOD_NAME(analyze) (ac, cm);
168  CHOLMOD_NAME(factorize) (ac, Lfactor, cm);
170 
171  is_pd = cm->status == CHOLMOD_OK;
172  info = (is_pd ? 0 : cm->status);
173 
174  if (is_pd || force)
175  {
177  cond = CHOLMOD_NAME(rcond) (Lfactor, cm);
179 
180  minor_p = Lfactor->minor;
181 
183  Lsparse = CHOLMOD_NAME(factor_to_sparse) (Lfactor, cm);
185 
186  if (minor_p > 0 && minor_p < a_nr)
187  {
188  size_t n1 = a_nr + 1;
189  Lsparse->p = CHOLMOD_NAME(realloc) (minor_p+1,
190  sizeof(octave_idx_type),
191  Lsparse->p, &n1, cm);
193  CHOLMOD_NAME(reallocate_sparse)
194  (static_cast<octave_idx_type *>(Lsparse->p)[minor_p], Lsparse, cm);
196  Lsparse->ncol = minor_p;
197  }
198 
199  drop_zeros (Lsparse);
200 
201  if (! natural)
202  {
203  perms.resize (a_nr);
204  for (octave_idx_type i = 0; i < a_nr; i++)
205  perms(i) = static_cast<octave_idx_type *>(Lfactor->Perm)[i];
206  }
207 
208  static char tmp[] = " ";
209 
211  CHOLMOD_NAME(free_factor) (&Lfactor, cm);
212  CHOLMOD_NAME(finish) (cm);
213  CHOLMOD_NAME(print_common) (tmp, cm);
215  }
216 #else
217  (*current_liboctave_error_handler)
218  ("Missing CHOLMOD. Sparse cholesky factorization disabled");
219 #endif
220  return info;
221 }
222 
223 template <class chol_type, class chol_elt, class p_type>
224 chol_type
226 {
227 #ifdef HAVE_CHOLMOD
228  cholmod_sparse *m = rep->L ();
229  octave_idx_type nc = m->ncol;
230  octave_idx_type nnz = m->nzmax;
231  chol_type ret (m->nrow, nc, nnz);
232  for (octave_idx_type j = 0; j < nc+1; j++)
233  ret.xcidx (j) = static_cast<octave_idx_type *>(m->p)[j];
234  for (octave_idx_type i = 0; i < nnz; i++)
235  {
236  ret.xridx (i) = static_cast<octave_idx_type *>(m->i)[i];
237  ret.xdata (i) = static_cast<chol_elt *>(m->x)[i];
238  }
239  return ret;
240 #else
241  return chol_type ();
242 #endif
243 }
244 
245 template <class chol_type, class chol_elt, class p_type>
246 p_type
249 {
250 #ifdef HAVE_CHOLMOD
251  octave_idx_type n = Lsparse->nrow;
252  p_type p (n, n, n);
253 
254  for (octave_idx_type i = 0; i < n; i++)
255  {
256  p.xcidx (i) = i;
257  p.xridx (i) = static_cast<octave_idx_type>(perms (i));
258  p.xdata (i) = 1;
259  }
260  p.xcidx (n) = n;
261 
262  return p;
263 #else
264  return p_type ();
265 #endif
266 }
267 
268 template <class chol_type, class chol_elt, class p_type>
269 chol_type
271 {
272  chol_type retval;
273 #ifdef HAVE_CHOLMOD
274  cholmod_sparse *m = rep->L ();
275  octave_idx_type n = m->ncol;
276  ColumnVector perms = rep->perm ();
277  chol_type ret;
278  double rcond2;
279  octave_idx_type info;
280  MatrixType mattype (MatrixType::Upper);
281  chol_type linv = L ().hermitian ().inverse (mattype, info, rcond2, 1, 0);
282 
283  if (perms.length () == n)
284  {
285  p_type Qc = Q ();
286  retval = Qc * linv * linv.hermitian () * Qc.transpose ();
287  }
288  else
289  retval = linv * linv.hermitian ();
290 #endif
291  return retval;
292 }