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
SparsedbleLU.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2004-2013 David Bateman
4 Copyright (C) 1998-2004 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 <vector>
29 
30 #include "lo-error.h"
31 #include "oct-locbuf.h"
32 
33 #include "SparsedbleLU.h"
34 #include "oct-spparms.h"
35 
36 // Instantiate the base LU class for the types we need.
37 
38 #include "sparse-base-lu.h"
39 #include "sparse-base-lu.cc"
40 
41 template class sparse_base_lu <SparseMatrix, double, SparseMatrix, double>;
42 
43 #include "oct-sparse.h"
44 
45 SparseLU::SparseLU (const SparseMatrix& a, const Matrix& piv_thres, bool scale)
46 {
47 #ifdef HAVE_UMFPACK
48  octave_idx_type nr = a.rows ();
49  octave_idx_type nc = a.cols ();
50 
51  // Setup the control parameters
52  Matrix Control (UMFPACK_CONTROL, 1);
53  double *control = Control.fortran_vec ();
54  UMFPACK_DNAME (defaults) (control);
55 
56  double tmp = octave_sparse_params::get_key ("spumoni");
57  if (!xisnan (tmp))
58  Control (UMFPACK_PRL) = tmp;
59 
60  if (piv_thres.nelem () == 2)
61  {
62  tmp = (piv_thres (0) > 1. ? 1. : piv_thres (0));
63  if (!xisnan (tmp))
64  Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
65  tmp = (piv_thres (1) > 1. ? 1. : piv_thres (1));
66  if (!xisnan (tmp))
67  Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
68  }
69  else
70  {
71  tmp = octave_sparse_params::get_key ("piv_tol");
72  if (!xisnan (tmp))
73  Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
74 
75  tmp = octave_sparse_params::get_key ("sym_tol");
76  if (!xisnan (tmp))
77  Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
78  }
79 
80  // Set whether we are allowed to modify Q or not
81  tmp = octave_sparse_params::get_key ("autoamd");
82  if (!xisnan (tmp))
83  Control (UMFPACK_FIXQ) = tmp;
84 
85  if (scale)
86  Control (UMFPACK_SCALE) = UMFPACK_SCALE_SUM;
87  else
88  Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
89 
90  UMFPACK_DNAME (report_control) (control);
91 
92  const octave_idx_type *Ap = a.cidx ();
93  const octave_idx_type *Ai = a.ridx ();
94  const double *Ax = a.data ();
95 
96  UMFPACK_DNAME (report_matrix) (nr, nc, Ap, Ai, Ax, 1, control);
97 
98  void *Symbolic;
99  Matrix Info (1, UMFPACK_INFO);
100  double *info = Info.fortran_vec ();
101  int status = UMFPACK_DNAME (qsymbolic) (nr, nc, Ap, Ai, Ax, 0,
102  &Symbolic, control, info);
103 
104  if (status < 0)
105  {
106  (*current_liboctave_error_handler)
107  ("SparseLU::SparseLU symbolic factorization failed");
108 
109  UMFPACK_DNAME (report_status) (control, status);
110  UMFPACK_DNAME (report_info) (control, info);
111 
112  UMFPACK_DNAME (free_symbolic) (&Symbolic) ;
113  }
114  else
115  {
116  UMFPACK_DNAME (report_symbolic) (Symbolic, control);
117 
118  void *Numeric;
119  status = UMFPACK_DNAME (numeric) (Ap, Ai, Ax, Symbolic,
120  &Numeric, control, info) ;
121  UMFPACK_DNAME (free_symbolic) (&Symbolic) ;
122 
123  cond = Info (UMFPACK_RCOND);
124 
125  if (status < 0)
126  {
127  (*current_liboctave_error_handler)
128  ("SparseLU::SparseLU numeric factorization failed");
129 
130  UMFPACK_DNAME (report_status) (control, status);
131  UMFPACK_DNAME (report_info) (control, info);
132 
133  UMFPACK_DNAME (free_numeric) (&Numeric);
134  }
135  else
136  {
137  UMFPACK_DNAME (report_numeric) (Numeric, control);
138 
139  octave_idx_type lnz, unz, ignore1, ignore2, ignore3;
140  status = UMFPACK_DNAME (get_lunz) (&lnz, &unz, &ignore1,
141  &ignore2, &ignore3, Numeric) ;
142 
143  if (status < 0)
144  {
145  (*current_liboctave_error_handler)
146  ("SparseLU::SparseLU extracting LU factors failed");
147 
148  UMFPACK_DNAME (report_status) (control, status);
149  UMFPACK_DNAME (report_info) (control, info);
150 
151  UMFPACK_DNAME (free_numeric) (&Numeric);
152  }
153  else
154  {
155  octave_idx_type n_inner = (nr < nc ? nr : nc);
156 
157  if (lnz < 1)
158  Lfact = SparseMatrix (n_inner, nr,
159  static_cast<octave_idx_type> (1));
160  else
161  Lfact = SparseMatrix (n_inner, nr, lnz);
162 
163  octave_idx_type *Ltp = Lfact.cidx ();
164  octave_idx_type *Ltj = Lfact.ridx ();
165  double *Ltx = Lfact.data ();
166 
167  if (unz < 1)
168  Ufact = SparseMatrix (n_inner, nc,
169  static_cast<octave_idx_type> (1));
170  else
171  Ufact = SparseMatrix (n_inner, nc, unz);
172 
173  octave_idx_type *Up = Ufact.cidx ();
174  octave_idx_type *Uj = Ufact.ridx ();
175  double *Ux = Ufact.data ();
176 
177  Rfact = SparseMatrix (nr, nr, nr);
178  for (octave_idx_type i = 0; i < nr; i++)
179  {
180  Rfact.xridx (i) = i;
181  Rfact.xcidx (i) = i;
182  }
183  Rfact.xcidx (nr) = nr;
184  double *Rx = Rfact.data ();
185 
186  P.resize (dim_vector (nr, 1));
187  octave_idx_type *p = P.fortran_vec ();
188 
189  Q.resize (dim_vector (nc, 1));
190  octave_idx_type *q = Q.fortran_vec ();
191 
192  octave_idx_type do_recip;
193  status = UMFPACK_DNAME (get_numeric) (Ltp, Ltj, Ltx,
194  Up, Uj, Ux, p, q, 0,
195  &do_recip, Rx,
196  Numeric) ;
197 
198  UMFPACK_DNAME (free_numeric) (&Numeric) ;
199 
200  if (status < 0)
201  {
202  (*current_liboctave_error_handler)
203  ("SparseLU::SparseLU extracting LU factors failed");
204 
205  UMFPACK_DNAME (report_status) (control, status);
206  }
207  else
208  {
209  Lfact = Lfact.transpose ();
210 
211  if (do_recip)
212  for (octave_idx_type i = 0; i < nr; i++)
213  Rx[i] = 1.0 / Rx[i];
214 
215  UMFPACK_DNAME (report_matrix) (nr, n_inner,
216  Lfact.cidx (), Lfact.ridx (),
217  Lfact.data (), 1, control);
218  UMFPACK_DNAME (report_matrix) (n_inner, nc,
219  Ufact.cidx (), Ufact.ridx (),
220  Ufact.data (), 1, control);
221  UMFPACK_DNAME (report_perm) (nr, p, control);
222  UMFPACK_DNAME (report_perm) (nc, q, control);
223  }
224 
225  UMFPACK_DNAME (report_info) (control, info);
226  }
227  }
228  }
229 #else
230  (*current_liboctave_error_handler) ("UMFPACK not installed");
231 #endif
232 }
233 
234 SparseLU::SparseLU (const SparseMatrix& a, const ColumnVector& Qinit,
235  const Matrix& piv_thres, bool scale, bool FixedQ,
236  double droptol, bool milu, bool udiag)
237 {
238 #ifdef HAVE_UMFPACK
239  if (milu)
240  (*current_liboctave_error_handler)
241  ("Modified incomplete LU not implemented");
242  else
243  {
244  octave_idx_type nr = a.rows ();
245  octave_idx_type nc = a.cols ();
246 
247  // Setup the control parameters
248  Matrix Control (UMFPACK_CONTROL, 1);
249  double *control = Control.fortran_vec ();
250  UMFPACK_DNAME (defaults) (control);
251 
252  double tmp = octave_sparse_params::get_key ("spumoni");
253  if (!xisnan (tmp))
254  Control (UMFPACK_PRL) = tmp;
255 
256  if (piv_thres.nelem () == 2)
257  {
258  tmp = (piv_thres (0) > 1. ? 1. : piv_thres (0));
259  if (!xisnan (tmp))
260  Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
261  tmp = (piv_thres (1) > 1. ? 1. : piv_thres (1));
262  if (!xisnan (tmp))
263  Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
264  }
265  else
266  {
267  tmp = octave_sparse_params::get_key ("piv_tol");
268  if (!xisnan (tmp))
269  Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
270 
271  tmp = octave_sparse_params::get_key ("sym_tol");
272  if (!xisnan (tmp))
273  Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
274  }
275 
276  if (droptol >= 0.)
277  Control (UMFPACK_DROPTOL) = droptol;
278 
279 
280  // Set whether we are allowed to modify Q or not
281  if (FixedQ)
282  Control (UMFPACK_FIXQ) = 1.0;
283  else
284  {
285  tmp = octave_sparse_params::get_key ("autoamd");
286  if (!xisnan (tmp))
287  Control (UMFPACK_FIXQ) = tmp;
288  }
289 
290  if (scale)
291  Control (UMFPACK_SCALE) = UMFPACK_SCALE_SUM;
292  else
293  Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
294 
295  UMFPACK_DNAME (report_control) (control);
296 
297  const octave_idx_type *Ap = a.cidx ();
298  const octave_idx_type *Ai = a.ridx ();
299  const double *Ax = a.data ();
300 
301  UMFPACK_DNAME (report_matrix) (nr, nc, Ap, Ai, Ax, 1,
302  control);
303 
304  void *Symbolic;
305  Matrix Info (1, UMFPACK_INFO);
306  double *info = Info.fortran_vec ();
307  int status;
308 
309  // Null loop so that qinit is imediately deallocated when not needed
310  do
311  {
312  OCTAVE_LOCAL_BUFFER (octave_idx_type, qinit, nc);
313 
314  for (octave_idx_type i = 0; i < nc; i++)
315  qinit[i] = static_cast<octave_idx_type> (Qinit (i));
316 
317  status = UMFPACK_DNAME (qsymbolic) (nr, nc, Ap, Ai, Ax,
318  qinit, &Symbolic, control, info);
319  }
320  while (0);
321 
322  if (status < 0)
323  {
324  (*current_liboctave_error_handler)
325  ("SparseLU::SparseLU symbolic factorization failed");
326 
327  UMFPACK_DNAME (report_status) (control, status);
328  UMFPACK_DNAME (report_info) (control, info);
329 
330  UMFPACK_DNAME (free_symbolic) (&Symbolic) ;
331  }
332  else
333  {
334  UMFPACK_DNAME (report_symbolic) (Symbolic, control);
335 
336  void *Numeric;
337  status = UMFPACK_DNAME (numeric) (Ap, Ai, Ax, Symbolic,
338  &Numeric, control, info) ;
339  UMFPACK_DNAME (free_symbolic) (&Symbolic) ;
340 
341  cond = Info (UMFPACK_RCOND);
342 
343  if (status < 0)
344  {
345  (*current_liboctave_error_handler)
346  ("SparseLU::SparseLU numeric factorization failed");
347 
348  UMFPACK_DNAME (report_status) (control, status);
349  UMFPACK_DNAME (report_info) (control, info);
350 
351  UMFPACK_DNAME (free_numeric) (&Numeric);
352  }
353  else
354  {
355  UMFPACK_DNAME (report_numeric) (Numeric, control);
356 
357  octave_idx_type lnz, unz, ignore1, ignore2, ignore3;
358  status = UMFPACK_DNAME (get_lunz) (&lnz, &unz,
359  &ignore1, &ignore2, &ignore3,
360  Numeric) ;
361 
362  if (status < 0)
363  {
364  (*current_liboctave_error_handler)
365  ("SparseLU::SparseLU extracting LU factors failed");
366 
367  UMFPACK_DNAME (report_status) (control, status);
368  UMFPACK_DNAME (report_info) (control, info);
369 
370  UMFPACK_DNAME (free_numeric) (&Numeric);
371  }
372  else
373  {
374  octave_idx_type n_inner = (nr < nc ? nr : nc);
375 
376  if (lnz < 1)
377  Lfact = SparseMatrix (n_inner, nr,
378  static_cast<octave_idx_type> (1));
379  else
380  Lfact = SparseMatrix (n_inner, nr, lnz);
381 
382  octave_idx_type *Ltp = Lfact.cidx ();
383  octave_idx_type *Ltj = Lfact.ridx ();
384  double *Ltx = Lfact.data ();
385 
386  if (unz < 1)
387  Ufact = SparseMatrix (n_inner, nc,
388  static_cast<octave_idx_type> (1));
389  else
390  Ufact = SparseMatrix (n_inner, nc, unz);
391 
392  octave_idx_type *Up = Ufact.cidx ();
393  octave_idx_type *Uj = Ufact.ridx ();
394  double *Ux = Ufact.data ();
395 
396  Rfact = SparseMatrix (nr, nr, nr);
397  for (octave_idx_type i = 0; i < nr; i++)
398  {
399  Rfact.xridx (i) = i;
400  Rfact.xcidx (i) = i;
401  }
402  Rfact.xcidx (nr) = nr;
403  double *Rx = Rfact.data ();
404 
405  P.resize (dim_vector (nr, 1));
406  octave_idx_type *p = P.fortran_vec ();
407 
408  Q.resize (dim_vector (nc, 1));
409  octave_idx_type *q = Q.fortran_vec ();
410 
411  octave_idx_type do_recip;
412  status = UMFPACK_DNAME (get_numeric) (Ltp, Ltj,
413  Ltx, Up, Uj, Ux, p, q,
414  0, &do_recip,
415  Rx, Numeric) ;
416 
417  UMFPACK_DNAME (free_numeric) (&Numeric) ;
418 
419  if (status < 0)
420  {
421  (*current_liboctave_error_handler)
422  ("SparseLU::SparseLU extracting LU factors failed");
423 
424  UMFPACK_DNAME (report_status) (control, status);
425  }
426  else
427  {
428  Lfact = Lfact.transpose ();
429 
430  if (do_recip)
431  for (octave_idx_type i = 0; i < nr; i++)
432  Rx[i] = 1.0 / Rx[i];
433 
434  UMFPACK_DNAME (report_matrix) (nr, n_inner,
435  Lfact.cidx (),
436  Lfact.ridx (),
437  Lfact.data (),
438  1, control);
439  UMFPACK_DNAME (report_matrix) (n_inner, nc,
440  Ufact.cidx (),
441  Ufact.ridx (),
442  Ufact.data (),
443  1, control);
444  UMFPACK_DNAME (report_perm) (nr, p, control);
445  UMFPACK_DNAME (report_perm) (nc, q, control);
446  }
447 
448  UMFPACK_DNAME (report_info) (control, info);
449  }
450  }
451  }
452 
453  if (udiag)
454  (*current_liboctave_error_handler)
455  ("Option udiag of incomplete LU not implemented");
456  }
457 #else
458  (*current_liboctave_error_handler) ("UMFPACK not installed");
459 #endif
460 }
Generated on Mon Dec 30 2013 03:04:50 for GNU Octave by   doxygen 1.8.1.2