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