Mercurial > octave-nkf

comparison liboctave/SparsedbleLU.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 dbeafbc0ff64
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 <vector>
27
28 #include "lo-error.h"
29
30 #include "SparsedbleLU.h"
31 #include "oct-spparms.h"
32
33 // Instantiate the base LU class for the types we need.
34
35 #include "sparse-base-lu.h"
36 #include "sparse-base-lu.cc"
37
38 template class sparse_base_lu <SparseMatrix, double, SparseMatrix, double>;
39
40 // Include the UMFPACK functions
41 extern "C" {
42 #include "umfpack.h"
43 }
44
45 SparseLU::SparseLU (const SparseMatrix& a, double piv_thres)
46 {
47 int nr = a.rows ();
48 int nc = a.cols ();
49
50 // Setup the control parameters
51 Matrix Control (UMFPACK_CONTROL, 1);
52 double *control = Control.fortran_vec ();
53 umfpack_di_defaults (control);
54
55 double tmp = Voctave_sparse_controls.get_key ("spumoni");
56 if (!xisnan (tmp))
57 Control (UMFPACK_PRL) = tmp;
58
59 if (piv_thres >= 0.)
60 {
61 piv_thres = (piv_thres > 1. ? 1. : piv_thres);
62 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = piv_thres;
63 Control (UMFPACK_PIVOT_TOLERANCE) = piv_thres;
64 }
65 else
66 {
67 tmp = Voctave_sparse_controls.get_key ("piv_tol");
68 if (!xisnan (tmp))
69 {
70 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
71 Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
72 }
73 }
74
75 // Set whether we are allowed to modify Q or not
76 tmp = Voctave_sparse_controls.get_key ("autoamd");
77 if (!xisnan (tmp))
78 Control (UMFPACK_FIXQ) = tmp;
79
80 // Turn-off UMFPACK scaling for LU
81 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
82
83 umfpack_di_report_control (control);
84
85 const int *Ap = a.cidx ();
86 const int *Ai = a.ridx ();
87 const double *Ax = a.data ();
88
89 umfpack_di_report_matrix (nr, nc, Ap, Ai, Ax, 1, control);
90
91 void *Symbolic;
92 Matrix Info (1, UMFPACK_INFO);
93 double *info = Info.fortran_vec ();
94 int status = umfpack_di_qsymbolic (nr, nc, Ap, Ai, Ax, NULL,
95 &Symbolic, control, info);
96
97 if (status < 0)
98 {
99 (*current_liboctave_error_handler)
100 ("SparseLU::SparseLU symbolic factorization failed");
101
102 umfpack_di_report_status (control, status);
103 umfpack_di_report_info (control, info);
104
105 umfpack_di_free_symbolic (&Symbolic) ;
106 }
107 else
108 {
109 umfpack_di_report_symbolic (Symbolic, control);
110
111 void *Numeric;
112 status = umfpack_di_numeric (Ap, Ai, Ax, Symbolic, &Numeric,
113 control, info) ;
114 umfpack_di_free_symbolic (&Symbolic) ;
115
116 cond = Info (UMFPACK_RCOND);
117
118 if (status < 0)
119 {
120 (*current_liboctave_error_handler)
121 ("SparseLU::SparseLU numeric factorization failed");
122
123 umfpack_di_report_status (control, status);
124 umfpack_di_report_info (control, info);
125
126 umfpack_di_free_numeric (&Numeric);
127 }
128 else
129 {
130 umfpack_di_report_numeric (Numeric, control);
131
132 int lnz, unz, ignore1, ignore2, ignore3;
133 status = umfpack_di_get_lunz (&lnz, &unz, &ignore1, &ignore2,
134 &ignore3, Numeric) ;
135
136 if (status < 0)
137 {
138 (*current_liboctave_error_handler)
139 ("SparseLU::SparseLU extracting LU factors failed");
140
141 umfpack_di_report_status (control, status);
142 umfpack_di_report_info (control, info);
143
144 umfpack_di_free_numeric (&Numeric);
145 }
146 else
147 {
148 int n_inner = (nr < nc ? nr : nc);
149
150 if (lnz < 1)
151 Lfact = SparseMatrix (n_inner, nr, 1);
152 else
153 Lfact = SparseMatrix (n_inner, nr, lnz);
154
155 int *Ltp = Lfact.cidx ();
156 int *Ltj = Lfact.ridx ();
157 double *Ltx = Lfact.data ();
158
159 if (unz < 1)
160 Ufact = SparseMatrix (n_inner, nc, 1);
161 else
162 Ufact = SparseMatrix (n_inner, nc, unz);
163
164 int *Up = Ufact.cidx ();
165 int *Uj = Ufact.ridx ();
166 double *Ux = Ufact.data ();
167
168 P.resize (nr);
169 int *p = P.fortran_vec ();
170
171 Q.resize (nc);
172 int *q = Q.fortran_vec ();
173
174 int do_recip;
175 status = umfpack_di_get_numeric (Ltp, Ltj, Ltx, Up, Uj,
176 Ux, p, q, (double *) NULL,
177 &do_recip, (double *) NULL,
178 Numeric) ;
179
180 umfpack_di_free_numeric (&Numeric) ;
181
182 if (status < 0 || do_recip)
183 {
184 (*current_liboctave_error_handler)
185 ("SparseLU::SparseLU extracting LU factors failed");
186
187 umfpack_di_report_status (control, status);
188 }
189 else
190 {
191 Lfact = Lfact.transpose ();
192
193 umfpack_di_report_matrix (nr, n_inner, Lfact.cidx (),
194 Lfact.ridx (), Lfact.data (),
195 1, control);
196 umfpack_di_report_matrix (n_inner, nc, Ufact.cidx (),
197 Ufact.ridx (), Ufact.data (),
198 1, control);
199 umfpack_di_report_perm (nr, p, control);
200 umfpack_di_report_perm (nc, q, control);
201 }
202
203 umfpack_di_report_info (control, info);
204 }
205 }
206 }
207 }
208
209 SparseLU::SparseLU (const SparseMatrix& a, const ColumnVector& Qinit,
210 double piv_thres, bool FixedQ)
211 {
212 int nr = a.rows ();
213 int nc = a.cols ();
214
215 // Setup the control parameters
216 Matrix Control (UMFPACK_CONTROL, 1);
217 double *control = Control.fortran_vec ();
218 umfpack_di_defaults (control);
219
220 double tmp = Voctave_sparse_controls.get_key ("spumoni");
221 if (!xisnan (tmp))
222 Control (UMFPACK_PRL) = tmp;
223 if (piv_thres >= 0.)
224 {
225 piv_thres = (piv_thres > 1. ? 1. : piv_thres);
226 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = piv_thres;
227 Control (UMFPACK_PIVOT_TOLERANCE) = piv_thres;
228 }
229 else
230 {
231 tmp = Voctave_sparse_controls.get_key ("piv_tol");
232 if (!xisnan (tmp))
233 {
234 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
235 Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
236 }
237 }
238
239 // Set whether we are allowed to modify Q or not
240 if (FixedQ)
241 Control (UMFPACK_FIXQ) = 1.0;
242 else
243 {
244 tmp = Voctave_sparse_controls.get_key ("autoamd");
245 if (!xisnan (tmp))
246 Control (UMFPACK_FIXQ) = tmp;
247 }
248
249 // Turn-off UMFPACK scaling for LU
250 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
251
252 umfpack_di_report_control (control);
253
254 const int *Ap = a.cidx ();
255 const int *Ai = a.ridx ();
256 const double *Ax = a.data ();
257
258 umfpack_di_report_matrix (nr, nc, Ap, Ai, Ax, 1, control);
259
260 void *Symbolic;
261 Matrix Info (1, UMFPACK_INFO);
262 double *info = Info.fortran_vec ();
263 int status;
264
265 // Null loop so that qinit is imediately deallocated when not needed
266 do {
267 OCTAVE_LOCAL_BUFFER (int, qinit, nc);
268
269 for (int i = 0; i < nc; i++)
270 qinit [i] = static_cast<int> (Qinit (i));
271
272 status = umfpack_di_qsymbolic (nr, nc, Ap, Ai, Ax, qinit,
273 &Symbolic, control, info);
274 } while (0);
275
276 if (status < 0)
277 {
278 (*current_liboctave_error_handler)
279 ("SparseLU::SparseLU symbolic factorization failed");
280
281 umfpack_di_report_status (control, status);
282 umfpack_di_report_info (control, info);
283
284 umfpack_di_free_symbolic (&Symbolic) ;
285 }
286 else
287 {
288 umfpack_di_report_symbolic (Symbolic, control);
289
290 void *Numeric;
291 status = umfpack_di_numeric (Ap, Ai, Ax, Symbolic, &Numeric,
292 control, info) ;
293 umfpack_di_free_symbolic (&Symbolic) ;
294
295 cond = Info (UMFPACK_RCOND);
296
297 if (status < 0)
298 {
299 (*current_liboctave_error_handler)
300 ("SparseLU::SparseLU numeric factorization failed");
301
302 umfpack_di_report_status (control, status);
303 umfpack_di_report_info (control, info);
304
305 umfpack_di_free_numeric (&Numeric);
306 }
307 else
308 {
309 umfpack_di_report_numeric (Numeric, control);
310
311 int lnz, unz, ignore1, ignore2, ignore3;
312 status = umfpack_di_get_lunz (&lnz, &unz, &ignore1, &ignore2,
313 &ignore3, Numeric) ;
314
315 if (status < 0)
316 {
317 (*current_liboctave_error_handler)
318 ("SparseLU::SparseLU extracting LU factors failed");
319
320 umfpack_di_report_status (control, status);
321 umfpack_di_report_info (control, info);
322
323 umfpack_di_free_numeric (&Numeric);
324 }
325 else
326 {
327 int n_inner = (nr < nc ? nr : nc);
328
329 if (lnz < 1)
330 Lfact = SparseMatrix (n_inner, nr, 1);
331 else
332 Lfact = SparseMatrix (n_inner, nr, lnz);
333
334 int *Ltp = Lfact.cidx ();
335 int *Ltj = Lfact.ridx ();
336 double *Ltx = Lfact.data ();
337
338 if (unz < 1)
339 Ufact = SparseMatrix (n_inner, nc, 1);
340 else
341 Ufact = SparseMatrix (n_inner, nc, unz);
342
343 int *Up = Ufact.cidx ();
344 int *Uj = Ufact.ridx ();
345 double *Ux = Ufact.data ();
346
347 P.resize (nr);
348 int *p = P.fortran_vec ();
349
350 Q.resize (nc);
351 int *q = Q.fortran_vec ();
352
353 int do_recip;
354 status = umfpack_di_get_numeric (Ltp, Ltj, Ltx, Up, Uj,
355 Ux, p, q, (double *) NULL,
356 &do_recip, (double *) NULL,
357 Numeric) ;
358
359 umfpack_di_free_numeric (&Numeric) ;
360
361 if (status < 0 || do_recip)
362 {
363 (*current_liboctave_error_handler)
364 ("SparseLU::SparseLU extracting LU factors failed");
365
366 umfpack_di_report_status (control, status);
367 }
368 else
369 {
370 Lfact = Lfact.transpose ();
371 umfpack_di_report_matrix (nr, n_inner, Lfact.cidx (),
372 Lfact.ridx (), Lfact.data (),
373 1, control);
374 umfpack_di_report_matrix (n_inner, nc, Ufact.cidx (),
375 Ufact.ridx (), Ufact.data (),
376 1, control);
377 umfpack_di_report_perm (nr, p, control);
378 umfpack_di_report_perm (nc, q, control);
379 }
380
381 umfpack_di_report_info (control, info);
382 }
383 }
384 }
385 }
386
387 /*
388 ;;; Local Variables: ***
389 ;;; mode: C++ ***
390 ;;; End: ***
391 */
392