Mercurial > octave-nkf

comparison src/DLD-FUNCTIONS/splu.cc @ 5164:57077d0ddc8e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2005-02-25 19:55:24 by jwe]
author jwe
date Fri, 25 Feb 2005 19:55:28 +0000
parents
children 5bdc3f24cd5f
comparison
equal deleted inserted replaced
5163:9f3299378193 5164:57077d0ddc8e
1 /*
2
3 Copyright (C) 2004 David Bateman
4 Copyright (C) 1998-2004 Andy Adler
5
6 Octave is free software; you can redistribute it and/or modify it
7 under the terms of the GNU General Public License as published by the
8 Free Software Foundation; either version 2, or (at your option) any
9 later version.
10
11 Octave is distributed in the hope that it will be useful, but WITHOUT
12 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
14 for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with this program; see the file COPYING. If not, write to the Free
18 Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
19
20 */
21
22 #ifdef HAVE_CONFIG_H
23 #include <config.h>
24 #endif
25
26 #include "defun-dld.h"
27 #include "error.h"
28 #include "gripes.h"
29 #include "oct-obj.h"
30 #include "utils.h"
31
32 #include "SparseCmplxLU.h"
33 #include "SparsedbleLU.h"
34 #include "ov-re-sparse.h"
35 #include "ov-cx-sparse.h"
36
37 // PKG_ADD: dispatch ("lu", "splu", "sparse matrix")
38 // PKG_ADD: dispatch ("lu", "splu", "sparse complex matrix")
39 // PKG_ADD: dispatch ("lu", "splu", "sparse bool matrix")
40 DEFUN_DLD (splu, args, nargout,
41 "-*- texinfo -*-\n\
42 @deftypefn {Loadable Function} {[@var{l}, @var{u}] =} splu (@var{a})\n\
43 @deftypefnx {Loadable Function} {[@var{l}, @var{u}, @var{P}] =} splu (@var{a})\n\
44 @deftypefnx {Loadable Function} {[@var{l}, @var{u}, @var{P}, @var{Q}] =} splu (@var{a})\n\
45 @deftypefnx {Loadable Function} {[@var{l}, @var{u}, @var{P}, @var{Q}] =} splu (@dots{}, @var{thres})\n\
46 @deftypefnx {Loadable Function} {[@var{l}, @var{u}, @var{P}] =} splu (@dots{}, @var{Q})\n\
47 @cindex LU decomposition\n\
48 Compute the LU decomposition of the sparse matrix @var{a}, using\n\
49 subroutines from UMFPACK. The result is returned in a permuted\n\
50 form, according to the optional return values @var{P} and @var{Q}.\n\
51 \n\
52 Called with two or three output arguments and a single input argument,\n\
53 @dfn{splu} is a replacement for @dfn{lu}, and therefore the sparsity\n\
54 preserving column permutations @var{Q} are not performed. Called with\n\
55 a fourth output argument, the sparsity preserving column transformation\n\
56 @var{Q} is returned, such that @code{@var{P} * @var{a} * @var{Q} =\n\
57 @var{l} * @var{u}}.\n\
58 \n\
59 An additional input argument @var{thres}, that defines the pivoting\n\
60 threshold can be given. Alternatively, the desired sparsity preserving\n\
61 column permutations @var{Q} can be passed. Note that @var{Q} is assumed\n\
62 to be fixed if three are fewer than four output arguments. Otherwise,\n\
63 the updated column permutations are returned as the fourth argument.\n\
64 \n\
65 With two output arguments, returns the permuted forms of the upper and\n\
66 lower triangular matrices, such that @code{@var{a} = @var{l} * @var{u}}.\n\
67 With two or three output arguments, if a user-defined @var{Q} is given,\n\
68 then @code{@var{u} * @var{Q}'} is returned. The matrix is not required to\n\
69 be square.\n\
70 @end deftypefn\n\
71 @seealso{sparse, spinv, colamd, symamd}")
72 {
73 octave_value_list retval;
74
75 int nargin = args.length ();
76
77 if (nargin < 1 || nargin > 3 || nargout > 4)
78 {
79 print_usage ("splu");
80 return retval;
81 }
82
83 octave_value arg = args(0);
84
85 int nr = arg.rows ();
86 int nc = arg.columns ();
87
88 int arg_is_empty = empty_arg ("splu", nr, nc);
89
90 if (arg_is_empty < 0)
91 return retval;
92 else if (arg_is_empty > 0)
93 return octave_value_list (3, SparseMatrix ());
94
95 ColumnVector Qinit;
96 bool have_Qinit = false;
97 double thres = -1.;
98
99 for (int k = 1; k < nargin; k++)
100 {
101 if (args(k).class_name () == "sparse")
102 {
103 SparseMatrix tmp = args (k).sparse_matrix_value ();
104
105 if (error_state)
106 {
107 error ("splu: Not a valid permutation/threshold");
108 return retval;
109 }
110
111 dim_vector dv = tmp.dims ();
112
113 if (dv(0) == 1 && dv(1) == 1)
114 thres = tmp (0);
115 else if (dv(0) == 1 || dv(1) == 1)
116 {
117 int nel = tmp.numel ();
118 Qinit.resize (nel);
119 for (int i = 0; i < nel; i++)
120 Qinit (i) = tmp (i) - 1;
121 have_Qinit = true;
122 }
123 else
124 {
125 int t_nc = tmp.cols ();
126
127 if (tmp.nnz () != t_nc)
128 error ("splu: Not a valid permutation matrix");
129 else
130 {
131 for (int i = 0; i < t_nc + 1; i++)
132 if (tmp.cidx(i) != i)
133 {
134 error ("splu: Not a valid permutation matrix");
135 break;
136 }
137 }
138
139 if (!error_state)
140 {
141 for (int i = 0; i < t_nc; i++)
142 if (tmp.data (i) != 1.)
143 {
144 error ("splu: Not a valid permutation matrix");
145 break;
146 }
147 else
148 Qinit (i) = tmp.ridx (i) - 1;
149 }
150
151 if (! error_state)
152 have_Qinit = true;
153 }
154 }
155 else
156 {
157 NDArray tmp = args(k).array_value ();
158
159 if (error_state)
160 return retval;
161
162 dim_vector dv = tmp.dims ();
163 if (dv.length () > 2)
164 {
165 error ("splu: second argument must be a vector/matrix or a scalar");
166 }
167 else if (dv(0) == 1 && dv(1) == 1)
168 thres = tmp (0);
169 else if (dv(0) == 1 || dv(1) == 1)
170 {
171 int nel = tmp.numel ();
172 Qinit.resize (nel);
173 for (int i = 0; i < nel; i++)
174 Qinit (i) = tmp (i) - 1;
175 have_Qinit = true;
176 }
177 else
178 {
179 SparseMatrix tmp2 (tmp);
180
181 int t_nc = tmp2.cols ();
182
183 if (tmp2.nnz () != t_nc)
184 error ("splu: Not a valid permutation matrix");
185 else
186 {
187 for (int i = 0; i < t_nc + 1; i++)
188 if (tmp2.cidx(i) != i)
189 {
190 error ("splu: Not a valid permutation matrix");
191 break;
192 }
193 }
194
195 if (!error_state)
196 {
197 for (int i = 0; i < t_nc; i++)
198 if (tmp2.data (i) != 1.)
199 {
200 error ("splu: Not a valid permutation matrix");
201 break;
202 }
203 else
204 Qinit (i) = tmp2.ridx (i) - 1;
205 }
206
207 if (! error_state)
208 have_Qinit = true;
209 }
210 }
211 }
212
213 if (error_state)
214 return retval;
215
216 if (arg.is_real_type ())
217 {
218 SparseMatrix m = arg.sparse_matrix_value ();
219
220 if (nargout < 4 && ! have_Qinit)
221 {
222 int m_nc = m.cols ();
223 Qinit.resize (m_nc);
224 for (int i = 0; i < m_nc; i++)
225 Qinit (i) = i;
226 }
227
228 if (! error_state)
229 {
230 switch (nargout)
231 {
232 case 0:
233 case 1:
234 case 2:
235 {
236 SparseLU fact (m, Qinit, thres, true);
237
238 SparseMatrix P = fact.Pr ();
239 SparseMatrix L = P.transpose () * fact.L ();
240 if (have_Qinit)
241 retval(1) = fact.U () * fact.Pc ().transpose ();
242 else
243 retval(1) = fact.U ();
244
245 retval(0) = L;
246 }
247 break;
248
249 case 3:
250 {
251 SparseLU fact (m, Qinit, thres, true);
252
253 retval(2) = fact.Pr ();
254 if (have_Qinit)
255 retval(1) = fact.U () * fact.Pc ().transpose ();
256 else
257 retval(1) = fact.U ();
258 retval(0) = fact.L ();
259 }
260 break;
261
262 case 4:
263 default:
264 {
265 if (have_Qinit)
266 {
267 SparseLU fact (m, Qinit, thres, false);
268
269 retval(3) = fact.Pc ();
270 retval(2) = fact.Pr ();
271 retval(1) = fact.U ();
272 retval(0) = fact.L ();
273 }
274 else
275 {
276 SparseLU fact (m, thres);
277
278 retval(3) = fact.Pc ();
279 retval(2) = fact.Pr ();
280 retval(1) = fact.U ();
281 retval(0) = fact.L ();
282 }
283 }
284 break;
285 }
286 }
287 }
288 else if (arg.is_complex_type ())
289 {
290 SparseComplexMatrix m = arg.sparse_complex_matrix_value ();
291
292 if (nargout < 4 && ! have_Qinit)
293 {
294 int m_nc = m.cols ();
295 Qinit.resize (m_nc);
296 for (int i = 0; i < m_nc; i++)
297 Qinit (i) = i;
298 }
299
300 if (! error_state)
301 {
302 switch (nargout)
303 {
304 case 0:
305 case 1:
306 case 2:
307 {
308 SparseComplexLU fact (m, Qinit, thres, true);
309
310 SparseMatrix P = fact.Pr ();
311 SparseComplexMatrix L = P.transpose () * fact.L ();
312
313 if (have_Qinit)
314 retval(1) = fact.U () * fact.Pc ().transpose ();
315 else
316 retval(1) = fact.U ();
317 retval(0) = L;
318 }
319 break;
320
321 case 3:
322 {
323 SparseComplexLU fact (m, Qinit, thres, true);
324
325 retval(2) = fact.Pr ();
326 if (have_Qinit)
327 retval(1) = fact.U () * fact.Pc ().transpose ();
328 else
329 retval(1) = fact.U ();
330 retval(0) = fact.L ();
331 }
332 break;
333
334 case 4:
335 default:
336 {
337 if (have_Qinit)
338 {
339 SparseComplexLU fact (m, Qinit, thres, false);
340
341 retval(3) = fact.Pc ();
342 retval(2) = fact.Pr ();
343 retval(1) = fact.U ();
344 retval(0) = fact.L ();
345 }
346 else
347 {
348 SparseComplexLU fact (m, thres);
349
350 retval(3) = fact.Pc ();
351 retval(2) = fact.Pr ();
352 retval(1) = fact.U ();
353 retval(0) = fact.L ();
354 }
355 }
356 break;
357 }
358 }
359 }
360 else
361 {
362 gripe_wrong_type_arg ("splu", arg);
363 }
364
365 return retval;
366 }
367
368 // PKG_ADD: dispatch ("inv", "spinv", "sparse matrix")
369 // PKG_ADD: dispatch ("inv", "spinv", "sparse complex matrix")
370 // PKG_ADD: dispatch ("inv", "spinv", "sparse bool matrix")
371 // PKG_ADD: dispatch ("inverse", "spinv", "sparse matrix")
372 // PKG_ADD: dispatch ("inverse", "spinv", "sparse complex matrix")
373 // PKG_ADD: dispatch ("inverse", "spinv", "sparse bool matrix")
374 DEFUN_DLD (spinv, args, nargout,
375 "-*- texinfo -*-\n\
376 @deftypefn {Loadable Function} {[@var{x}, @var{rcond}] = } spinv (@var{a}, @var{Q})\n\
377 @deftypefnx {Loadable Function} {[@var{x}, @var{rcond}, @var{Q}] = } spinv (@var{a}, @var{Q})\n\
378 Compute the inverse of the square matrix @var{a}. Return an estimate\n\
379 of the reciprocal condition number if requested, otherwise warn of an\n\
380 ill-conditioned matrix if the reciprocal condition number is small.\n\
381 \n\
382 An optional second input argument @var{Q} is the optional pre-ordering of\n\
383 the matrix, such that @code{@var{x} = inv (@var{a} (:, @var{Q}))}. @var{Q}\n\
384 can equally be a matrix, in which case @code{@var{x} = inv (@var{a} *\n\
385 @var{Q}))}.\n\
386 \n\
387 If a third output argument is given then the permuations to achieve a sparse\n\
388 inverse are returned. It is not required that the return column permutations\n\
389 @var{Q} and the same as the user supplied permutations\n\
390 @end deftypefn")
391 {
392 error ("spinv: not implemented yet");
393 return octave_value ();
394 }
395
396 /*
397 ;;; Local Variables: ***
398 ;;; mode: C++ ***
399 ;;; End: ***
400 */