GNU Octave  4.0.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
luinc.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2005-2015 David Bateman
4 
5 This file is part of Octave.
6 
7 Octave is free software; you can redistribute it and/or modify it
8 under the terms of the GNU General Public License as published by the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11 
12 Octave is distributed in the hope that it will be useful, but WITHOUT
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 for more details.
16 
17 You should have received a copy of the GNU General Public License
18 along with Octave; see the file COPYING. If not, see
19 <http://www.gnu.org/licenses/>.
20 
21 */
22 
23 #ifdef HAVE_CONFIG_H
24 #include <config.h>
25 #endif
26 
27 #include "defun.h"
28 #include "error.h"
29 #include "gripes.h"
30 #include "oct-obj.h"
31 #include "utils.h"
32 #include "oct-map.h"
33 
34 #include "MatrixType.h"
35 #include "SparseCmplxLU.h"
36 #include "SparsedbleLU.h"
37 #include "ov-re-sparse.h"
38 #include "ov-cx-sparse.h"
39 
40 // FIXME: Deprecated in 4.0 and should be removed in 4.4.
41 DEFUN (__luinc__, args, nargout,
42  "-*- texinfo -*-\n\
43 @deftypefn {Built-in Function} {[@var{L}, @var{U}, @var{P}, @var{Q}] =} luinc (@var{A}, '0')\n\
44 @deftypefnx {Built-in Function} {[@var{L}, @var{U}, @var{P}, @var{Q}] =} luinc (@var{A}, @var{droptol})\n\
45 @deftypefnx {Built-in Function} {[@var{L}, @var{U}, @var{P}, @var{Q}] =} luinc (@var{A}, @var{opts})\n\
46 @cindex LU decomposition\n\
47 Produce the incomplete LU@tie{}factorization of the sparse matrix @var{A}.\n\
48 \n\
49 Two types of incomplete factorization are possible, and the type\n\
50 is determined by the second argument to @code{luinc}.\n\
51 \n\
52 Called with a second argument of @qcode{'0'}, the zero-level incomplete\n\
53 LU@tie{}factorization is produced. This creates a factorization of @var{A}\n\
54 where the position of the nonzero arguments correspond to the same\n\
55 positions as in the matrix @var{A}.\n\
56 \n\
57 Alternatively, the fill-in of the incomplete LU@tie{}factorization can\n\
58 be controlled through the variable @var{droptol} or the structure\n\
59 @var{opts}. The @sc{umfpack} multifrontal factorization code by Tim A.\n\
60 Davis is used for the incomplete LU@tie{}factorization, (availability\n\
61 @url{http://www.cise.ufl.edu/research/sparse/umfpack/})\n\
62 \n\
63 @var{droptol} determines the values below which the values in the\n\
64 LU@tie{} factorization are dropped and replaced by zero. It must be a\n\
65 positive scalar, and any values in the factorization whose absolute value\n\
66 are less than this value are dropped, expect if leaving them increase the\n\
67 sparsity of the matrix. Setting @var{droptol} to zero results in a complete\n\
68 LU@tie{}factorization which is the default.\n\
69 \n\
70 @var{opts} is a structure containing one or more of the fields\n\
71 \n\
72 @table @code\n\
73 @item droptol\n\
74 The drop tolerance as above. If @var{opts} only contains @code{droptol}\n\
75 then this is equivalent to using the variable @var{droptol}.\n\
76 \n\
77 @item milu\n\
78 A logical variable flagging whether to use the modified incomplete\n\
79 LU@tie{} factorization. In the case that @code{milu} is true, the dropped\n\
80 values are subtracted from the diagonal of the matrix @var{U} of the\n\
81 factorization. The default is @code{false}.\n\
82 \n\
83 @item udiag\n\
84 A logical variable that flags whether zero elements on the diagonal of\n\
85 @var{U} should be replaced with @var{droptol} to attempt to avoid singular\n\
86 factors. The default is @code{false}.\n\
87 \n\
88 @item thresh\n\
89 Defines the pivot threshold in the interval [0,1]. Values outside that\n\
90 range are ignored.\n\
91 @end table\n\
92 \n\
93 All other fields in @var{opts} are ignored. The outputs from @code{luinc}\n\
94 are the same as for @code{lu}.\n\
95 \n\
96 Given the string argument @qcode{\"vector\"}, @code{luinc} returns the\n\
97 values of @var{p} @var{q} as vector values.\n\
98 @seealso{sparse, lu, ilu, ichol}\n\
99 @end deftypefn")
100 {
101  int nargin = args.length ();
102  octave_value_list retval;
103 
104  if (nargin == 0)
105  print_usage ();
106  else if (nargin < 2 || nargin > 3)
107  error ("luinc: incorrect number of arguments");
108  else
109  {
110  bool zero_level = false;
111  bool milu = false;
112  bool udiag = false;
113  Matrix thresh;
114  double droptol = -1.;
115  bool vecout = false;
116 
117  if (args(1).is_string ())
118  {
119  if (args(1).string_value () == "0")
120  zero_level = true;
121  else
122  error ("luinc: unrecognized string argument");
123  }
124  else if (args(1).is_map ())
125  {
126  octave_scalar_map map = args(1).scalar_map_value ();
127 
128  if (! error_state)
129  {
130  octave_value tmp;
131 
132  tmp = map.getfield ("droptol");
133  if (tmp.is_defined ())
134  droptol = tmp.double_value ();
135 
136  tmp = map.getfield ("milu");
137  if (tmp.is_defined ())
138  {
139  double val = tmp.double_value ();
140 
141  milu = (val == 0. ? false : true);
142  }
143 
144  tmp = map.getfield ("udiag");
145  if (tmp.is_defined ())
146  {
147  double val = tmp.double_value ();
148 
149  udiag = (val == 0. ? false : true);
150  }
151 
152  tmp = map.getfield ("thresh");
153  if (tmp.is_defined ())
154  {
155  thresh = tmp.matrix_value ();
156 
157  if (thresh.nelem () == 1)
158  {
159  thresh.resize (1,2);
160  thresh(1) = thresh(0);
161  }
162  else if (thresh.nelem () != 2)
163  {
164  error ("luinc: expecting 2-element vector for thresh");
165  return retval;
166  }
167  }
168  }
169  else
170  {
171  error ("luinc: OPTS must be a scalar structure");
172  return retval;
173  }
174  }
175  else
176  droptol = args(1).double_value ();
177 
178  if (nargin == 3)
179  {
180  std::string tmp = args(2).string_value ();
181 
182  if (! error_state)
183  {
184  if (tmp.compare ("vector") == 0)
185  vecout = true;
186  else
187  error ("luinc: unrecognized string argument");
188  }
189  }
190 
191  // FIXME: Add code for zero-level factorization
192  if (zero_level)
193  error ("luinc: zero-level factorization not implemented");
194 
195  if (!error_state)
196  {
197  if (args(0).type_name () == "sparse matrix")
198  {
199  SparseMatrix sm = args(0).sparse_matrix_value ();
200  octave_idx_type sm_nr = sm.rows ();
201  octave_idx_type sm_nc = sm.cols ();
202  ColumnVector Qinit (sm_nc);
203 
204  for (octave_idx_type i = 0; i < sm_nc; i++)
205  Qinit (i) = i;
206 
207  if (! error_state)
208  {
209  switch (nargout)
210  {
211  case 0:
212  case 1:
213  case 2:
214  {
215  SparseLU fact (sm, Qinit, thresh, false, true, droptol,
216  milu, udiag);
217 
218  if (! error_state)
219  {
220  SparseMatrix P = fact.Pr ();
221  SparseMatrix L = P.transpose () * fact.L ();
222  retval(1)
223  = octave_value (fact.U (),
225  retval(0)
228  sm_nr, fact.row_perm ()));
229  }
230  }
231  break;
232 
233  case 3:
234  {
235  SparseLU fact (sm, Qinit, thresh, false, true, droptol,
236  milu, udiag);
237 
238  if (! error_state)
239  {
240  if (vecout)
241  retval(2) = fact.Pr_vec ();
242  else
243  retval(2) = fact.Pr_mat ();
244  retval(1)
245  = octave_value (fact.U (),
247  retval(0)
248  = octave_value (fact.L (),
250  }
251  }
252  break;
253 
254  case 4:
255  default:
256  {
257  SparseLU fact (sm, Qinit, thresh, false, false, droptol,
258  milu, udiag);
259 
260  if (! error_state)
261  {
262  if (vecout)
263  {
264  retval(3) = fact.Pc_vec ();
265  retval(2) = fact.Pr_vec ();
266  }
267  else
268  {
269  retval(3) = fact.Pc_mat ();
270  retval(2) = fact.Pr_mat ();
271  }
272  retval(1)
273  = octave_value (fact.U (),
275  retval(0)
276  = octave_value (fact.L (),
278  }
279  }
280  break;
281  }
282  }
283  }
284  else if (args(0).type_name () == "sparse complex matrix")
285  {
287  args(0).sparse_complex_matrix_value ();
288  octave_idx_type sm_nr = sm.rows ();
289  octave_idx_type sm_nc = sm.cols ();
290  ColumnVector Qinit (sm_nc);
291 
292  for (octave_idx_type i = 0; i < sm_nc; i++)
293  Qinit (i) = i;
294 
295  if (! error_state)
296  {
297  switch (nargout)
298  {
299  case 0:
300  case 1:
301  case 2:
302  {
303  SparseComplexLU fact (sm, Qinit, thresh, false, true,
304  droptol, milu, udiag);
305 
306 
307  if (! error_state)
308  {
309  SparseMatrix P = fact.Pr ();
310  SparseComplexMatrix L = P.transpose () * fact.L ();
311  retval(1)
312  = octave_value (fact.U (),
314  retval(0)
317  sm_nr, fact.row_perm ()));
318  }
319  }
320  break;
321 
322  case 3:
323  {
324  SparseComplexLU fact (sm, Qinit, thresh, false, true,
325  droptol, milu, udiag);
326 
327  if (! error_state)
328  {
329  if (vecout)
330  retval(2) = fact.Pr_vec ();
331  else
332  retval(2) = fact.Pr_mat ();
333  retval(1)
334  = octave_value (fact.U (),
336  retval(0)
337  = octave_value (fact.L (),
339  }
340  }
341  break;
342 
343  case 4:
344  default:
345  {
346  SparseComplexLU fact (sm, Qinit, thresh, false, false,
347  droptol, milu, udiag);
348 
349  if (! error_state)
350  {
351  if (vecout)
352  {
353  retval(3) = fact.Pc_vec ();
354  retval(2) = fact.Pr_vec ();
355  }
356  else
357  {
358  retval(3) = fact.Pc_mat ();
359  retval(2) = fact.Pr_mat ();
360  }
361  retval(1)
362  = octave_value (fact.U (),
364  retval(0)
365  = octave_value (fact.L (),
367  }
368  }
369  break;
370  }
371  }
372  }
373  else
374  error ("luinc: matrix A must be sparse");
375  }
376  }
377 
378  return retval;
379 }
380 
381 /*
382 %!testif HAVE_UMFPACK
383 %! a = sparse ([1,2,0,0;0,1,2,0;1e-14,0,3,0;0,0,0,1]);
384 %! [l,u] = luinc (a, 1e-10);
385 %! assert (l*u, sparse ([1,2,0,0;0,1,2,0;0,0,3,0;0,0,0,1]), 1e-10);
386 %! opts.droptol = 1e-10;
387 %! [l,u] = luinc (a, opts);
388 %! assert (l*u, sparse ([1,2,0,0;0,1,2,0;0,0,3,0;0,0,0,1]), 1e-10);
389 
390 %!testif HAVE_UMFPACK
391 %! a = sparse ([1i,2,0,0;0,1,2,0;1e-14,0,3,0;0,0,0,1]);
392 %! [l,u] = luinc (a, 1e-10);
393 %! assert (l*u, sparse ([1i,2,0,0;0,1,2,0;0,0,3,0;0,0,0,1]), 1e-10);
394 %! opts.droptol = 1e-10;
395 %! [l,u] = luinc (a, opts);
396 %! assert (l*u, sparse ([1i,2,0,0;0,1,2,0;0,0,3,0;0,0,0,1]), 1e-10);
397 */
octave_idx_type cols(void) const
Definition: Sparse.h:264
octave_idx_type rows(void) const
Definition: Sparse.h:263
lu_type U(void) const
void resize(octave_idx_type nr, octave_idx_type nc, double rfv=0)
Definition: dMatrix.h:130
SparseMatrix transpose(void) const
Definition: dSparse.h:133
const octave_idx_type * row_perm(void) const
OCTINTERP_API void print_usage(void)
Definition: defun.cc:51
bool is_defined(void) const
Definition: ov.h:520
#define DEFUN(name, args_name, nargout_name, doc)
Definition: defun.h:44
void error(const char *fmt,...)
Definition: error.cc:476
PermMatrix Pr_mat(void) const
octave_idx_type nelem(void) const
Number of elements in the array.
Definition: Array.h:271
p_type Pr(void) const
ColumnVector Pr_vec(void) const
lu_type L(void) const
int error_state
Definition: error.cc:101
Definition: dMatrix.h:35
Matrix matrix_value(bool frc_str_conv=false) const
Definition: ov.h:773
PermMatrix Pc_mat(void) const
ColumnVector Pc_vec(void) const
octave_value getfield(const std::string &key) const
Definition: oct-map.cc:164
double double_value(bool frc_str_conv=false) const
Definition: ov.h:759
return octave_value(v1.char_array_value().concat(v2.char_array_value(), ra_idx),((a1.is_sq_string()||a2.is_sq_string())? '\'': '"'))