5164
|
1 --- UMFPACKv4.4.orig/AMD/Makefile 2004-01-29 20:40:44.000000000 +0100 |
|
2 +++ UMFPACK/AMD/Makefile 2005-02-01 22:00:38.917972334 +0100 |
|
3 @@ -28,11 +28,18 @@ |
|
4 ( cd Demo ; make cross ) |
|
5 - cat Doc/License |
|
6 |
|
7 +# compile a Octave version |
|
8 +# (not compiled by "make all") |
|
9 +octave: |
|
10 + ( cd OCTAVE ; make ) |
|
11 + - cat Doc/License |
|
12 + |
|
13 # remove object files, but keep the compiled programs and library archives |
|
14 clean: |
|
15 ( cd Source ; make clean ) |
|
16 ( cd Demo ; make clean ) |
|
17 ( cd MATLAB ; make clean ) |
|
18 + ( cd OCTAVE ; make clean ) |
|
19 ( cd Doc ; make clean ) |
|
20 |
|
21 # clean, and then remove compiled programs and library archives |
|
22 @@ -40,6 +47,7 @@ |
|
23 ( cd Source ; make purge ) |
|
24 ( cd Demo ; make purge ) |
|
25 ( cd MATLAB ; make purge ) |
|
26 + ( cd OCTAVE ; make purge ) |
|
27 ( cd Doc ; make purge ) |
|
28 |
|
29 # create PDF documents for the original distribution |
|
30 --- UMFPACKv4.4.orig/UMFPACK/Makefile 2004-01-29 20:40:48.000000000 +0100 |
|
31 +++ UMFPACK/UMFPACK/Makefile 2005-02-01 22:00:38.916972398 +0100 |
|
32 @@ -32,12 +32,19 @@ |
|
33 hb: |
|
34 ( cd Demo ; make hb ) |
|
35 |
|
36 +# compile a Octave version |
|
37 +# (not compiled by "make all") |
|
38 +octave: |
|
39 + ( cd OCTAVE ; make ) |
|
40 + - cat Doc/License |
|
41 + |
|
42 # remove object files, but keep the compiled programs and library archives |
|
43 clean: |
|
44 ( cd ../AMD ; make clean ) |
|
45 ( cd Source ; make clean ) |
|
46 ( cd Demo ; make clean ) |
|
47 ( cd MATLAB ; make clean ) |
|
48 + ( cd OCTAVE ; make clean ) |
|
49 ( cd Doc ; make clean ) |
|
50 |
|
51 # clean, and then remove compiled programs and library archives |
|
52 @@ -46,6 +53,7 @@ |
|
53 ( cd Source ; make purge ) |
|
54 ( cd Demo ; make purge ) |
|
55 ( cd MATLAB ; make purge ) |
|
56 + ( cd OCTAVE ; make purge ) |
|
57 ( cd Doc ; make purge ) |
|
58 |
|
59 # create PDF documents for the original distribution |
|
60 diff -uNr UMFPACKv4.4.orig/UMFPACK/OCTAVE/Contents.m UMFPACK/UMFPACK/OCTAVE/Contents.m |
|
61 --- UMFPACKv4.4.orig/UMFPACK/OCTAVE/Contents.m 1970-01-01 01:00:00.000000000 +0100 |
|
62 +++ UMFPACK/UMFPACK/OCTAVE/Contents.m 2004-12-30 01:58:46.000000000 +0100 |
|
63 @@ -0,1 +1,22 @@ |
|
64 +%Contents of the UMFPACK sparse matrix toolbox: |
|
65 +% |
|
66 +% umfpack computes x=A\b, x=A/b, or lu (A) for a sparse matrix A |
|
67 +% umfpack_make to compile umfpack for use in MATLAB |
|
68 +% umfpack_report prints optional control settings and statistics |
|
69 +% umfpack_demo a long demo |
|
70 +% umfpack_simple a simple demo |
|
71 +% umfpack_btf factorize A using a block triangular form |
|
72 +% umfpack_solve x = A\b or x = b/A |
|
73 +% lu_normest estimates norm (L*U-A, 1) without forming L*U-A |
|
74 +% luflop given L and U, computes # of flops required to compute them |
|
75 +% |
|
76 +% See also: |
|
77 +% amd symmetric minimum degree ordering |
|
78 +% colamd unsymmetric column approx minimum degree ordering |
|
79 +% symamd symmetric approx minimum degree ordering, based on colamd |
|
80 +% |
|
81 +% UMFPACK Version 4.3 (Jan. 16, 2004), Copyright (c) 2004 by Timothy A. |
|
82 +% Davis. All Rights Reserved. Type umfpack_details for License. |
|
83 + |
|
84 +help Contents |
|
85 --- UMFPACKv4.4.orig/UMFPACK/OCTAVE/GNUmakefile 1970-01-01 01:00:00.000000000 +0100 |
|
86 +++ UMFPACK/UMFPACK/OCTAVE/GNUmakefile 2004-12-30 01:58:46.000000000 +0100 |
|
87 @@ -0,0 +1,248 @@ |
|
88 +#------------------------------------------------------------------------------- |
|
89 +# UMFPACK GNUmakefile for the UMFPACK OCTAVE oct-file (GNU "make" only) |
|
90 +#------------------------------------------------------------------------------- |
|
91 + |
|
92 +.PRECIOUS: %.o |
|
93 + |
|
94 +# UMFPACK Version 4.3 (Jan. 16, 2004), Copyright (c) 2004 by Timothy A. |
|
95 +# Davis. All Rights Reserved. See ../README for License. |
|
96 + |
|
97 +all: umfpack luflop |
|
98 + |
|
99 +include ../Make/Make.include |
|
100 + |
|
101 +MKOCT = mkoctfile $(CONFIG) -DNRECIPROCAL -I/usr/include/atlas -I../Include -I../Source -I../../AMD/Include -I../../AMD/Source |
|
102 + |
|
103 +OCT_SPARSE_INC = -I../../../ |
|
104 + |
|
105 + |
|
106 +#------------------------------------------------------------------------------- |
|
107 +# source files |
|
108 +#------------------------------------------------------------------------------- |
|
109 + |
|
110 +# non-user-callable umf_*.[ch] files: |
|
111 +UMFCH = umf_assemble umf_blas3_update \ |
|
112 + umf_build_tuples umf_create_element \ |
|
113 + umf_dump umf_extend_front umf_garbage_collection \ |
|
114 + umf_get_memory umf_init_front umf_kernel \ |
|
115 + umf_kernel_init umf_kernel_wrapup \ |
|
116 + umf_local_search umf_lsolve umf_ltsolve \ |
|
117 + umf_mem_alloc_element umf_mem_alloc_head_block \ |
|
118 + umf_mem_alloc_tail_block umf_mem_free_tail_block \ |
|
119 + umf_mem_init_memoryspace \ |
|
120 + umf_report_vector umf_row_search umf_scale_column \ |
|
121 + umf_set_stats umf_solve umf_symbolic_usage umf_transpose \ |
|
122 + umf_tuple_lengths umf_usolve umf_utsolve umf_valid_numeric \ |
|
123 + umf_valid_symbolic umf_grow_front umf_start_front umf_2by2 \ |
|
124 + umf_store_lu umf_scale |
|
125 + |
|
126 +# non-user-callable umf_*.[ch] files, int/long versions only (no real/complex): |
|
127 +UMFINT = umf_analyze umf_apply_order umf_colamd umf_free umf_fsize \ |
|
128 + umf_is_permutation umf_malloc umf_realloc umf_report_perm \ |
|
129 + umf_singletons |
|
130 + |
|
131 +# non-user-callable and user-callable amd_*.[ch] files (int/long versions only): |
|
132 +AMD = amd_aat amd_1 amd_2 amd_dump amd_postorder amd_post_tree amd_defaults \ |
|
133 + amd_order amd_control amd_info amd_valid |
|
134 + |
|
135 +# non-user-callable, created from umf_ltsolve.c, umf_utsolve.c, |
|
136 +# umf_triplet.c, and umf_assemble.c , with int/long and real/complex versions: |
|
137 +UMF_CREATED = umf_lhsolve umf_uhsolve umf_triplet_map_nox \ |
|
138 + umf_triplet_nomap_x umf_triplet_nomap_nox umf_triplet_map_x \ |
|
139 + umf_assemble_fixq umf_store_lu_drop |
|
140 + |
|
141 +# non-user-callable, int/long and real/complex versions: |
|
142 +UMF = $(UMF_CREATED) $(UMFCH) |
|
143 + |
|
144 +# user-callable umfpack_*.[ch] files (int/long and real/complex): |
|
145 +UMFPACK = umfpack_col_to_triplet umfpack_defaults umfpack_free_numeric \ |
|
146 + umfpack_free_symbolic umfpack_get_numeric umfpack_get_lunz \ |
|
147 + umfpack_get_symbolic umfpack_get_determinant umfpack_numeric \ |
|
148 + umfpack_qsymbolic umfpack_report_control umfpack_report_info \ |
|
149 + umfpack_report_matrix umfpack_report_numeric umfpack_report_perm \ |
|
150 + umfpack_report_status umfpack_report_symbolic umfpack_report_triplet \ |
|
151 + umfpack_report_vector umfpack_solve umfpack_symbolic \ |
|
152 + umfpack_transpose umfpack_triplet_to_col umfpack_scale \ |
|
153 + umfpack_load_numeric umfpack_save_numeric \ |
|
154 + umfpack_load_symbolic umfpack_save_symbolic |
|
155 + |
|
156 +# user-callable, created from umfpack_solve.c (umfpack_wsolve.h exists, though): |
|
157 +# with int/long and real/complex versions: |
|
158 +UMFPACKW = umfpack_wsolve |
|
159 + |
|
160 +USER = $(UMFPACKW) $(UMFPACK) |
|
161 + |
|
162 +# user-callable, only one version for int/long, real/complex, *.[ch] files: |
|
163 +GENERIC = umfpack_timer umfpack_tictoc |
|
164 + |
|
165 +#------------------------------------------------------------------------------- |
|
166 +# include files: |
|
167 +#------------------------------------------------------------------------------- |
|
168 + |
|
169 +AMDH = ../../AMD/Source/amd_internal.h ../../AMD/Include/amd.h |
|
170 + |
|
171 +INC1 = umf_config.h umf_version.h umf_internal.h umf_triplet.h |
|
172 + |
|
173 +INC = ../Include/umfpack.h \ |
|
174 + $(addprefix ../Source/, $(INC1)) \ |
|
175 + $(addprefix ../Source/, $(addsuffix .h,$(UMFCH))) \ |
|
176 + $(addprefix ../Source/, $(addsuffix .h,$(UMFINT))) \ |
|
177 + $(addprefix ../Include/, $(addsuffix .h,$(USER))) \ |
|
178 + $(addprefix ../Include/, $(addsuffix .h,$(GENERIC))) \ |
|
179 + $(AMDH) |
|
180 + |
|
181 +#------------------------------------------------------------------------------- |
|
182 +# Create the umfpack and amd oct-file for OCTAVE (int versions only) |
|
183 +#------------------------------------------------------------------------------- |
|
184 + |
|
185 +OCTI = $(addsuffix .o, $(subst umf_,umf_o_,$(UMFINT))) |
|
186 +OCTDI = $(addsuffix .o, $(subst umf_,umf_od_,$(UMF)) $(subst umfpack_,umfpack_od_,$(USER))) |
|
187 +OCTZI = $(addsuffix .o, $(subst umf_,umf_oz_,$(UMF)) $(subst umfpack_,umfpack_oz_,$(USER)) ) |
|
188 +OCTAMD = $(addsuffix .o, $(subst amd_,amd_o_,$(AMD))) |
|
189 +OCTGN = $(addsuffix .o, $(subst umfpack_,umfpack_o_,$(GENERIC))) |
|
190 + |
|
191 +OCTUMFPACK = $(OCTI) $(OCTDI) $(OCTZI) $(OCTGN) |
|
192 + |
|
193 +OCTUMFPACK_LIB = umfpack_octave.o |
|
194 + |
|
195 +# Note that mkoctfile has an "-o" option, but it doesn't work in conjunction |
|
196 +# with the "-c" option, thus the need for $(MV) commands. |
|
197 +# If it did, then the rules would be much simpler: |
|
198 +# $(MKOCT) -DDINT -c $< -o $@ |
|
199 + |
|
200 +#---------------------------------------- |
|
201 +# integer-only routines (no real/complex): |
|
202 +#---------------------------------------- |
|
203 + |
|
204 +amd_o_%.o: ../../AMD/Source/amd_%.c $(AMDH) |
|
205 + $(MKOCT) -DDINT -c $< |
|
206 + - $(MV) ../../AMD/Source/amd_$*.o $@ |
|
207 + |
|
208 + |
|
209 +umf_o_%.o: ../Source/umf_%.c $(INC) |
|
210 + $(MKOCT) -DDINT -c $< |
|
211 + - $(MV) ../Source/umf_$*.o $@ |
|
212 + |
|
213 +#---------------------------------------- |
|
214 +# Double precision, int version, for OCTAVE |
|
215 +#---------------------------------------- |
|
216 + |
|
217 +umf_od_%.o: ../Source/umf_%.c $(INC) |
|
218 + $(MKOCT) -DDINT -c $< |
|
219 + - $(MV) ../Source/umf_$*.o $@ |
|
220 + |
|
221 +umf_od_%hsolve.o: ../Source/umf_%tsolve.c $(INC) |
|
222 + $(MKOCT) -DDINT -DCONJUGATE_SOLVE -c $< |
|
223 + - $(MV) ../Source/umf_$*tsolve.o $@ |
|
224 + |
|
225 +umf_od_triplet_map_x.o: ../Source/umf_triplet.c $(INC) |
|
226 + $(MKOCT) -DDINT -DDO_MAP -DDO_VALUES -c $< |
|
227 + - $(MV) ../Source/umf_triplet.o $@ |
|
228 + |
|
229 +umf_od_triplet_map_nox.o: ../Source/umf_triplet.c $(INC) |
|
230 + $(MKOCT) -DDINT -DDO_MAP -c $< |
|
231 + - $(MV) ../Source/umf_triplet.o $@ |
|
232 + |
|
233 +umf_od_triplet_nomap_x.o: ../Source/umf_triplet.c $(INC) |
|
234 + $(MKOCT) -DDINT -DDO_VALUES -c $< |
|
235 + - $(MV) ../Source/umf_triplet.o $@ |
|
236 + |
|
237 +umf_od_triplet_nomap_nox.o: ../Source/umf_triplet.c $(INC) |
|
238 + $(MKOCT) -DDINT -c $< |
|
239 + - $(MV) ../Source/umf_triplet.o $@ |
|
240 + |
|
241 +umf_od_assemble_fixq.o: ../Source/umf_assemble.c $(INC) |
|
242 + $(MKOCT) -DDINT -DFIXQ -c $< |
|
243 + - $(MV) ../Source/umf_assemble.o $@ |
|
244 + |
|
245 +umf_od_store_lu_drop.o: ../Source/umf_store_lu.c $(INC) |
|
246 + $(MKOCT) -DDINT -DDROP -c $< |
|
247 + - $(MV) ../Source/umf_store_lu.o $@ |
|
248 + |
|
249 +umfpack_od_wsolve.o: ../Source/umfpack_solve.c $(INC) |
|
250 + $(MKOCT) -DDINT -DWSOLVE -c $< |
|
251 + - $(MV) ../Source/umfpack_solve.o $@ |
|
252 + |
|
253 +umfpack_od_%.o: ../Source/umfpack_%.c $(INC) |
|
254 + $(MKOCT) -DDINT -c $< |
|
255 + - $(MV) ../Source/umfpack_$*.o $@ |
|
256 + |
|
257 +#---------------------------------------- |
|
258 +# Complex double precision, int version, for OCTAVE |
|
259 +#---------------------------------------- |
|
260 + |
|
261 +umf_oz_%.o: ../Source/umf_%.c $(INC) |
|
262 + $(MKOCT) -DZINT -c $< |
|
263 + - $(MV) ../Source/umf_$*.o $@ |
|
264 + |
|
265 +umf_oz_%hsolve.o: ../Source/umf_%tsolve.c $(INC) |
|
266 + $(MKOCT) -DZINT -DCONJUGATE_SOLVE -c $< |
|
267 + - $(MV) ../Source/umf_$*tsolve.o $@ |
|
268 + |
|
269 +umf_oz_triplet_map_x.o: ../Source/umf_triplet.c $(INC) |
|
270 + $(MKOCT) -DZINT -DDO_MAP -DDO_VALUES -c $< |
|
271 + - $(MV) ../Source/umf_triplet.o $@ |
|
272 + |
|
273 +umf_oz_triplet_map_nox.o: ../Source/umf_triplet.c $(INC) |
|
274 + $(MKOCT) -DZINT -DDO_MAP -c $< |
|
275 + - $(MV) ../Source/umf_triplet.o $@ |
|
276 + |
|
277 +umf_oz_triplet_nomap_x.o: ../Source/umf_triplet.c $(INC) |
|
278 + $(MKOCT) -DZINT -DDO_VALUES -c $< |
|
279 + - $(MV) ../Source/umf_triplet.o $@ |
|
280 + |
|
281 +umf_oz_triplet_nomap_nox.o: ../Source/umf_triplet.c $(INC) |
|
282 + $(MKOCT) -DZINT -c $< |
|
283 + - $(MV) ../Source/umf_triplet.o $@ |
|
284 + |
|
285 +umf_oz_assemble_fixq.o: ../Source/umf_assemble.c $(INC) |
|
286 + $(MKOCT) -DZINT -DFIXQ -c $< |
|
287 + - $(MV) ../Source/umf_assemble.o $@ |
|
288 + |
|
289 +umf_oz_store_lu_drop.o: ../Source/umf_store_lu.c $(INC) |
|
290 + $(MKOCT) -DZINT -DDROP -c $< |
|
291 + - $(MV) ../Source/umf_store_lu.o $@ |
|
292 + |
|
293 +umfpack_oz_wsolve.o: ../Source/umfpack_solve.c $(INC) |
|
294 + $(MKOCT) -DZINT -DWSOLVE -c $< |
|
295 + - $(MV) ../Source/umfpack_solve.o $@ |
|
296 + |
|
297 +umfpack_oz_%.o: ../Source/umfpack_%.c $(INC) |
|
298 + $(MKOCT) -DZINT -c $< |
|
299 + - $(MV) ../Source/umfpack_$*.o $@ |
|
300 + |
|
301 +#---------------------------------------- |
|
302 +# Generic routines for OCTAVE |
|
303 +#---------------------------------------- |
|
304 + |
|
305 +umfpack_o_timer.o: ../Source/umfpack_timer.c $(INC) |
|
306 + $(MKOCT) -c $< |
|
307 + - $(MV) ../Source/umfpack_timer.o $@ |
|
308 + |
|
309 +umfpack_o_tictoc.o: ../Source/umfpack_tictoc.c $(INC) |
|
310 + $(MKOCT) -c $< |
|
311 + - $(MV) ../Source/umfpack_tictoc.o $@ |
|
312 + |
|
313 +#---------------------------------------- |
|
314 +# umfpack oct-files |
|
315 +#---------------------------------------- |
|
316 + |
|
317 +umfpack: umfpack.cc $(OCTUMFPACK) $(OCTAMD) |
|
318 + $(MKOCT) $(OCT_SPARSE_INC) umfpack.cc $(OCTUMFPACK) $(OCTAMD) -o umfpack.oct |
|
319 + |
|
320 +luflop: luflop.cc |
|
321 + $(MKOCT) $(OCT_SPARSE_INC) luflop.cc -o luflop.oct |
|
322 + |
|
323 +#---------------------------------------- |
|
324 +# umfpack library to link with octave |
|
325 +#---------------------------------------- |
|
326 + |
|
327 +octave: $(OCTUMFPACK) $(OCTAMD) |
|
328 + ld -r $(OCTUMFPACK) $(OCTAMD) -o ../../../$(OCTUMFPACK_LIB) |
|
329 + |
|
330 +#------------------------------------------------------------------------------- |
|
331 +# Remove all but the files in the original distribution |
|
332 +#------------------------------------------------------------------------------- |
|
333 + |
|
334 +purge: clean |
|
335 + - $(RM) *.oct *.so |
|
336 diff -uNr UMFPACKv4.4.orig/UMFPACK/OCTAVE/luflop.cc UMFPACK/UMFPACK/OCTAVE/luflop.cc |
|
337 --- UMFPACKv4.4.orig/UMFPACK/OCTAVE/luflop.cc 1970-01-01 01:00:00.000000000 +0100 |
|
338 +++ UMFPACK/UMFPACK/OCTAVE/luflop.cc 2004-12-30 01:58:46.000000000 +0100 |
|
339 @@ -0,0 +1,223 @@ |
|
340 +/* |
|
341 + |
|
342 +Copyright (C) 2004 David Bateman |
|
343 + |
|
344 +This program is free software; you can redistribute it and/or modify it |
|
345 +under the terms of the GNU General Public License as published by the |
|
346 +Free Software Foundation; either version 2, or (at your option) any |
|
347 +later version. |
|
348 + |
|
349 +This program is distributed in the hope that it will be useful, but WITHOUT |
|
350 +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
|
351 +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
|
352 +for more details. |
|
353 + |
|
354 +You should have received a copy of the GNU General Public License |
|
355 +along with this program; see the file COPYING. If not, write to the Free |
|
356 +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. |
|
357 + |
|
358 +In addition to the terms of the GPL, you are permitted to link |
|
359 +this program with any Open Source program, as defined by the |
|
360 +Open Source Initiative (www.opensource.org) |
|
361 + |
|
362 +*/ |
|
363 + |
|
364 +/* |
|
365 + |
|
366 +This is the Octave interface to the UMFPACK code, which bore the following |
|
367 +copyright |
|
368 + |
|
369 + UMFPACK Version 4.3 (Jan. 16, 2004), Copyright (c) 2004 by Timothy A. |
|
370 + Davis. All Rights Reserved. See ../README for License. |
|
371 + email: davis@cise.ufl.edu CISE Department, Univ. of Florida. |
|
372 + web: http://www.cise.ufl.edu/research/sparse/umfpack |
|
373 + |
|
374 +*/ |
|
375 + |
|
376 + |
|
377 +#include <cstdlib> |
|
378 +#include <string> |
|
379 + |
|
380 +#include <octave/config.h> |
|
381 +#include <octave/ov.h> |
|
382 +#include <octave/defun-dld.h> |
|
383 +#include <octave/pager.h> |
|
384 +#include <octave/ov-re-mat.h> |
|
385 + |
|
386 +#include "ov-re-sparse.h" |
|
387 +#include "ov-cx-sparse.h" |
|
388 + |
|
389 +DEFUN_DLD (luflop, args, nargout, |
|
390 + "-*- texinfo -*-\n\ |
|
391 +@deftypefn {Loadable Function} {@var{f} =} luflup (@var{l}, @var{u})\n\ |
|
392 +\n\ |
|
393 +Given an LU factorization, compute how many flops took to compute it. This\n\ |
|
394 +is the same as (assuming U has a zero-free diagonal):\n\ |
|
395 +\n\ |
|
396 +@example\n\ |
|
397 +@group\n\ |
|
398 + Lnz = full (sum (spones (L))) - 1 ;\n\ |
|
399 + Unz = full (sum (spones (U')))' - 1 ;\n\ |
|
400 + f = 2*Lnz*Unz + sum (Lnz) ;\n\ |
|
401 +@end group\n\ |
|
402 +@end example\n\ |
|
403 +\n\ |
|
404 +except that no extra workspace is allocated for spones (L) and spones (U).\n\ |
|
405 +L and U must be sparse.\n\ |
|
406 +\n\ |
|
407 +Note: the above expression has a subtle undercount when exact numerical\n\ |
|
408 +cancelation occurs. Try [L,U,P] = lu (sparse (ones (10))) and then\n\ |
|
409 +luflop (L,U).\n\ |
|
410 +@end deftypefn") |
|
411 +{ |
|
412 + int nargin = args.length (); |
|
413 + octave_value retval; |
|
414 + |
|
415 + // These are here only so that the C++ destructors don't prematurally |
|
416 + // remove the underlying data we are interested in |
|
417 + SparseMatrix Lmatrix, Umatrix; |
|
418 + SparseComplexMatrix CLmatrix, CUmatrix; |
|
419 + int *Lp, *Li, *Up, *Ui; |
|
420 + int Lm, Ln, Um, Un; |
|
421 + |
|
422 + if (nargin != 2) |
|
423 + { |
|
424 + usage ("f = luflop (L, U)"); |
|
425 + return retval; |
|
426 + } |
|
427 + |
|
428 + if (args(0).class_name () == "sparse") |
|
429 + { |
|
430 + if (args(0).is_complex_type ()) |
|
431 + { |
|
432 + CLmatrix = |
|
433 + (((const octave_sparse_complex_matrix&) args(0).get_rep ()) |
|
434 + .sparse_complex_matrix_value ()); |
|
435 + Lp = CLmatrix.cidx (); |
|
436 + Li = CLmatrix.ridx (); |
|
437 + Lm = CLmatrix.rows (); |
|
438 + Ln = CLmatrix.cols (); |
|
439 + } |
|
440 + else |
|
441 + { |
|
442 + Lmatrix = (((const octave_sparse_matrix&) args(0).get_rep ()) |
|
443 + .sparse_matrix_value ()); |
|
444 + Lp = Lmatrix.cidx (); |
|
445 + Li = Lmatrix.ridx (); |
|
446 + Lm = Lmatrix.rows (); |
|
447 + Ln = Lmatrix.cols (); |
|
448 + } |
|
449 + } |
|
450 + else |
|
451 + { |
|
452 + if (args(0).is_complex_type ()) |
|
453 + { |
|
454 + CLmatrix = SparseComplexMatrix (args(0).complex_matrix_value ()); |
|
455 + Lp = CLmatrix.cidx (); |
|
456 + Li = CLmatrix.ridx (); |
|
457 + Lm = CLmatrix.rows (); |
|
458 + Ln = CLmatrix.cols (); |
|
459 + } |
|
460 + else |
|
461 + { |
|
462 + Lmatrix = SparseMatrix (args(0).matrix_value ()); |
|
463 + Lp = Lmatrix.cidx (); |
|
464 + Li = Lmatrix.ridx (); |
|
465 + Lm = Lmatrix.rows (); |
|
466 + Ln = Lmatrix.cols (); |
|
467 + } |
|
468 + } |
|
469 + |
|
470 + |
|
471 + if (args(0).class_name () == "sparse") |
|
472 + { |
|
473 + if (args(1).is_complex_type ()) |
|
474 + { |
|
475 + CUmatrix = |
|
476 + (((const octave_sparse_complex_matrix&) args(1).get_rep ()) |
|
477 + .sparse_complex_matrix_value ()); |
|
478 + Up = CUmatrix.cidx (); |
|
479 + Ui = CUmatrix.ridx (); |
|
480 + Um = CUmatrix.rows (); |
|
481 + Un = CUmatrix.cols (); |
|
482 + } |
|
483 + else |
|
484 + { |
|
485 + Umatrix = |
|
486 + (((const octave_sparse_matrix&) args(1).get_rep ()) |
|
487 + .sparse_matrix_value ()); |
|
488 + Up = Umatrix.cidx (); |
|
489 + Ui = Umatrix.ridx (); |
|
490 + Um = Umatrix.rows (); |
|
491 + Un = Umatrix.cols (); |
|
492 + } |
|
493 + } |
|
494 + else |
|
495 + { |
|
496 + if (args(1).is_complex_type ()) |
|
497 + { |
|
498 + CUmatrix = SparseComplexMatrix (args(1).complex_matrix_value ()); |
|
499 + Up = CUmatrix.cidx (); |
|
500 + Ui = CUmatrix.ridx (); |
|
501 + Um = CUmatrix.rows (); |
|
502 + Un = CUmatrix.cols (); |
|
503 + } |
|
504 + else |
|
505 + { |
|
506 + Umatrix = SparseMatrix (args(1).matrix_value ()); |
|
507 + Up = Umatrix.cidx (); |
|
508 + Ui = Umatrix.ridx (); |
|
509 + Um = Umatrix.rows (); |
|
510 + Un = Umatrix.cols (); |
|
511 + } |
|
512 + } |
|
513 + |
|
514 + |
|
515 + if (error_state) |
|
516 + return retval; |
|
517 + |
|
518 + |
|
519 + int n = Lm; |
|
520 + |
|
521 + if (n != Ln || n != Um || n != Un) |
|
522 + error ("luflop: L and U must be square"); |
|
523 + else |
|
524 + { |
|
525 + |
|
526 + OCTAVE_LOCAL_BUFFER (int, Unz, n); |
|
527 + |
|
528 + // count the nonzeros in each row of U |
|
529 + for (int row = 0 ; row < n ; row++) |
|
530 + Unz [row] = 0 ; |
|
531 + |
|
532 + for (int col = 0 ; col < n ; col++) |
|
533 + { |
|
534 + for (int p = Up [col] ; p < Up [col+1] ; p++) |
|
535 + { |
|
536 + int row = Ui [p] ; |
|
537 + Unz [row]++ ; |
|
538 + } |
|
539 + } |
|
540 + |
|
541 + // count the flops |
|
542 + double flop_count = 0.0 ; |
|
543 + for (int k = 0 ; k < n ; k++) |
|
544 + { |
|
545 + /* off-diagonal nonzeros in column k of L: */ |
|
546 + int Lnz_k = Lp [k+1] - Lp [k] - 1 ; |
|
547 + int Unz_k = Unz [k] - 1 ; |
|
548 + flop_count += (2 * Lnz_k * Unz_k) + Lnz_k ; |
|
549 + } |
|
550 + |
|
551 + // return the result |
|
552 + retval = flop_count; |
|
553 + } |
|
554 + |
|
555 + return retval; |
|
556 +} |
|
557 + |
|
558 +/* |
|
559 +;;; Local Variables: *** |
|
560 +;;; mode: C++ *** |
|
561 +;;; End: *** |
|
562 +*/ |
|
563 diff -uNr UMFPACKv4.4.orig/UMFPACK/OCTAVE/lu_normest.m UMFPACK/UMFPACK/OCTAVE/lu_normest.m |
|
564 --- UMFPACKv4.4.orig/UMFPACK/OCTAVE/lu_normest.m 1970-01-01 01:00:00.000000000 +0100 |
|
565 +++ UMFPACK/UMFPACK/OCTAVE/lu_normest.m 2004-12-30 01:58:46.000000000 +0100 |
|
566 @@ -0,0 +1,106 @@ |
|
567 +function rho = lu_normest (A, L, U) |
|
568 +% LU_NORMEST: estimate the 1-norm of A-L*U without computing L*U |
|
569 +% |
|
570 +% Usage: |
|
571 +% |
|
572 +% rho = lu_normest (A, L, U) |
|
573 +% |
|
574 +% which estimates the computation of the 1-norm: |
|
575 +% |
|
576 +% rho = norm (A-L*U, 1) |
|
577 +% |
|
578 +% Authors: William W. Hager, Math Dept., Univ. of Florida |
|
579 +% Timothy A. Davis, CISE Dept., Univ. of Florida |
|
580 +% Gainesville, FL, 32611, USA. |
|
581 +% based on normest1, contributed on November, 1997 |
|
582 +% |
|
583 +% This code can be quite easily adapted to estimate the 1-norm of any |
|
584 +% matrix E, where E itself is dense or not explicitly represented, but the |
|
585 +% computation of E (and E') times a vector is easy. In this case, our matrix |
|
586 +% of interest is: |
|
587 +% |
|
588 +% E = A-L*U |
|
589 +% |
|
590 +% That is, L*U is the LU factorization of A, where A, L and U |
|
591 +% are sparse. This code works for dense matrices A and L too, |
|
592 +% but it would not be needed in that case, since E is easy to compute |
|
593 +% explicitly. For sparse A, L, and U, computing E explicitly would be quite |
|
594 +% expensive, and thus normest (A-L*U) would be prohibitive. |
|
595 +% |
|
596 +% For a detailed description, see Davis, T. A. and Hager, W. W., |
|
597 +% Modifying a sparse Cholesky factorization, SIAM J. Matrix Analysis and |
|
598 +% Applications, 1999, vol. 20, no. 3, 606-627. |
|
599 + |
|
600 +% The three places that the matrix-vector multiply E*x is used are highlighted. |
|
601 +% Note that E is never formed explicity. |
|
602 + |
|
603 +[m n] = size (A) ; |
|
604 + |
|
605 +if (m ~= n) |
|
606 + % pad A, L, and U with zeros so that they are all square |
|
607 + if (m < n) |
|
608 + U = [ U ; (sparse (n-m,n)) ] ; |
|
609 + L = [ L , (sparse (m,n-m)) ; (sparse (n-m,n)) ] ; |
|
610 + A = [ A ; (sparse (n-m,n)) ] ; |
|
611 + else |
|
612 + U = [ U , (sparse (n,m-n)) ; (sparse (m-n,m)) ] ; |
|
613 + L = [ L , (sparse (m,m-n)) ] ; |
|
614 + A = [ A , (sparse (m,m-n)) ] ; |
|
615 + end |
|
616 +end |
|
617 + |
|
618 +[m n] = size (A) ; |
|
619 + |
|
620 +notvisited = ones (m, 1) ; % nonvisited(j) is zero if j is visited, 1 otherwise |
|
621 +rho = 0 ; % the global rho |
|
622 + |
|
623 +At = A' ; |
|
624 +Lt = L' ; |
|
625 + |
|
626 +for trial = 1:3 % { |
|
627 + |
|
628 + x = notvisited ./ sum (notvisited) ; |
|
629 + rho1 = 0 ; % the current rho for this trial |
|
630 + |
|
631 + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
|
632 + %%% COMPUTE Ex1 = E*x EFFICIENTLY: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
|
633 + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
|
634 + Ex1 = (A*x) - L*(U*x) ; |
|
635 + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
|
636 + |
|
637 + rho2 = norm (Ex1, 1) ; |
|
638 + |
|
639 + while rho2 > rho1 % { |
|
640 + |
|
641 + rho1 = rho2 ; |
|
642 + y = 2*(Ex1 >= 0) - 1 ; |
|
643 + |
|
644 + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
|
645 + %%% COMPUTE z = E'*y EFFICIENTLY: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
|
646 + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
|
647 + z = (A'*y) - U'*(L'*y) ; |
|
648 + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
|
649 + |
|
650 + [zj, j] = max (abs (z .* notvisited)) ; |
|
651 + j = j (1) ; |
|
652 + if (abs (z (j)) > z'*x) % { |
|
653 + x = zeros (m, 1) ; |
|
654 + x (j) = 1 ; |
|
655 + notvisited (j) = 0 ; |
|
656 + else % } { |
|
657 + break ; |
|
658 + end % } |
|
659 + |
|
660 + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
|
661 + %%% COMPUTE Ex1 = E*x EFFICIENTLY: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
|
662 + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
|
663 + Ex1 = (A*x) - L*(U*x) ; |
|
664 + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
|
665 + |
|
666 + rho2 = norm (Ex1, 1) ; |
|
667 + |
|
668 + end % } |
|
669 + |
|
670 + rho = max (rho, rho1) ; |
|
671 + |
|
672 +end % } |
|
673 diff -uNr UMFPACKv4.4.orig/UMFPACK/OCTAVE/Makefile UMFPACK/UMFPACK/OCTAVE/Makefile |
|
674 --- UMFPACKv4.4.orig/UMFPACK/OCTAVE/Makefile 1970-01-01 01:00:00.000000000 +0100 |
|
675 +++ UMFPACK/UMFPACK/OCTAVE/Makefile 2004-12-30 01:58:46.000000000 +0100 |
|
676 @@ -0,0 +1,517 @@ |
|
677 +#------------------------------------------------------------------------------- |
|
678 +# UMFPACK Makefile for the UMFPACK MATLAB mexFunction (old "make" only) |
|
679 +#------------------------------------------------------------------------------- |
|
680 + |
|
681 +# UMFPACK Version 4.3 (Jan. 16, 2004), Copyright (c) 2004 by Timothy A. |
|
682 +# Davis. All Rights Reserved. See ../README for License. |
|
683 + |
|
684 +# This is a very ugly Makefile, and is only provided for those who do not |
|
685 +# have GNU make. Note that it is not used if you have GNU make. It ignores |
|
686 +# dependency checking and just compiles everything. It was created |
|
687 +# automatically, via make -n using the GNUmakefile. That way, I don't have |
|
688 +# maintain two Makefiles. |
|
689 + |
|
690 +all: umfpack luflop |
|
691 + |
|
692 +include ../Make/Make.include |
|
693 + |
|
694 +MKOCT = mkoctfile $(CONFIG) -DNRECIPROCAL -I/usr/include/atlas -I../Include -I../Source -I../../AMD/Include -I../../AMD/Source |
|
695 + |
|
696 +OCT_SPARSE_INC = -I../../../ |
|
697 +OCTUMFPACK_LIB = umfpack_octave.o |
|
698 + |
|
699 +umfpack: |
|
700 + $(MKOCT) -DDINT -c ../Source/umf_analyze.c |
|
701 + $(MV) -f ../Source/umf_analyze.o umf_m_analyze.o |
|
702 + $(MKOCT) -DDINT -c ../Source/umf_apply_order.c |
|
703 + $(MV) -f ../Source/umf_apply_order.o umf_m_apply_order.o |
|
704 + $(MKOCT) -DDINT -c ../Source/umf_colamd.c |
|
705 + $(MV) -f ../Source/umf_colamd.o umf_m_colamd.o |
|
706 + $(MKOCT) -DDINT -c ../Source/umf_free.c |
|
707 + $(MV) -f ../Source/umf_free.o umf_m_free.o |
|
708 + $(MKOCT) -DDINT -c ../Source/umf_fsize.c |
|
709 + $(MV) -f ../Source/umf_fsize.o umf_m_fsize.o |
|
710 + $(MKOCT) -DDINT -c ../Source/umf_is_permutation.c |
|
711 + $(MV) -f ../Source/umf_is_permutation.o umf_m_is_permutation.o |
|
712 + $(MKOCT) -DDINT -c ../Source/umf_malloc.c |
|
713 + $(MV) -f ../Source/umf_malloc.o umf_m_malloc.o |
|
714 + $(MKOCT) -DDINT -c ../Source/umf_realloc.c |
|
715 + $(MV) -f ../Source/umf_realloc.o umf_m_realloc.o |
|
716 + $(MKOCT) -DDINT -c ../Source/umf_report_perm.c |
|
717 + $(MV) -f ../Source/umf_report_perm.o umf_m_report_perm.o |
|
718 + $(MKOCT) -DDINT -c ../Source/umf_singletons.c |
|
719 + $(MV) -f ../Source/umf_singletons.o umf_m_singletons.o |
|
720 + $(MKOCT) -DDINT -DCONJUGATE_SOLVE -c ../Source/umf_ltsolve.c |
|
721 + $(MV) -f ../Source/umf_ltsolve.o umf_od_lhsolve.o |
|
722 + $(MKOCT) -DDINT -DCONJUGATE_SOLVE -c ../Source/umf_utsolve.c |
|
723 + $(MV) -f ../Source/umf_utsolve.o umf_od_uhsolve.o |
|
724 + $(MKOCT) -DDINT -DDO_MAP -c ../Source/umf_triplet.c |
|
725 + $(MV) -f ../Source/umf_triplet.o umf_od_triplet_map_nox.o |
|
726 + $(MKOCT) -DDINT -DDO_VALUES -c ../Source/umf_triplet.c |
|
727 + $(MV) -f ../Source/umf_triplet.o umf_od_triplet_nomap_x.o |
|
728 + $(MKOCT) -DDINT -c ../Source/umf_triplet.c |
|
729 + $(MV) -f ../Source/umf_triplet.o umf_od_triplet_nomap_nox.o |
|
730 + $(MKOCT) -DDINT -DDO_MAP -DDO_VALUES -c ../Source/umf_triplet.c |
|
731 + $(MV) -f ../Source/umf_triplet.o umf_od_triplet_map_x.o |
|
732 + $(MKOCT) -DDINT -DFIXQ -c ../Source/umf_assemble.c |
|
733 + $(MV) -f ../Source/umf_assemble.o umf_od_assemble_fixq.o |
|
734 + $(MKOCT) -DDINT -DDROP -c ../Source/umf_store_lu.c |
|
735 + $(MV) -f ../Source/umf_store_lu.o umf_od_store_lu_drop.o |
|
736 + $(MKOCT) -DDINT -c ../Source/umf_assemble.c |
|
737 + $(MV) -f ../Source/umf_assemble.o umf_od_assemble.o |
|
738 + $(MKOCT) -DDINT -c ../Source/umf_blas3_update.c |
|
739 + $(MV) -f ../Source/umf_blas3_update.o umf_od_blas3_update.o |
|
740 + $(MKOCT) -DDINT -c ../Source/umf_build_tuples.c |
|
741 + $(MV) -f ../Source/umf_build_tuples.o umf_od_build_tuples.o |
|
742 + $(MKOCT) -DDINT -c ../Source/umf_create_element.c |
|
743 + $(MV) -f ../Source/umf_create_element.o umf_od_create_element.o |
|
744 + $(MKOCT) -DDINT -c ../Source/umf_dump.c |
|
745 + $(MV) -f ../Source/umf_dump.o umf_od_dump.o |
|
746 + $(MKOCT) -DDINT -c ../Source/umf_extend_front.c |
|
747 + $(MV) -f ../Source/umf_extend_front.o umf_od_extend_front.o |
|
748 + $(MKOCT) -DDINT -c ../Source/umf_garbage_collection.c |
|
749 + $(MV) -f ../Source/umf_garbage_collection.o umf_od_garbage_collection.o |
|
750 + $(MKOCT) -DDINT -c ../Source/umf_get_memory.c |
|
751 + $(MV) -f ../Source/umf_get_memory.o umf_od_get_memory.o |
|
752 + $(MKOCT) -DDINT -c ../Source/umf_init_front.c |
|
753 + $(MV) -f ../Source/umf_init_front.o umf_od_init_front.o |
|
754 + $(MKOCT) -DDINT -c ../Source/umf_kernel.c |
|
755 + $(MV) -f ../Source/umf_kernel.o umf_od_kernel.o |
|
756 + $(MKOCT) -DDINT -c ../Source/umf_kernel_init.c |
|
757 + $(MV) -f ../Source/umf_kernel_init.o umf_od_kernel_init.o |
|
758 + $(MKOCT) -DDINT -c ../Source/umf_kernel_wrapup.c |
|
759 + $(MV) -f ../Source/umf_kernel_wrapup.o umf_od_kernel_wrapup.o |
|
760 + $(MKOCT) -DDINT -c ../Source/umf_local_search.c |
|
761 + $(MV) -f ../Source/umf_local_search.o umf_od_local_search.o |
|
762 + $(MKOCT) -DDINT -c ../Source/umf_lsolve.c |
|
763 + $(MV) -f ../Source/umf_lsolve.o umf_od_lsolve.o |
|
764 + $(MKOCT) -DDINT -c ../Source/umf_ltsolve.c |
|
765 + $(MV) -f ../Source/umf_ltsolve.o umf_od_ltsolve.o |
|
766 + $(MKOCT) -DDINT -c ../Source/umf_mem_alloc_element.c |
|
767 + $(MV) -f ../Source/umf_mem_alloc_element.o umf_od_mem_alloc_element.o |
|
768 + $(MKOCT) -DDINT -c ../Source/umf_mem_alloc_head_block.c |
|
769 + $(MV) -f ../Source/umf_mem_alloc_head_block.o umf_od_mem_alloc_head_block.o |
|
770 + $(MKOCT) -DDINT -c ../Source/umf_mem_alloc_tail_block.c |
|
771 + $(MV) -f ../Source/umf_mem_alloc_tail_block.o umf_od_mem_alloc_tail_block.o |
|
772 + $(MKOCT) -DDINT -c ../Source/umf_mem_free_tail_block.c |
|
773 + $(MV) -f ../Source/umf_mem_free_tail_block.o umf_od_mem_free_tail_block.o |
|
774 + $(MKOCT) -DDINT -c ../Source/umf_mem_init_memoryspace.c |
|
775 + $(MV) -f ../Source/umf_mem_init_memoryspace.o umf_od_mem_init_memoryspace.o |
|
776 + $(MKOCT) -DDINT -c ../Source/umf_report_vector.c |
|
777 + $(MV) -f ../Source/umf_report_vector.o umf_od_report_vector.o |
|
778 + $(MKOCT) -DDINT -c ../Source/umf_row_search.c |
|
779 + $(MV) -f ../Source/umf_row_search.o umf_od_row_search.o |
|
780 + $(MKOCT) -DDINT -c ../Source/umf_scale_column.c |
|
781 + $(MV) -f ../Source/umf_scale_column.o umf_od_scale_column.o |
|
782 + $(MKOCT) -DDINT -c ../Source/umf_set_stats.c |
|
783 + $(MV) -f ../Source/umf_set_stats.o umf_od_set_stats.o |
|
784 + $(MKOCT) -DDINT -c ../Source/umf_solve.c |
|
785 + $(MV) -f ../Source/umf_solve.o umf_od_solve.o |
|
786 + $(MKOCT) -DDINT -c ../Source/umf_symbolic_usage.c |
|
787 + $(MV) -f ../Source/umf_symbolic_usage.o umf_od_symbolic_usage.o |
|
788 + $(MKOCT) -DDINT -c ../Source/umf_transpose.c |
|
789 + $(MV) -f ../Source/umf_transpose.o umf_od_transpose.o |
|
790 + $(MKOCT) -DDINT -c ../Source/umf_tuple_lengths.c |
|
791 + $(MV) -f ../Source/umf_tuple_lengths.o umf_od_tuple_lengths.o |
|
792 + $(MKOCT) -DDINT -c ../Source/umf_usolve.c |
|
793 + $(MV) -f ../Source/umf_usolve.o umf_od_usolve.o |
|
794 + $(MKOCT) -DDINT -c ../Source/umf_utsolve.c |
|
795 + $(MV) -f ../Source/umf_utsolve.o umf_od_utsolve.o |
|
796 + $(MKOCT) -DDINT -c ../Source/umf_valid_numeric.c |
|
797 + $(MV) -f ../Source/umf_valid_numeric.o umf_od_valid_numeric.o |
|
798 + $(MKOCT) -DDINT -c ../Source/umf_valid_symbolic.c |
|
799 + $(MV) -f ../Source/umf_valid_symbolic.o umf_od_valid_symbolic.o |
|
800 + $(MKOCT) -DDINT -c ../Source/umf_grow_front.c |
|
801 + $(MV) -f ../Source/umf_grow_front.o umf_od_grow_front.o |
|
802 + $(MKOCT) -DDINT -c ../Source/umf_start_front.c |
|
803 + $(MV) -f ../Source/umf_start_front.o umf_od_start_front.o |
|
804 + $(MKOCT) -DDINT -c ../Source/umf_2by2.c |
|
805 + $(MV) -f ../Source/umf_2by2.o umf_od_2by2.o |
|
806 + $(MKOCT) -DDINT -c ../Source/umf_store_lu.c |
|
807 + $(MV) -f ../Source/umf_store_lu.o umf_od_store_lu.o |
|
808 + $(MKOCT) -DDINT -c ../Source/umf_scale.c |
|
809 + $(MV) -f ../Source/umf_scale.o umf_od_scale.o |
|
810 + $(MKOCT) -DDINT -DWSOLVE -c ../Source/umfpack_solve.c |
|
811 + $(MV) -f ../Source/umfpack_solve.o umfpack_od_wsolve.o |
|
812 + $(MKOCT) -DDINT -c ../Source/umfpack_col_to_triplet.c |
|
813 + $(MV) -f ../Source/umfpack_col_to_triplet.o umfpack_od_col_to_triplet.o |
|
814 + $(MKOCT) -DDINT -c ../Source/umfpack_defaults.c |
|
815 + $(MV) -f ../Source/umfpack_defaults.o umfpack_od_defaults.o |
|
816 + $(MKOCT) -DDINT -c ../Source/umfpack_free_numeric.c |
|
817 + $(MV) -f ../Source/umfpack_free_numeric.o umfpack_od_free_numeric.o |
|
818 + $(MKOCT) -DDINT -c ../Source/umfpack_free_symbolic.c |
|
819 + $(MV) -f ../Source/umfpack_free_symbolic.o umfpack_od_free_symbolic.o |
|
820 + $(MKOCT) -DDINT -c ../Source/umfpack_get_numeric.c |
|
821 + $(MV) -f ../Source/umfpack_get_numeric.o umfpack_od_get_numeric.o |
|
822 + $(MKOCT) -DDINT -c ../Source/umfpack_get_lunz.c |
|
823 + $(MV) -f ../Source/umfpack_get_lunz.o umfpack_od_get_lunz.o |
|
824 + $(MKOCT) -DDINT -c ../Source/umfpack_get_symbolic.c |
|
825 + $(MV) -f ../Source/umfpack_get_symbolic.o umfpack_od_get_symbolic.o |
|
826 + $(MKOCT) -DDINT -c ../Source/umfpack_get_determinant.c |
|
827 + $(MV) -f ../Source/umfpack_get_determinant.o umfpack_od_get_determinant.o |
|
828 + $(MKOCT) -DDINT -c ../Source/umfpack_numeric.c |
|
829 + $(MV) -f ../Source/umfpack_numeric.o umfpack_od_numeric.o |
|
830 + $(MKOCT) -DDINT -c ../Source/umfpack_qsymbolic.c |
|
831 + $(MV) -f ../Source/umfpack_qsymbolic.o umfpack_od_qsymbolic.o |
|
832 + $(MKOCT) -DDINT -c ../Source/umfpack_report_control.c |
|
833 + $(MV) -f ../Source/umfpack_report_control.o umfpack_od_report_control.o |
|
834 + $(MKOCT) -DDINT -c ../Source/umfpack_report_info.c |
|
835 + $(MV) -f ../Source/umfpack_report_info.o umfpack_od_report_info.o |
|
836 + $(MKOCT) -DDINT -c ../Source/umfpack_report_matrix.c |
|
837 + $(MV) -f ../Source/umfpack_report_matrix.o umfpack_od_report_matrix.o |
|
838 + $(MKOCT) -DDINT -c ../Source/umfpack_report_numeric.c |
|
839 + $(MV) -f ../Source/umfpack_report_numeric.o umfpack_od_report_numeric.o |
|
840 + $(MKOCT) -DDINT -c ../Source/umfpack_report_perm.c |
|
841 + $(MV) -f ../Source/umfpack_report_perm.o umfpack_od_report_perm.o |
|
842 + $(MKOCT) -DDINT -c ../Source/umfpack_report_status.c |
|
843 + $(MV) -f ../Source/umfpack_report_status.o umfpack_od_report_status.o |
|
844 + $(MKOCT) -DDINT -c ../Source/umfpack_report_symbolic.c |
|
845 + $(MV) -f ../Source/umfpack_report_symbolic.o umfpack_od_report_symbolic.o |
|
846 + $(MKOCT) -DDINT -c ../Source/umfpack_report_triplet.c |
|
847 + $(MV) -f ../Source/umfpack_report_triplet.o umfpack_od_report_triplet.o |
|
848 + $(MKOCT) -DDINT -c ../Source/umfpack_report_vector.c |
|
849 + $(MV) -f ../Source/umfpack_report_vector.o umfpack_od_report_vector.o |
|
850 + $(MKOCT) -DDINT -c ../Source/umfpack_solve.c |
|
851 + $(MV) -f ../Source/umfpack_solve.o umfpack_od_solve.o |
|
852 + $(MKOCT) -DDINT -c ../Source/umfpack_symbolic.c |
|
853 + $(MV) -f ../Source/umfpack_symbolic.o umfpack_od_symbolic.o |
|
854 + $(MKOCT) -DDINT -c ../Source/umfpack_transpose.c |
|
855 + $(MV) -f ../Source/umfpack_transpose.o umfpack_od_transpose.o |
|
856 + $(MKOCT) -DDINT -c ../Source/umfpack_triplet_to_col.c |
|
857 + $(MV) -f ../Source/umfpack_triplet_to_col.o umfpack_od_triplet_to_col.o |
|
858 + $(MKOCT) -DDINT -c ../Source/umfpack_scale.c |
|
859 + $(MV) -f ../Source/umfpack_scale.o umfpack_od_scale.o |
|
860 + $(MKOCT) -DDINT -c ../Source/umfpack_load_numeric.c |
|
861 + $(MV) -f ../Source/umfpack_load_numeric.o umfpack_od_load_numeric.o |
|
862 + $(MKOCT) -DDINT -c ../Source/umfpack_save_numeric.c |
|
863 + $(MV) -f ../Source/umfpack_save_numeric.o umfpack_od_save_numeric.o |
|
864 + $(MKOCT) -DDINT -c ../Source/umfpack_load_symbolic.c |
|
865 + $(MV) -f ../Source/umfpack_load_symbolic.o umfpack_od_load_symbolic.o |
|
866 + $(MKOCT) -DDINT -c ../Source/umfpack_save_symbolic.c |
|
867 + $(MV) -f ../Source/umfpack_save_symbolic.o umfpack_od_save_symbolic.o |
|
868 + $(MKOCT) -DZINT -DCONJUGATE_SOLVE -c ../Source/umf_ltsolve.c |
|
869 + $(MV) -f ../Source/umf_ltsolve.o umf_oz_lhsolve.o |
|
870 + $(MKOCT) -DZINT -DCONJUGATE_SOLVE -c ../Source/umf_utsolve.c |
|
871 + $(MV) -f ../Source/umf_utsolve.o umf_oz_uhsolve.o |
|
872 + $(MKOCT) -DZINT -DDO_MAP -c ../Source/umf_triplet.c |
|
873 + $(MV) -f ../Source/umf_triplet.o umf_oz_triplet_map_nox.o |
|
874 + $(MKOCT) -DZINT -DDO_VALUES -c ../Source/umf_triplet.c |
|
875 + $(MV) -f ../Source/umf_triplet.o umf_oz_triplet_nomap_x.o |
|
876 + $(MKOCT) -DZINT -c ../Source/umf_triplet.c |
|
877 + $(MV) -f ../Source/umf_triplet.o umf_oz_triplet_nomap_nox.o |
|
878 + $(MKOCT) -DZINT -DDO_MAP -DDO_VALUES -c ../Source/umf_triplet.c |
|
879 + $(MV) -f ../Source/umf_triplet.o umf_oz_triplet_map_x.o |
|
880 + $(MKOCT) -DZINT -DFIXQ -c ../Source/umf_assemble.c |
|
881 + $(MV) -f ../Source/umf_assemble.o umf_oz_assemble_fixq.o |
|
882 + $(MKOCT) -DZINT -DDROP -c ../Source/umf_store_lu.c |
|
883 + $(MV) -f ../Source/umf_store_lu.o umf_oz_store_lu_drop.o |
|
884 + $(MKOCT) -DZINT -c ../Source/umf_assemble.c |
|
885 + $(MV) -f ../Source/umf_assemble.o umf_oz_assemble.o |
|
886 + $(MKOCT) -DZINT -c ../Source/umf_blas3_update.c |
|
887 + $(MV) -f ../Source/umf_blas3_update.o umf_oz_blas3_update.o |
|
888 + $(MKOCT) -DZINT -c ../Source/umf_build_tuples.c |
|
889 + $(MV) -f ../Source/umf_build_tuples.o umf_oz_build_tuples.o |
|
890 + $(MKOCT) -DZINT -c ../Source/umf_create_element.c |
|
891 + $(MV) -f ../Source/umf_create_element.o umf_oz_create_element.o |
|
892 + $(MKOCT) -DZINT -c ../Source/umf_dump.c |
|
893 + $(MV) -f ../Source/umf_dump.o umf_oz_dump.o |
|
894 + $(MKOCT) -DZINT -c ../Source/umf_extend_front.c |
|
895 + $(MV) -f ../Source/umf_extend_front.o umf_oz_extend_front.o |
|
896 + $(MKOCT) -DZINT -c ../Source/umf_garbage_collection.c |
|
897 + $(MV) -f ../Source/umf_garbage_collection.o umf_oz_garbage_collection.o |
|
898 + $(MKOCT) -DZINT -c ../Source/umf_get_memory.c |
|
899 + $(MV) -f ../Source/umf_get_memory.o umf_oz_get_memory.o |
|
900 + $(MKOCT) -DZINT -c ../Source/umf_init_front.c |
|
901 + $(MV) -f ../Source/umf_init_front.o umf_oz_init_front.o |
|
902 + $(MKOCT) -DZINT -c ../Source/umf_kernel.c |
|
903 + $(MV) -f ../Source/umf_kernel.o umf_oz_kernel.o |
|
904 + $(MKOCT) -DZINT -c ../Source/umf_kernel_init.c |
|
905 + $(MV) -f ../Source/umf_kernel_init.o umf_oz_kernel_init.o |
|
906 + $(MKOCT) -DZINT -c ../Source/umf_kernel_wrapup.c |
|
907 + $(MV) -f ../Source/umf_kernel_wrapup.o umf_oz_kernel_wrapup.o |
|
908 + $(MKOCT) -DZINT -c ../Source/umf_local_search.c |
|
909 + $(MV) -f ../Source/umf_local_search.o umf_oz_local_search.o |
|
910 + $(MKOCT) -DZINT -c ../Source/umf_lsolve.c |
|
911 + $(MV) -f ../Source/umf_lsolve.o umf_oz_lsolve.o |
|
912 + $(MKOCT) -DZINT -c ../Source/umf_ltsolve.c |
|
913 + $(MV) -f ../Source/umf_ltsolve.o umf_oz_ltsolve.o |
|
914 + $(MKOCT) -DZINT -c ../Source/umf_mem_alloc_element.c |
|
915 + $(MV) -f ../Source/umf_mem_alloc_element.o umf_oz_mem_alloc_element.o |
|
916 + $(MKOCT) -DZINT -c ../Source/umf_mem_alloc_head_block.c |
|
917 + $(MV) -f ../Source/umf_mem_alloc_head_block.o umf_oz_mem_alloc_head_block.o |
|
918 + $(MKOCT) -DZINT -c ../Source/umf_mem_alloc_tail_block.c |
|
919 + $(MV) -f ../Source/umf_mem_alloc_tail_block.o umf_oz_mem_alloc_tail_block.o |
|
920 + $(MKOCT) -DZINT -c ../Source/umf_mem_free_tail_block.c |
|
921 + $(MV) -f ../Source/umf_mem_free_tail_block.o umf_oz_mem_free_tail_block.o |
|
922 + $(MKOCT) -DZINT -c ../Source/umf_mem_init_memoryspace.c |
|
923 + $(MV) -f ../Source/umf_mem_init_memoryspace.o umf_oz_mem_init_memoryspace.o |
|
924 + $(MKOCT) -DZINT -c ../Source/umf_report_vector.c |
|
925 + $(MV) -f ../Source/umf_report_vector.o umf_oz_report_vector.o |
|
926 + $(MKOCT) -DZINT -c ../Source/umf_row_search.c |
|
927 + $(MV) -f ../Source/umf_row_search.o umf_oz_row_search.o |
|
928 + $(MKOCT) -DZINT -c ../Source/umf_scale_column.c |
|
929 + $(MV) -f ../Source/umf_scale_column.o umf_oz_scale_column.o |
|
930 + $(MKOCT) -DZINT -c ../Source/umf_set_stats.c |
|
931 + $(MV) -f ../Source/umf_set_stats.o umf_oz_set_stats.o |
|
932 + $(MKOCT) -DZINT -c ../Source/umf_solve.c |
|
933 + $(MV) -f ../Source/umf_solve.o umf_oz_solve.o |
|
934 + $(MKOCT) -DZINT -c ../Source/umf_symbolic_usage.c |
|
935 + $(MV) -f ../Source/umf_symbolic_usage.o umf_oz_symbolic_usage.o |
|
936 + $(MKOCT) -DZINT -c ../Source/umf_transpose.c |
|
937 + $(MV) -f ../Source/umf_transpose.o umf_oz_transpose.o |
|
938 + $(MKOCT) -DZINT -c ../Source/umf_tuple_lengths.c |
|
939 + $(MV) -f ../Source/umf_tuple_lengths.o umf_oz_tuple_lengths.o |
|
940 + $(MKOCT) -DZINT -c ../Source/umf_usolve.c |
|
941 + $(MV) -f ../Source/umf_usolve.o umf_oz_usolve.o |
|
942 + $(MKOCT) -DZINT -c ../Source/umf_utsolve.c |
|
943 + $(MV) -f ../Source/umf_utsolve.o umf_oz_utsolve.o |
|
944 + $(MKOCT) -DZINT -c ../Source/umf_valid_numeric.c |
|
945 + $(MV) -f ../Source/umf_valid_numeric.o umf_oz_valid_numeric.o |
|
946 + $(MKOCT) -DZINT -c ../Source/umf_valid_symbolic.c |
|
947 + $(MV) -f ../Source/umf_valid_symbolic.o umf_oz_valid_symbolic.o |
|
948 + $(MKOCT) -DZINT -c ../Source/umf_grow_front.c |
|
949 + $(MV) -f ../Source/umf_grow_front.o umf_oz_grow_front.o |
|
950 + $(MKOCT) -DZINT -c ../Source/umf_start_front.c |
|
951 + $(MV) -f ../Source/umf_start_front.o umf_oz_start_front.o |
|
952 + $(MKOCT) -DZINT -c ../Source/umf_2by2.c |
|
953 + $(MV) -f ../Source/umf_2by2.o umf_oz_2by2.o |
|
954 + $(MKOCT) -DZINT -c ../Source/umf_store_lu.c |
|
955 + $(MV) -f ../Source/umf_store_lu.o umf_oz_store_lu.o |
|
956 + $(MKOCT) -DZINT -c ../Source/umf_scale.c |
|
957 + $(MV) -f ../Source/umf_scale.o umf_oz_scale.o |
|
958 + $(MKOCT) -DZINT -DWSOLVE -c ../Source/umfpack_solve.c |
|
959 + $(MV) -f ../Source/umfpack_solve.o umfpack_oz_wsolve.o |
|
960 + $(MKOCT) -DZINT -c ../Source/umfpack_col_to_triplet.c |
|
961 + $(MV) -f ../Source/umfpack_col_to_triplet.o umfpack_oz_col_to_triplet.o |
|
962 + $(MKOCT) -DZINT -c ../Source/umfpack_defaults.c |
|
963 + $(MV) -f ../Source/umfpack_defaults.o umfpack_oz_defaults.o |
|
964 + $(MKOCT) -DZINT -c ../Source/umfpack_free_numeric.c |
|
965 + $(MV) -f ../Source/umfpack_free_numeric.o umfpack_oz_free_numeric.o |
|
966 + $(MKOCT) -DZINT -c ../Source/umfpack_free_symbolic.c |
|
967 + $(MV) -f ../Source/umfpack_free_symbolic.o umfpack_oz_free_symbolic.o |
|
968 + $(MKOCT) -DZINT -c ../Source/umfpack_get_numeric.c |
|
969 + $(MV) -f ../Source/umfpack_get_numeric.o umfpack_oz_get_numeric.o |
|
970 + $(MKOCT) -DZINT -c ../Source/umfpack_get_lunz.c |
|
971 + $(MV) -f ../Source/umfpack_get_lunz.o umfpack_oz_get_lunz.o |
|
972 + $(MKOCT) -DZINT -c ../Source/umfpack_get_symbolic.c |
|
973 + $(MV) -f ../Source/umfpack_get_symbolic.o umfpack_oz_get_symbolic.o |
|
974 + $(MKOCT) -DZINT -c ../Source/umfpack_get_determinant.c |
|
975 + $(MV) -f ../Source/umfpack_get_determinant.o umfpack_oz_get_determinant.o |
|
976 + $(MKOCT) -DZINT -c ../Source/umfpack_numeric.c |
|
977 + $(MV) -f ../Source/umfpack_numeric.o umfpack_oz_numeric.o |
|
978 + $(MKOCT) -DZINT -c ../Source/umfpack_qsymbolic.c |
|
979 + $(MV) -f ../Source/umfpack_qsymbolic.o umfpack_oz_qsymbolic.o |
|
980 + $(MKOCT) -DZINT -c ../Source/umfpack_report_control.c |
|
981 + $(MV) -f ../Source/umfpack_report_control.o umfpack_oz_report_control.o |
|
982 + $(MKOCT) -DZINT -c ../Source/umfpack_report_info.c |
|
983 + $(MV) -f ../Source/umfpack_report_info.o umfpack_oz_report_info.o |
|
984 + $(MKOCT) -DZINT -c ../Source/umfpack_report_matrix.c |
|
985 + $(MV) -f ../Source/umfpack_report_matrix.o umfpack_oz_report_matrix.o |
|
986 + $(MKOCT) -DZINT -c ../Source/umfpack_report_numeric.c |
|
987 + $(MV) -f ../Source/umfpack_report_numeric.o umfpack_oz_report_numeric.o |
|
988 + $(MKOCT) -DZINT -c ../Source/umfpack_report_perm.c |
|
989 + $(MV) -f ../Source/umfpack_report_perm.o umfpack_oz_report_perm.o |
|
990 + $(MKOCT) -DZINT -c ../Source/umfpack_report_status.c |
|
991 + $(MV) -f ../Source/umfpack_report_status.o umfpack_oz_report_status.o |
|
992 + $(MKOCT) -DZINT -c ../Source/umfpack_report_symbolic.c |
|
993 + $(MV) -f ../Source/umfpack_report_symbolic.o umfpack_oz_report_symbolic.o |
|
994 + $(MKOCT) -DZINT -c ../Source/umfpack_report_triplet.c |
|
995 + $(MV) -f ../Source/umfpack_report_triplet.o umfpack_oz_report_triplet.o |
|
996 + $(MKOCT) -DZINT -c ../Source/umfpack_report_vector.c |
|
997 + $(MV) -f ../Source/umfpack_report_vector.o umfpack_oz_report_vector.o |
|
998 + $(MKOCT) -DZINT -c ../Source/umfpack_solve.c |
|
999 + $(MV) -f ../Source/umfpack_solve.o umfpack_oz_solve.o |
|
1000 + $(MKOCT) -DZINT -c ../Source/umfpack_symbolic.c |
|
1001 + $(MV) -f ../Source/umfpack_symbolic.o umfpack_oz_symbolic.o |
|
1002 + $(MKOCT) -DZINT -c ../Source/umfpack_transpose.c |
|
1003 + $(MV) -f ../Source/umfpack_transpose.o umfpack_oz_transpose.o |
|
1004 + $(MKOCT) -DZINT -c ../Source/umfpack_triplet_to_col.c |
|
1005 + $(MV) -f ../Source/umfpack_triplet_to_col.o umfpack_oz_triplet_to_col.o |
|
1006 + $(MKOCT) -DZINT -c ../Source/umfpack_scale.c |
|
1007 + $(MV) -f ../Source/umfpack_scale.o umfpack_oz_scale.o |
|
1008 + $(MKOCT) -DZINT -c ../Source/umfpack_load_numeric.c |
|
1009 + $(MV) -f ../Source/umfpack_load_numeric.o umfpack_oz_load_numeric.o |
|
1010 + $(MKOCT) -DZINT -c ../Source/umfpack_save_numeric.c |
|
1011 + $(MV) -f ../Source/umfpack_save_numeric.o umfpack_oz_save_numeric.o |
|
1012 + $(MKOCT) -DZINT -c ../Source/umfpack_load_symbolic.c |
|
1013 + $(MV) -f ../Source/umfpack_load_symbolic.o umfpack_oz_load_symbolic.o |
|
1014 + $(MKOCT) -DZINT -c ../Source/umfpack_save_symbolic.c |
|
1015 + $(MV) -f ../Source/umfpack_save_symbolic.o umfpack_oz_save_symbolic.o |
|
1016 + $(MKOCT) -c ../Source/umfpack_timer.c |
|
1017 + $(MV) -f ../Source/umfpack_timer.o umfpack_m_timer.o |
|
1018 + $(MKOCT) -c ../Source/umfpack_tictoc.c |
|
1019 + $(MV) -f ../Source/umfpack_tictoc.o umfpack_m_tictoc.o |
|
1020 + $(MKOCT) -DDINT -c ../../AMD/Source/amd_aat.c |
|
1021 + $(MV) -f ../../AMD/Source/amd_aat.o amd_m_aat.o |
|
1022 + $(MKOCT) -DDINT -c ../../AMD/Source/amd_1.c |
|
1023 + $(MV) -f ../../AMD/Source/amd_1.o amd_m_1.o |
|
1024 + $(MKOCT) -DDINT -c ../../AMD/Source/amd_2.c |
|
1025 + $(MV) -f ../../AMD/Source/amd_2.o amd_m_2.o |
|
1026 + $(MKOCT) -DDINT -c ../../AMD/Source/amd_dump.c |
|
1027 + $(MV) -f ../../AMD/Source/amd_dump.o amd_m_dump.o |
|
1028 + $(MKOCT) -DDINT -c ../../AMD/Source/amd_postorder.c |
|
1029 + $(MV) -f ../../AMD/Source/amd_postorder.o amd_m_postorder.o |
|
1030 + $(MKOCT) -DDINT -c ../../AMD/Source/amd_post_tree.c |
|
1031 + $(MV) -f ../../AMD/Source/amd_post_tree.o amd_m_post_tree.o |
|
1032 + $(MKOCT) -DDINT -c ../../AMD/Source/amd_defaults.c |
|
1033 + $(MV) -f ../../AMD/Source/amd_defaults.o amd_m_defaults.o |
|
1034 + $(MKOCT) -DDINT -c ../../AMD/Source/amd_order.c |
|
1035 + $(MV) -f ../../AMD/Source/amd_order.o amd_m_order.o |
|
1036 + $(MKOCT) -DDINT -c ../../AMD/Source/amd_control.c |
|
1037 + $(MV) -f ../../AMD/Source/amd_control.o amd_m_control.o |
|
1038 + $(MKOCT) -DDINT -c ../../AMD/Source/amd_info.c |
|
1039 + $(MV) -f ../../AMD/Source/amd_info.o amd_m_info.o |
|
1040 + $(MKOCT) -DDINT -c ../../AMD/Source/amd_valid.c |
|
1041 + $(MV) -f ../../AMD/Source/amd_valid.o amd_m_valid.o |
|
1042 + $(MKOCT) -o umfpack.oct $(OCT_SPARSE_INC) umfpack.cc \ |
|
1043 + umf_o_analyze.o umf_o_apply_order.o umf_o_colamd.o umf_o_free.o \ |
|
1044 + umf_o_fsize.o umf_o_is_permutation.o umf_o_malloc.o \ |
|
1045 + umf_o_realloc.o umf_o_report_perm.o umf_o_singletons.o \ |
|
1046 + umf_od_lhsolve.o umf_od_uhsolve.o umf_od_triplet_map_nox.o \ |
|
1047 + umf_od_triplet_nomap_x.o umf_od_triplet_nomap_nox.o \ |
|
1048 + umf_od_triplet_map_x.o umf_od_assemble_fixq.o \ |
|
1049 + umf_od_store_lu_drop.o umf_od_assemble.o umf_od_blas3_update.o \ |
|
1050 + umf_od_build_tuples.o umf_od_create_element.o umf_od_dump.o \ |
|
1051 + umf_od_extend_front.o umf_od_garbage_collection.o \ |
|
1052 + umf_od_get_memory.o umf_od_init_front.o umf_od_kernel.o \ |
|
1053 + umf_od_kernel_init.o umf_od_kernel_wrapup.o umf_od_local_search.o \ |
|
1054 + umf_od_lsolve.o umf_od_ltsolve.o umf_od_mem_alloc_element.o \ |
|
1055 + umf_od_mem_alloc_head_block.o umf_od_mem_alloc_tail_block.o \ |
|
1056 + umf_od_mem_free_tail_block.o umf_od_mem_init_memoryspace.o \ |
|
1057 + umf_od_report_vector.o umf_od_row_search.o umf_od_scale_column.o \ |
|
1058 + umf_od_set_stats.o umf_od_solve.o umf_od_symbolic_usage.o \ |
|
1059 + umf_od_transpose.o umf_od_tuple_lengths.o umf_od_usolve.o \ |
|
1060 + umf_od_utsolve.o umf_od_valid_numeric.o umf_od_valid_symbolic.o \ |
|
1061 + umf_od_grow_front.o umf_od_start_front.o umf_od_2by2.o \ |
|
1062 + umf_od_store_lu.o umf_od_scale.o umfpack_od_wsolve.o \ |
|
1063 + umfpack_od_col_to_triplet.o umfpack_od_defaults.o \ |
|
1064 + umfpack_od_free_numeric.o umfpack_od_free_symbolic.o \ |
|
1065 + umfpack_od_get_numeric.o umfpack_od_get_lunz.o \ |
|
1066 + umfpack_od_get_symbolic.o umfpack_od_numeric.o \ |
|
1067 + umfpack_od_qsymbolic.o umfpack_od_report_control.o \ |
|
1068 + umfpack_od_report_info.o umfpack_od_report_matrix.o \ |
|
1069 + umfpack_od_report_numeric.o umfpack_od_report_perm.o \ |
|
1070 + umfpack_od_report_status.o umfpack_od_report_symbolic.o \ |
|
1071 + umfpack_od_report_triplet.o umfpack_od_report_vector.o \ |
|
1072 + umfpack_od_solve.o umfpack_od_symbolic.o umfpack_od_transpose.o \ |
|
1073 + umfpack_od_triplet_to_col.o umfpack_od_scale.o \ |
|
1074 + umfpack_od_load_numeric.o umfpack_od_save_numeric.o \ |
|
1075 + umfpack_od_load_symbolic.o umfpack_od_save_symbolic.o \ |
|
1076 + umf_oz_lhsolve.o umf_oz_uhsolve.o umf_oz_triplet_map_nox.o \ |
|
1077 + umf_oz_triplet_nomap_x.o umf_oz_triplet_nomap_nox.o \ |
|
1078 + umf_oz_triplet_map_x.o umf_oz_assemble_fixq.o \ |
|
1079 + umf_oz_store_lu_drop.o umf_oz_assemble.o umf_oz_blas3_update.o \ |
|
1080 + umf_oz_build_tuples.o umf_oz_create_element.o umf_oz_dump.o \ |
|
1081 + umf_oz_extend_front.o umf_oz_garbage_collection.o \ |
|
1082 + umf_oz_get_memory.o umf_oz_init_front.o umf_oz_kernel.o \ |
|
1083 + umf_oz_kernel_init.o umf_oz_kernel_wrapup.o umf_oz_local_search.o \ |
|
1084 + umf_oz_lsolve.o umf_oz_ltsolve.o umf_oz_mem_alloc_element.o \ |
|
1085 + umf_oz_mem_alloc_head_block.o umf_oz_mem_alloc_tail_block.o \ |
|
1086 + umf_oz_mem_free_tail_block.o umf_oz_mem_init_memoryspace.o \ |
|
1087 + umf_oz_report_vector.o umf_oz_row_search.o umf_oz_scale_column.o \ |
|
1088 + umf_oz_set_stats.o umf_oz_solve.o umf_oz_symbolic_usage.o \ |
|
1089 + umf_oz_transpose.o umf_oz_tuple_lengths.o umf_oz_usolve.o \ |
|
1090 + umf_oz_utsolve.o umf_oz_valid_numeric.o umf_oz_valid_symbolic.o \ |
|
1091 + umf_oz_grow_front.o umf_oz_start_front.o umf_oz_2by2.o \ |
|
1092 + umf_oz_store_lu.o umf_oz_scale.o umfpack_oz_wsolve.o \ |
|
1093 + umfpack_oz_col_to_triplet.o umfpack_oz_defaults.o \ |
|
1094 + umfpack_oz_free_numeric.o umfpack_oz_free_symbolic.o \ |
|
1095 + umfpack_oz_get_numeric.o umfpack_oz_get_lunz.o \ |
|
1096 + umfpack_oz_get_symbolic.o umfpack_oz_numeric.o \ |
|
1097 + umfpack_oz_qsymbolic.o umfpack_oz_report_control.o \ |
|
1098 + umfpack_oz_report_info.o umfpack_oz_report_matrix.o \ |
|
1099 + umfpack_oz_report_numeric.o umfpack_oz_report_perm.o \ |
|
1100 + umfpack_oz_report_status.o umfpack_oz_report_symbolic.o \ |
|
1101 + umfpack_oz_report_triplet.o umfpack_oz_report_vector.o \ |
|
1102 + umfpack_oz_solve.o umfpack_oz_symbolic.o umfpack_oz_transpose.o \ |
|
1103 + umfpack_oz_triplet_to_col.o umfpack_oz_scale.o \ |
|
1104 + umfpack_oz_load_numeric.o umfpack_oz_save_numeric.o \ |
|
1105 + umfpack_oz_load_symbolic.o umfpack_oz_save_symbolic.o \ |
|
1106 + umfpack_o_timer.o umfpack_o_tictoc.o \ |
|
1107 + amd_o_aat.o amd_o_1.o amd_o_2.o amd_o_dump.o \ |
|
1108 + amd_o_postorder.o amd_o_post_tree.o amd_o_defaults.o amd_o_order.o \ |
|
1109 + amd_o_control.o amd_o_info.o amd_o_valid.o |
|
1110 + |
|
1111 +luflop: luflop.cc |
|
1112 + $(MKOCT) luflop.cc -I$(OCT_SPARSE_INC) -o luflop.oct |
|
1113 + |
|
1114 +#---------------------------------------- |
|
1115 +# umfpack library to link with octave |
|
1116 +#---------------------------------------- |
|
1117 + |
|
1118 +octave: umfpack |
|
1119 + ld -r -o ../../../$(OCTUMFPACK_LIB) \ |
|
1120 + umf_o_analyze.o umf_o_apply_order.o umf_o_colamd.o umf_o_free.o \ |
|
1121 + umf_o_fsize.o umf_o_is_permutation.o umf_o_malloc.o \ |
|
1122 + umf_o_realloc.o umf_o_report_perm.o umf_o_singletons.o \ |
|
1123 + umf_od_lhsolve.o umf_od_uhsolve.o umf_od_triplet_map_nox.o \ |
|
1124 + umf_od_triplet_nomap_x.o umf_od_triplet_nomap_nox.o \ |
|
1125 + umf_od_triplet_map_x.o umf_od_assemble_fixq.o \ |
|
1126 + umf_od_store_lu_drop.o umf_od_assemble.o umf_od_blas3_update.o \ |
|
1127 + umf_od_build_tuples.o umf_od_create_element.o umf_od_dump.o \ |
|
1128 + umf_od_extend_front.o umf_od_garbage_collection.o \ |
|
1129 + umf_od_get_memory.o umf_od_init_front.o umf_od_kernel.o \ |
|
1130 + umf_od_kernel_init.o umf_od_kernel_wrapup.o umf_od_local_search.o \ |
|
1131 + umf_od_lsolve.o umf_od_ltsolve.o umf_od_mem_alloc_element.o \ |
|
1132 + umf_od_mem_alloc_head_block.o umf_od_mem_alloc_tail_block.o \ |
|
1133 + umf_od_mem_free_tail_block.o umf_od_mem_init_memoryspace.o \ |
|
1134 + umf_od_report_vector.o umf_od_row_search.o umf_od_scale_column.o \ |
|
1135 + umf_od_set_stats.o umf_od_solve.o umf_od_symbolic_usage.o \ |
|
1136 + umf_od_transpose.o umf_od_tuple_lengths.o umf_od_usolve.o \ |
|
1137 + umf_od_utsolve.o umf_od_valid_numeric.o umf_od_valid_symbolic.o \ |
|
1138 + umf_od_grow_front.o umf_od_start_front.o umf_od_2by2.o \ |
|
1139 + umf_od_store_lu.o umf_od_scale.o umfpack_od_wsolve.o \ |
|
1140 + umfpack_od_col_to_triplet.o umfpack_od_defaults.o \ |
|
1141 + umfpack_od_free_numeric.o umfpack_od_free_symbolic.o \ |
|
1142 + umfpack_od_get_numeric.o umfpack_od_get_lunz.o \ |
|
1143 + umfpack_od_get_symbolic.o umfpack_od_numeric.o \ |
|
1144 + umfpack_od_qsymbolic.o umfpack_od_report_control.o \ |
|
1145 + umfpack_od_report_info.o umfpack_od_report_matrix.o \ |
|
1146 + umfpack_od_report_numeric.o umfpack_od_report_perm.o \ |
|
1147 + umfpack_od_report_status.o umfpack_od_report_symbolic.o \ |
|
1148 + umfpack_od_report_triplet.o umfpack_od_report_vector.o \ |
|
1149 + umfpack_od_solve.o umfpack_od_symbolic.o umfpack_od_transpose.o \ |
|
1150 + umfpack_od_triplet_to_col.o umfpack_od_scale.o \ |
|
1151 + umfpack_od_load_numeric.o umfpack_od_save_numeric.o \ |
|
1152 + umfpack_od_load_symbolic.o umfpack_od_save_symbolic.o \ |
|
1153 + umf_oz_lhsolve.o umf_oz_uhsolve.o umf_oz_triplet_map_nox.o \ |
|
1154 + umf_oz_triplet_nomap_x.o umf_oz_triplet_nomap_nox.o \ |
|
1155 + umf_oz_triplet_map_x.o umf_oz_assemble_fixq.o \ |
|
1156 + umf_oz_store_lu_drop.o umf_oz_assemble.o umf_oz_blas3_update.o \ |
|
1157 + umf_oz_build_tuples.o umf_oz_create_element.o umf_oz_dump.o \ |
|
1158 + umf_oz_extend_front.o umf_oz_garbage_collection.o \ |
|
1159 + umf_oz_get_memory.o umf_oz_init_front.o umf_oz_kernel.o \ |
|
1160 + umf_oz_kernel_init.o umf_oz_kernel_wrapup.o umf_oz_local_search.o \ |
|
1161 + umf_oz_lsolve.o umf_oz_ltsolve.o umf_oz_mem_alloc_element.o \ |
|
1162 + umf_oz_mem_alloc_head_block.o umf_oz_mem_alloc_tail_block.o \ |
|
1163 + umf_oz_mem_free_tail_block.o umf_oz_mem_init_memoryspace.o \ |
|
1164 + umf_oz_report_vector.o umf_oz_row_search.o umf_oz_scale_column.o \ |
|
1165 + umf_oz_set_stats.o umf_oz_solve.o umf_oz_symbolic_usage.o \ |
|
1166 + umf_oz_transpose.o umf_oz_tuple_lengths.o umf_oz_usolve.o \ |
|
1167 + umf_oz_utsolve.o umf_oz_valid_numeric.o umf_oz_valid_symbolic.o \ |
|
1168 + umf_oz_grow_front.o umf_oz_start_front.o umf_oz_2by2.o \ |
|
1169 + umf_oz_store_lu.o umf_oz_scale.o umfpack_oz_wsolve.o \ |
|
1170 + umfpack_oz_col_to_triplet.o umfpack_oz_defaults.o \ |
|
1171 + umfpack_oz_free_numeric.o umfpack_oz_free_symbolic.o \ |
|
1172 + umfpack_oz_get_numeric.o umfpack_oz_get_lunz.o \ |
|
1173 + umfpack_oz_get_symbolic.o umfpack_oz_numeric.o \ |
|
1174 + umfpack_oz_qsymbolic.o umfpack_oz_report_control.o \ |
|
1175 + umfpack_oz_report_info.o umfpack_oz_report_matrix.o \ |
|
1176 + umfpack_oz_report_numeric.o umfpack_oz_report_perm.o \ |
|
1177 + umfpack_oz_report_status.o umfpack_oz_report_symbolic.o \ |
|
1178 + umfpack_oz_report_triplet.o umfpack_oz_report_vector.o \ |
|
1179 + umfpack_oz_solve.o umfpack_oz_symbolic.o umfpack_oz_transpose.o \ |
|
1180 + umfpack_oz_triplet_to_col.o umfpack_oz_scale.o \ |
|
1181 + umfpack_oz_load_numeric.o umfpack_oz_save_numeric.o \ |
|
1182 + umfpack_oz_load_symbolic.o umfpack_oz_save_symbolic.o \ |
|
1183 + umfpack_o_timer.o umfpack_o_tictoc.o \ |
|
1184 + amd_o_aat.o amd_o_1.o amd_o_2.o amd_o_dump.o \ |
|
1185 + amd_o_postorder.o amd_o_post_tree.o amd_o_defaults.o amd_o_order.o \ |
|
1186 + amd_o_control.o amd_o_info.o amd_o_valid.o |
|
1187 + |
|
1188 +#------------------------------------------------------------------------------- |
|
1189 +# Remove all but the files in the original distribution |
|
1190 +#------------------------------------------------------------------------------- |
|
1191 + |
|
1192 +purge: clean |
|
1193 + - $(RM) *.oct* *.dll |
|
1194 diff -uNr UMFPACKv4.4.orig/UMFPACK/OCTAVE/umfpack_btf.m UMFPACK/UMFPACK/OCTAVE/umfpack_btf.m |
|
1195 --- UMFPACKv4.4.orig/UMFPACK/OCTAVE/umfpack_btf.m 1970-01-01 01:00:00.000000000 +0100 |
|
1196 +++ UMFPACK/UMFPACK/OCTAVE/umfpack_btf.m 2004-12-30 01:58:46.000000000 +0100 |
|
1197 @@ -0,0 +1,129 @@ |
|
1198 +function x = umfpack_btf (A, b, Control) |
|
1199 +% UMFPACK_BTF |
|
1200 +% |
|
1201 +% x = umfpack_btf (A, b, Control) |
|
1202 +% |
|
1203 +% solve Ax=b by first permuting the matrix A to block triangular form via dmperm |
|
1204 +% and then using UMFPACK to factorize each diagonal block. Adjacent 1-by-1 |
|
1205 +% blocks are merged into a single upper triangular block, and solved via |
|
1206 +% MATLAB's \ operator. The Control parameter is optional (Type umfpack_details |
|
1207 +% and umfpack_report for details on its use). A must be square. |
|
1208 +% |
|
1209 +% See also: umfpack, umfpack_factorize, umfpack_details, dmperm |
|
1210 + |
|
1211 +% UMFPACK Version 4.3 (Jan. 16, 2004), Copyright (c) 2004 by Timothy A. |
|
1212 +% Davis. All Rights Reserved. Type umfpack_details for License. |
|
1213 + |
|
1214 +if (nargin < 2) |
|
1215 + help umfpack_btf |
|
1216 + error ('Usage: x = umfpack_btf (A, b, Control)') ; |
|
1217 +end |
|
1218 + |
|
1219 +[m n] = size (A) ; |
|
1220 +if (m ~= n) |
|
1221 + help umfpack_btf |
|
1222 + error ('umfpack_btf: A must be square') ; |
|
1223 +end |
|
1224 +[m1 n1] = size (b) ; |
|
1225 +if (m1 ~= n) |
|
1226 + help umfpack_btf |
|
1227 + error ('umfpack_btf: b has the wrong dimensions') ; |
|
1228 +end |
|
1229 + |
|
1230 +if (nargin < 3) |
|
1231 + Control = umfpack ; |
|
1232 +end |
|
1233 + |
|
1234 +%------------------------------------------------------------------------------- |
|
1235 +% find the block triangular form |
|
1236 +%------------------------------------------------------------------------------- |
|
1237 + |
|
1238 +[p,q,r] = dmperm (A) ; |
|
1239 +nblocks = length (r) - 1 ; |
|
1240 + |
|
1241 +%------------------------------------------------------------------------------- |
|
1242 +% solve the system |
|
1243 +%------------------------------------------------------------------------------- |
|
1244 + |
|
1245 +if (nblocks == 1 | sprank (A) < n) |
|
1246 + |
|
1247 + %--------------------------------------------------------------------------- |
|
1248 + % matrix is irreducible or structurally singular |
|
1249 + %--------------------------------------------------------------------------- |
|
1250 + |
|
1251 + x = umfpack_solve (A, '\\', b, Control) ; |
|
1252 + |
|
1253 +else |
|
1254 + |
|
1255 + %--------------------------------------------------------------------------- |
|
1256 + % A (p,q) is in block triangular form |
|
1257 + %--------------------------------------------------------------------------- |
|
1258 + |
|
1259 + b = b (p,:) ; |
|
1260 + A = A (p,q) ; |
|
1261 + x = zeros (size (b)) ; |
|
1262 + |
|
1263 + %--------------------------------------------------------------------------- |
|
1264 + % merge adjacent singletons into a single upper triangular block |
|
1265 + %--------------------------------------------------------------------------- |
|
1266 + |
|
1267 + [r, nblocks, is_triangular] = merge_singletons (r) ; |
|
1268 + |
|
1269 + %--------------------------------------------------------------------------- |
|
1270 + % solve the system: x (q) = A\b |
|
1271 + %--------------------------------------------------------------------------- |
|
1272 + |
|
1273 + for k = nblocks:-1:1 |
|
1274 + |
|
1275 + % get the kth block |
|
1276 + k1 = r (k) ; |
|
1277 + k2 = r (k+1) - 1 ; |
|
1278 + |
|
1279 + % solve the system |
|
1280 + x (k1:k2,:) = solver (A (k1:k2, k1:k2), b (k1:k2,:), ... |
|
1281 + is_triangular (k), Control) ; |
|
1282 + |
|
1283 + % off-diagonal block back substitution |
|
1284 + b (1:k1-1,:) = b (1:k1-1,:) - A (1:k1-1, k1:k2) * x (k1:k2,:) ; |
|
1285 + |
|
1286 + end |
|
1287 + |
|
1288 + x (q,:) = x ; |
|
1289 + |
|
1290 +end |
|
1291 + |
|
1292 +%------------------------------------------------------------------------------- |
|
1293 +% merge_singletons |
|
1294 +%------------------------------------------------------------------------------- |
|
1295 + |
|
1296 +function [r, nblocks, is_triangular] = merge_singletons (r) |
|
1297 +% |
|
1298 +% Given r from [p,q,r] = dmperm (A), where A is square, return a modified r that |
|
1299 +% reflects the merger of adjacent singletons into a single upper triangular |
|
1300 +% block. is_triangular (k) is 1 if the kth block is upper triangular. nblocks |
|
1301 +% is the number of new blocks. |
|
1302 + |
|
1303 +nblocks = length (r) - 1 ; |
|
1304 +bsize = r (2:nblocks+1) - r (1:nblocks) ; |
|
1305 +t = [0 (bsize == 1)] ; |
|
1306 +z = (t (1:nblocks) == 0 & t (2:nblocks+1) == 1) | t (2:nblocks+1) == 0 ; |
|
1307 +y = [(find (z)) nblocks+1] ; |
|
1308 +r = r (y) ; |
|
1309 +nblocks = length (y) - 1 ; |
|
1310 +is_triangular = y (2:nblocks+1) - y (1:nblocks) > 1 ; |
|
1311 + |
|
1312 +%------------------------------------------------------------------------------- |
|
1313 +% solve Ax=b, but check for small and/or triangular systems |
|
1314 +%------------------------------------------------------------------------------- |
|
1315 + |
|
1316 +function x = solver (A, b, is_triangular, Control) |
|
1317 +if (is_triangular) |
|
1318 + % back substitution only |
|
1319 + x = A \ b ; |
|
1320 +elseif (size (A,1) < 4) |
|
1321 + % a very small matrix, solve it as a dense linear system |
|
1322 + x = full (A) \ b ; |
|
1323 +else |
|
1324 + % solve it as a sparse linear system |
|
1325 + x = umfpack_solve (A, '\\', b, Control) ; |
|
1326 +end |
|
1327 diff -uNr UMFPACKv4.4.orig/UMFPACK/OCTAVE/umfpack.cc UMFPACK/UMFPACK/OCTAVE/umfpack.cc |
|
1328 --- UMFPACKv4.4.orig/UMFPACK/OCTAVE/umfpack.cc 1970-01-01 01:00:00.000000000 +0100 |
|
1329 +++ UMFPACK/UMFPACK/OCTAVE/umfpack.cc 2005-02-01 21:58:15.885225363 +0100 |
|
1330 @@ -0,0 +1,1399 @@ |
|
1331 +/* |
|
1332 + |
|
1333 +Copyright (C) 2004 David Bateman |
|
1334 + |
|
1335 +This program is free software; you can redistribute it and/or modify it |
|
1336 +under the terms of the GNU General Public License as published by the |
|
1337 +Free Software Foundation; either version 2, or (at your option) any |
|
1338 +later version. |
|
1339 + |
|
1340 +This program is distributed in the hope that it will be useful, but WITHOUT |
|
1341 +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
|
1342 +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
|
1343 +for more details. |
|
1344 + |
|
1345 +You should have received a copy of the GNU General Public License |
|
1346 +along with this program; see the file COPYING. If not, write to the Free |
|
1347 +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. |
|
1348 + |
|
1349 +In addition to the terms of the GPL, you are permitted to link |
|
1350 +this program with any Open Source program, as defined by the |
|
1351 +Open Source Initiative (www.opensource.org) |
|
1352 + |
|
1353 +*/ |
|
1354 + |
|
1355 +/* |
|
1356 + |
|
1357 +This is the Octave interface to the UMFPACK code, which bore the following |
|
1358 +copyright |
|
1359 + |
|
1360 + UMFPACK Version 4.3 (Jan. 16, 2004), Copyright (c) 2004 by Timothy A. |
|
1361 + Davis. All Rights Reserved. See ../README for License. |
|
1362 + email: davis@cise.ufl.edu CISE Department, Univ. of Florida. |
|
1363 + web: http://www.cise.ufl.edu/research/sparse/umfpack |
|
1364 + |
|
1365 +*/ |
|
1366 + |
|
1367 +#include <cfloat> |
|
1368 +#include <cstdlib> |
|
1369 +#include <string> |
|
1370 + |
|
1371 +#include <octave/config.h> |
|
1372 +#include <octave/ov.h> |
|
1373 +#include <octave/defun-dld.h> |
|
1374 +#include <octave/pager.h> |
|
1375 +#include <octave/ov-re-mat.h> |
|
1376 + |
|
1377 +#include "ov-re-sparse.h" |
|
1378 +#include "ov-cx-sparse.h" |
|
1379 + |
|
1380 +// External UMFPACK functions in C |
|
1381 +extern "C" { |
|
1382 +#include "umfpack.h" |
|
1383 +} |
|
1384 + |
|
1385 +#define MIN(a,b) (((a) < (b)) ? (a) : (b)) |
|
1386 +#define MAX(a,b) (((a) > (b)) ? (a) : (b)) |
|
1387 + |
|
1388 +// Return an error message |
|
1389 +static |
|
1390 +void umfpack_error (const char *s, int A_is_complex, int nargout, |
|
1391 + octave_value_list retval, NDArray Control, |
|
1392 + NDArray Info, int status, int do_info) |
|
1393 +{ |
|
1394 + if (A_is_complex) |
|
1395 + { |
|
1396 + umfpack_zi_report_status (Control.fortran_vec (), status) ; |
|
1397 + umfpack_zi_report_info (Control.fortran_vec (), Info.fortran_vec ()) ; |
|
1398 + } |
|
1399 + else |
|
1400 + { |
|
1401 + umfpack_di_report_status (Control.fortran_vec (), status) ; |
|
1402 + umfpack_di_report_info (Control.fortran_vec (), Info.fortran_vec ()) ; |
|
1403 + } |
|
1404 + if (do_info > 0) |
|
1405 + // return Info |
|
1406 + retval (do_info) = octave_value (Info); |
|
1407 + |
|
1408 + error (s); |
|
1409 +} |
|
1410 + |
|
1411 +DEFUN_DLD (umfpack, args, nargout, |
|
1412 + "-*- texinfo -*-\n\ |
|
1413 +@deftypefn {Loadable Function} {[@var{x}, @var{Info}] =} umfpack (@var{a}, '\\', @var{b})\n\ |
|
1414 +@deftypefnx {Loadable Function} {[@var{x}, @var{Info}] =} umfpack (@var{a}, '\\', @var{b}, @var{Control})\n\ |
|
1415 +@deftypefnx {Loadable Function} {[@var{x}, @var{Info}] =} umfpack (@var{a}, @var{Qinit}, '\\', @var{b}, @var{Control})\n\ |
|
1416 +@deftypefnx {Loadable Function} {[@var{x}, @var{Info}] =} umfpack (@var{a}, @var{Qinit}, '\\', b)\n\ |
|
1417 +@deftypefnx {Loadable Function} {[@var{x}, @var{Info}] =} umfpack (@var{b}, '/', A) ;\n\ |
|
1418 +@deftypefnx {Loadable Function} {[@var{x}, @var{Info}] =} umfpack (@var{b}, '/', @var{a}, @var{Control}) ;\n\ |
|
1419 +@deftypefnx {Loadable Function} {[@var{x}, @var{Info}] =} umfpack (@var{b}, '/', @var{a}, @var{Qinit}) ;\n\ |
|
1420 +@deftypefnx {Loadable Function} {[@var{x}, @var{Info}] =} umfpack (@var{b}, '/', @var{a}, @var{Qinit}, @var{Control}) ;\n\ |
|
1421 +@deftypefnx {Loadable Function} {[@var{l}, @var{u}, @var{p}, @var{q}, @var{r}, @var{Info}] =} umfpack (@var{a}) ;\n\ |
|
1422 +@deftypefnx {Loadable Function} {[@var{l}, @var{u}, @var{p}, @var{q}, @var{r}, @var{Info}] =} umfpack (@var{a}, @var{Control}) ;\n\ |
|
1423 +@deftypefnx {Loadable Function} {[@var{l}, @var{u}, @var{p}, @var{q}, @var{r}, @var{Info}] =} umfpack (@var{a}, @var{Qinit}) ;\n\ |
|
1424 +@deftypefnx {Loadable Function} {[@var{l}, @var{u}, @var{p}, @var{q}, @var{r}, @var{Info}] =} umfpack (@var{a}, @var{Qinit}, @var{Control}) ;\n\ |
|
1425 +@deftypefnx {Loadable Function} {[@var{l}, @var{u}, @var{p}, @var{q}] =} umfpack (@var{a}) ;\n\ |
|
1426 +@deftypefnx {Loadable Function} {[@var{l}, @var{u}, @var{p}, @var{q}] =} umfpack (@var{a}, @var{Control}) ;\n\ |
|
1427 +@deftypefnx {Loadable Function} {[@var{l}, @var{u}, @var{p}, @var{q}] =} umfpack (@var{a}, @var{Qinit}) ;\n\ |
|
1428 +@deftypefnx {Loadable Function} {[@var{l}, @var{u}, @var{p}, @var{q}] =} umfpack (@var{a}, @var{Qinit}, @var{Control}) ;\n\ |
|
1429 +@deftypefnx {Loadable Function} {[P1, Q1, Fr, Ch, Info] =} umfpack (@var{a}, 'symbolic') ;\n\ |
|
1430 +@deftypefnx {Loadable Function} {[P1, Q1, Fr, Ch, Info] =} umfpack (@var{a}, 'symbolic', @var{Control}) ;\n\ |
|
1431 +@deftypefnx {Loadable Function} {[P1, Q1, Fr, Ch, Info] =} umfpack (@var{a}, @var{Qinit}, 'symbolic') ;\n\ |
|
1432 +@deftypefnx {Loadable Function} {[P1, Q1, Fr, Ch, Info] =} umfpack (@var{a}, @var{Qinit}, 'symbolic', @var{Control});\n\ |
|
1433 +@deftypefnx {Loadable Function} {@var{Control} =} umfpack ;\n\ |
|
1434 +\n\ |
|
1435 +UMFPACK v4.3 is a Octave oct-file for solving sparse linear systems.\n\ |
|
1436 +\n\ |
|
1437 +@iftex\n\ |
|
1438 +@tex\n\ |
|
1439 +\\vskip 2ex\n\ |
|
1440 +\\hfil\\vbox{\n\ |
|
1441 +\\offinterlineskip\n\ |
|
1442 +\\tabskip=0pt\n\ |
|
1443 +\\halign{\n\ |
|
1444 +\\vrule height1.75ex depth1.25ex width 0.6pt #\\tabskip=1em &\n\ |
|
1445 +\\hfil #\\hfil &\\vrule # & \n\ |
|
1446 +\\hfil #\\hfil &\\vrule # width 0.6pt \\tabskip=0pt\\cr\n\ |
|
1447 +\\noalign{\\hrule height 0.6pt}\n\ |
|
1448 +& UMFPACK v4.3 && OCTAVE approximate equivalent &\\cr\n\ |
|
1449 +\\noalign{\\hrule} \n\ |
|
1450 +& x = umfpack (A, '\\', b) ; && x = A \\ b &\\cr\n\ |
|
1451 +& && &\\cr\n\ |
|
1452 +&x = umfpack (b, '/', A) ; && x = b / A &\\cr\n\ |
|
1453 +& && &\\cr\n\ |
|
1454 +&[L,U,P,Q] = umfpack (A) ; && [m,n] = size (A) ; &\\cr\n\ |
|
1455 +& && I = speye (n) ; &\\cr\n\ |
|
1456 +& && Q = I (:, colamd (A)) ; &\\cr\n\ |
|
1457 +& && [L,U,P] = lu (A*Q) ; &\\cr\n\ |
|
1458 +& && &\\cr\n\ |
|
1459 +&[L,U,P,Q,R] = umfpack (A) ; && [m,n] = size (A) ; &\\cr\n\ |
|
1460 +& && I = speye (n) ; &\\cr\n\ |
|
1461 +& && Q = I (:, colamd (A)) ; &\\cr\n\ |
|
1462 +& && r = full (sum (abs (A), 2)) ; &\\cr\n\ |
|
1463 +& && r (find (r == 0)) = 1 ; &\\cr\n\ |
|
1464 +& && R = spdiags (r, 0, m, m) ; &\\cr\n\ |
|
1465 +& && [L,U,P] = lu ((R\\A)*Q) ; &\\cr\n\ |
|
1466 +& && &\\cr\n\ |
|
1467 +&[P,Q,F,C] = umfpack (A, 'symbolic')&& [m,n] = size (A) ; &\\cr\n\ |
|
1468 +& && I = speye (n) ; &\\cr\n\ |
|
1469 +& && Q = I (:, colamd (A)) ; &\\cr\n\ |
|
1470 +& && [count,h,parent,post] = ... &\\cr\n\ |
|
1471 +& && symbfact (A*Q, 'col') ; &\\cr\n\ |
|
1472 +\\noalign{\\hrule height 0.6pt}\n\ |
|
1473 +}}\\hfil\n\ |
|
1474 +\\vskip 1ex\n\ |
|
1475 +@end tex\n\ |
|
1476 +@end iftex\n\ |
|
1477 +@ifinfo\n\ |
|
1478 +@multitable @columnfractions 0.43 .02 .43\n\ |
|
1479 +@item UMFPACK v4.3: @tab | \n\ |
|
1480 +@tab OCTAVE approx. equivalent\n\ |
|
1481 +@item ------------------------------- @tab | \n\ |
|
1482 +@tab --------------------------------\n\ |
|
1483 +@item x = umfpack (A, '\\', b) ; @tab | \n\ |
|
1484 +@tab x = A \\ b\n\ |
|
1485 +@item @tab | \n\ |
|
1486 +@tab\n\ |
|
1487 +@item x = umfpack (b, '/', A) ; @tab | \n\ |
|
1488 +@tab x = b / A\n\ |
|
1489 +@item @tab | \n\ |
|
1490 +@tab\n\ |
|
1491 +@item [L,U,P,Q] = umfpack (A) ; @tab | \n\ |
|
1492 +@tab [m,n] = size (A) ;\n\ |
|
1493 +@item @tab | \n\ |
|
1494 +@tab I = speye (n) ;\n\ |
|
1495 +@item @tab | \n\ |
|
1496 +@tab Q = I (:, colamd (A)) ;\n\ |
|
1497 +@item @tab | \n\ |
|
1498 +@tab [L,U,P] = lu (A*Q) ;\n\ |
|
1499 +@item @tab | \n\ |
|
1500 +@tab\n\ |
|
1501 +@item [L,U,P,Q,R] = umfpack (A) ; @tab | \n\ |
|
1502 +@tab [m,n] = size (A) ;\n\ |
|
1503 +@item @tab | \n\ |
|
1504 +@tab I = speye (n) ;\n\ |
|
1505 +@item @tab | \n\ |
|
1506 +@tab Q = I (:, colamd (A)) ;\n\ |
|
1507 +@item @tab | \n\ |
|
1508 +@tab r = full (sum (abs (A), 2)) ;\n\ |
|
1509 +@item @tab | \n\ |
|
1510 +@tab r (find (r == 0)) = 1 ;\n\ |
|
1511 +@item @tab | \n\ |
|
1512 +@tab R = spdiags (r, 0, m, m) ;\n\ |
|
1513 +@item @tab | \n\ |
|
1514 +@tab [L,U,P] = lu ((R\\A)*Q) ;\n\ |
|
1515 +@item @tab | \n\ |
|
1516 +@tab\n\ |
|
1517 +@item [P,Q,F,C] = umfpack (A, 'symbolic') @tab | \n\ |
|
1518 +@tab [m,n] = size (A) ; \n\ |
|
1519 +@item @tab | \n\ |
|
1520 +@tab I = speye (n) ;\n\ |
|
1521 +@item @tab | \n\ |
|
1522 +@tab Q = I (:, colamd (A)) ;\n\ |
|
1523 +@item @tab | \n\ |
|
1524 +@tab [count,h,parent,post] = ...\n\ |
|
1525 +@item @tab | \n\ |
|
1526 +@tab symbfact (A*Q, 'col') ;\n\ |
|
1527 +@end multitable\n\ |
|
1528 +@end ifinfo\n\ |
|
1529 +\n\ |
|
1530 +A must be sparse. It can be complex, singular, and/or rectangular.\n\ |
|
1531 +A must be square for '/' or '\\'. b must be a full real or complex\n\ |
|
1532 +vector. For @code{[@var{l}, @var{u}, @var{p}, @var{q}, @var{r}] =\n\ |
|
1533 +umfpack (@var{a})}, the factorization is @code{@var{l} * @var{u} =\n\ |
|
1534 +@var{p} * (@var{r} \\ @var{a}) * @var{q}}. If @var{a} has a mostly\n\ |
|
1535 +symmetric nonzero pattern, then replace @dfn{colamd} with @dfn{amd}\n\ |
|
1536 +in the OCTAVE-equivalent column in the table above.\n\ |
|
1537 +\n\ |
|
1538 +Factor or solve a sparse linear system, returning either the solution\n\ |
|
1539 +@var{x} to @code{@var{A} * @var{x} = @var{b}} or @code{@var{A}' * @var{x}'\n\ |
|
1540 += @var{b}'}, the factorization LU=PAQ, or LU=P(R\\A)Q. A must be sparse.\n\ |
|
1541 +For the solve, A must be square and b must be a dense n-by-1 vector. For LU\n\ |
|
1542 +factorization, A can be rectangular. In both cases, A and/or b can be real\n\ |
|
1543 +or complex.\n\ |
|
1544 +\n\ |
|
1545 +UMFPACK analyzes the matrix and selects one of three strategies to factorize\n\ |
|
1546 +the matrix. It first finds a set of k initial pivot entries of zero\n\ |
|
1547 +Markowitz cost. This forms the first k rows and columns of L and U. The\n\ |
|
1548 +remaining submatrix S is then analyzed, based on the symmetry of the nonzero\n\ |
|
1549 +pattern of the submatrix and the values on the diagaonal. The strategies\n\ |
|
1550 +include:\n\ |
|
1551 +\n\ |
|
1552 +@table @asis\n\ |
|
1553 +@item unsymmetric\n\ |
|
1554 +Use a COLAMD pre-ordering, a column elimination tree\n\ |
|
1555 +post-ordering, refine the column ordering during factorization,\n\ |
|
1556 +and make no effort at selecting pivots on the diagonal.\n\ |
|
1557 +@item 2-by-2\n\ |
|
1558 +Like the symmetric strategy (see below), except that local\n\ |
|
1559 +row permutations are first made to attempt to place large entries\n\ |
|
1560 +on the diagonal.\n\ |
|
1561 +@item symmetric\n\ |
|
1562 +Use an AMD pre-ordering on the matrix @code{@var{s} + @var{s}'}, an\n\ |
|
1563 +elimination tree post-ordering, do not refine the column ordering during\n\ |
|
1564 +factorization, and attempt to select pivots on the diagonal.\n\ |
|
1565 +@end table\n\ |
|
1566 +\n\ |
|
1567 +Each of the following uses of umfpack (except for 'Control = umfpack') is\n\ |
|
1568 +stand-alone. That is, no call to umfpack is required for any subsequent\n\ |
|
1569 +call. In each usage, the Info output argument is optional.\n\ |
|
1570 +\n\ |
|
1571 +Usage:\n\ |
|
1572 +\n\ |
|
1573 +[x, Info] = umfpack (A, '\\', b) ;\n\ |
|
1574 +[x, Info] = umfpack (A, '\\', b, Control) ;\n\ |
|
1575 +[x, Info] = umfpack (A, Qinit, '\\', b, Control) ;\n\ |
|
1576 +[x, Info] = umfpack (A, Qinit, '\\', b) ;\n\ |
|
1577 +\n\ |
|
1578 + Solves Ax=b (similar to x = A\\b in OCTAVE).\n\ |
|
1579 +\n\ |
|
1580 +[x, Info] = umfpack (b, '/', A) ;\n\ |
|
1581 +[x, Info] = umfpack (b, '/', A, Control) ;\n\ |
|
1582 +[x, Info] = umfpack (b, '/', A, Qinit) ;\n\ |
|
1583 +[x, Info] = umfpack (b, '/', A, Qinit, Control) ;\n\ |
|
1584 +\n\ |
|
1585 + Solves A'x'=b' (similar to x = b/A in OCTAVE).\n\ |
|
1586 +\n\ |
|
1587 +[L, U, P, Q, R, Info] = umfpack (A) ;\n\ |
|
1588 +[L, U, P, Q, R, Info] = umfpack (A, Control) ;\n\ |
|
1589 +[L, U, P, Q, R, Info] = umfpack (A, Qinit) ;\n\ |
|
1590 +[L, U, P, Q, R, Info] = umfpack (A, Qinit, Control) ;\n\ |
|
1591 +\n\ |
|
1592 + Returns the LU factorization of A. P and Q are returned as permutation\n\ |
|
1593 + matrices. R is a diagonal sparse matrix of scale factors for the rows\n\ |
|
1594 + of A, L is lower triangular, and U is upper triangular. The\n\ |
|
1595 + factorization is L*U = P*(R\\A)*Q. You can turn off scaling by setting\n\ |
|
1596 + Control (17) to zero (in which case R = speye (m)), or by using the\n\ |
|
1597 + following syntaxes (in which case Control (17) is ignored):\n\ |
|
1598 +\n\ |
|
1599 +[L, U, P, Q] = umfpack (A) ;\n\ |
|
1600 +[L, U, P, Q] = umfpack (A, Control) ;\n\ |
|
1601 +[L, U, P, Q] = umfpack (A, Qinit) ;\n\ |
|
1602 +[L, U, P, Q] = umfpack (A, Qinit, Control) ;\n\ |
|
1603 +\n\ |
|
1604 + Same as above, except that no row scaling is performed. The Info array\n\ |
|
1605 + is not returned, either.\n\ |
|
1606 +\n\ |
|
1607 +[P1, Q1, Fr, Ch, Info] = umfpack (A, 'symbolic') ;\n\ |
|
1608 +[P1, Q1, Fr, Ch, Info] = umfpack (A, 'symbolic', Control) ;\n\ |
|
1609 +[P1, Q1, Fr, Ch, Info] = umfpack (A, Qinit, 'symbolic') ;\n\ |
|
1610 +[P1, Q1, Fr, Ch, Info] = umfpack (A, Qinit, 'symbolic', Control);\n\ |
|
1611 +\n\ |
|
1612 + Performs only the fill-reducing column pre-ordering (including the\n\ |
|
1613 + elimination tree post-ordering) and symbolic factorization. Q1 is the\n\ |
|
1614 + initial column permutation (either from colamd, amd, or the input\n\ |
|
1615 + ordering Qinit), possibly followed by a column elimination tree post-\n\ |
|
1616 + ordering or a symmetric elimination tree post-ordering, depending on\n\ |
|
1617 + the strategy used.\n\ |
|
1618 +\n\ |
|
1619 + For the unsymmetric strategy, P1 is the row ordering induced by Q1\n\ |
|
1620 + (row-merge order). For the 2-by-2 strategy, P1 is the row ordering that\n\ |
|
1621 + places large entries on the diagonal of P1*A*Q1. For the symmetric\n\ |
|
1622 + strategy, P1 = Q1.\n\ |
|
1623 +\n\ |
|
1624 + Fr is a (nfr+1)-by-4 array containing information about each frontal\n\ |
|
1625 + matrix, where nfr <= n is the number of frontal matrices. Fr (:,1) is\n\ |
|
1626 + the number of pivot columns in each front, and Fr (:,2) is the parent\n\ |
|
1627 + of each front in the supercolumn elimination tree. Fr (k,2) is zero if\n\ |
|
1628 + k is a root. The first Fr (1,1) columns of P1*A*Q1 are the pivot\n\ |
|
1629 + columns for the first front, the next Fr (2,1) columns of P1*A*Q1\n\ |
|
1630 + are the pivot columns for the second front, and so on.\n\ |
|
1631 +\n\ |
|
1632 + For the unsymmetric strategy, Fr (:,3) is the row index of the first\n\ |
|
1633 + row in P1*A*Q1 whose leftmost nonzero entry is in a pivot column for\n\ |
|
1634 + the kth front. Fr (:,4) is the leftmost descendent of the kth front.\n\ |
|
1635 + Rows in the range Fr (Fr (k,4),3) to Fr (k+1,3)-1 form the entire set\n\ |
|
1636 + of candidate pivot rows for the kth front (some of these will typically\n\ |
|
1637 + have been selected as pivot rows of fronts Fr (k,3) to k-1, before the\n\ |
|
1638 + factorization reaches the kth front. If front k is a leaf node, then\n\ |
|
1639 + Fr (k,4) is k.\n\ |
|
1640 +\n\ |
|
1641 + Ch is a (nchains+1)-by-3 array containing information about each\n\ |
|
1642 + 'chain' (unifrontal sequence) of frontal matrices, and where\n\ |
|
1643 + nchains <= nfr is the number of chains. The ith chain consists of\n\ |
|
1644 + frontal matrices. Chain (i,1) to Chain (i+1,1)-1, and the largest\n\ |
|
1645 + front in chain i is Chain (i,2)-by-Chain (i,3).\n\ |
|
1646 +\n\ |
|
1647 + This use of umfpack is not required to factor or solve a linear system\n\ |
|
1648 + in OCTAVE. It analyzes the matrix A and provides information only.\n\ |
|
1649 + The OCTAVE statement @code{treeplot (Fr (:,2)')} plots the column\n\ |
|
1650 + elimination tree.\n\ |
|
1651 +\n\ |
|
1652 +Control = umfpack ;\n\ |
|
1653 +\n\ |
|
1654 + Returns a 20-by-1 vector of default parameter settings for umfpack.\n\ |
|
1655 +\n\ |
|
1656 +umfpack_report (Control, Info) ;\n\ |
|
1657 +\n\ |
|
1658 + Prints the current Control settings, and Info\n\ |
|
1659 +\n\ |
|
1660 +If present, Qinit is a user-supplied 1-by-n permutation vector. It is an\n\ |
|
1661 +initial fill-reducing column pre-ordering for A; if not present, then colamd\n\ |
|
1662 +or amd are used instead. If present, Control is a user-supplied 20-by-1\n\ |
|
1663 +array. Control and Info are optional; if Control is not present, defaults\n\ |
|
1664 +are used. If a Control entry is NaN, then the default is used for that entry.\n\ |
|
1665 +\n\ |
|
1666 +UMFPACK Version 4.3 (Jan. 16, 2004), Copyright @copyright{} 2004 by\n\ |
|
1667 +Timothy A. Davis. All Rights Reserved.\n\ |
|
1668 +\n\ |
|
1669 +UMFPACK License:\n\ |
|
1670 +\n\ |
|
1671 +@example\n\ |
|
1672 +Your use or distribution of UMFPACK or any modified version of\n\ |
|
1673 +UMFPACK implies that you agree to this License.\n\ |
|
1674 +\n\ |
|
1675 +THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY\n\ |
|
1676 +EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK.\n\ |
|
1677 +\n\ |
|
1678 +Permission is hereby granted to use or copy this program, provided\n\ |
|
1679 +that the Copyright, this License, and the Availability of the original\n\ |
|
1680 +version is retained on all copies. User documentation of any code that\n\ |
|
1681 +uses UMFPACK or any modified version of UMFPACK code must cite the\n\ |
|
1682 +Copyright, this License, the Availability note, and 'Used by permission.'\n\ |
|
1683 +Permission to modify the code and to distribute modified code is granted,\n\ |
|
1684 +provided the Copyright, this License, and the Availability note are\n\ |
|
1685 +retained, and a notice that the code was modified is included. This\n\ |
|
1686 +software was developed with support from the National Science Foundation,\n\ |
|
1687 +and is provided to you free of charge.\n\ |
|
1688 +@end example\n\ |
|
1689 +\n\ |
|
1690 +Availability: http://www.cise.ufl.edu/research/sparse/umfpack\n\ |
|
1691 +@end deftypefn\n\ |
|
1692 +@seealso{lu_normtest, colamd, amd, umfpack_solve}") |
|
1693 +{ |
|
1694 + octave_value_list retval; |
|
1695 + int nargin = args.length (); |
|
1696 + int op = 0; |
|
1697 + std::string operation; |
|
1698 + bool do_solve = false; |
|
1699 + int do_info = 0; |
|
1700 + |
|
1701 + ColumnVector User_Qinit; |
|
1702 + SparseComplexMatrix CAmatrix; |
|
1703 + ComplexMatrix CBmatrix; |
|
1704 + SparseMatrix Amatrix; |
|
1705 + Matrix Bmatrix; |
|
1706 + NDArray User_Control_matrix; |
|
1707 + |
|
1708 + bool A_is_complex = false; |
|
1709 + bool B_is_complex = false; |
|
1710 + bool X_is_complex = false; |
|
1711 + bool transpose = false; |
|
1712 + bool have_User_Qinit = false; |
|
1713 + bool have_User_Control_matrix = false; |
|
1714 + bool do_numeric = true; |
|
1715 + bool no_scale = false; |
|
1716 + |
|
1717 + // find the operator |
|
1718 + for (int i = 0 ; i < nargin ; i++) |
|
1719 + { |
|
1720 + if (args(i).is_string ()) |
|
1721 + { |
|
1722 + op = i; |
|
1723 + break; |
|
1724 + } |
|
1725 + } |
|
1726 + |
|
1727 + if (op > 0) |
|
1728 + { |
|
1729 + std::string op_type = args (op).string_value (); |
|
1730 + |
|
1731 + if (op_type == "\\") |
|
1732 + { |
|
1733 + |
|
1734 + // matrix left divide, x = A\b |
|
1735 + |
|
1736 + // [x, Info] = umfpack (A, '\', b) ; |
|
1737 + // [x, Info] = umfpack (A, '\', b, Control) ; |
|
1738 + // [x, Info] = umfpack (A, Qinit, '\', b, Control) ; |
|
1739 + // [x, Info] = umfpack (A, Qinit, '\', b) ; |
|
1740 + |
|
1741 + do_solve = true; |
|
1742 + operation = "x = A\\b"; |
|
1743 + |
|
1744 + if (args(0).class_name () != "sparse") |
|
1745 + { |
|
1746 + error ("umfpack: input matrix A must be sparse"); |
|
1747 + return retval; |
|
1748 + } |
|
1749 + |
|
1750 + if (args(0).is_complex_type ()) |
|
1751 + { |
|
1752 + CAmatrix = args(0).sparse_complex_matrix_value (); |
|
1753 + A_is_complex = true; |
|
1754 + } |
|
1755 + else |
|
1756 + Amatrix = args(0).sparse_matrix_value (); |
|
1757 + |
|
1758 + |
|
1759 + if (args(op+1).is_complex_type ()) |
|
1760 + { |
|
1761 + CBmatrix = args(op+1).complex_matrix_value (); |
|
1762 + B_is_complex = true; |
|
1763 + } |
|
1764 + else |
|
1765 + Bmatrix = args(op+1).matrix_value (); |
|
1766 + |
|
1767 + if (nargout == 2) |
|
1768 + do_info = 1; |
|
1769 + |
|
1770 + if (op == 2) |
|
1771 + { |
|
1772 + User_Qinit = args(1).column_vector_value (); |
|
1773 + have_User_Qinit = true; |
|
1774 + } |
|
1775 + |
|
1776 + if ((op == 1 && nargin == 4) || (op == 2 && nargin == 5)) |
|
1777 + { |
|
1778 + User_Control_matrix = args(nargin-1).array_value (); |
|
1779 + have_User_Control_matrix = true; |
|
1780 + } |
|
1781 + |
|
1782 + if (error_state) |
|
1783 + { |
|
1784 + error ("umfpack: incorrect argument type"); |
|
1785 + return retval; |
|
1786 + } |
|
1787 + |
|
1788 + if (nargin < 3 || nargin > 5 || nargout > 2) |
|
1789 + { |
|
1790 + error ("umfpack: wrong number of arguments"); |
|
1791 + return retval; |
|
1792 + } |
|
1793 + |
|
1794 + } |
|
1795 + else if (op_type == "/") |
|
1796 + { |
|
1797 + // matrix right divide, x = b/A |
|
1798 + |
|
1799 + // [x, Info] = umfpack (b, '/', A) ; |
|
1800 + // [x, Info] = umfpack (b, '/', A, Control) ; |
|
1801 + // [x, Info] = umfpack (b, '/', A, Qinit) ; |
|
1802 + // [x, Info] = umfpack (b, '/', A, Qinit, Control) ; |
|
1803 + |
|
1804 + do_solve = true; |
|
1805 + operation = "x = b/A" ; |
|
1806 + |
|
1807 + transpose = true; |
|
1808 + |
|
1809 + if (args(2).class_name () != "sparse") |
|
1810 + { |
|
1811 + error ("umfpack: input matrix A must be sparse"); |
|
1812 + return retval; |
|
1813 + } |
|
1814 + |
|
1815 + if (args(2).is_complex_type ()) |
|
1816 + { |
|
1817 + CAmatrix = args(2).sparse_complex_matrix_value (); |
|
1818 + A_is_complex = true; |
|
1819 + } |
|
1820 + else |
|
1821 + Amatrix = args(2).sparse_matrix_value (); |
|
1822 + |
|
1823 + if (args(0).is_complex_type ()) |
|
1824 + { |
|
1825 + CBmatrix = args(0).complex_matrix_value (); |
|
1826 + B_is_complex = true; |
|
1827 + } |
|
1828 + else |
|
1829 + Bmatrix = args(0).matrix_value (); |
|
1830 + |
|
1831 + if (nargout == 2) |
|
1832 + do_info = 1; |
|
1833 + |
|
1834 + if (nargin == 5) |
|
1835 + { |
|
1836 + User_Qinit = args(3).column_vector_value (); |
|
1837 + User_Control_matrix = args(4).array_value (); |
|
1838 + have_User_Qinit = true; |
|
1839 + have_User_Control_matrix = true; |
|
1840 + } |
|
1841 + else if (nargin == 4) |
|
1842 + { |
|
1843 + User_Control_matrix = args(3).array_value (); |
|
1844 + |
|
1845 + if (User_Control_matrix.rows () == 1) |
|
1846 + { |
|
1847 + User_Qinit = args(3).column_vector_value (); |
|
1848 + have_User_Qinit = true; |
|
1849 + } |
|
1850 + else |
|
1851 + have_User_Control_matrix = true; |
|
1852 + } |
|
1853 + else if (nargin < 3 || nargin > 5 || nargout > 2) |
|
1854 + { |
|
1855 + error ("umfpack: wrong number of arguments"); |
|
1856 + return retval; |
|
1857 + } |
|
1858 + |
|
1859 + if (error_state) |
|
1860 + { |
|
1861 + error ("umfpack: incorrect argument type"); |
|
1862 + return retval; |
|
1863 + } |
|
1864 + } |
|
1865 + else if (op_type == "symbolic") |
|
1866 + { |
|
1867 + // symbolic factorization only |
|
1868 + |
|
1869 + // [P Q Fr Ch Info] = umfpack (A, 'symbolic') ; |
|
1870 + // [P Q Fr Ch Info] = umfpack (A, 'symbolic', Control) ; |
|
1871 + // [P Q Fr Ch Info] = umfpack (A, Qinit, 'symbolic') ; |
|
1872 + // [P Q Fr Ch Info] = umfpack (A, Qinit, 'symbolic', Control) ; |
|
1873 + |
|
1874 + operation = "symbolic factorization"; |
|
1875 + do_numeric = false; |
|
1876 + |
|
1877 + if (args(0).class_name () != "sparse") |
|
1878 + { |
|
1879 + error ("umfpack: input matrix A must be sparse"); |
|
1880 + return retval; |
|
1881 + } |
|
1882 + |
|
1883 + if (args(0).is_complex_type ()) |
|
1884 + { |
|
1885 + CAmatrix = args(0).sparse_complex_matrix_value (); |
|
1886 + A_is_complex = true; |
|
1887 + } |
|
1888 + else |
|
1889 + Amatrix = args(0).sparse_matrix_value (); |
|
1890 + |
|
1891 + if (nargout == 5) |
|
1892 + do_info = 4 ; |
|
1893 + |
|
1894 + if (op == 2) |
|
1895 + { |
|
1896 + User_Qinit = args(1).column_vector_value (); |
|
1897 + have_User_Qinit = true; |
|
1898 + } |
|
1899 + if ((op == 1 && nargin == 3) || (op == 2 && nargin == 4)) |
|
1900 + { |
|
1901 + User_Control_matrix = args(nargin-1).array_value (); |
|
1902 + have_User_Control_matrix = true; |
|
1903 + } |
|
1904 + |
|
1905 + if (error_state) |
|
1906 + { |
|
1907 + error ("umfpack: incorrect argument type"); |
|
1908 + return retval; |
|
1909 + } |
|
1910 + |
|
1911 + if (nargin < 2 || nargin > 4 || nargout > 5 || nargout < 4) |
|
1912 + { |
|
1913 + error ("umfpack: wrong number of arguments") ; |
|
1914 + return retval; |
|
1915 + } |
|
1916 + } |
|
1917 + else |
|
1918 + { |
|
1919 + error ("operator must be '/', '\\', or 'symbolic'") ; |
|
1920 + return retval; |
|
1921 + } |
|
1922 + } |
|
1923 + else if (nargin > 0) |
|
1924 + { |
|
1925 + // LU factorization |
|
1926 + |
|
1927 + // with scaling: |
|
1928 + // [L, U, P, Q, R, Info] = umfpack (A) ; |
|
1929 + // [L, U, P, Q, R, Info] = umfpack (A, Qinit) ; |
|
1930 + // |
|
1931 + // scaling determined by Control settings: |
|
1932 + // [L, U, P, Q, R, Info] = umfpack (A, Control) ; |
|
1933 + // [L, U, P, Q, R, Info] = umfpack (A, Qinit, Control) ; |
|
1934 + // |
|
1935 + // with no scaling: |
|
1936 + // [L, U, P, Q] = umfpack (A) ; |
|
1937 + // [L, U, P, Q] = umfpack (A, Control) ; |
|
1938 + // [L, U, P, Q] = umfpack (A, Qinit) ; |
|
1939 + // [L, U, P, Q] = umfpack (A, Qinit, Control) ; |
|
1940 + |
|
1941 + operation = "numeric factorization" ; |
|
1942 + |
|
1943 + if (args(0).is_complex_type ()) |
|
1944 + { |
|
1945 + CAmatrix = args(0).sparse_complex_matrix_value (); |
|
1946 + A_is_complex = true; |
|
1947 + } |
|
1948 + else |
|
1949 + Amatrix = args(0).sparse_matrix_value (); |
|
1950 + |
|
1951 + no_scale = nargout <= 4 ; |
|
1952 + |
|
1953 + if (nargout == 6) |
|
1954 + do_info = 5 ; |
|
1955 + |
|
1956 + if (nargin == 3) |
|
1957 + { |
|
1958 + User_Qinit = args(1).column_vector_value (); |
|
1959 + User_Control_matrix = args(2).array_value (); |
|
1960 + have_User_Qinit = true; |
|
1961 + have_User_Control_matrix = true; |
|
1962 + } |
|
1963 + else if (nargin == 2) |
|
1964 + { |
|
1965 + User_Control_matrix = args(1).array_value (); |
|
1966 + |
|
1967 + if (User_Control_matrix.rows () == 1) |
|
1968 + { |
|
1969 + User_Qinit = args(1).column_vector_value (); |
|
1970 + have_User_Qinit = true; |
|
1971 + } |
|
1972 + else |
|
1973 + have_User_Control_matrix = true; |
|
1974 + } |
|
1975 + else if (nargin > 3 || nargout > 6 || nargout < 4) |
|
1976 + { |
|
1977 + error ("umfpack: wrong number of arguments") ; |
|
1978 + return retval; |
|
1979 + } |
|
1980 + } |
|
1981 + else |
|
1982 + { |
|
1983 + // return default control settings |
|
1984 + |
|
1985 + // Control = umfpack ; |
|
1986 + // umfpack ; |
|
1987 + |
|
1988 + if (nargout > 1) |
|
1989 + { |
|
1990 + error ("umfpack: wrong number of arguments") ; |
|
1991 + return retval; |
|
1992 + } |
|
1993 + |
|
1994 + NDArray user_control (dim_vector (UMFPACK_CONTROL, 1)); |
|
1995 + double *user_control_ptr = user_control.fortran_vec (); |
|
1996 + umfpack_di_defaults (user_control_ptr); |
|
1997 + retval(0) = user_control; |
|
1998 + return retval; |
|
1999 + } |
|
2000 + |
|
2001 + // check inputs |
|
2002 + |
|
2003 + int n_row = Amatrix.rows (); |
|
2004 + int n_col = Amatrix.cols (); |
|
2005 + int nn = MAX (n_row, n_col) ; |
|
2006 + int n_inner = MIN (n_row, n_col) ; |
|
2007 + if (do_solve && n_row != n_col) |
|
2008 + { |
|
2009 + error ("umfpack: input matrix A must square for '\\' or '/'") ; |
|
2010 + return retval; |
|
2011 + } |
|
2012 + if (n_row == 0 || n_col == 0) |
|
2013 + { |
|
2014 + error ("umfpack: input matrix A cannot have zero rows or zero columns") ; |
|
2015 + return retval; |
|
2016 + } |
|
2017 + |
|
2018 + /* The real/complex status of A determines which version to use, */ |
|
2019 + /* (umfpack_di_* or umfpack_zi_*). */ |
|
2020 + const int *Ap; |
|
2021 + const int *Ai; |
|
2022 + const double *Ax; |
|
2023 + const double *Bx; |
|
2024 + |
|
2025 + if (A_is_complex) |
|
2026 + { |
|
2027 + Ap = CAmatrix.cidx (); |
|
2028 + Ai = CAmatrix.ridx (); |
|
2029 + Ax = X_CAST (const double *, CAmatrix.data ()); |
|
2030 + } |
|
2031 + else |
|
2032 + { |
|
2033 + Ap = Amatrix.cidx (); |
|
2034 + Ai = Amatrix.ridx (); |
|
2035 + Ax = Amatrix.data (); |
|
2036 + } |
|
2037 + |
|
2038 + if (B_is_complex) |
|
2039 + Bx = X_CAST (const double *, CBmatrix.fortran_vec ()); |
|
2040 + else |
|
2041 + Bx = Bmatrix.fortran_vec (); |
|
2042 + |
|
2043 + if (do_solve) |
|
2044 + { |
|
2045 + int b_row = Bmatrix.rows (); |
|
2046 + int b_col = Bmatrix.cols (); |
|
2047 + |
|
2048 + if (!((transpose && b_row == 1 && b_col == nn) || |
|
2049 + (!transpose && b_row == nn && b_col == 1))) |
|
2050 + { |
|
2051 + error ("umfpack: b has the wrong dimensions") ; |
|
2052 + return retval; |
|
2053 + } |
|
2054 + |
|
2055 + X_is_complex = A_is_complex || B_is_complex ; |
|
2056 + } |
|
2057 + |
|
2058 + // set the Control parameters |
|
2059 + NDArray Control (dim_vector (UMFPACK_CONTROL, 1)); |
|
2060 + double *Control_ptr = Control.fortran_vec (); |
|
2061 + if (A_is_complex) |
|
2062 + umfpack_zi_defaults (Control_ptr) ; |
|
2063 + else |
|
2064 + umfpack_di_defaults (Control_ptr) ; |
|
2065 + |
|
2066 + if (have_User_Control_matrix) |
|
2067 + { |
|
2068 + int size = MIN (UMFPACK_CONTROL, User_Control_matrix.length ()); |
|
2069 + for (int i = 0 ; i < size ; i++) |
|
2070 + Control (i) = User_Control_matrix (i) ; |
|
2071 + } |
|
2072 + |
|
2073 + if (no_scale) |
|
2074 + { |
|
2075 + // turn off scaling for [L, U, P, Q] = umfpack (A) ; |
|
2076 + // ignoring the input value of Control (24) for the usage |
|
2077 + // [L, U, P, Q] = umfpack (A, Control) ; |
|
2078 + Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE ; |
|
2079 + } |
|
2080 + |
|
2081 + int print_level; |
|
2082 + if (xisnan (Control (UMFPACK_PRL))) |
|
2083 + print_level = UMFPACK_DEFAULT_PRL ; |
|
2084 + else |
|
2085 + print_level = int (Control (UMFPACK_PRL)) ; |
|
2086 + |
|
2087 + Control (UMFPACK_PRL) = print_level ; |
|
2088 + |
|
2089 + // check Qinit, if present |
|
2090 + int *Qinit = NULL; |
|
2091 + if (have_User_Qinit) |
|
2092 + { |
|
2093 + if(User_Qinit.rows () != 1 || User_Qinit.cols () != n_col) |
|
2094 + { |
|
2095 + error ("umfpack: Qinit must be 1-by-n_col") ; |
|
2096 + return retval; |
|
2097 + } |
|
2098 + |
|
2099 + Qinit = new int [n_col]; |
|
2100 + for (int i = 0; i < n_col; i++) |
|
2101 + Qinit[i] = static_cast<int> (User_Qinit (i)); |
|
2102 + } |
|
2103 + |
|
2104 + // report the inputs A and Qinit |
|
2105 + |
|
2106 + if (print_level >= 2) |
|
2107 + // print the operation |
|
2108 + octave_stdout << "\numfpack: " << operation; |
|
2109 + |
|
2110 + if (A_is_complex) |
|
2111 + { |
|
2112 + umfpack_zi_report_control (Control_ptr) ; |
|
2113 + |
|
2114 + if (print_level >= 3) |
|
2115 + octave_stdout << "\nA: " ; |
|
2116 + |
|
2117 + umfpack_zi_report_matrix (n_row, n_col, Ap, Ai, Ax, NULL, |
|
2118 + 1, Control_ptr) ; |
|
2119 + if (have_User_Qinit) |
|
2120 + { |
|
2121 + if (print_level >= 3) |
|
2122 + octave_stdout << "\nQinit: " ; |
|
2123 + |
|
2124 + umfpack_zi_report_perm (n_col, Qinit, Control_ptr) ; |
|
2125 + } |
|
2126 + } |
|
2127 + else |
|
2128 + { |
|
2129 + umfpack_di_report_control (Control_ptr) ; |
|
2130 + |
|
2131 + if (print_level >= 3) |
|
2132 + octave_stdout << "\nA: " ; |
|
2133 + |
|
2134 + umfpack_di_report_matrix (n_row, n_col, Ap, Ai, Ax, |
|
2135 + 1, Control_ptr) ; |
|
2136 + if (have_User_Qinit) |
|
2137 + { |
|
2138 + if (print_level >= 3) |
|
2139 + octave_stdout << "\nQinit: " ; |
|
2140 + |
|
2141 + umfpack_di_report_perm (n_col, Qinit, Control_ptr) ; |
|
2142 + } |
|
2143 + } |
|
2144 + |
|
2145 +#ifndef NO_TRANSPOSE_FORWARD_SLASH |
|
2146 + // create the array transpose for x = b/A |
|
2147 + if (transpose) |
|
2148 + if (A_is_complex) |
|
2149 + { |
|
2150 + CAmatrix = CAmatrix.transpose (); |
|
2151 + Ap = Amatrix.cidx (); |
|
2152 + Ai = Amatrix.ridx (); |
|
2153 + Ax = X_CAST (const double *, CAmatrix.data ()); |
|
2154 + } |
|
2155 + else |
|
2156 + { |
|
2157 + Amatrix = Amatrix.transpose (); |
|
2158 + Ap = Amatrix.cidx (); |
|
2159 + Ai = Amatrix.ridx (); |
|
2160 + Ax = Amatrix.data (); |
|
2161 + } |
|
2162 +#endif |
|
2163 + |
|
2164 + // perform the symbolic factorization |
|
2165 + |
|
2166 + NDArray InfoOut (dim_vector (1, UMFPACK_INFO)); |
|
2167 + double * Info = InfoOut.fortran_vec (); |
|
2168 + void *Symbolic; |
|
2169 + int status, status2; |
|
2170 + if (A_is_complex) |
|
2171 + status = umfpack_zi_qsymbolic (n_row, n_col, Ap, Ai, Ax, NULL, |
|
2172 + Qinit, &Symbolic, Control_ptr, |
|
2173 + Info); |
|
2174 + else |
|
2175 + status = umfpack_di_qsymbolic (n_row, n_col, Ap, Ai, Ax, |
|
2176 + Qinit, &Symbolic, Control_ptr, |
|
2177 + Info); |
|
2178 + |
|
2179 + if (status < 0) |
|
2180 + { |
|
2181 + umfpack_error ("symbolic factorization failed", A_is_complex, |
|
2182 + nargout, retval, Control, InfoOut, status, do_info) ; |
|
2183 + return retval; |
|
2184 + } |
|
2185 + |
|
2186 + if (have_User_Qinit) |
|
2187 + delete [] Qinit; |
|
2188 + |
|
2189 + // report the Symbolic object |
|
2190 + |
|
2191 + if (A_is_complex) |
|
2192 + umfpack_zi_report_symbolic (Symbolic, Control_ptr) ; |
|
2193 + else |
|
2194 + umfpack_di_report_symbolic (Symbolic, Control_ptr) ; |
|
2195 + |
|
2196 + // perform numeric factorization, or just return symbolic factorization |
|
2197 + |
|
2198 + if (do_numeric) |
|
2199 + { |
|
2200 + // perform the numeric factorization |
|
2201 + void *Numeric; |
|
2202 + |
|
2203 + if (A_is_complex) |
|
2204 + status = umfpack_zi_numeric (Ap, Ai, Ax, NULL, Symbolic, &Numeric, |
|
2205 + Control_ptr, Info) ; |
|
2206 + else |
|
2207 + status = umfpack_di_numeric (Ap, Ai, Ax, Symbolic, &Numeric, |
|
2208 + Control_ptr, Info) ; |
|
2209 + |
|
2210 + // free the symbolic factorization |
|
2211 + if (A_is_complex) |
|
2212 + umfpack_zi_free_symbolic (&Symbolic) ; |
|
2213 + else |
|
2214 + umfpack_di_free_symbolic (&Symbolic) ; |
|
2215 + |
|
2216 + // report the Numeric object |
|
2217 + if (status < 0) |
|
2218 + { |
|
2219 + umfpack_error ("numeric factorization failed", A_is_complex, |
|
2220 + nargout, retval, Control, InfoOut, status, do_info); |
|
2221 + return retval; |
|
2222 + } |
|
2223 + |
|
2224 + if (A_is_complex) |
|
2225 + (void) umfpack_zi_report_numeric (Numeric, Control_ptr) ; |
|
2226 + else |
|
2227 + (void) umfpack_di_report_numeric (Numeric, Control_ptr) ; |
|
2228 + |
|
2229 + // return the solution or the factorization |
|
2230 + |
|
2231 + if (do_solve) |
|
2232 + { |
|
2233 + int sys; |
|
2234 + ComplexNDArray Xcmplx; |
|
2235 + NDArray Xreal; |
|
2236 + |
|
2237 + // solve Ax=b or A'x'=b', and return just the solution x |
|
2238 + |
|
2239 +#ifndef NO_TRANSPOSE_FORWARD_SLASH |
|
2240 + if (transpose) |
|
2241 + { |
|
2242 + // A.'x.'=b.' gives the same x=b/A as solving A'x'=b' |
|
2243 + // since C=A.' was factorized, solve with sys = UMFPACK_A |
|
2244 + // since x and b are vectors, x.' and b.' are implicit |
|
2245 + if (X_is_complex) |
|
2246 + Xcmplx.resize (dim_vector (1, nn)); |
|
2247 + else |
|
2248 + Xreal.resize (dim_vector (1, nn)); |
|
2249 + } |
|
2250 + else |
|
2251 + { |
|
2252 + if (X_is_complex) |
|
2253 + Xcmplx.resize (dim_vector (nn, 1)); |
|
2254 + else |
|
2255 + Xreal.resize (dim_vector (nn, 1)); |
|
2256 + } |
|
2257 + |
|
2258 + sys = UMFPACK_A ; |
|
2259 +#else |
|
2260 + if (transpose) |
|
2261 + { |
|
2262 + // If A is real, A'x=b is the same as A.'x=b. |
|
2263 + // x and b are vectors, so x and b are the same as x' and b'. |
|
2264 + // If A is complex, then A.'x.'=b.' gives the same solution x |
|
2265 + // as the complex conjugate transpose. If we used the A'x=b |
|
2266 + // option in umfpack_*_solve, we would have to form b' on |
|
2267 + // input and x' on output (negating the imaginary part). |
|
2268 + // We can save this work by just using the A.'x=b option in |
|
2269 + // umfpack_*_solve. Then, forming x.' and b.' is implicit, |
|
2270 + // since x and b are just vectors anyway. |
|
2271 + // In both cases, the system to solve is A.'x=b |
|
2272 + if (X_is_complex) |
|
2273 + Xcmplx.resize (dim_vector (1, nn)); |
|
2274 + else |
|
2275 + Xreal.resize (dim_vector (1, nn)); |
|
2276 + |
|
2277 + sys = UMFPACK_Aat ; |
|
2278 + } |
|
2279 + else |
|
2280 + { |
|
2281 + if (X_is_complex) |
|
2282 + Xcmplx.resize (dim_vector (nn, 1)); |
|
2283 + else |
|
2284 + Xreal.resize (dim_vector (nn, 1)); |
|
2285 + sys = UMFPACK_A ; |
|
2286 + } |
|
2287 +#endif |
|
2288 + |
|
2289 + // print the right-hand-side, B |
|
2290 + if (print_level >= 3) |
|
2291 + octave_stdout << "\nright-hand side, b: "; |
|
2292 + |
|
2293 + if (B_is_complex) |
|
2294 + (void) umfpack_zi_report_vector (nn, Bx, NULL, Control_ptr) ; |
|
2295 + else |
|
2296 + (void) umfpack_di_report_vector (nn, Bx, Control_ptr) ; |
|
2297 + |
|
2298 + // solve the system |
|
2299 + double * Xx; |
|
2300 + if (X_is_complex) |
|
2301 + Xx = X_CAST (double *, Xcmplx.fortran_vec ()); |
|
2302 + else |
|
2303 + Xx = Xreal.fortran_vec (); |
|
2304 + status2 = UMFPACK_OK ; |
|
2305 + |
|
2306 + if (A_is_complex) |
|
2307 + { |
|
2308 + if (!B_is_complex) |
|
2309 + { |
|
2310 + OCTAVE_LOCAL_BUFFER (double, Bz, nn); |
|
2311 + for (int i = 0; i < nn; i++) |
|
2312 + Bz[i] = 0.; |
|
2313 + |
|
2314 + status = umfpack_zi_solve (sys, Ap, Ai, Ax, NULL, Xx, NULL, |
|
2315 + Bx, Bz, Numeric, Control_ptr, |
|
2316 + Info); |
|
2317 + } |
|
2318 + else |
|
2319 + status = umfpack_zi_solve (sys, Ap, Ai, Ax, NULL, Xx, NULL, |
|
2320 + Bx, NULL, Numeric, Control_ptr, |
|
2321 + Info); |
|
2322 + } |
|
2323 + else |
|
2324 + { |
|
2325 + if (B_is_complex) |
|
2326 + { |
|
2327 + // Ax=b when b is complex and A is sparse can be split |
|
2328 + // into two systems, A*xr=br and A*xi=bi, where r denotes |
|
2329 + // the real part and i the imaginary part of x and b. |
|
2330 + OCTAVE_LOCAL_BUFFER (double, Tx, nn); |
|
2331 + OCTAVE_LOCAL_BUFFER (double, Tz, nn); |
|
2332 + |
|
2333 + status = umfpack_di_solve (sys, Ap, Ai, Ax, Tx, Bx, |
|
2334 + Numeric, Control_ptr, Info); |
|
2335 + status2 = umfpack_di_solve (sys, Ap, Ai, Ax, Tz, Bx, |
|
2336 + Numeric, Control_ptr, Info) ; |
|
2337 + |
|
2338 + for (int i = 0; i < nn; i++) |
|
2339 + Xcmplx (i) = Complex (Tx[i], Tz[i]); |
|
2340 + } |
|
2341 + else |
|
2342 + status = umfpack_di_solve (sys, Ap, Ai, Ax, Xx, Bx, |
|
2343 + Numeric, Control_ptr, Info); |
|
2344 + } |
|
2345 + |
|
2346 + // free the Numeric object |
|
2347 + if (A_is_complex) |
|
2348 + umfpack_zi_free_numeric (&Numeric) ; |
|
2349 + else |
|
2350 + umfpack_di_free_numeric (&Numeric) ; |
|
2351 + |
|
2352 + // check error status |
|
2353 + if (status < 0 || status2 < 0) |
|
2354 + { |
|
2355 + umfpack_error ("solve failed", A_is_complex, nargout, |
|
2356 + retval, Control, InfoOut, status, do_info) ; |
|
2357 + return retval; |
|
2358 + } |
|
2359 + |
|
2360 + // print the solution, X |
|
2361 + if (print_level >= 3) |
|
2362 + octave_stdout << "\nsolution, x: "; |
|
2363 + |
|
2364 + if (X_is_complex) |
|
2365 + (void) umfpack_zi_report_vector (nn, Xx, NULL, Control_ptr); |
|
2366 + else |
|
2367 + (void) umfpack_di_report_vector (nn, Xx, Control_ptr); |
|
2368 + |
|
2369 + // warn about singular or near-singular matrices |
|
2370 + // no warning is given if Control (1) is zero |
|
2371 + if (Control (UMFPACK_PRL) >= 1) |
|
2372 + { |
|
2373 + if (status == UMFPACK_WARNING_singular_matrix) |
|
2374 + { |
|
2375 + warning ("matrix is singular"); |
|
2376 + warning ("Try increasing Control (%d) and Control (%d).", |
|
2377 + 1+UMFPACK_PIVOT_TOLERANCE, |
|
2378 + 1+UMFPACK_SYM_PIVOT_TOLERANCE); |
|
2379 + warning ("(Suppress this warning with Control (%d) = 0.)", |
|
2380 + 1+UMFPACK_PRL); |
|
2381 + } |
|
2382 + else if (InfoOut (UMFPACK_RCOND) < DBL_EPSILON) |
|
2383 + { |
|
2384 + warning ("matrix is nearly singular, rcond = %g", |
|
2385 + InfoOut (UMFPACK_RCOND)); |
|
2386 + warning ("Try increasing Control (%d) and Control (%d).", |
|
2387 + 1+UMFPACK_PIVOT_TOLERANCE, |
|
2388 + 1+UMFPACK_SYM_PIVOT_TOLERANCE); |
|
2389 + warning ("(Suppress this warning with Control (%d) = 0.)", |
|
2390 + 1+UMFPACK_PRL); |
|
2391 + } |
|
2392 + } |
|
2393 + |
|
2394 + // Setup the return value |
|
2395 + if (X_is_complex) |
|
2396 + retval (0) = octave_value (Xcmplx); |
|
2397 + else |
|
2398 + retval (0) = octave_value (Xreal); |
|
2399 + } |
|
2400 + else |
|
2401 + { |
|
2402 + // get L, U, P, Q, and r |
|
2403 + int lnz, unz, ignore1, ignore2, ignore3; |
|
2404 + |
|
2405 + if (A_is_complex) |
|
2406 + status = umfpack_zi_get_lunz (&lnz, &unz, &ignore1, &ignore2, |
|
2407 + &ignore3, Numeric) ; |
|
2408 + else |
|
2409 + status = umfpack_di_get_lunz (&lnz, &unz, &ignore1, &ignore2, |
|
2410 + &ignore3, Numeric) ; |
|
2411 + |
|
2412 + if (status < 0) |
|
2413 + { |
|
2414 + if (A_is_complex) |
|
2415 + umfpack_zi_free_numeric (&Numeric) ; |
|
2416 + else |
|
2417 + umfpack_di_free_numeric (&Numeric) ; |
|
2418 + |
|
2419 + umfpack_error ("extracting LU factors failed", A_is_complex, |
|
2420 + nargout, retval, Control, InfoOut, status, |
|
2421 + do_info); |
|
2422 + return retval; |
|
2423 + } |
|
2424 + |
|
2425 + // avoid malloc of zero-sized arrays |
|
2426 + lnz = MAX (lnz, 1) ; |
|
2427 + unz = MAX (unz, 1) ; |
|
2428 + |
|
2429 + // get space for the *** ROW *** form of L |
|
2430 + SparseMatrix Lreal; |
|
2431 + SparseComplexMatrix Limag; |
|
2432 + int *Ltp, *Ltj; |
|
2433 + double *Ltx; |
|
2434 + if (A_is_complex) |
|
2435 + { |
|
2436 + Limag = SparseComplexMatrix (n_inner, n_row, lnz); |
|
2437 + Ltp = Limag.cidx (); |
|
2438 + Ltj = Limag.ridx (); |
|
2439 + Ltx = X_CAST (double *, Limag.data ()); |
|
2440 + } |
|
2441 + else |
|
2442 + { |
|
2443 + Lreal = SparseMatrix (n_inner, n_row, lnz); |
|
2444 + Ltp = Lreal.cidx (); |
|
2445 + Ltj = Lreal.ridx (); |
|
2446 + Ltx = Lreal.data (); |
|
2447 + } |
|
2448 + |
|
2449 + // create permanent copy of the output matrix U |
|
2450 + int *Up, *Ui; |
|
2451 + double *Ux; |
|
2452 + SparseMatrix Ureal; |
|
2453 + SparseComplexMatrix Uimag; |
|
2454 + |
|
2455 + if (A_is_complex) |
|
2456 + { |
|
2457 + Uimag = SparseComplexMatrix (n_inner, n_col, unz); |
|
2458 + Up = Uimag.cidx (); |
|
2459 + Ui = Uimag.ridx (); |
|
2460 + Ux = X_CAST (double *, Uimag.data ()); |
|
2461 + } |
|
2462 + else |
|
2463 + { |
|
2464 + Ureal = SparseMatrix (n_inner, n_col, unz); |
|
2465 + Up = Ureal.cidx (); |
|
2466 + Ui = Ureal.ridx (); |
|
2467 + Ux = Ureal.data (); |
|
2468 + } |
|
2469 + |
|
2470 + // temporary space for the integer permutation vectors |
|
2471 + OCTAVE_LOCAL_BUFFER (int, P, n_row); |
|
2472 + OCTAVE_LOCAL_BUFFER (int, Q, n_col); |
|
2473 + |
|
2474 + // get scale factors, if requested |
|
2475 + status2 = UMFPACK_OK ; |
|
2476 + SparseMatrix Rsout; |
|
2477 + double * Rs; |
|
2478 + if (!no_scale) |
|
2479 + { |
|
2480 + // create a diagonal sparse matrix for the scale factors |
|
2481 + Rsout = SparseMatrix (n_row, n_row, n_row); |
|
2482 + for (int i = 0 ; i < n_row ; i++) |
|
2483 + { |
|
2484 + Rsout.cidx (i) = i; |
|
2485 + Rsout.ridx (i) = i; |
|
2486 + } |
|
2487 + Rsout.cidx (n_row) = n_row; |
|
2488 + Rs = Rsout.data (); |
|
2489 + } |
|
2490 + else |
|
2491 + Rs = (double *) NULL ; |
|
2492 + |
|
2493 + // get Lt, U, P, Q, and Rs from the numeric object |
|
2494 + int do_recip; |
|
2495 + if (A_is_complex) |
|
2496 + { |
|
2497 + status = umfpack_zi_get_numeric (Ltp, Ltj, Ltx, NULL, Up, Ui, |
|
2498 + Ux, NULL, P, Q, |
|
2499 + (double *) NULL, |
|
2500 + (double *) NULL, |
|
2501 + &do_recip, Rs, Numeric) ; |
|
2502 + umfpack_zi_free_numeric (&Numeric) ; |
|
2503 + } |
|
2504 + else |
|
2505 + { |
|
2506 + status = umfpack_di_get_numeric (Ltp, Ltj, Ltx, Up, Ui, |
|
2507 + Ux, P, Q, (double *) NULL, |
|
2508 + &do_recip, Rs, Numeric) ; |
|
2509 + umfpack_di_free_numeric (&Numeric) ; |
|
2510 + } |
|
2511 + |
|
2512 + if (!no_scale) |
|
2513 + retval (4) = octave_vale (Rsout); |
|
2514 + |
|
2515 + // for the oct-file, -DNRECIPROCAL must be set, |
|
2516 + // so do_recip must be FALSE |
|
2517 + |
|
2518 + if (status < 0 || status2 < 0 || do_recip) |
|
2519 + { |
|
2520 + umfpack_error ("extracting LU factors failed", A_is_complex, |
|
2521 + nargout, retval, Control, InfoOut, status, |
|
2522 + do_info); |
|
2523 + return retval; |
|
2524 + } |
|
2525 + |
|
2526 + if (A_is_complex) |
|
2527 + retval (1) = octave_value (Uimag); |
|
2528 + else |
|
2529 + retval (1) = octave_valye (Ureal); |
|
2530 + |
|
2531 + // create sparse permutation matrix for P |
|
2532 + SparseMatrix Pout (n_row, n_row, n_row); |
|
2533 + for (int k = 0 ; k < n_row ; k++) |
|
2534 + { |
|
2535 + Pout.cidx (k) = k ; |
|
2536 + Pout.ridx (P [k]) = k; |
|
2537 + Pout.data (k) = 1; |
|
2538 + } |
|
2539 + Pout.cidx (n_row) = n_row; |
|
2540 + retval (2) = octave_value (Pout); |
|
2541 + |
|
2542 + // create sparse permutation matrix for Q |
|
2543 + SparseMatrix Qout (n_col, n_col, n_col); |
|
2544 + for (int k = 0 ; k < n_col ; k++) |
|
2545 + { |
|
2546 + Qout.cidx (k) = k ; |
|
2547 + Qout.ridx (k) = Q[k]; |
|
2548 + Qout.data (k) = 1; |
|
2549 + } |
|
2550 + Qout.cidx (n_col) = n_col; |
|
2551 + retval (3) = octave_value (Qout); |
|
2552 + |
|
2553 + // permanent copy of L |
|
2554 + if (A_is_complex) |
|
2555 + retval (0) = octave_value (Limag.transpose()); |
|
2556 + else |
|
2557 + retval (0) = octave_value (Lreal.transpose()); |
|
2558 + |
|
2559 + if (status < 0) |
|
2560 + { |
|
2561 + umfpack_error ("constructing L failed", A_is_complex, |
|
2562 + nargout, retval, Control, InfoOut, status, |
|
2563 + do_info) ; |
|
2564 + return octave_value (); |
|
2565 + } |
|
2566 + |
|
2567 + // print L, U, P, and Q |
|
2568 + if (A_is_complex) |
|
2569 + { |
|
2570 + if (print_level >= 3) |
|
2571 + { |
|
2572 + octave_stdout << "\nL: "; |
|
2573 + int *Lp = Limag.cidx (); |
|
2574 + int *Li = Limag.ridx (); |
|
2575 + double *Lx = X_CAST (double *, Limag.data ()); |
|
2576 + |
|
2577 + (void) umfpack_zi_report_matrix (n_row, n_inner, Lp, Li, |
|
2578 + Lx, NULL, 1, Control_ptr) ; |
|
2579 + } |
|
2580 + |
|
2581 + if (print_level >= 3) |
|
2582 + octave_stdout << "\nU: "; |
|
2583 + (void) umfpack_zi_report_matrix (n_inner, n_col, Up, Ui, |
|
2584 + Ux, NULL, 1, Control_ptr) ; |
|
2585 + if (print_level >= 3) |
|
2586 + octave_stdout << "\nP: "; |
|
2587 + (void) umfpack_zi_report_perm (n_row, P, Control_ptr); |
|
2588 + if (print_level >= 3) |
|
2589 + octave_stdout << "\nQ: "; |
|
2590 + (void) umfpack_zi_report_perm (n_col, Q, Control_ptr); |
|
2591 + } |
|
2592 + else |
|
2593 + { |
|
2594 + if (print_level >= 3) |
|
2595 + { |
|
2596 + int *Lp = Lreal.cidx (); |
|
2597 + int *Li = Lreal.ridx (); |
|
2598 + double *Lx = Lreal.data (); |
|
2599 + octave_stdout << "\nL: "; |
|
2600 + (void) umfpack_di_report_matrix (n_row, n_inner, Lp, Li, |
|
2601 + Lx, 1, Control_ptr); |
|
2602 + } |
|
2603 + |
|
2604 + if (print_level >= 3) |
|
2605 + octave_stdout << "\nU: "; |
|
2606 + (void) umfpack_di_report_matrix (n_inner, n_col, Up, Ui, |
|
2607 + Ux, 1, Control_ptr); |
|
2608 + if (print_level >= 3) |
|
2609 + octave_stdout << "\nP: "; |
|
2610 + (void) umfpack_di_report_perm (n_row, P, Control_ptr); |
|
2611 + if (print_level >= 3) |
|
2612 + octave_stdout << "\nQ: "; |
|
2613 + (void) umfpack_di_report_perm (n_col, Q, Control_ptr); |
|
2614 + } |
|
2615 + } |
|
2616 + } |
|
2617 + else |
|
2618 + { |
|
2619 + // return the symbolic factorization |
|
2620 + int ignore1, ignore2, ignore3; |
|
2621 + OCTAVE_LOCAL_BUFFER (int, Q, n_col); |
|
2622 + OCTAVE_LOCAL_BUFFER (int, P, n_row); |
|
2623 + OCTAVE_LOCAL_BUFFER (int, Front_npivcol, nn + 1); |
|
2624 + OCTAVE_LOCAL_BUFFER (int, Front_parent, nn + 1); |
|
2625 + OCTAVE_LOCAL_BUFFER (int, Front_1strow, nn + 1); |
|
2626 + OCTAVE_LOCAL_BUFFER (int, Front_leftmostdesc, nn + 1); |
|
2627 + OCTAVE_LOCAL_BUFFER (int, Chain_start, nn + 1); |
|
2628 + OCTAVE_LOCAL_BUFFER (int, Chain_maxrows, nn + 1); |
|
2629 + OCTAVE_LOCAL_BUFFER (int, Chain_maxcols, nn + 1); |
|
2630 + |
|
2631 + int nz, nfronts, nchains; |
|
2632 + |
|
2633 + if (A_is_complex) |
|
2634 + { |
|
2635 + status = umfpack_zi_get_symbolic (&ignore1, &ignore2, &ignore3, |
|
2636 + &nz, &nfronts, &nchains, P, Q, |
|
2637 + Front_npivcol, Front_parent, |
|
2638 + Front_1strow, |
|
2639 + Front_leftmostdesc, |
|
2640 + Chain_start, Chain_maxrows, |
|
2641 + Chain_maxcols, Symbolic) ; |
|
2642 + umfpack_zi_free_symbolic (&Symbolic) ; |
|
2643 + } |
|
2644 + else |
|
2645 + { |
|
2646 + status = umfpack_di_get_symbolic (&ignore1, &ignore2, &ignore3, |
|
2647 + &nz, &nfronts, &nchains, P, Q, |
|
2648 + Front_npivcol, Front_parent, |
|
2649 + Front_1strow, |
|
2650 + Front_leftmostdesc, |
|
2651 + Chain_start, Chain_maxrows, |
|
2652 + Chain_maxcols, Symbolic) ; |
|
2653 + umfpack_di_free_symbolic (&Symbolic) ; |
|
2654 + } |
|
2655 + |
|
2656 + if (status < 0) |
|
2657 + { |
|
2658 + umfpack_error ("extracting symbolic factors failed", |
|
2659 + A_is_complex, nargout, retval, Control, |
|
2660 + InfoOut, status, do_info) ; |
|
2661 + return retval; |
|
2662 + } |
|
2663 + |
|
2664 + // create sparse permutation matrix for P |
|
2665 + SparseMatrix Pout (n_row, n_row, n_row); |
|
2666 + for (int k = 0 ; k < n_row ; k++) |
|
2667 + { |
|
2668 + Pout.cidx (k) = k; |
|
2669 + Pout.ridx (P [k]) = k; |
|
2670 + Pout.data (k) = 1; |
|
2671 + } |
|
2672 + Pout.cidx (n_row) = n_row; |
|
2673 + retval (0) = octave_value (Pout); |
|
2674 + |
|
2675 + // create sparse permutation matrix for Q |
|
2676 + SparseMatrix Qout (n_col, n_col, n_col); |
|
2677 + for (int k = 0 ; k < n_col ; k++) |
|
2678 + { |
|
2679 + Qout.cidx (k) = k; |
|
2680 + Qout.ridx (k) = Q[k]; |
|
2681 + Qout.data (k) = 1; |
|
2682 + } |
|
2683 + Qout.cidx (n_col) = n_col; |
|
2684 + retval (1) = octave_value (Qout); |
|
2685 + |
|
2686 + // create Fr |
|
2687 + Matrix Frout (nfronts + 1, 4); |
|
2688 + for (int i = 0 ; i <= nfronts ; i++) |
|
2689 + { |
|
2690 + // convert parent, 1strow, and leftmostdesc to 1-based |
|
2691 + Frout (i, 0) = (double) (Front_npivcol [i]) ; |
|
2692 + Frout (i, 1) = (double) (Front_parent [i] + 1) ; |
|
2693 + Frout (i, 2) = (double) (Front_1strow [i] + 1) ; |
|
2694 + Frout (i, 3) = (double) (Front_leftmostdesc [i] + 1) ; |
|
2695 + } |
|
2696 + retval (2) = octave_value (Frout); |
|
2697 + |
|
2698 + // create Ch |
|
2699 + Matrix Chout (nchains + 1, 3); |
|
2700 + for (int i = 0 ; i <= nchains ; i++) |
|
2701 + { |
|
2702 + // convert to 1-based |
|
2703 + Chout (i, 0) = (double) (Chain_start [i] + 1) ; |
|
2704 + Chout (i, 1) = (double) (Chain_maxrows [i]) ; |
|
2705 + Chout (i, 2) = (double) (Chain_maxcols [i]) ; |
|
2706 + } |
|
2707 + Chout (0, nchains) = (double) Chain_start [nchains] + 1 ; |
|
2708 + Chout (1, nchains) = 0.; |
|
2709 + Chout (2, nchains) = 0.; |
|
2710 + retval (3) = octave_value (Chout); |
|
2711 + } |
|
2712 + |
|
2713 + // report Info |
|
2714 + if (A_is_complex) |
|
2715 + umfpack_zi_report_info (Control_ptr, Info); |
|
2716 + else |
|
2717 + umfpack_di_report_info (Control_ptr, Info); |
|
2718 + |
|
2719 + if (do_info > 0) |
|
2720 + retval (do_info) = InfoOut; |
|
2721 + |
|
2722 + return retval; |
|
2723 +} |
|
2724 + |
|
2725 +/* |
|
2726 +;;; Local Variables: *** |
|
2727 +;;; mode: C++ *** |
|
2728 +;;; End: *** |
|
2729 +*/ |
|
2730 diff -uNr UMFPACKv4.4.orig/UMFPACK/OCTAVE/umfpack_demo.m UMFPACK/UMFPACK/OCTAVE/umfpack_demo.m |
|
2731 --- UMFPACKv4.4.orig/UMFPACK/OCTAVE/umfpack_demo.m 1970-01-01 01:00:00.000000000 +0100 |
|
2732 +++ UMFPACK/UMFPACK/OCTAVE/umfpack_demo.m 2004-12-30 01:58:46.000000000 +0100 |
|
2733 @@ -0,0 +1,191 @@ |
|
2734 +function umfpack_demo |
|
2735 +% UMFPACK DEMO |
|
2736 +% |
|
2737 +% A demo of UMFPACK for OCTAVE. |
|
2738 +% |
|
2739 +% See also umfpack, umfpack_make, umfpack_details, umfpack_report, |
|
2740 +% and umfpack_simple. |
|
2741 + |
|
2742 +% UMFPACK Version 4.3 (Jan. 16, 2004), Copyright (c) 2004 by Timothy A. |
|
2743 +% Davis. All Rights Reserved. Type umfpack_details for License. |
|
2744 + |
|
2745 +%------------------------------------------------------------------------------- |
|
2746 +% get default control parameters |
|
2747 +%------------------------------------------------------------------------------- |
|
2748 + |
|
2749 +control = umfpack ; |
|
2750 +fprintf ('\nEnter the printing level for UMFPACK''s output statistics:\n') ; |
|
2751 +fprintf ('0: none, 1: errors only, 2: statistics, 4: print some of outputs\n') ; |
|
2752 +c = input ('5: print all output [default is 1]: ') ; |
|
2753 +if (isempty (c)) |
|
2754 + c = 1 ; |
|
2755 +end |
|
2756 +control (1) = c ; |
|
2757 + |
|
2758 +%------------------------------------------------------------------------------- |
|
2759 +% solve a simple system |
|
2760 +%------------------------------------------------------------------------------- |
|
2761 + |
|
2762 +fprintf ('\n--------------------------------------------------------------\n') ; |
|
2763 +fprintf ('Factor and solve a small system, Ax=b, using default parameters\n') ; |
|
2764 +if (control (1) > 1) |
|
2765 + fprintf ('(except for verbose printing enabled)\n') ; |
|
2766 +end |
|
2767 + |
|
2768 +load west0067 ; |
|
2769 +A = Problem.A ; |
|
2770 +n = size (A, 1) ; |
|
2771 +b = rand (n, 1) ; |
|
2772 + |
|
2773 +fprintf ('Solving Ax=b via UMFPACK:\n') ; |
|
2774 +[xu, info] = umfpack (A, '\\', b, control) ; |
|
2775 +x = xu ; |
|
2776 + |
|
2777 +fprintf ('Solving Ax=b via OCTAVE:\n') ; |
|
2778 +xm = A\b ; |
|
2779 +x = xm ; |
|
2780 + |
|
2781 +fprintf ('Difference between UMFPACK and OCTAVE solution: %g\n', ... |
|
2782 + norm (xu - xm, Inf)) ; |
|
2783 + |
|
2784 +%------------------------------------------------------------------------------- |
|
2785 +% spy the results |
|
2786 +%------------------------------------------------------------------------------- |
|
2787 + |
|
2788 +figure (1) ; |
|
2789 +clf |
|
2790 + |
|
2791 +subplot (2,3,1) |
|
2792 +title ('The matrix A') ; |
|
2793 +spy (A) |
|
2794 + |
|
2795 +subplot (2,3,2) |
|
2796 +[P1, Q1, Fr, Ch, Info] = umfpack (A, 'symbolic') ; |
|
2797 +title ('Supernodal column elimination tree') ; |
|
2798 +%% Disable for now !! |
|
2799 +%% treeplot (Fr (1:end-1,2)') ; |
|
2800 + |
|
2801 +subplot (2,3,3) |
|
2802 +title ('A, with initial row and column order') ; |
|
2803 +spy (P1 * A * Q1) |
|
2804 + |
|
2805 +subplot (2,3,4) |
|
2806 +fprintf ('\n--------------------------------------------------------------\n') ; |
|
2807 +fprintf ('\nFactorizing [L, U, P, Q, R] = umfpack (A)\n') ; |
|
2808 +[L, U, P, Q, R] = umfpack (A) ; |
|
2809 +title ('A, with final row/column order') ; |
|
2810 +spy (P*A*Q) |
|
2811 + |
|
2812 +fprintf ('\nP * (R\\A) * Q - L*U should be zero:\n') ; |
|
2813 +fprintf ('norm (P*(R\\A)*Q - L*U, 1) = %g (exact) %g (estimated)\n', ... |
|
2814 + norm (P * (R\A) * Q - L*U, 1), lu_normest (P * (R\A) * Q, L, U)) ; |
|
2815 + |
|
2816 +fprintf ('\nSolution to Ax=b via UMFPACK factorization:\n') ; |
|
2817 +fprintf ('x = Q * (U \\ (L \\ (P * (R \\ b))))\n') ; |
|
2818 +xu = Q * (U \ (L \ (P * (R \ b)))) ; |
|
2819 +x = xu ; |
|
2820 + |
|
2821 +fprintf ('\nUMFPACK flop count: %d\n', luflop (L, U)) ; |
|
2822 + |
|
2823 +subplot (2,3,5) |
|
2824 +title ('UMFPACK LU factors') ; |
|
2825 +spy (spones (L) + spones (U)) |
|
2826 + |
|
2827 +subplot (2,3,6) |
|
2828 +fprintf ('\nFactorizing [L, U, P] = lu (A (:, q))\n') ; |
|
2829 +try |
|
2830 + q = colamd (A) ; |
|
2831 +catch |
|
2832 + fprintf ('\n *** colamd not found, using colmmd instead *** \n') ; |
|
2833 + q = colmmd (A) ; |
|
2834 +end |
|
2835 +[L, U, P] = lu (A (:,q)) ; |
|
2836 +title ('OCTAVE LU factors') ; |
|
2837 +spy (spones (L) + spones (U)) |
|
2838 + |
|
2839 +fprintf ('\nSolution to Ax=b via OCTAVE factorization:\n') ; |
|
2840 +fprintf ('x = U \\ (L \\ (P * b)) ; x (q) = x ;\n') ; |
|
2841 +xm = U \ (L \ (P * b)) ; |
|
2842 +xm (q) = xm ; |
|
2843 + |
|
2844 +fprintf ('Difference between UMFPACK and OCTAVE solution: %g\n', ... |
|
2845 + norm (xu - xm, Inf)) ; |
|
2846 + |
|
2847 +fprintf ('\nOCTAVE LU flop count: %d\n', luflop (L, U)) ; |
|
2848 + |
|
2849 +%------------------------------------------------------------------------------- |
|
2850 +% solve A'x=b |
|
2851 +%------------------------------------------------------------------------------- |
|
2852 + |
|
2853 +fprintf ('\n--------------------------------------------------------------\n') ; |
|
2854 +fprintf ('Solve A''x=b:\n') ; |
|
2855 + |
|
2856 +fprintf ('Solving A''x=b via UMFPACK:\n') ; |
|
2857 +[xu, info] = umfpack (b', '/', A, control) ; |
|
2858 +xu = xu' ; |
|
2859 + |
|
2860 +fprintf ('Solving A''x=b via OCTAVE:\n') ; |
|
2861 +xm = (b'/A)' ; |
|
2862 +x = xm ; |
|
2863 + |
|
2864 +fprintf ('Difference between UMFPACK and OCTAVE solution: %g\n', ... |
|
2865 + norm (xu - xm, Inf)) ; |
|
2866 + |
|
2867 +%------------------------------------------------------------------------------- |
|
2868 +% factor A' and then solve Ax=b using the factors of A' |
|
2869 +%------------------------------------------------------------------------------- |
|
2870 + |
|
2871 +fprintf ('\n--------------------------------------------------------------\n') ; |
|
2872 +fprintf ('Compute C = A'', and compute the LU factorization of C.\n') ; |
|
2873 +fprintf ('Factorizing A'' can sometimes be better than factorizing A itself\n'); |
|
2874 +fprintf ('(less work and memory usage). Solve C''x=b; the solution is the\n') ; |
|
2875 +fprintf ('same as the solution to Ax=b for the original A.\n'); |
|
2876 + |
|
2877 +C = A' ; |
|
2878 + |
|
2879 +% factorize C (P,Q) = L*U |
|
2880 +[L, U, P, Q, R, info] = umfpack (C, control) ; |
|
2881 + |
|
2882 +fprintf ('\nP * (R\\C) * Q - L*U should be zero:\n') ; |
|
2883 +fprintf ('norm (P*(R\\C)*Q - L*U, 1) = %g (exact) %g (estimated)\n', ... |
|
2884 + norm (P * (R\C) * Q - L*U, 1), lu_normest (P * (R\C) * Q, L, U)) ; |
|
2885 + |
|
2886 +fprintf ('\nSolution to Ax=b via UMFPACK, using the factors of C:\n') ; |
|
2887 +fprintf ('x = R \\ (P'' * (L'' \\ (U'' \\ (Q'' * b)))) ;\n') ; |
|
2888 +xu = R \ (P' * (L' \ (U' \ (Q' * b)))) ; |
|
2889 +x = xu ; |
|
2890 + |
|
2891 +fprintf ('Solution to Ax=b via OCTAVE:\n') ; |
|
2892 +xm = A\b ; |
|
2893 +x = xm ; |
|
2894 + |
|
2895 +fprintf ('Difference between UMFPACK and OCTAVE solution: %g\n', ... |
|
2896 + norm (xu - xm, Inf)) ; |
|
2897 + |
|
2898 +%------------------------------------------------------------------------------- |
|
2899 +% solve Ax=B |
|
2900 +%------------------------------------------------------------------------------- |
|
2901 + |
|
2902 +fprintf ('\n--------------------------------------------------------------\n') ; |
|
2903 +fprintf ('\nSolve AX=B, where B is n-by-10, and sparse\n') ; |
|
2904 +B = sprandn (n, 10, 0.05) ; |
|
2905 +XU = umfpack_solve (A, '\\', B, control) ; |
|
2906 +XM = A\B ; |
|
2907 + |
|
2908 +fprintf ('Difference between UMFPACK and OCTAVE solution: %g\n', ... |
|
2909 + norm (XU - XM, Inf)) ; |
|
2910 + |
|
2911 +fprintf ('\n--------------------------------------------------------------\n') ; |
|
2912 +fprintf ('\nSolve AX=B, where B is n-by-10, and sparse, using umfpack_btf\n') ; |
|
2913 +XU = umfpack_btf (A, B, control) ; |
|
2914 + |
|
2915 +fprintf ('Difference between UMFPACK and OCTAVE solution: %g\n', ... |
|
2916 + norm (XU - XM, Inf)) ; |
|
2917 + |
|
2918 +fprintf ('\n--------------------------------------------------------------\n') ; |
|
2919 +fprintf ('\nSolve A''X=B, where B is n-by-10, and sparse\n') ; |
|
2920 +XU = umfpack_solve (B', '/', A, control) ; |
|
2921 +XM = B'/A ; |
|
2922 + |
|
2923 +fprintf ('Difference between UMFPACK and OCTAVE solution: %g\n', ... |
|
2924 + norm (XU - XM, Inf)) ; |
|
2925 diff -uNr UMFPACKv4.4.orig/UMFPACK/OCTAVE/umfpack_demo.m.out UMFPACK/UMFPACK/OCTAVE/umfpack_demo.m.out |
|
2926 --- UMFPACKv4.4.orig/UMFPACK/OCTAVE/umfpack_demo.m.out 1970-01-01 01:00:00.000000000 +0100 |
|
2927 +++ UMFPACK/UMFPACK/OCTAVE/umfpack_demo.m.out 2004-12-30 01:58:46.000000000 +0100 |
|
2928 @@ -0,0 +1,72 @@ |
|
2929 +>> umfpack_demo |
|
2930 + |
|
2931 +Enter the printing level for UMFPACK's output statistics: |
|
2932 +0: none, 1: errors only, 2: statistics, 4: print some of outputs |
|
2933 +5: print all output [default is 1]: |
|
2934 + |
|
2935 +-------------------------------------------------------------- |
|
2936 +Factor and solve a small system, Ax=b, using default parameters |
|
2937 +Solving Ax=b via UMFPACK: |
|
2938 +Solving Ax=b via MATLAB: |
|
2939 +Difference between UMFPACK and MATLAB solution: 1.24345e-14 |
|
2940 + |
|
2941 +-------------------------------------------------------------- |
|
2942 + |
|
2943 +Factorizing [L, U, P, Q, R] = umfpack (A) |
|
2944 + |
|
2945 +P * (R\A) * Q - L*U should be zero: |
|
2946 +norm (P*(R\A)*Q - L*U, 1) = 4.2068e-16 (exact) 3.74627e-16 (estimated) |
|
2947 + |
|
2948 +Solution to Ax=b via UMFPACK factorization: |
|
2949 +x = Q * (U \ (L \ (P * (R \ b)))) |
|
2950 + |
|
2951 +UMFPACK flop count: 2362 |
|
2952 + |
|
2953 +Factorizing [L, U, P] = lu (A (:, q)) |
|
2954 +If you are using a version of MATLAB prior to V6.0, then the |
|
2955 +following statement (q = colamd (A)) may fail. Either download |
|
2956 +colamd from http://www.cise.ufl.edu/research/sparse, upgrade to |
|
2957 +MATLAB V6.0 or later, or replace the statement with |
|
2958 +q = colmmd (A) ; |
|
2959 + |
|
2960 +Solution to Ax=b via MATLAB factorization: |
|
2961 +x = U \ (L \ (P * b)) ; x (q) = x ; |
|
2962 +Difference between UMFPACK and MATLAB solution: 1.37668e-14 |
|
2963 + |
|
2964 +MATLAB LU flop count: 3164 |
|
2965 + |
|
2966 +-------------------------------------------------------------- |
|
2967 +Solve A'x=b: |
|
2968 +Solving A'x=b via UMFPACK: |
|
2969 +Solving A'x=b via MATLAB: |
|
2970 +Difference between UMFPACK and MATLAB solution: 3.10862e-15 |
|
2971 + |
|
2972 +-------------------------------------------------------------- |
|
2973 +Compute C = A', and compute the LU factorization of C. |
|
2974 +Factorizing A' can sometimes be better than factorizing A itself |
|
2975 +(less work and memory usage). Solve C'x=b; the solution is the |
|
2976 +same as the solution to Ax=b for the original A. |
|
2977 + |
|
2978 +P * (R\C) * Q - L*U should be zero: |
|
2979 +norm (P*(R\C)*Q - L*U, 1) = 1.31839e-16 (exact) 6.41848e-17 (estimated) |
|
2980 + |
|
2981 +Solution to Ax=b via UMFPACK, using the factors of C: |
|
2982 +x = R \ (P' * (L' \ (U' \ (Q' * b)))) ; |
|
2983 +Solution to Ax=b via MATLAB: |
|
2984 +Difference between UMFPACK and MATLAB solution: 1.77636e-14 |
|
2985 + |
|
2986 +-------------------------------------------------------------- |
|
2987 + |
|
2988 +Solve AX=B, where B is n-by-10, and sparse |
|
2989 +Difference between UMFPACK and MATLAB solution: 2.88198e-14 |
|
2990 + |
|
2991 +-------------------------------------------------------------- |
|
2992 + |
|
2993 +Solve AX=B, where B is n-by-10, and sparse, using umfpack_btf |
|
2994 +Difference between UMFPACK and MATLAB solution: 9.79736e-14 |
|
2995 + |
|
2996 +-------------------------------------------------------------- |
|
2997 + |
|
2998 +Solve A'X=B, where B is n-by-10, and sparse |
|
2999 +Difference between UMFPACK and MATLAB solution: 1.05244e-13 |
|
3000 +>> diary off |
|
3001 diff -uNr UMFPACKv4.4.orig/UMFPACK/OCTAVE/umfpack_make.m UMFPACK/UMFPACK/OCTAVE/umfpack_make.m |
|
3002 --- UMFPACKv4.4.orig/UMFPACK/OCTAVE/umfpack_make.m 1970-01-01 01:00:00.000000000 +0100 |
|
3003 +++ UMFPACK/UMFPACK/OCTAVE/umfpack_make.m 2004-12-30 01:58:46.000000000 +0100 |
|
3004 @@ -0,0 +1,356 @@ |
|
3005 +function umfpack_make |
|
3006 +% UMFPACK_MAKE |
|
3007 +% |
|
3008 +% Compiles the UMFPACK mexFunction and then runs a simple demo. |
|
3009 +% |
|
3010 +% UMFPACK Version 4.3 (Jan. 16, 2004), Copyright (c) 2004 by Timothy A. |
|
3011 +% Davis. All Rights Reserved. Type umfpack_details for License. |
|
3012 +% |
|
3013 +% See also: umfpack, umfpack_details, umfpack_report, umfpack_demo, and |
|
3014 +% umfpack_simple. |
|
3015 + |
|
3016 +help umfpack_make |
|
3017 + |
|
3018 +fprintf ('\n--------------------------------------------------------------\n') ; |
|
3019 +fprintf ('Now compiling the UMFPACK and AMD mexFunctions.\n') ; |
|
3020 +fprintf ('--------------------------------------------------------------\n') ; |
|
3021 + |
|
3022 +try |
|
3023 + % ispc does not appear in MATLAB 5.3 |
|
3024 + pc = ispc ; |
|
3025 +catch |
|
3026 + % if ispc fails, assume we aren't on a Windows PC. |
|
3027 + pc = 0 ; |
|
3028 +end |
|
3029 + |
|
3030 +obj = 'o' ; |
|
3031 +blas_lib = '' ; |
|
3032 +if (pc) |
|
3033 + obj = 'obj' ; |
|
3034 +end |
|
3035 + |
|
3036 +%------------------------------------------------------------------------------- |
|
3037 +% BLAS option |
|
3038 +%------------------------------------------------------------------------------- |
|
3039 + |
|
3040 +msg = [ ... |
|
3041 + '\nUsing the BLAS is faster, but might not compile correctly.\n', ... |
|
3042 + 'If you get an error stating that dgemm, dgemv, dger, zgemm,\n', ... |
|
3043 + 'zgemv, and/or zger are not defined, then recompile without the\n', ... |
|
3044 + 'BLAS. You can ignore warnings that these routines are implicitly\n', ... |
|
3045 + 'declared.\n\nPlease select one of the following options: \n', ... |
|
3046 + ' 1: attempt to compile with the BLAS (default)\n', ... |
|
3047 + ' 2: do not use the BLAS\n'] ; |
|
3048 +fprintf (msg) ; |
|
3049 +blas = input (': ') ; |
|
3050 +if (isempty (blas)) |
|
3051 + blas = 1 ; |
|
3052 +end |
|
3053 +if (blas == 1) |
|
3054 + % try to link to MATLAB's built-in BLAS |
|
3055 + blas = '' ; |
|
3056 + if (pc) |
|
3057 + % the default lcc compiler needs this library to access the BLAS |
|
3058 + blas_lib = ' libmwlapack.lib' ; |
|
3059 + msg = [ ... |
|
3060 + '\nCheck to see if you have a file called libmwlapack.lib in the\n', ... |
|
3061 + '<matlab>\\extern\\lib\\win32\\lcc\\ directory, where <matlab> is ', ... |
|
3062 + 'the\ndirectory where MATLAB is installed. If a file of that ', ... |
|
3063 + 'name is already\nthere, then you don''t have to do anything. ', ... |
|
3064 + 'Otherwise, you must first\ncopy the libmwlapack.lib file from ', ... |
|
3065 + 'the umfpack\\lcc_lib\\ directory to the\n', ... |
|
3066 + '<matlab>\\extern\\lib\\win32\\lcc\\ directory. Next, type\n\n', ... |
|
3067 + ' mex -setup\n\n', ... |
|
3068 + 'at the MATLAB prompt, and ask MATLAB to select the lcc compiler. ',... |
|
3069 + 'You can skip\nall of this if you have already done it, or have ', ... |
|
3070 + 'configured mex to use\na different compiler. If you are using ', ... |
|
3071 + 'Norton anti-virus software on Windows\n98SE, then you need to ', ... |
|
3072 + 'exit MATLAB, turn off virus checking, and restart MATLAB\n', ... |
|
3073 + 'before you can use the mex command or compile UMFPACK.\n', ... |
|
3074 + 'You may also need to turn off virus checking in other cases.\n', ... |
|
3075 + '\nHit enter to continue, or type control-C if you do not wish to '] ; |
|
3076 + fprintf (msg) ; |
|
3077 + input ('proceed: ') ; |
|
3078 + end |
|
3079 + fprintf ('\nUsing the BLAS (recommended).\n') ; |
|
3080 +else |
|
3081 + % No BLAS |
|
3082 + fprintf ('\nNot using the BLAS. UMFPACK will be slow.\n') ; |
|
3083 + blas = ' -DNBLAS' ; |
|
3084 +end |
|
3085 + |
|
3086 +%------------------------------------------------------------------------------- |
|
3087 +% -DNUTIL option (using utMalloc or mxMalloc) |
|
3088 +%------------------------------------------------------------------------------- |
|
3089 + |
|
3090 +utils = '' ; |
|
3091 + |
|
3092 +if (~pc) |
|
3093 + msg = [ ... |
|
3094 + '--------------------------------------------------------------\n', ... |
|
3095 + '\nUMFPACK uses MATLAB''s memory allocation routines. The internal', ... |
|
3096 + '\nutMalloc, utFree, and utRealloc allow for better use of memory,', ... |
|
3097 + '\nbut they are internal utility routines that are not documented.\n', ... |
|
3098 + 'Thus, they might not always work. Using mxMalloc, mxFree, and\n', ... |
|
3099 + 'mxRealloc works, but UMFPACK might run out of memory when solving\n', ... |
|
3100 + 'problems that it could otherwise solve. Try using the default.\n', ... |
|
3101 + 'If you get an error stating that utMalloc, utFree, and/or\n', ... |
|
3102 + 'utRealloc are not defined, then recompile with the mx* routines.\n', ... |
|
3103 + '\nPlease select one of the following options:\n', ... |
|
3104 + ' 1: attempt to use the ut* routines (default)\n', ... |
|
3105 + ' 2: use the standard mx* routines\n'] ; |
|
3106 + fprintf (msg) ; |
|
3107 + utils = input (': ') ; |
|
3108 + if (isempty (utils)) |
|
3109 + utils = 1 ; |
|
3110 + end |
|
3111 + if (utils == 2) |
|
3112 + fprintf ('\nNot using utMalloc, utFree, or utRealloc\n') ; |
|
3113 + utils = ' -DNUTIL' ; |
|
3114 + else |
|
3115 + fprintf ('\nUsing utMalloc, utFree, and utRealloc\n') ; |
|
3116 + utils = '' ; |
|
3117 + end |
|
3118 +end |
|
3119 + |
|
3120 +%------------------------------------------------------------------------------- |
|
3121 +% -DNPOSIX option (for sysconf and times timer routines) |
|
3122 +%------------------------------------------------------------------------------- |
|
3123 + |
|
3124 +posix = '' ; |
|
3125 + |
|
3126 +if (~pc) |
|
3127 + msg = [ ... |
|
3128 + '--------------------------------------------------------------\n', ... |
|
3129 + '\nUMFPACK can use the POSIX routines sysconf () and times ()\n', ... |
|
3130 + 'to provide CPU time and wallclock time statistics. If you do not\n', ... |
|
3131 + 'have a POSIX-compliant operating system, then UMFPACK won''t\n', ... |
|
3132 + 'compile. If you don''t know which option to pick, try the\n', ... |
|
3133 + 'default. If you get an error saying that sysconf and/or times\n', ... |
|
3134 + 'are not defined, then recompile with the non-POSIX option.\n', ... |
|
3135 + '\nPlease select one of the following options:\n', ... |
|
3136 + ' 1: use POSIX sysconf and times routines (default)\n', ... |
|
3137 + ' 2: do not use POSIX routines\n'] ; |
|
3138 + fprintf (msg) ; |
|
3139 + posix = input (': ') ; |
|
3140 + if (isempty (posix)) |
|
3141 + posix = 1 ; |
|
3142 + end |
|
3143 + if (posix == 2) |
|
3144 + fprintf ('\nNot using POSIX sysconf and times routines.\n') ; |
|
3145 + posix = ' -DNPOSIX' ; |
|
3146 + else |
|
3147 + fprintf ('\nUsing POSIX sysconf and times routines.\n') ; |
|
3148 + posix = '' ; |
|
3149 + end |
|
3150 +end |
|
3151 + |
|
3152 +%------------------------------------------------------------------------------- |
|
3153 +% mex command |
|
3154 +%------------------------------------------------------------------------------- |
|
3155 + |
|
3156 +umfdir = sprintf ('..%sSource%s', filesep, filesep) ; |
|
3157 +amddir = sprintf ('..%s..%sAMD%sSource%s', filesep, filesep, filesep, filesep) ; |
|
3158 +incdir = sprintf ( ... |
|
3159 +' -I..%sInclude -I..%sSource -I..%s..%sAMD%sInclude -I..%s..%sAMD%sSource', ... |
|
3160 +filesep,filesep, filesep, filesep, filesep, filesep, filesep, filesep) ; |
|
3161 + |
|
3162 +mx = sprintf ('mex -inline -O%s%s%s%s', blas, utils, posix, incdir) ; |
|
3163 +msg = [ ... |
|
3164 + '--------------------------------------------------------------\n', ... |
|
3165 + '\nCompile options:\n%s\nNow compiling. Please wait.\n'] ; |
|
3166 +fprintf (msg, mx) ; |
|
3167 + |
|
3168 +% The following is adapted from GNUmakefile |
|
3169 + |
|
3170 +%------------------------------------------------------------------------------- |
|
3171 +% source files |
|
3172 +%------------------------------------------------------------------------------- |
|
3173 + |
|
3174 +% non-user-callable umf_*.[ch] files: |
|
3175 +umfch = { 'assemble', 'blas3_update', ... |
|
3176 + 'build_tuples', 'create_element', ... |
|
3177 + 'dump', 'extend_front', 'garbage_collection', ... |
|
3178 + 'get_memory', 'init_front', 'kernel', ... |
|
3179 + 'kernel_init', 'kernel_wrapup', ... |
|
3180 + 'local_search', 'lsolve', 'ltsolve', ... |
|
3181 + 'mem_alloc_element', 'mem_alloc_head_block', ... |
|
3182 + 'mem_alloc_tail_block', 'mem_free_tail_block', ... |
|
3183 + 'mem_init_memoryspace', ... |
|
3184 + 'report_vector', 'row_search', 'scale_column', ... |
|
3185 + 'set_stats', 'solve', 'symbolic_usage', 'transpose', ... |
|
3186 + 'tuple_lengths', 'usolve', 'utsolve', 'valid_numeric', ... |
|
3187 + 'valid_symbolic', 'grow_front', 'start_front', '2by2', ... |
|
3188 + 'store_lu', 'scale' } ; |
|
3189 + |
|
3190 +% non-user-callable umf_*.[ch] files, int versions only (no real/complex): |
|
3191 +umfint = { 'analyze', 'apply_order', 'colamd', 'free', 'fsize', ... |
|
3192 + 'is_permutation', 'malloc', 'realloc', 'report_perm', ... |
|
3193 + 'singletons' } ; |
|
3194 + |
|
3195 +% non-user-callable and user-callable amd_*.[ch] files (int versions only): |
|
3196 +amd = { 'aat', '1', '2', 'dump', 'postorder', 'post_tree', 'defaults', ... |
|
3197 + 'order', 'control', 'info', 'valid' } ; |
|
3198 + |
|
3199 +% user-callable umfpack_*.[ch] files (real/complex): |
|
3200 +user = { 'col_to_triplet', 'defaults', 'free_numeric', ... |
|
3201 + 'free_symbolic', 'get_numeric', 'get_lunz', ... |
|
3202 + 'get_symbolic', 'numeric', 'qsymbolic', ... |
|
3203 + 'report_control', 'report_info', 'report_matrix', ... |
|
3204 + 'report_numeric', 'report_perm', 'report_status', ... |
|
3205 + 'report_symbolic', 'report_triplet', ... |
|
3206 + 'report_vector', 'solve', 'symbolic', ... |
|
3207 + 'transpose', 'triplet_to_col', 'scale' ... |
|
3208 + 'load_numeric', 'save_numeric', 'load_symbolic', 'save_symbolic' } ; |
|
3209 + |
|
3210 +% user-callable umfpack_*.[ch], only one version |
|
3211 +generic = { 'timer', 'tictoc' } ; |
|
3212 + |
|
3213 +M = cell (0) ; |
|
3214 + |
|
3215 +%------------------------------------------------------------------------------- |
|
3216 +% Create the umfpack and amd mexFunctions for MATLAB (int versions only) |
|
3217 +%------------------------------------------------------------------------------- |
|
3218 + |
|
3219 +for k = 1:length(umfint) |
|
3220 + M = make (M, '%s -DDINT -c %sumf_%s.c', 'umf_%s.%s', 'umf_%s_%s.%s', ... |
|
3221 + mx, umfint {k}, umfint {k}, 'm', obj, umfdir) ; |
|
3222 +end |
|
3223 + |
|
3224 +rules = { [mx ' -DDINT'] , [mx ' -DZINT'] } ; |
|
3225 +kinds = { 'md', 'mz' } ; |
|
3226 + |
|
3227 +for what = 1:2 |
|
3228 + |
|
3229 + rule = rules {what} ; |
|
3230 + kind = kinds {what} ; |
|
3231 + |
|
3232 + M = make (M, '%s -DCONJUGATE_SOLVE -c %sumf_%s.c', 'umf_%s.%s', ... |
|
3233 + 'umf_%s_%s.%s', rule, 'ltsolve', 'lhsolve', kind, obj, umfdir) ; |
|
3234 + |
|
3235 + M = make (M, '%s -DCONJUGATE_SOLVE -c %sumf_%s.c', 'umf_%s.%s', ... |
|
3236 + 'umf_%s_%s.%s', rule, 'utsolve', 'uhsolve', kind, obj, umfdir) ; |
|
3237 + |
|
3238 + M = make (M, '%s -DDO_MAP -c %sumf_%s.c', 'umf_%s.%s', ... |
|
3239 + 'umf_%s_%s_map_nox.%s', rule, 'triplet', 'triplet', kind, obj, umfdir) ; |
|
3240 + |
|
3241 + M = make (M, '%s -DDO_VALUES -c %sumf_%s.c', 'umf_%s.%s', ... |
|
3242 + 'umf_%s_%s_nomap_x.%s', rule, 'triplet', 'triplet', kind, obj, umfdir) ; |
|
3243 + |
|
3244 + M = make (M, '%s -c %sumf_%s.c', 'umf_%s.%s', ... |
|
3245 + 'umf_%s_%s_nomap_nox.%s', rule, 'triplet', 'triplet', kind, obj, ... |
|
3246 + umfdir) ; |
|
3247 + |
|
3248 + M = make (M, '%s -DDO_MAP -DDO_VALUES -c %sumf_%s.c', 'umf_%s.%s', ... |
|
3249 + 'umf_%s_%s_map_x.%s', rule, 'triplet', 'triplet', kind, obj, umfdir) ; |
|
3250 + |
|
3251 + M = make (M, '%s -DFIXQ -c %sumf_%s.c', 'umf_%s.%s', ... |
|
3252 + 'umf_%s_%s_fixq.%s', rule, 'assemble', 'assemble', kind, obj, umfdir) ; |
|
3253 + |
|
3254 + M = make (M, '%s -DDROP -c %sumf_%s.c', 'umf_%s.%s', ... |
|
3255 + 'umf_%s_%s_drop.%s', rule, 'store_lu', 'store_lu', kind, obj, umfdir) ; |
|
3256 + |
|
3257 + for k = 1:length(umfch) |
|
3258 + M = make (M, '%s -c %sumf_%s.c', 'umf_%s.%s', 'umf_%s_%s.%s', ... |
|
3259 + rule, umfch {k}, umfch {k}, kind, obj, umfdir) ; |
|
3260 + end |
|
3261 + |
|
3262 + M = make (M, '%s -DWSOLVE -c %sumfpack_%s.c', 'umfpack_%s.%s', ... |
|
3263 + 'umfpack_%s_w%s.%s', rule, 'solve', 'solve', kind, obj, umfdir) ; |
|
3264 + |
|
3265 + for k = 1:length(user) |
|
3266 + M = make (M, '%s -c %sumfpack_%s.c', 'umfpack_%s.%s', ... |
|
3267 + 'umfpack_%s_%s.%s', rule, user {k}, user {k}, kind, obj, umfdir) ; |
|
3268 + end |
|
3269 +end |
|
3270 + |
|
3271 +for k = 1:length(generic) |
|
3272 + M = make (M, '%s -c %sumfpack_%s.c', 'umfpack_%s.%s', ... |
|
3273 + 'umfpack_%s_%s.%s', mx, generic {k}, generic {k}, 'm', obj, umfdir) ; |
|
3274 +end |
|
3275 + |
|
3276 +%---------------------------------------- |
|
3277 +% AMD routines (int only) |
|
3278 +%---------------------------------------- |
|
3279 + |
|
3280 +for k = 1:length(amd) |
|
3281 + M = make (M, '%s -DDINT -c %samd_%s.c', 'amd_%s.%s', 'amd_%s_%s.%s', ... |
|
3282 + mx, amd {k}, amd {k}, 'm', obj, amddir) ; |
|
3283 +end |
|
3284 + |
|
3285 +%---------------------------------------- |
|
3286 +% compile the umfpack mexFunction |
|
3287 +%---------------------------------------- |
|
3288 + |
|
3289 +C = sprintf ('%s -output umfpack umfpackmex.c', mx) ; |
|
3290 +for i = 1:length (M) |
|
3291 + C = [C ' ' (M {i})] ; |
|
3292 +end |
|
3293 +C = [C ' ' blas_lib] ; |
|
3294 +cmd (C) ; |
|
3295 + |
|
3296 +%---------------------------------------- |
|
3297 +% delete the object files |
|
3298 +%---------------------------------------- |
|
3299 + |
|
3300 +for i = 1:length (M) |
|
3301 + rmfile (M {i}) ; |
|
3302 +end |
|
3303 + |
|
3304 +%---------------------------------------- |
|
3305 +% compile the luflop mexFunction |
|
3306 +%---------------------------------------- |
|
3307 + |
|
3308 +cmd (sprintf ('%s -output luflop luflopmex.c', mx)) ; |
|
3309 + |
|
3310 +fprintf ('\n\nCompilation has completed. Now trying the umfpack_simple demo.\n'); |
|
3311 +umfpack_simple |
|
3312 + |
|
3313 +%------------------------------------------------------------------------------- |
|
3314 +% rmfile: delete a file, but only if it exists |
|
3315 +%------------------------------------------------------------------------------- |
|
3316 + |
|
3317 +function rmfile (file) |
|
3318 +if (length (dir (file)) > 0) |
|
3319 + delete (file) ; |
|
3320 +end |
|
3321 + |
|
3322 +%------------------------------------------------------------------------------- |
|
3323 +% cpfile: copy the src file to the filename dst, overwriting dst if it exists |
|
3324 +%------------------------------------------------------------------------------- |
|
3325 + |
|
3326 +function cpfile (src, dst) |
|
3327 +rmfile (dst) |
|
3328 +if (length (dir (src)) == 0) |
|
3329 + help umfpack_make |
|
3330 + error (sprintf ('File does not exist: %s\n', src)) ; |
|
3331 +end |
|
3332 +copyfile (src, dst) ; |
|
3333 + |
|
3334 +%------------------------------------------------------------------------------- |
|
3335 +% mvfile: move the src file to the filename dst, overwriting dst if it exists |
|
3336 +%------------------------------------------------------------------------------- |
|
3337 + |
|
3338 +function mvfile (src, dst) |
|
3339 +cpfile (src, dst) ; |
|
3340 +rmfile (src) ; |
|
3341 + |
|
3342 +%------------------------------------------------------------------------------- |
|
3343 +% cmd: display and execute a command |
|
3344 +%------------------------------------------------------------------------------- |
|
3345 + |
|
3346 +function cmd (s) |
|
3347 +fprintf ('.') ; |
|
3348 +eval (s) ; |
|
3349 + |
|
3350 +%------------------------------------------------------------------------------- |
|
3351 +% make: execute a "make" command for a source file |
|
3352 +%------------------------------------------------------------------------------- |
|
3353 + |
|
3354 +function M = make (M, s, src, dst, rule, file1, file2, kind, obj, srcdir) |
|
3355 +cmd (sprintf (s, rule, srcdir, file1)) ; |
|
3356 +src = sprintf (src, file1, obj) ; |
|
3357 +dst = sprintf (dst, kind, file2, obj) ; |
|
3358 +mvfile (src, dst) ; |
|
3359 +M {end + 1} = dst ; |
|
3360 + |
|
3361 diff -uNr UMFPACKv4.4.orig/UMFPACK/OCTAVE/umfpack_report.m UMFPACK/UMFPACK/OCTAVE/umfpack_report.m |
|
3362 --- UMFPACKv4.4.orig/UMFPACK/OCTAVE/umfpack_report.m 1970-01-01 01:00:00.000000000 +0100 |
|
3363 +++ UMFPACK/UMFPACK/OCTAVE/umfpack_report.m 2004-12-30 01:58:46.000000000 +0100 |
|
3364 @@ -0,0 +1,346 @@ |
|
3365 +function umfpack_report (Control, Info) |
|
3366 +% UMFPACK_REPORT |
|
3367 +% |
|
3368 +% umfpack_report (Control, Info) ; |
|
3369 +% |
|
3370 +% Prints the current Control settings for umfpack, and the statistical |
|
3371 +% information returned by umfpack in the Info array. If Control is |
|
3372 +% an empty matrix, then the default control settings are printed. |
|
3373 +% |
|
3374 +% Control is 20-by-1, and Info is 90-by-1. Not all entries are used. |
|
3375 +% |
|
3376 +% Alternative usages: |
|
3377 +% |
|
3378 +% umfpack_report ([ ], Info) ; print the default control parameters |
|
3379 +% and the Info array. |
|
3380 +% umfpack_report (Control) ; print the control parameters only. |
|
3381 +% umfpack_report ; print the default control parameters |
|
3382 +% and an empty Info array. |
|
3383 +% |
|
3384 +% See also umfpack, umfpack_make, umfpack_details, |
|
3385 +% umfpack_demo, and umfpack_simple. |
|
3386 + |
|
3387 +% UMFPACK Version 4.3 (Jan. 16, 2004), Copyright (c) 2004 by Timothy A. |
|
3388 +% Davis. All Rights Reserved. See ../README for License. |
|
3389 + |
|
3390 +%------------------------------------------------------------------------------- |
|
3391 +% get inputs, use defaults if input arguments not present |
|
3392 +%------------------------------------------------------------------------------- |
|
3393 + |
|
3394 +% The contents of Control and Info are defined in umfpack.h |
|
3395 +if (nargin < 1) |
|
3396 + Control = [] ; |
|
3397 +end |
|
3398 +if (nargin < 2) |
|
3399 + Info = [] ; |
|
3400 +end |
|
3401 +if (isempty (Control)) |
|
3402 + Control = umfpack ; |
|
3403 +end |
|
3404 +if (isempty (Info)) |
|
3405 + Info = [ 0 (-ones (1, 89)) ] ; |
|
3406 +end |
|
3407 + |
|
3408 +%------------------------------------------------------------------------------- |
|
3409 +% control settings |
|
3410 +%------------------------------------------------------------------------------- |
|
3411 + |
|
3412 +fprintf ('\nUMFPACK Version 4.3: Control settings:\n\n') ; |
|
3413 +fprintf (' Control (1): print level: %d\n', Control (1)) ; |
|
3414 +fprintf (' Control (2): dense row parameter: %g\n', Control (2)) ; |
|
3415 +fprintf (' "dense" rows have > max (16, (%g)*16*sqrt(n_col)) entries\n', Control (2)) ; |
|
3416 +fprintf (' Control (3): dense column parameter: %g\n', Control (3)) ; |
|
3417 +fprintf (' "dense" columns have > max (16, (%g)*16*sqrt(n_row)) entries\n', Control (3)) ; |
|
3418 +fprintf (' Control (4): pivot tolerance: %g\n', Control (4)) ; |
|
3419 +fprintf (' Control (5): max block size for dense matrix kernels: %d\n', Control (5)) ; |
|
3420 +prstrat (' Control (6): strategy: %g ', Control (6)) ; |
|
3421 +fprintf (' Control (7): initial allocation ratio: %g\n', Control (7)) ; |
|
3422 +fprintf (' Control (8): max iterative refinement steps: %d\n', Control (8)) ; |
|
3423 +fprintf (' Control (13): 2-by-2 pivot tolerance: %g\n', Control (13)) ; |
|
3424 +fprintf (' Control (14): Q fixed during numeric factorization: %g ', Control (14)) ; |
|
3425 +if (Control (14) > 0) |
|
3426 + fprintf ('(yes)\n') ; |
|
3427 +elseif (Control (14) < 0) |
|
3428 + fprintf ('(no)\n') ; |
|
3429 +else |
|
3430 + fprintf ('(auto)\n') ; |
|
3431 +end |
|
3432 +fprintf (' Control (15): AMD dense row/column parameter: %g\n', Control (15)) ; |
|
3433 +fprintf (' "dense" rows/columns in A+A'' have > max (16, (%g)*sqrt(n)) entries.\n', Control (15)) ; |
|
3434 +fprintf (' Only used if the AMD ordering is used.\n') ; |
|
3435 +fprintf (' Control (16): diagonal pivot tolerance: %g\n', Control (16)) ; |
|
3436 +fprintf (' Only used if diagonal pivoting is attempted.\n') ; |
|
3437 + |
|
3438 +fprintf (' Control (17): scaling option: %g ', Control (17)) ; |
|
3439 +if (Control (17) == 0) |
|
3440 + fprintf ('(none)\n') ; |
|
3441 +elseif (Control (17) == 2) |
|
3442 + fprintf ('(scale the matrix by\n') ; |
|
3443 + fprintf (' dividing each row by max. abs. value in each row)\n') ; |
|
3444 +else |
|
3445 + fprintf ('(scale the matrix by\n') ; |
|
3446 + fprintf (' dividing each row by sum of abs. values in each row)\n') ; |
|
3447 +end |
|
3448 + |
|
3449 +fprintf (' Control (18): frontal matrix allocation ratio: %g\n', Control (18)) ; |
|
3450 +fprintf (' Control (19): drop tolerance: %g\n', Control (19)) ; |
|
3451 +fprintf (' Control (20): AMD and COLAMD aggressive absorption: %g ', Control (20)) ; |
|
3452 +yes_no (Control (20)) ; |
|
3453 + |
|
3454 +% compile-time options: |
|
3455 + |
|
3456 +fprintf ('\n The following options can only be changed at compile-time:\n') ; |
|
3457 + |
|
3458 +if (Control (9) == 1) |
|
3459 + fprintf (' Control (9): compiled to use the BLAS\n') ; |
|
3460 +else |
|
3461 + fprintf (' Control (9): compiled without the BLAS\n') ; |
|
3462 + fprintf (' (you will not get the best possible performance)\n') ; |
|
3463 +end |
|
3464 + |
|
3465 +if (Control (10) == 1) |
|
3466 + fprintf (' Control (10): compiled for MATLAB\n') ; |
|
3467 +elseif (Control (10) == 2) |
|
3468 + fprintf (' Control (10): compiled for MATLAB\n') ; |
|
3469 + fprintf (' Uses internal utMalloc, utFree, utRealloc, utPrintf\n') ; |
|
3470 + fprintf (' utDivideComplex, and utFdlibm_hypot routines.\n') ; |
|
3471 +else |
|
3472 + fprintf (' Control (10): not compiled for MATLAB\n') ; |
|
3473 + fprintf (' Uses ANSI C malloc, free, realloc, and printf\n') ; |
|
3474 + fprintf (' instead of mxMalloc, mxFree, mxRealloc, and mexPrintf.\n') ; |
|
3475 + fprintf (' Printing will be in terms of 0-based matrix indexing,\n') ; |
|
3476 + fprintf (' not 1-based as is expected in MATLAB. Diary output may\n') ; |
|
3477 + fprintf (' not be properly recorded.\n') ; |
|
3478 +end |
|
3479 + |
|
3480 +if (Control (11) == 2) |
|
3481 + fprintf (' Control (11): uses POSIX times ( ) to get CPU time and wallclock time.\n') ; |
|
3482 +elseif (Control (11) == 1) |
|
3483 + fprintf (' Control (11): uses getrusage to get CPU time.\n') ; |
|
3484 +else |
|
3485 + fprintf (' Control (11): uses ANSI C clock to get CPU time.\n') ; |
|
3486 + fprintf (' The CPU time may wrap around, type "help cputime".\n') ; |
|
3487 +end |
|
3488 + |
|
3489 +if (Control (12) == 1) |
|
3490 + fprintf (' Control (12): compiled with debugging enabled\n') ; |
|
3491 + fprintf (' ###########################################\n') ; |
|
3492 + fprintf (' ### This will be exceedingly slow! ########\n') ; |
|
3493 + fprintf (' ###########################################\n') ; |
|
3494 + if (Control (10) == 1) |
|
3495 + fprintf (' Uses mxAssert.\n') ; |
|
3496 + elseif (Control (10) == 2) |
|
3497 + fprintf (' Uses utAssert.\n') ; |
|
3498 + else |
|
3499 + fprintf (' Uses ANSI C assert instead of mxAssert.\n') ; |
|
3500 + end |
|
3501 +else |
|
3502 + fprintf (' Control (12): compiled for normal operation (no debugging)\n') ; |
|
3503 +end |
|
3504 + |
|
3505 +%------------------------------------------------------------------------------- |
|
3506 +% Info: |
|
3507 +%------------------------------------------------------------------------------- |
|
3508 + |
|
3509 +if (nargin == 1) |
|
3510 + return |
|
3511 +end |
|
3512 + |
|
3513 +status = Info (1) ; |
|
3514 +fprintf ('\nUMFPACK status: Info (1): %d, ', status) ; |
|
3515 + |
|
3516 +if (status == 0) |
|
3517 + fprintf ('OK\n') ; |
|
3518 +elseif (status == 1) |
|
3519 + fprintf ('WARNING matrix is singular\n') ; |
|
3520 +elseif (status == -1) |
|
3521 + fprintf ('ERROR out of memory\n') ; |
|
3522 +elseif (status == -3) |
|
3523 + fprintf ('ERROR numeric LU factorization is invalid\n') ; |
|
3524 +elseif (status == -4) |
|
3525 + fprintf ('ERROR symbolic LU factorization is invalid\n') ; |
|
3526 +elseif (status == -5) |
|
3527 + fprintf ('ERROR required argument is missing\n') ; |
|
3528 +elseif (status == -6) |
|
3529 + fprintf ('ERROR n <= 0\n') ; |
|
3530 +elseif (status <= -7 & status >= -12 | status == -14) |
|
3531 + fprintf ('ERROR matrix A is corrupted\n') ; |
|
3532 +elseif (status == -13) |
|
3533 + fprintf ('ERROR invalid system\n') ; |
|
3534 +elseif (status == -15) |
|
3535 + fprintf ('ERROR invalid permutation\n') ; |
|
3536 +elseif (status == -911) |
|
3537 + fprintf ('ERROR internal error!\n') ; |
|
3538 + fprintf ('Please report this error to Tim Davis (davis@cise.ufl.edu)\n') ; |
|
3539 +else |
|
3540 + fprintf ('ERROR unrecognized error. Info array corrupted\n') ; |
|
3541 +end |
|
3542 + |
|
3543 +fprintf (' (a -1 means the entry has not been computed):\n') ; |
|
3544 + |
|
3545 +fprintf ('\n Basic statistics:\n') ; |
|
3546 +fprintf (' Info (2): %d, # of rows of A\n', Info (2)) ; |
|
3547 +fprintf (' Info (17): %d, # of columns of A\n', Info (17)) ; |
|
3548 +fprintf (' Info (3): %d, nnz (A)\n', Info (3)) ; |
|
3549 +fprintf (' Info (4): %d, Unit size, in bytes, for memory usage reported below\n', Info (4)) ; |
|
3550 +fprintf (' Info (5): %d, size of int (in bytes)\n', Info (5)) ; |
|
3551 +fprintf (' Info (6): %d, size of long (in bytes)\n', Info (6)) ; |
|
3552 +fprintf (' Info (7): %d, size of pointer (in bytes)\n', Info (7)) ; |
|
3553 +fprintf (' Info (8): %d, size of numerical entry (in bytes)\n', Info (8)) ; |
|
3554 + |
|
3555 +fprintf ('\n Pivots with zero Markowitz cost removed to obtain submatrix S:\n') ; |
|
3556 +fprintf (' Info (57): %d, # of pivots with one entry in pivot column\n', Info (57)) ; |
|
3557 +fprintf (' Info (58): %d, # of pivots with one entry in pivot row\n', Info (58)) ; |
|
3558 +fprintf (' Info (59): %d, # of rows/columns in submatrix S (if square)\n', Info (59)) ; |
|
3559 +fprintf (' Info (60): %d ') ; |
|
3560 +if (Info (60) > 0) |
|
3561 + fprintf ('submatrix S square and diagonal preserved\n') ; |
|
3562 +elseif (Info (60) == 0) |
|
3563 + fprintf ('submatrix S not square or diagonal not preserved\n') ; |
|
3564 +else |
|
3565 + fprintf ('\n') ; |
|
3566 +end |
|
3567 +fprintf (' Info (9): %d, # of "dense" rows in S\n', Info (9)) ; |
|
3568 +fprintf (' Info (10): %d, # of empty rows in S\n', Info (10)) ; |
|
3569 +fprintf (' Info (11): %d, # of "dense" columns in S\n', Info (11)) ; |
|
3570 +fprintf (' Info (12): %d, # of empty columns in S\n', Info (12)) ; |
|
3571 +fprintf (' Info (34): %g, symmetry of pattern of S\n', Info (34)) ; |
|
3572 +fprintf (' Info (35): %d, # of off-diagonal nonzeros in S+S''\n', Info (35)) ; |
|
3573 +fprintf (' Info (36): %d, nnz (diag (S))\n', Info (36)) ; |
|
3574 + |
|
3575 +fprintf ('\n 2-by-2 pivoting to place large entries on diagonal:\n') ; |
|
3576 +fprintf (' Info (52): %d, # of small diagonal entries of S\n', Info (52)) ; |
|
3577 +fprintf (' Info (53): %d, # of unmatched small diagonal entries\n', Info (53)) ; |
|
3578 +fprintf (' Info (54): %g, symmetry of P2*S\n', Info (54)) ; |
|
3579 +fprintf (' Info (55): %d, # of off-diagonal entries in (P2*S)+(P2*S)''\n', Info (55)) ; |
|
3580 +fprintf (' Info (56): %d, nnz (diag (P2*S))\n', Info (56)) ; |
|
3581 + |
|
3582 +fprintf ('\n AMD results, for strict diagonal pivoting:\n') ; |
|
3583 +fprintf (' Info (37): %d, est. nz in L and U\n', Info (37)) ; |
|
3584 +fprintf (' Info (38): %g, est. flop count\n', Info (38)) ; |
|
3585 +fprintf (' Info (39): %g, # of "dense" rows in S+S''\n', Info (39)) ; |
|
3586 +fprintf (' Info (40): %g, est. max. nz in any column of L\n', Info (40)) ; |
|
3587 + |
|
3588 +fprintf ('\n Final strategy selection, based on the analysis above:\n') ; |
|
3589 +prstrat (' Info (19): %d, strategy used ', Info (19)) ; |
|
3590 +fprintf (' Info (20): %d, ordering used ', Info (20)) ; |
|
3591 +if (Info (20) == 0) |
|
3592 + fprintf ('(COLAMD on A)\n') ; |
|
3593 +elseif (Info (20) == 1) |
|
3594 + fprintf ('(AMD on A+A'')\n') ; |
|
3595 +elseif (Info (20) == 2) |
|
3596 + fprintf ('(provided by user)\n') ; |
|
3597 +else |
|
3598 + fprintf ('(undefined ordering option)\n') ; |
|
3599 +end |
|
3600 +fprintf (' Info (32): %d, Q fixed during numeric factorization: ', Info (32)) ; |
|
3601 +yes_no (Info (32)) ; |
|
3602 +fprintf (' Info (33): %d, prefer diagonal pivoting: ', Info (33)) ; |
|
3603 +yes_no (Info (33)) ; |
|
3604 + |
|
3605 +fprintf ('\n symbolic analysis time and memory usage:\n') ; |
|
3606 +fprintf (' Info (13): %d, defragmentations during symbolic analysis\n', Info (13)) ; |
|
3607 +fprintf (' Info (14): %d, memory used during symbolic analysis (Units)\n', Info (14)) ; |
|
3608 +fprintf (' Info (15): %d, final size of symbolic factors (Units)\n', Info (15)) ; |
|
3609 +fprintf (' Info (16): %.2f, symbolic analysis CPU time (seconds)\n', Info (16)) ; |
|
3610 +fprintf (' Info (18): %.2f, symbolic analysis wall clock time (seconds)\n', Info (18)) ; |
|
3611 + |
|
3612 +fprintf ('\n Estimates computed in the symbolic analysis:\n') ; |
|
3613 +fprintf (' Info (21): %d, est. size of LU factors (Units)\n', Info (21)) ; |
|
3614 +fprintf (' Info (22): %d, est. total peak memory usage (Units)\n', Info (22)) ; |
|
3615 +fprintf (' Info (23): %d, est. factorization flop count\n', Info (23)) ; |
|
3616 +fprintf (' Info (24): %d, est. nnz (L)\n', Info (24)) ; |
|
3617 +fprintf (' Info (25): %d, est. nnz (U)\n', Info (25)) ; |
|
3618 +fprintf (' Info (26): %d, est. initial size, variable-part of LU (Units)\n', Info (26)) ; |
|
3619 +fprintf (' Info (27): %d, est. peak size, of variable-part of LU (Units)\n', Info (27)) ; |
|
3620 +fprintf (' Info (28): %d, est. final size, of variable-part of LU (Units)\n', Info (28)) ; |
|
3621 +fprintf (' Info (29): %d, est. max frontal matrix size (# of entries)\n', Info (29)) ; |
|
3622 +fprintf (' Info (30): %d, est. max # of rows in frontal matrix\n', Info (30)) ; |
|
3623 +fprintf (' Info (31): %d, est. max # of columns in frontal matrix\n', Info (31)) ; |
|
3624 + |
|
3625 +fprintf ('\n Computed in the numeric factorization (estimates shown above):\n') ; |
|
3626 +fprintf (' Info (41): %d, size of LU factors (Units)\n', Info (41)) ; |
|
3627 +fprintf (' Info (42): %d, total peak memory usage (Units)\n', Info (42)) ; |
|
3628 +fprintf (' Info (43): %d, factorization flop count\n', Info (43)) ; |
|
3629 +fprintf (' Info (44): %d, nnz (L)\n', Info (44)) ; |
|
3630 +fprintf (' Info (45): %d, nnz (U)\n', Info (45)) ; |
|
3631 +fprintf (' Info (46): %d, initial size of variable-part of LU (Units)\n', Info (46)) ; |
|
3632 +fprintf (' Info (47): %d, peak size of variable-part of LU (Units)\n', Info (47)) ; |
|
3633 +fprintf (' Info (48): %d, final size of variable-part of LU (Units)\n', Info (48)) ; |
|
3634 +fprintf (' Info (49): %d, max frontal matrix size (# of numerical entries)\n', Info (49)) ; |
|
3635 +fprintf (' Info (50): %d, max # of rows in frontal matrix\n', Info (50)) ; |
|
3636 +fprintf (' Info (51): %d, max # of columns in frontal matrix\n', Info (51)) ; |
|
3637 + |
|
3638 +fprintf ('\n Computed in the numeric factorization (no estimates computed a priori):\n') ; |
|
3639 +fprintf (' Info (61): %d, defragmentations during numeric factorization\n', Info (61)) ; |
|
3640 +fprintf (' Info (62): %d, reallocations during numeric factorization\n', Info (62)) ; |
|
3641 +fprintf (' Info (63): %d, costly reallocations during numeric factorization\n', Info (63)) ; |
|
3642 +fprintf (' Info (64): %d, integer indices in compressed pattern of L and U\n', Info (64)) ; |
|
3643 +fprintf (' Info (65): %d, numerical values stored in L and U\n', Info (65)) ; |
|
3644 +fprintf (' Info (66): %.2f, numeric factorization CPU time (seconds)\n', Info (66)) ; |
|
3645 +fprintf (' Info (76): %.2f, numeric factorization wall clock time (seconds)\n', Info (76)) ; |
|
3646 +if (Info (66) > 0.05 & Info (43) > 0) |
|
3647 +fprintf (' mflops in numeric factorization phase: %.2f\n', 1e-6 * Info (43) / Info (66)) ; |
|
3648 +end |
|
3649 +fprintf (' Info (67): %d, nnz (diag (U))\n', Info (67)) ; |
|
3650 +fprintf (' Info (68): %g, reciprocal condition number estimate\n', Info (68)) ; |
|
3651 +fprintf (' Info (69): %g, matrix was ', Info (69)) ; |
|
3652 +if (Info (69) == 0) |
|
3653 + fprintf ('not scaled\n') ; |
|
3654 +elseif (Info (69) == 2) |
|
3655 + fprintf ('scaled (row max)\n') ; |
|
3656 +else |
|
3657 + fprintf ('scaled (row sum)\n') ; |
|
3658 +end |
|
3659 +fprintf (' Info (70): %g, min. scale factor of rows of A\n', Info (70)) ; |
|
3660 +fprintf (' Info (71): %g, max. scale factor of rows of A\n', Info (71)) ; |
|
3661 +fprintf (' Info (72): %g, min. abs. on diagonal of U\n', Info (72)) ; |
|
3662 +fprintf (' Info (73): %g, max. abs. on diagonal of U\n', Info (73)) ; |
|
3663 +fprintf (' Info (74): %g, initial allocation parameter used\n', Info (74)) ; |
|
3664 +fprintf (' Info (75): %g, # of forced updates due to frontal growth\n', Info (75)) ; |
|
3665 +fprintf (' Info (77): %d, # of off-diaogonal pivots\n', Info (77)) ; |
|
3666 +fprintf (' Info (78): %d, nnz (L), if no small entries dropped\n', Info (78)) ; |
|
3667 +fprintf (' Info (79): %d, nnz (U), if no small entries dropped\n', Info (79)) ; |
|
3668 +fprintf (' Info (80): %d, # of small entries dropped\n', Info (80)) ; |
|
3669 + |
|
3670 +fprintf ('\n Computed in the solve step:\n') ; |
|
3671 +fprintf (' Info (81): %d, iterative refinement steps taken\n', Info (81)) ; |
|
3672 +fprintf (' Info (82): %d, iterative refinement steps attempted\n', Info (82)) ; |
|
3673 +fprintf (' Info (83): %g, omega(1), sparse-backward error estimate\n', Info (83)) ; |
|
3674 +fprintf (' Info (84): %g, omega(2), sparse-backward error estimate\n', Info (84)) ; |
|
3675 +fprintf (' Info (85): %d, solve flop count\n', Info (85)) ; |
|
3676 +fprintf (' Info (86): %.2f, solve CPU time (seconds)\n', Info (86)) ; |
|
3677 +fprintf (' Info (87): %.2f, solve wall clock time (seconds)\n', Info (87)) ; |
|
3678 + |
|
3679 +fprintf ('\n Info (88:90): unused\n\n') ; |
|
3680 + |
|
3681 +%------------------------------------------------------------------------------- |
|
3682 + |
|
3683 +function prstrat (fmt, strategy) |
|
3684 +fprintf (fmt, strategy) ; |
|
3685 +if (strategy == 1) |
|
3686 + fprintf ('(unsymmetric)\n') ; |
|
3687 + fprintf (' Q = COLAMD (A), Q refined during numerical\n') ; |
|
3688 + fprintf (' factorization, and no attempt at diagonal pivoting.\n') ; |
|
3689 +elseif (strategy == 2) |
|
3690 + fprintf ('(symmetric, with 2-by-2 pivoting)\n') ; |
|
3691 + fprintf (' P2 = row permutation to place large values on the diagonal\n') ; |
|
3692 + fprintf (' Q = AMD (P2*A+(P2*A)''), Q not refined during numeric factorization,\n') ; |
|
3693 + fprintf (' and diagonal pivoting attempted.\n') ; |
|
3694 +elseif (strategy == 3) |
|
3695 + fprintf ('(symmetric)\n') ; |
|
3696 + fprintf (' Q = AMD (A+A''), Q not refined during numeric factorization,\n') ; |
|
3697 + fprintf (' and diagonal pivoting (P=Q'') attempted.\n') ; |
|
3698 +else |
|
3699 + strategy = 0 ; |
|
3700 + fprintf ('(auto)\n') ; |
|
3701 +end |
|
3702 + |
|
3703 +%------------------------------------------------------------------------------- |
|
3704 + |
|
3705 +function yes_no (s) |
|
3706 +if (s == 0) |
|
3707 + fprintf ('(no)\n') ; |
|
3708 +else |
|
3709 + fprintf ('(yes)\n') ; |
|
3710 +end |
|
3711 diff -uNr UMFPACKv4.4.orig/UMFPACK/OCTAVE/umfpack_simple.m UMFPACK/UMFPACK/OCTAVE/umfpack_simple.m |
|
3712 --- UMFPACKv4.4.orig/UMFPACK/OCTAVE/umfpack_simple.m 1970-01-01 01:00:00.000000000 +0100 |
|
3713 +++ UMFPACK/UMFPACK/OCTAVE/umfpack_simple.m 2004-12-30 01:58:46.000000000 +0100 |
|
3714 @@ -0,0 +1,61 @@ |
|
3715 +% umfpack_simple: a simple demo of UMFPACK |
|
3716 +% |
|
3717 +% UMFPACK Version 4.3 (Jan. 16, 2004), Copyright (c) 2004 by Timothy A. |
|
3718 +% Davis. All Rights Reserved. |
|
3719 +% |
|
3720 +% UMFPACK License: |
|
3721 +% |
|
3722 +% Your use or distribution of UMFPACK or any modified version of |
|
3723 +% UMFPACK implies that you agree to this License. |
|
3724 +% |
|
3725 +% THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY |
|
3726 +% EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK. |
|
3727 +% |
|
3728 +% Permission is hereby granted to use or copy this program, provided |
|
3729 +% that the Copyright, this License, and the Availability of the original |
|
3730 +% version is retained on all copies. User documentation of any code that |
|
3731 +% uses UMFPACK or any modified version of UMFPACK code must cite the |
|
3732 +% Copyright, this License, the Availability note, and "Used by permission." |
|
3733 +% Permission to modify the code and to distribute modified code is granted, |
|
3734 +% provided the Copyright, this License, and the Availability note are |
|
3735 +% retained, and a notice that the code was modified is included. This |
|
3736 +% software was developed with support from the National Science Foundation, |
|
3737 +% and is provided to you free of charge. |
|
3738 +% |
|
3739 +% Availability: http://www.cise.ufl.edu/research/sparse/umfpack |
|
3740 +% |
|
3741 +% See also: umfpack, umfpack_details |
|
3742 + |
|
3743 +help umfpack_simple |
|
3744 +i = input ('Hit enter to agree to the above License: ', 's') ; |
|
3745 +if (~isempty (i)) |
|
3746 + error ('terminating') ; |
|
3747 +end |
|
3748 + |
|
3749 +format short |
|
3750 + |
|
3751 +A = [ |
|
3752 + 2 3 0 0 0 |
|
3753 + 3 0 4 0 6 |
|
3754 + 0 -1 -3 2 0 |
|
3755 + 0 0 1 0 0 |
|
3756 + 0 4 2 0 1 |
|
3757 +] |
|
3758 + |
|
3759 +A = sparse (A) ; |
|
3760 + |
|
3761 +b = [8 45 -3 3 19]' |
|
3762 + |
|
3763 +fprintf ('Solution to Ax=b via UMFPACK:\n') ; |
|
3764 +fprintf ('x1 = umfpack (A, ''\\'', b)\n') ; |
|
3765 + |
|
3766 +x1 = umfpack (A, '\\', b) |
|
3767 + |
|
3768 +fprintf ('Solution to Ax=b via OCTAVE:\n') ; |
|
3769 +fprintf ('x2 = A\\b\n') ; |
|
3770 + |
|
3771 +x2 = A\b |
|
3772 + |
|
3773 +fprintf ('norm (x1-x2) should be small: %g\n', norm (x1-x2)) ; |
|
3774 + |
|
3775 +fprintf ('Type ''umfpack_demo'' for a full demo of UMFPACK\n') ; |
|
3776 diff -uNr UMFPACKv4.4.orig/UMFPACK/OCTAVE/umfpack_simple.m.out UMFPACK/UMFPACK/OCTAVE/umfpack_simple.m.out |
|
3777 --- UMFPACKv4.4.orig/UMFPACK/OCTAVE/umfpack_simple.m.out 1970-01-01 01:00:00.000000000 +0100 |
|
3778 +++ UMFPACK/UMFPACK/OCTAVE/umfpack_simple.m.out 2004-12-30 01:58:46.000000000 +0100 |
|
3779 @@ -0,0 +1,79 @@ |
|
3780 +octave:4> umfpack_simple |
|
3781 +umfpack_simple is the file: /home/dbateman/octave/devel/octave-sparse/UMFPACKv4.3/UMFPACK/OCTAVE/umfpack_simple.m |
|
3782 + |
|
3783 +umfpack_simple: a simple demo of UMFPACK |
|
3784 + |
|
3785 +UMFPACK Version 4.3 (Jan. 16, 2004), Copyright (c) 2004 by Timothy A. |
|
3786 +Davis. All Rights Reserved. |
|
3787 + |
|
3788 +UMFPACK License: |
|
3789 + |
|
3790 + Your use or distribution of UMFPACK or any modified version of |
|
3791 + UMFPACK implies that you agree to this License. |
|
3792 + |
|
3793 + THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY |
|
3794 + EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK. |
|
3795 + |
|
3796 + Permission is hereby granted to use or copy this program, provided |
|
3797 + that the Copyright, this License, and the Availability of the original |
|
3798 + version is retained on all copies. User documentation of any code that |
|
3799 + uses UMFPACK or any modified version of UMFPACK code must cite the |
|
3800 + Copyright, this License, the Availability note, and "Used by permission." |
|
3801 + Permission to modify the code and to distribute modified code is granted, |
|
3802 + provided the Copyright, this License, and the Availability note are |
|
3803 + retained, and a notice that the code was modified is included. This |
|
3804 + software was developed with support from the National Science Foundation, |
|
3805 + and is provided to you free of charge. |
|
3806 + |
|
3807 +Availability: http://www.cise.ufl.edu/research/sparse/umfpack |
|
3808 + |
|
3809 +See also: umfpack, umfpack_details |
|
3810 + |
|
3811 + |
|
3812 +Additional help for built-in functions, operators, and variables |
|
3813 +is available in the on-line version of the manual. Use the command |
|
3814 +`help -i <topic>' to search the manual index. |
|
3815 + |
|
3816 +Help and information about Octave is also available on the WWW |
|
3817 +at http://www.octave.org and via the help@octave.org |
|
3818 +mailing list. |
|
3819 +Hit enter to agree to the above License: |
|
3820 +A = |
|
3821 + |
|
3822 + 2 3 0 0 0 |
|
3823 + 3 0 4 0 6 |
|
3824 + 0 -1 -3 2 0 |
|
3825 + 0 0 1 0 0 |
|
3826 + 0 4 2 0 1 |
|
3827 + |
|
3828 +b = |
|
3829 + |
|
3830 + 8 |
|
3831 + 45 |
|
3832 + -3 |
|
3833 + 3 |
|
3834 + 19 |
|
3835 + |
|
3836 +Solution to Ax=b via UMFPACK: |
|
3837 +x1 = umfpack (A, '\', b) |
|
3838 +x1 = |
|
3839 + |
|
3840 + 1.00000 |
|
3841 + 2.00000 |
|
3842 + 3.00000 |
|
3843 + 4.00000 |
|
3844 + 5.00000 |
|
3845 + |
|
3846 +Solution to Ax=b via OCTAVE: |
|
3847 +x2 = A\b |
|
3848 +x2 = |
|
3849 + |
|
3850 + 1.00000 |
|
3851 + 2.00000 |
|
3852 + 3.00000 |
|
3853 + 4.00000 |
|
3854 + 5.00000 |
|
3855 + |
|
3856 +norm (x1-x2) should be small: 0 |
|
3857 +Type 'umfpack_demo' for a full demo of UMFPACK |
|
3858 +octave:5> diary off |
|
3859 diff -uNr UMFPACKv4.4.orig/UMFPACK/OCTAVE/umfpack_solve.m UMFPACK/UMFPACK/OCTAVE/umfpack_solve.m |
|
3860 --- UMFPACKv4.4.orig/UMFPACK/OCTAVE/umfpack_solve.m 1970-01-01 01:00:00.000000000 +0100 |
|
3861 +++ UMFPACK/UMFPACK/OCTAVE/umfpack_solve.m 2004-12-30 01:58:46.000000000 +0100 |
|
3862 @@ -0,0 +1,97 @@ |
|
3863 +function x = umfpack_solve (arg1, op, arg2, Control) |
|
3864 +% UMFPACK_SOLVE |
|
3865 +% |
|
3866 +% x = umfpack_solve (A, '\', b, Control) |
|
3867 +% x = umfpack_solve (b, '/', A, Control) |
|
3868 +% |
|
3869 +% Computes x = A\b, or b/A, where A is square. Uses UMFPACK if A is sparse. |
|
3870 +% The Control argument is optional. |
|
3871 +% |
|
3872 +% See also umfpack, umfpack_make, umfpack_details, umfpack_report, |
|
3873 +% and umfpack_simple. |
|
3874 + |
|
3875 +% UMFPACK Version 4.3 (Jan. 16, 2004), Copyright (c) 2004 by Timothy A. |
|
3876 +% Davis. All Rights Reserved. Type umfpack_details for License. |
|
3877 + |
|
3878 +%------------------------------------------------------------------------------- |
|
3879 +% check inputs and get default control parameters |
|
3880 +%------------------------------------------------------------------------------- |
|
3881 + |
|
3882 +if (op == '\\') |
|
3883 + A = arg1 ; |
|
3884 + b = arg2 ; |
|
3885 +elseif (op == '/') |
|
3886 + A = arg2 ; |
|
3887 + b = arg1 ; |
|
3888 +else |
|
3889 + help umfack_solve |
|
3890 + error ('umfpack_solve: unrecognized operator') ; |
|
3891 +end |
|
3892 + |
|
3893 +[m n] = size (A) ; |
|
3894 +if (m ~= n) |
|
3895 + help umfpack_solve |
|
3896 + error ('umfpack_solve: A must be square') ; |
|
3897 +end |
|
3898 + |
|
3899 +[m1 n1] = size (b) ; |
|
3900 +if ((op == '\\' & n ~= m1) | (op == '/' & n1 ~= m)) |
|
3901 + help umfpack_solve |
|
3902 + error ('umfpack_solve: b has the wrong dimensions') ; |
|
3903 +end |
|
3904 + |
|
3905 +if (nargin < 4) |
|
3906 + Control = umfpack ; |
|
3907 +end |
|
3908 + |
|
3909 +%------------------------------------------------------------------------------- |
|
3910 +% solve the system |
|
3911 +%------------------------------------------------------------------------------- |
|
3912 + |
|
3913 +if (op == '\\') |
|
3914 + |
|
3915 + if (~issparse (A)) |
|
3916 + |
|
3917 + % A is not sparse, so just use MATLAB |
|
3918 + x = A\b ; |
|
3919 + |
|
3920 + elseif (n1 == 1 & ~issparse (b)) |
|
3921 + |
|
3922 + % the UMFPACK '\' requires b to be a dense column vector |
|
3923 + x = umfpack (A, '\\', b, Control) ; |
|
3924 + |
|
3925 + else |
|
3926 + |
|
3927 + % factorize with UMFPACK and do the forward/back solves in MATLAB |
|
3928 + [L, U, P, Q, R] = umfpack (A, Control) ; |
|
3929 + keyboard |
|
3930 + x = Q * (U \ (L \ (P * (R \ b)))) ; |
|
3931 + |
|
3932 + end |
|
3933 + |
|
3934 +else |
|
3935 + |
|
3936 + if (~issparse (A)) |
|
3937 + |
|
3938 + % A is not sparse, so just use MATLAB |
|
3939 + x = b/A ; |
|
3940 + |
|
3941 + elseif (m1 == 1 & ~issparse (b)) |
|
3942 + |
|
3943 + % the UMFPACK '/' requires b to be a dense column vector |
|
3944 + x = umfpack (b, '/', A, Control) ; |
|
3945 + |
|
3946 + else |
|
3947 + |
|
3948 + % factorize with UMFPACK and do the forward/back solves in MATLAB |
|
3949 + % this mimics the behavior of x = b/A, except for the row scaling |
|
3950 + [L, U, P, Q, R] = umfpack (A.', Control) ; |
|
3951 + x = (Q * (U \ (L \ (P * (R \ (b.')))))).' ; |
|
3952 + |
|
3953 + % an alternative method: |
|
3954 + % [L, U, P, Q, r] = umfpack (A, Control) ; |
|
3955 + % x = (R \ (P' * (L.' \ (U.' \ (Q' * b.'))))).' ; |
|
3956 + |
|
3957 + end |
|
3958 + |
|
3959 +end |
|
3960 diff -uNr UMFPACKv4.4.orig/AMD/OCTAVE/amd.cc UMFPACK/AMD/OCTAVE/amd.cc |
|
3961 --- UMFPACKv4.4.orig/AMD/OCTAVE/amd.cc 1970-01-01 01:00:00.000000000 +0100 |
|
3962 +++ UMFPACK/AMD/OCTAVE/amd.cc 2005-02-01 21:51:23.944874478 +0100 |
|
3963 @@ -0,0 +1,299 @@ |
|
3964 +/* |
|
3965 + |
|
3966 +Copyright (C) 2004 David Bateman |
|
3967 + |
|
3968 +This program is free software; you can redistribute it and/or modify it |
|
3969 +under the terms of the GNU General Public License as published by the |
|
3970 +Free Software Foundation; either version 2, or (at your option) any |
|
3971 +later version. |
|
3972 + |
|
3973 +This program is distributed in the hope that it will be useful, but WITHOUT |
|
3974 +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
|
3975 +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
|
3976 +for more details. |
|
3977 + |
|
3978 +You should have received a copy of the GNU General Public License |
|
3979 +along with this program; see the file COPYING. If not, write to the Free |
|
3980 +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. |
|
3981 + |
|
3982 +In addition to the terms of the GPL, you are permitted to link |
|
3983 +this program with any Open Source program, as defined by the |
|
3984 +Open Source Initiative (www.opensource.org) |
|
3985 + |
|
3986 +*/ |
|
3987 + |
|
3988 +/* |
|
3989 + |
|
3990 +This is the Octave interface to the UMFPACK AMD code, which bore the following |
|
3991 +copyright |
|
3992 + |
|
3993 + AMD Version 1.1 (Jan. 21, 2004), Copyright (c) 2004 by Timothy A. Davis, |
|
3994 + Patrick R. Amestoy, and Iain S. Duff. See ../README for License. |
|
3995 + email: davis@cise.ufl.edu CISE Department, Univ. of Florida. |
|
3996 + web: http://www.cise.ufl.edu/research/sparse/amd |
|
3997 + -------------------------------------------------------------------------- |
|
3998 + |
|
3999 + Acknowledgements: This work was supported by the National Science |
|
4000 + Foundation, under grants ASC-9111263, DMS-9223088, and CCR-0203270. |
|
4001 + |
|
4002 +*/ |
|
4003 + |
|
4004 +#include <cstdlib> |
|
4005 +#include <string> |
|
4006 + |
|
4007 +#include <octave/config.h> |
|
4008 +#include <octave/ov.h> |
|
4009 +#include <octave/defun-dld.h> |
|
4010 +#include <octave/pager.h> |
|
4011 +#include <octave/ov-re-mat.h> |
|
4012 + |
|
4013 +#include "ov-re-sparse.h" |
|
4014 +#include "ov-cx-sparse.h" |
|
4015 + |
|
4016 +// External AMD functions in C |
|
4017 +extern "C" { |
|
4018 +#include "amd.h" |
|
4019 +} |
|
4020 + |
|
4021 +DEFUN_DLD (amd, args, nargout, |
|
4022 + "-*- texinfo -*-\n\ |
|
4023 +@deftypefn {Loadable Function} {@var{p} =} amd (@var{s})\n\ |
|
4024 +@deftypefnx {Loadable Function} {@var{Control} =} amd ()\n\ |
|
4025 +@deftypefnx {Loadable Function} {[@var{p}, @var{info}] =} amd (@var{s})\n\ |
|
4026 +\n\ |
|
4027 +AMD Approximate minimum degree permutation. Returns the approximate\n\ |
|
4028 +minimum degree permutation vector for the sparse matrix\n\ |
|
4029 +@code{@var{c} = @var{S} + @var{S}'}. The Cholesky factorization of\n\ |
|
4030 +@code{@var{c} (@var{p}, @var{p})}, or @code{@var{s} (@var{p}, @var{p})},\n\ |
|
4031 +tends to be sparser than that of @var{c} or @var{s}.\n\ |
|
4032 +@var{s} must be square. If @var{s} is full, @code{amd (@var{S})} is\n\ |
|
4033 +equivalent to @code{amd (sparse (@var{s}))}.\n\ |
|
4034 +\n\ |
|
4035 +@table @asis\n\ |
|
4036 +@item @var{Control} (1)\n\ |
|
4037 +If S is n-by-n, then rows/columns with more than\n\ |
|
4038 +@code{@dfn{max} (16, (@var{Control} (1)) * @dfn{sqrt} (@var{n}))} entries\n\ |
|
4039 +in @code{@var{s} + @var{S}'} are considered @emph{dense}, and ignored during\n\ |
|
4040 +ordering. They are placed last in the output permutation. The default is\n\ |
|
4041 +10.0 if @var{Control} is not present.\n\ |
|
4042 +@item @var{Control} (2)\n\ |
|
4043 +If nonzero, then aggressive absorption is performed. This is the default if\n\ |
|
4044 +@var{Control} is not present.\n\ |
|
4045 +@item @var{Control} (3)\n\ |
|
4046 +If nonzero, print statistics about the ordering.\n\ |
|
4047 +@end table\n\ |
|
4048 +\n\ |
|
4049 +@table @asis\n\ |
|
4050 +@item @var{Info} (1)\n\ |
|
4051 +status (0: ok, -1: out of memory, -2: matrix invalid)\n\ |
|
4052 +@item @var{Info} (2)\n\ |
|
4053 +@code{@var{n} = size (@var{a}, 1)}\n\ |
|
4054 +@item @var{Info} (3)\n\ |
|
4055 +@code{nnz (A)}\n\ |
|
4056 +@item @var{Info} (4)\n\ |
|
4057 +The symmetry of the matrix @var{s} (0.0 means purely unsymmetric, 1.0 means\n\ |
|
4058 +purely symmetric). Computed as: @code{@var{b} = tril (@var{s}, -1) +\n\ |
|
4059 +triu (@var{s}, 1); @var{symmetry} = nnz (@var{b} & @var{b}')\n\ |
|
4060 +/ nnz (@var{b});}\n\ |
|
4061 +@item @var{Info} (5)\n\ |
|
4062 +@code{nnz (diag (@var{s}))}\n\ |
|
4063 +@item @var{Info} (6)\n\ |
|
4064 +@dfn{nnz} in @code{@var{s} + @var{s}'}, excluding the diagonal\n\ |
|
4065 +(= @code{nnz (@var{b} + @var{b}')})\n\ |
|
4066 +@item @var{Info} (7)\n\ |
|
4067 +Number of @emph{dense} rows/columns in @code{@var{s} + @var{s}'}\n\ |
|
4068 +@item @var{Info} (8)\n\ |
|
4069 +The amount of memory used by AMD, in bytes\n\ |
|
4070 +@item @var{Info} (9)\n\ |
|
4071 +The number of memory compactions performed by AMD\n\ |
|
4072 +@end table\n\ |
|
4073 +\n\ |
|
4074 +The following statistics are slight upper bounds because of the\n\ |
|
4075 +approximate degree in AMD. The bounds are looser if @emph{dense}\n\ |
|
4076 +rows/columns are ignored during ordering @code{(@var{Info} (7) > 0)}.\n\ |
|
4077 +The statistics are for a subsequent factorization of the matrix\n\ |
|
4078 +@code{@var{c} (@var{p},@var{p})}. The LU factorization statistics assume\n \ |
|
4079 +no pivoting.\n\ |
|
4080 +\n\ |
|
4081 +@table @asis\n\ |
|
4082 +@item @var{Info} (10)\n\ |
|
4083 +The number of nonzeros in L, excluding the diagonal\n\ |
|
4084 +@item @var{Info} (11)\n\ |
|
4085 +The number of divide operations for LL', LDL', or LU\n\ |
|
4086 +@item @var{Info (12)}\n\ |
|
4087 +The number of multiply-subtract pairs for LL' or LDL'\n\ |
|
4088 +@item @var{Info} (13)\n\ |
|
4089 +The number of multiply-subtract pairs for LU\n\ |
|
4090 +@item @var{Info} (14)\n\ |
|
4091 +The max number of nonzeros in any column of L (incl. diagonal)\n\ |
|
4092 +@item @var{Info} (15:20)\n\ |
|
4093 +unused, reserved for future use\n\ |
|
4094 +@end table\n\ |
|
4095 +\n\ |
|
4096 +An assembly tree post-ordering is performed, which is typically the same\n\ |
|
4097 +as an elimination tree post-ordering. It is not always identical because\n\ |
|
4098 +of the approximate degree update used, and because @emph{dense} rows/columns\n\ |
|
4099 +do not take part in the post-order. It well-suited for a subsequent\n\ |
|
4100 +@dfn{chol}, however. If you require a precise elimination tree\n\ |
|
4101 +post-ordering, then do:\n\ |
|
4102 +\n\ |
|
4103 +@group\n\ |
|
4104 + @var{p} = @dfn{amd} (@var{s});\n\ |
|
4105 + % skip this if S already symmetric\n\ |
|
4106 + @var{c} = spones (@var{s}) + spones (@var{s}');\n\ |
|
4107 + [@var{ignore}, @var{q}] = sparsfun ('symetree', @var{c} (@var{p}, @var{p}));\n\ |
|
4108 + @var{p} = @var{p} (@var{q});\n\ |
|
4109 +@end group\n\ |
|
4110 +\n\ |
|
4111 +AMD Version 1.1 (Jan. 21, 2004), Copyright @copyright{} 2004 by\n\ |
|
4112 +Timothy A. Davis, Patrick R. Amestoy, and Iain S. Duff.\n\ |
|
4113 +\n\ |
|
4114 +email: davis@@cise.ufl.edu (CISE Department, Univ. of Florida).\n\ |
|
4115 +\n\ |
|
4116 +web: http://www.cise.ufl.edu/research/sparse/amd\n\ |
|
4117 +\n\ |
|
4118 +Acknowledgements: This work was supported by the National Science\n\ |
|
4119 +Foundation, under grants ASC-9111263, DMS-9223088, and CCR-0203270.\n\ |
|
4120 +@end deftypefn") |
|
4121 +{ |
|
4122 + int nargin = args.length (); |
|
4123 + octave_value_list retval; |
|
4124 + int spumoni = 0; |
|
4125 + |
|
4126 + if (nargin > 2 || nargout > 2) |
|
4127 + usage ("p = amd (A) or [p, Info] = amd (A, Control)"); |
|
4128 + else if (nargin == 0) |
|
4129 + { |
|
4130 + // Get the default control parameter, and return |
|
4131 + NDArray control (dim_vector (AMD_CONTROL, 1)); |
|
4132 + double *control_ptr = control.fortran_vec (); |
|
4133 + amd_defaults (control_ptr); |
|
4134 + if (nargout == 0) |
|
4135 + amd_control (control_ptr); |
|
4136 + else |
|
4137 + retval(0) = control; |
|
4138 + } |
|
4139 + else |
|
4140 + { |
|
4141 + NDArray control (dim_vector (AMD_CONTROL, 1)); |
|
4142 + double *control_ptr = control.fortran_vec (); |
|
4143 + amd_defaults (control_ptr); |
|
4144 + |
|
4145 + if (nargin > 1) |
|
4146 + { |
|
4147 + NDArray control_in = args(1).array_value(); |
|
4148 + |
|
4149 + if (error_state) |
|
4150 + { |
|
4151 + error ("amd: could not read control vector"); |
|
4152 + return retval; |
|
4153 + } |
|
4154 + |
|
4155 + dim_vector dv = control_in.dims (); |
|
4156 + if (dv.length() > 2 || (dv(0) != 1 && dv(1) != 1)) |
|
4157 + { |
|
4158 + error ("amd: control vector isn't a vector"); |
|
4159 + return retval; |
|
4160 + } |
|
4161 + |
|
4162 + int nc = dv.numel (); |
|
4163 + control (AMD_DENSE) = (nc > 0 ? control_in (AMD_DENSE) : |
|
4164 + AMD_DEFAULT_DENSE); |
|
4165 + control (AMD_AGGRESSIVE) = (nc > 1 ? control_in (AMD_AGGRESSIVE) : |
|
4166 + AMD_DEFAULT_AGGRESSIVE); |
|
4167 + spumoni = (nc > 2 ? (control_in (2) != 0) : 0); |
|
4168 + } |
|
4169 + |
|
4170 + if (spumoni > 0) |
|
4171 + amd_control (control_ptr); |
|
4172 + |
|
4173 + int *Ap, *Ai; |
|
4174 + int n, m, nz; |
|
4175 + |
|
4176 + // These are here only so that the C++ destructors don't prematurally |
|
4177 + // remove the underlying data we are interested in |
|
4178 + SparseMatrix sm; |
|
4179 + SparseComplexMatrix scm; |
|
4180 + Matrix mm; |
|
4181 + ComplexMatrix cm; |
|
4182 + |
|
4183 + if (args(0).class_name () != "sparse" && spumoni > 0) |
|
4184 + octave_stdout << " input matrix A is full (sparse copy" |
|
4185 + << " of A will be created)" << std::endl; |
|
4186 + |
|
4187 + if (args(0).is_complex_type ()) |
|
4188 + { |
|
4189 + scm = args(0).sparse_complex_matrix_value (); |
|
4190 + Ai = scm.ridx (); |
|
4191 + Ap = scm.cidx (); |
|
4192 + m = scm.rows (); |
|
4193 + n = scm.cols (); |
|
4194 + nz = scm.nnz (); |
|
4195 + } |
|
4196 + else |
|
4197 + { |
|
4198 + sm = args(0).sparse_matrix_value (); |
|
4199 + Ai = sm.ridx (); |
|
4200 + Ap = sm.cidx (); |
|
4201 + m = sm.rows (); |
|
4202 + n = sm.cols (); |
|
4203 + nz = sm.nnz (); |
|
4204 + } |
|
4205 + |
|
4206 + if (spumoni > 0) |
|
4207 + octave_stdout << " input matrix A is " << m << "-by-" << n |
|
4208 + << std::endl; |
|
4209 + |
|
4210 + if (m != n) |
|
4211 + { |
|
4212 + error ("amd: A must be square"); |
|
4213 + return retval; |
|
4214 + } |
|
4215 + |
|
4216 + if (spumoni > 0) |
|
4217 + octave_stdout << " input matrix A has " << nz << |
|
4218 + " nonzero entries" << std::endl; |
|
4219 + |
|
4220 + // allocate workspace for output permutation |
|
4221 + Array<int> P(n+1); |
|
4222 + int *P_ptr = P.fortran_vec (); |
|
4223 + NDArray info (dim_vector (AMD_INFO, 1)); |
|
4224 + double *info_ptr = info.fortran_vec (); |
|
4225 + int result; |
|
4226 + |
|
4227 + // order the matrix |
|
4228 + result = amd_order (n, Ap, Ai, P_ptr, control_ptr, info_ptr); |
|
4229 + |
|
4230 + // print results (including return value) |
|
4231 + if (spumoni > 0) |
|
4232 + amd_info (info_ptr); |
|
4233 + |
|
4234 + // check error conditions |
|
4235 + if (result == AMD_OUT_OF_MEMORY) |
|
4236 + error ("amd: out of memory"); |
|
4237 + else if (result == AMD_INVALID) |
|
4238 + error ("amd: input matrix A is corrupted"); |
|
4239 + else |
|
4240 + { |
|
4241 + // copy the outputs to Octave |
|
4242 + |
|
4243 + // output permutation, P |
|
4244 + NDArray perm (dim_vector (1, n)); |
|
4245 + for (int i = 0; i < n; i++) |
|
4246 + perm (i) = double (P(i) + 1); // 1-based indexing for Octave |
|
4247 + |
|
4248 + retval (0) = perm; |
|
4249 + |
|
4250 + // Info |
|
4251 + if (nargout > 1) |
|
4252 + retval (1) = info; |
|
4253 + } |
|
4254 + } |
|
4255 + return retval; |
|
4256 +} |
|
4257 + |
|
4258 +/* |
|
4259 +;;; Local Variables: *** |
|
4260 +;;; mode: C++ *** |
|
4261 +;;; End: *** |
|
4262 +*/ |
|
4263 diff -uNr UMFPACKv4.4.orig/AMD/OCTAVE/amd_demo.m UMFPACK/AMD/OCTAVE/amd_demo.m |
|
4264 --- UMFPACKv4.4.orig/AMD/OCTAVE/amd_demo.m 1970-01-01 01:00:00.000000000 +0100 |
|
4265 +++ UMFPACK/AMD/OCTAVE/amd_demo.m 2004-12-30 01:58:45.000000000 +0100 |
|
4266 @@ -0,0 +1,61 @@ |
|
4267 +function amd_demo |
|
4268 +% AMD DEMO |
|
4269 +% |
|
4270 +% A demo of AMD for OCTAVE. |
|
4271 +% |
|
4272 +% -------------------------------------------------------------------------- |
|
4273 +% AMD Version 1.1 (Jan. 21, 2004), Copyright (c) 2004 by Timothy A. Davis, |
|
4274 +% Patrick R. Amestoy, and Iain S. Duff. See ../README for License. |
|
4275 +% email: davis@cise.ufl.edu CISE Department, Univ. of Florida. |
|
4276 +% web: http://www.cise.ufl.edu/research/sparse/amd |
|
4277 +% -------------------------------------------------------------------------- |
|
4278 +% |
|
4279 +% See also: amd |
|
4280 + |
|
4281 +% Get the Harwell/Boeing can_24 matrix. This is an example matrix from the |
|
4282 +% MATLAB-accessible UF sparse matrix collection, and can be loaded into |
|
4283 +% MATLAB with the statment "Problem = UFget ('HB/can_24')", after obtaining |
|
4284 +% the UFget function and its supporting routines at |
|
4285 +% http://www.cise.ufl.edu/sparse/mat . |
|
4286 + |
|
4287 +load can_24.mat |
|
4288 +A = Problem.A ; |
|
4289 +n = size (A,1) ; |
|
4290 + |
|
4291 +figure (1) |
|
4292 +clf |
|
4293 +hold off |
|
4294 +subplot (2,1,1) ; |
|
4295 +% remove the "_" from the name before printing it in the plot title |
|
4296 +title (sprintf ('%s', strrep (Problem.name, '_', '-'))) ; |
|
4297 +fprintf ('Matrix name: %s\n', Problem.name) ; |
|
4298 +fprintf ('Matrix title: %s\n', Problem.title) ; |
|
4299 +spy (A) |
|
4300 + |
|
4301 +% print the details during AMD ordering and SYMBFACT |
|
4302 +control = amd (); |
|
4303 +control (3) = 1; |
|
4304 + |
|
4305 +% order the matrix. Note that the Info argument is optional. |
|
4306 +fprintf ('\nIf the next step fails, then you have\n') ; |
|
4307 +fprintf ('not yet compiled the AMD mexFunction.\n') ; |
|
4308 +[p, Info] = amd (A, control) ; |
|
4309 + |
|
4310 +% order again, but this time print some statistics |
|
4311 +[p, Info] = amd (A, [10 1 1]) ; |
|
4312 + |
|
4313 +fprintf ('Permutation vector:\n') ; |
|
4314 +fprintf (' %2d', p) ; |
|
4315 +fprintf ('\n\n') ; |
|
4316 + |
|
4317 +subplot (2,1,2) ; |
|
4318 +title ('Permuted matrix') ; |
|
4319 +spy (A (p,p)) |
|
4320 + |
|
4321 +% approximations from amd: |
|
4322 +lnz2 = n + Info (10) ; |
|
4323 +fl2 = n + Info (11) + 2 * Info (12) ; |
|
4324 +fprintf ('\nResults from AMD''s approximate analysis:\n') ; |
|
4325 +fprintf ('number of nonzeros in L (including diagonal): %d\n', lnz2) ; |
|
4326 +fprintf ('floating point operation count for chol (A (p,p)): %d\n\n', fl2) ; |
|
4327 + |
|
4328 diff -uNr UMFPACKv4.4.orig/AMD/OCTAVE/amd_demo.m.out UMFPACK/AMD/OCTAVE/amd_demo.m.out |
|
4329 --- UMFPACKv4.4.orig/AMD/OCTAVE/amd_demo.m.out 1970-01-01 01:00:00.000000000 +0100 |
|
4330 +++ UMFPACK/AMD/OCTAVE/amd_demo.m.out 2004-12-30 01:58:45.000000000 +0100 |
|
4331 @@ -0,1 +1,21 @@ |
|
4332 +octave:3> amd_demo |
|
4333 +ans = 1 |
|
4334 +Matrix name: HB/can_24 |
|
4335 +Matrix title: 1SYMMETRIC PATTERN FROM CANNES,LUCIEN MARRO,JUNE 1981. |
|
4336 + |
|
4337 +If the next step fails, then you have |
|
4338 +not yet compiled the AMD mexFunction. |
|
4339 + input matrix A is 24-by-24 |
|
4340 + input matrix A has 160 nonzero entries |
|
4341 + input matrix A is 24-by-24 |
|
4342 + input matrix A has 160 nonzero entries |
|
4343 +Permutation vector: |
|
4344 + 23 21 11 24 13 6 17 9 15 5 16 8 2 10 14 18 1 3 4 7 12 19 22 20 |
|
4345 + |
|
4346 + |
|
4347 +Results from AMD's approximate analysis: |
|
4348 +number of nonzeros in L (including diagonal): 121 |
|
4349 +floating point operation count for chol (A (p,p)): 671 |
|
4350 + |
|
4351 +octave:4> quit |
|
4352 --- UMFPACKv4.4.orig/AMD/OCTAVE/GNUmakefile 1970-01-01 01:00:00.000000000 +0100 |
|
4353 +++ UMFPACK/AMD/OCTAVE/GNUmakefile 2004-12-30 01:58:45.000000000 +0100 |
|
4354 @@ -0,0 +1,33 @@ |
|
4355 +#------------------------------------------------------------------------------ |
|
4356 +# GNUmakefile for the AMD MATLAB mexFunction |
|
4357 +#------------------------------------------------------------------------------ |
|
4358 + |
|
4359 +all: amd |
|
4360 + |
|
4361 +include ../Make/Make.include |
|
4362 + |
|
4363 +MKOCT = mkoctfile -I../Include |
|
4364 + |
|
4365 +OCT_SPARSE_INC = -I../../../ |
|
4366 + |
|
4367 +AMD = amd_aat amd_1 amd_2 amd_dump amd_postorder amd_post_tree amd_defaults \ |
|
4368 + amd_order amd_control amd_info amd_valid |
|
4369 + |
|
4370 +INC = ../Include/amd.h ../Source/amd_internal.h |
|
4371 + |
|
4372 +OCTAMD = $(addsuffix .o, $(subst amd_,amd_o_,$(AMD))) |
|
4373 + |
|
4374 +amd_o_%.o: ../Source/amd_%.c $(INC) |
|
4375 + $(MKOCT) -DDINT -c $< -o $@ |
|
4376 + - $(MV) ../Source/amd_$*.o |
|
4377 + |
|
4378 +# Note temporary addition of octave sparse path |
|
4379 +amd: amd.cc $(OCTAMD) $(INC) |
|
4380 + $(MKOCT) amd.cc $(OCTAMD) $(OCT_SPARSE_INC) -o amd.oct |
|
4381 + |
|
4382 +#------------------------------------------------------------------------------ |
|
4383 +# Remove all but the files in the original distribution |
|
4384 +#------------------------------------------------------------------------------ |
|
4385 + |
|
4386 +purge: clean |
|
4387 + - $(RM) amd.oct |
|
4388 diff -uNr UMFPACKv4.4.orig/AMD/OCTAVE/Makefile UMFPACK/AMD/OCTAVE/Makefile |
|
4389 --- UMFPACKv4.4.orig/AMD/OCTAVE/Makefile 1970-01-01 01:00:00.000000000 +0100 |
|
4390 +++ UMFPACK/AMD/OCTAVE/Makefile 2004-12-30 01:58:45.000000000 +0100 |
|
4391 @@ -0,0 +1,41 @@ |
|
4392 +#------------------------------------------------------------------------------ |
|
4393 +# compile the AMD mexFunction for MATLAB (original make only) |
|
4394 +#------------------------------------------------------------------------------ |
|
4395 + |
|
4396 +# This is a very ugly Makefile, and is only provided for those who do not |
|
4397 +# have GNU make. Note that it is not used if you have GNU make. It ignores |
|
4398 +# dependency checking and just compiles everything. It was created |
|
4399 +# automatically, via make -n using the GNUmakefile. That way, I don't have |
|
4400 +# maintain two Makefiles. |
|
4401 + |
|
4402 +all: amd |
|
4403 + |
|
4404 +include ../Make/Make.include |
|
4405 + |
|
4406 +MKOCT = mkoctfile -I../Include |
|
4407 + |
|
4408 +OCT_SPARSE_INC = ../../../ |
|
4409 + |
|
4410 +amd: |
|
4411 + $(MKOCT) -DDINT -o amd_o_aat.o -c ../Source/amd_aat.c |
|
4412 + $(MKOCT) -DDINT -o amd_o_1.o -c ../Source/amd_1.c |
|
4413 + $(MKOCT) -DDINT -o amd_o_2.o -c ../Source/amd_2.c |
|
4414 + $(MKOCT) -DDINT -o amd_o_dump.o -c ../Source/amd_dump.c |
|
4415 + $(MKOCT) -DDINT -o amd_o_postorder.o -c ../Source/amd_postorder.c |
|
4416 + $(MKOCT) -DDINT -o amd_o_post_tree.o -c ../Source/amd_post_tree.c |
|
4417 + $(MKOCT) -DDINT -o amd_o_defaults.o -c ../Source/amd_defaults.c |
|
4418 + $(MKOCT) -DDINT -o amd_o_order.o -c ../Source/amd_order.c |
|
4419 + $(MKOCT) -DDINT -o amd_o_control.o -c ../Source/amd_control.c |
|
4420 + $(MKOCT) -DDINT -o amd_o_info.o -c ../Source/amd_info.c |
|
4421 + $(MKOCT) -DDINT -o amd_o_valid.o -c ../Source/amd_valid.c |
|
4422 + $(MKOCT) -output amd.oct amd_mex.c amd_o_aat.o \ |
|
4423 + amd_o_1.o amd_o_2.o amd_o_dump.o amd_o_postorder.o \ |
|
4424 + amd_o_post_tree.o amd_o_defaults.o amd_o_order.o amd_o_control.o \ |
|
4425 + amd_o_info.o amd_o_valid.o $(OCT_SPARSE_INC) -o amd.oct |
|
4426 + |
|
4427 +#------------------------------------------------------------------------------ |
|
4428 +# Remove all but the files in the original distribution |
|
4429 +#------------------------------------------------------------------------------ |
|
4430 + |
|
4431 +purge: clean |
|
4432 + - $(RM) amd.oct |
|
4433 --- UMFPACKv4.4.orig/UMFPACK/Source/umf_solve.c 2005-01-17 17:21:29.000000000 +0100 |
|
4434 +++ UMFPACK/UMFPACK/Source/umf_solve.c 2005-02-02 00:01:22.904346651 +0100 |
|
4435 @@ -79,7 +79,8 @@ |
|
4436 double *Z2, *Y, *B2, *Rs ; |
|
4437 Int *Rperm, *Cperm, i, n, p, step, j, nz, status, p2, do_scale ; |
|
4438 #ifdef COMPLEX |
|
4439 - Int split ; |
|
4440 + Int AXsplit ; |
|
4441 + Int Bsplit ; |
|
4442 #endif |
|
4443 #ifndef NRECIPROCAL |
|
4444 Int do_recip = Numeric->do_recip ; |
|
4445 @@ -141,7 +142,7 @@ |
|
4446 return (UMFPACK_ERROR_argument_missing) ; |
|
4447 } |
|
4448 /* A, B, and X in split format if Az, Bz, and Xz present */ |
|
4449 - split = SPLIT (Az) && SPLIT (Bz) && SPLIT (Xz) ; |
|
4450 + AXsplit = SPLIT (Az) || SPLIT(Xz); |
|
4451 Z = (Entry *) (SolveWork + 4*n) ; /* Entry Z [0..n-1] */ |
|
4452 S = (Entry *) (SolveWork + 6*n) ; /* Entry S [0..n-1] */ |
|
4453 Y = (double *) (SolveWork + 8*n) ; /* double Y [0..n-1] */ |
|
4454 @@ -150,10 +151,12 @@ |
|
4455 } |
|
4456 else |
|
4457 { |
|
4458 - /* A is ignored, only look at X and B for split/packed cases */ |
|
4459 - split = SPLIT (Bz) && SPLIT (Xz) ; |
|
4460 + /* A is ignored, only look at X for split/packed cases */ |
|
4461 + AXsplit = SPLIT(Xz); |
|
4462 } |
|
4463 - if (split) |
|
4464 + Bsplit = SPLIT (Bz); |
|
4465 + |
|
4466 + if (AXsplit) |
|
4467 { |
|
4468 X = (Entry *) (SolveWork + 2*n) ; /* Entry X [0..n-1] */ |
|
4469 } |
|
4470 @@ -209,7 +212,7 @@ |
|
4471 for (p = 0 ; p < p2 ; p++) |
|
4472 { |
|
4473 /* Y [Ai [p]] += ABS (Ax [p]) ; */ |
|
4474 - ASSIGN (aij, Ax, Az, p, split) ; |
|
4475 + ASSIGN (aij, Ax, Az, p, AXsplit) ; |
|
4476 ABS (d, aij) ; |
|
4477 Y [Ai [p]] += d ; |
|
4478 } |
|
4479 @@ -219,7 +222,7 @@ |
|
4480 for (i = 0 ; i < n ; i++) |
|
4481 { |
|
4482 /* B2 [i] = ABS (B [i]) ; */ |
|
4483 - ASSIGN (bi, Bx, Bz, i, split) ; |
|
4484 + ASSIGN (bi, Bx, Bz, i, Bsplit) ; |
|
4485 ABS (B2 [i], bi) ; |
|
4486 } |
|
4487 |
|
4488 @@ -276,7 +279,7 @@ |
|
4489 /* multiply by the scale factors */ |
|
4490 for (i = 0 ; i < n ; i++) |
|
4491 { |
|
4492 - ASSIGN (X [i], Bx, Bz, i, split) ; |
|
4493 + ASSIGN (X [i], Bx, Bz, i, Bsplit) ; |
|
4494 SCALE (X [i], Rs [i]) ; |
|
4495 } |
|
4496 } |
|
4497 @@ -286,7 +289,7 @@ |
|
4498 /* divide by the scale factors */ |
|
4499 for (i = 0 ; i < n ; i++) |
|
4500 { |
|
4501 - ASSIGN (X [i], Bx, Bz, i, split) ; |
|
4502 + ASSIGN (X [i], Bx, Bz, i, Bsplit) ; |
|
4503 SCALE_DIV (X [i], Rs [i]) ; |
|
4504 } |
|
4505 } |
|
4506 @@ -302,7 +305,7 @@ |
|
4507 for (i = 0 ; i < n ; i++) |
|
4508 { |
|
4509 /* W [i] = B [Rperm [i]] ; */ |
|
4510 - ASSIGN (W [i], Bx, Bz, Rperm [i], split) ; |
|
4511 + ASSIGN (W [i], Bx, Bz, Rperm [i], Bsplit) ; |
|
4512 } |
|
4513 } |
|
4514 } |
|
4515 @@ -311,7 +314,7 @@ |
|
4516 for (i = 0 ; i < n ; i++) |
|
4517 { |
|
4518 /* Z [i] = B [i] ; */ |
|
4519 - ASSIGN (Z [i], Bx, Bz, i, split) ; |
|
4520 + ASSIGN (Z [i], Bx, Bz, i, Bsplit) ; |
|
4521 } |
|
4522 flops += MULTSUB_FLOPS * nz ; |
|
4523 for (i = 0 ; i < n ; i++) |
|
4524 @@ -321,7 +324,7 @@ |
|
4525 for (p = Ap [i] ; p < p2 ; p++) |
|
4526 { |
|
4527 /* Z [Ai [p]] -= Ax [p] * xi ; */ |
|
4528 - ASSIGN (aij, Ax, Az, p, split) ; |
|
4529 + ASSIGN (aij, Ax, Az, p, AXsplit) ; |
|
4530 MULT_SUB (Z [Ai [p]], aij, xi) ; |
|
4531 } |
|
4532 } |
|
4533 @@ -390,7 +393,7 @@ |
|
4534 for (i = 0 ; i < n ; i++) |
|
4535 { |
|
4536 /* W [i] = B [i] ; */ |
|
4537 - ASSIGN (W [i], Bx, Bz, i, split) ; |
|
4538 + ASSIGN (W [i], Bx, Bz, i, Bsplit) ; |
|
4539 Z2 [i] = 0. ; |
|
4540 } |
|
4541 flops += (MULT_FLOPS + DECREMENT_FLOPS + ABS_FLOPS + 1) * nz ; |
|
4542 @@ -403,7 +406,7 @@ |
|
4543 i = Ai [p] ; |
|
4544 |
|
4545 /* axx = Ax [p] * xj ; */ |
|
4546 - ASSIGN (aij, Ax, Az, p, split) ; |
|
4547 + ASSIGN (aij, Ax, Az, p, AXsplit) ; |
|
4548 MULT (axx, aij, xj) ; |
|
4549 |
|
4550 /* W [i] -= axx ; */ |
|
4551 @@ -493,7 +496,7 @@ |
|
4552 /* yi += ABS (Ax [p]) * Rs [Ai [p]] ; */ |
|
4553 /* note that abs (aij) is the same as |
|
4554 * abs (conj (aij)) */ |
|
4555 - ASSIGN (aij, Ax, Az, p, split) ; |
|
4556 + ASSIGN (aij, Ax, Az, p, AXsplit) ; |
|
4557 ABS (d, aij) ; |
|
4558 yi += (d * Rs [Ai [p]]) ; |
|
4559 } |
|
4560 @@ -513,7 +516,7 @@ |
|
4561 /* yi += ABS (Ax [p]) / Rs [Ai [p]] ; */ |
|
4562 /* note that abs (aij) is the same as |
|
4563 * abs (conj (aij)) */ |
|
4564 - ASSIGN (aij, Ax, Az, p, split) ; |
|
4565 + ASSIGN (aij, Ax, Az, p, AXsplit) ; |
|
4566 ABS (d, aij) ; |
|
4567 yi += (d / Rs [Ai [p]]) ; |
|
4568 } |
|
4569 @@ -534,7 +537,7 @@ |
|
4570 /* yi += ABS (Ax [p]) ; */ |
|
4571 /* note that abs (aij) is the same as |
|
4572 * abs (conj (aij)) */ |
|
4573 - ASSIGN (aij, Ax, Az, p, split) ; |
|
4574 + ASSIGN (aij, Ax, Az, p, AXsplit) ; |
|
4575 ABS (d, aij) ; |
|
4576 yi += d ; |
|
4577 } |
|
4578 @@ -546,7 +549,7 @@ |
|
4579 for (i = 0 ; i < n ; i++) |
|
4580 { |
|
4581 /* B2 [i] = ABS (B [i]) ; */ |
|
4582 - ASSIGN (bi, Bx, Bz, i, split) ; |
|
4583 + ASSIGN (bi, Bx, Bz, i, Bsplit) ; |
|
4584 ABS (B2 [i], bi) ; |
|
4585 } |
|
4586 |
|
4587 @@ -568,7 +571,7 @@ |
|
4588 for (i = 0 ; i < n ; i++) |
|
4589 { |
|
4590 /* W [i] = B [Cperm [i]] ; */ |
|
4591 - ASSIGN (W [i], Bx, Bz, Cperm [i], split) ; |
|
4592 + ASSIGN (W [i], Bx, Bz, Cperm [i], Bsplit) ; |
|
4593 } |
|
4594 } |
|
4595 else |
|
4596 @@ -577,7 +580,7 @@ |
|
4597 for (i = 0 ; i < n ; i++) |
|
4598 { |
|
4599 /* Z [i] = B [i] ; */ |
|
4600 - ASSIGN (Z [i], Bx, Bz, i, split) ; |
|
4601 + ASSIGN (Z [i], Bx, Bz, i, Bsplit) ; |
|
4602 } |
|
4603 flops += MULTSUB_FLOPS * nz ; |
|
4604 for (i = 0 ; i < n ; i++) |
|
4605 @@ -587,7 +590,7 @@ |
|
4606 for (p = Ap [i] ; p < p2 ; p++) |
|
4607 { |
|
4608 /* zi -= conjugate (Ax [p]) * X [Ai [p]] ; */ |
|
4609 - ASSIGN (aij, Ax, Az, p, split) ; |
|
4610 + ASSIGN (aij, Ax, Az, p, Bsplit) ; |
|
4611 MULT_SUB_CONJ (zi, X [Ai [p]], aij) ; |
|
4612 } |
|
4613 Z [i] = zi ; |
|
4614 @@ -696,13 +699,13 @@ |
|
4615 for (i = 0 ; i < n ; i++) |
|
4616 { |
|
4617 /* wi = B [i] ; */ |
|
4618 - ASSIGN (wi, Bx, Bz, i, split) ; |
|
4619 + ASSIGN (wi, Bx, Bz, i, Bsplit) ; |
|
4620 z2i = 0. ; |
|
4621 p2 = Ap [i+1] ; |
|
4622 for (p = Ap [i] ; p < p2 ; p++) |
|
4623 { |
|
4624 /* axx = conjugate (Ax [p]) * X [Ai [p]] ; */ |
|
4625 - ASSIGN (aij, Ax, Az, p, split) ; |
|
4626 + ASSIGN (aij, Ax, Az, p, AXsplit) ; |
|
4627 MULT_CONJ (axx, X [Ai [p]], aij) ; |
|
4628 |
|
4629 /* wi -= axx ; */ |
|
4630 @@ -766,7 +769,7 @@ |
|
4631 /* yi += ABS (Ax [p]) * Rs [Ai [p]] ; */ |
|
4632 /* note that A.' is the array transpose, |
|
4633 * so no conjugate */ |
|
4634 - ASSIGN (aij, Ax, Az, p, split) ; |
|
4635 + ASSIGN (aij, Ax, Az, p, AXsplit) ; |
|
4636 ABS (d, aij) ; |
|
4637 yi += (d * Rs [Ai [p]]) ; |
|
4638 } |
|
4639 @@ -786,7 +789,7 @@ |
|
4640 /* yi += ABS (Ax [p]) / Rs [Ai [p]] ; */ |
|
4641 /* note that A.' is the array transpose, |
|
4642 * so no conjugate */ |
|
4643 - ASSIGN (aij, Ax, Az, p, split) ; |
|
4644 + ASSIGN (aij, Ax, Az, p, AXsplit) ; |
|
4645 ABS (d, aij) ; |
|
4646 yi += (d / Rs [Ai [p]]) ; |
|
4647 } |
|
4648 @@ -807,7 +810,7 @@ |
|
4649 /* yi += ABS (Ax [p]) */ |
|
4650 /* note that A.' is the array transpose, |
|
4651 * so no conjugate */ |
|
4652 - ASSIGN (aij, Ax, Az, p, split) ; |
|
4653 + ASSIGN (aij, Ax, Az, p, AXsplit) ; |
|
4654 ABS (d, aij) ; |
|
4655 yi += d ; |
|
4656 } |
|
4657 @@ -819,7 +822,7 @@ |
|
4658 for (i = 0 ; i < n ; i++) |
|
4659 { |
|
4660 /* B2 [i] = ABS (B [i]) ; */ |
|
4661 - ASSIGN (bi, Bx, Bz, i, split) ; |
|
4662 + ASSIGN (bi, Bx, Bz, i, Bsplit) ; |
|
4663 ABS (B2 [i], bi) ; |
|
4664 } |
|
4665 |
|
4666 @@ -841,7 +844,7 @@ |
|
4667 for (i = 0 ; i < n ; i++) |
|
4668 { |
|
4669 /* W [i] = B [Cperm [i]] ; */ |
|
4670 - ASSIGN (W [i], Bx, Bz, Cperm [i], split) ; |
|
4671 + ASSIGN (W [i], Bx, Bz, Cperm [i], Bsplit) ; |
|
4672 } |
|
4673 } |
|
4674 else |
|
4675 @@ -850,7 +853,7 @@ |
|
4676 for (i = 0 ; i < n ; i++) |
|
4677 { |
|
4678 /* Z [i] = B [i] ; */ |
|
4679 - ASSIGN (Z [i], Bx, Bz, i, split) ; |
|
4680 + ASSIGN (Z [i], Bx, Bz, i, Bsplit) ; |
|
4681 } |
|
4682 flops += MULTSUB_FLOPS * nz ; |
|
4683 for (i = 0 ; i < n ; i++) |
|
4684 @@ -860,7 +863,7 @@ |
|
4685 for (p = Ap [i] ; p < p2 ; p++) |
|
4686 { |
|
4687 /* zi -= Ax [p] * X [Ai [p]] ; */ |
|
4688 - ASSIGN (aij, Ax, Az, p, split) ; |
|
4689 + ASSIGN (aij, Ax, Az, p, AXsplit) ; |
|
4690 MULT_SUB (zi, aij, X [Ai [p]]) ; |
|
4691 } |
|
4692 Z [i] = zi ; |
|
4693 @@ -969,13 +972,13 @@ |
|
4694 for (i = 0 ; i < n ; i++) |
|
4695 { |
|
4696 /* wi = B [i] ; */ |
|
4697 - ASSIGN (wi, Bx, Bz, i, split) ; |
|
4698 + ASSIGN (wi, Bx, Bz, i, Bsplit) ; |
|
4699 z2i = 0. ; |
|
4700 p2 = Ap [i+1] ; |
|
4701 for (p = Ap [i] ; p < p2 ; p++) |
|
4702 { |
|
4703 /* axx = Ax [p] * X [Ai [p]] ; */ |
|
4704 - ASSIGN (aij, Ax, Az, p, split) ; |
|
4705 + ASSIGN (aij, Ax, Az, p, AXsplit) ; |
|
4706 MULT (axx, aij, X [Ai [p]]) ; |
|
4707 |
|
4708 /* wi -= axx ; */ |
|
4709 @@ -1011,7 +1014,7 @@ |
|
4710 for (i = 0 ; i < n ; i++) |
|
4711 { |
|
4712 /* X [i] = B [Rperm [i]] ; */ |
|
4713 - ASSIGN (X [i], Bx, Bz, Rperm [i], split) ; |
|
4714 + ASSIGN (X [i], Bx, Bz, Rperm [i], Bsplit) ; |
|
4715 } |
|
4716 flops = UMF_lsolve (Numeric, X, Pattern) ; |
|
4717 status = UMFPACK_OK ; |
|
4718 @@ -1027,7 +1030,7 @@ |
|
4719 for (i = 0 ; i < n ; i++) |
|
4720 { |
|
4721 /* X [i] = B [i] ; */ |
|
4722 - ASSIGN (X [i], Bx, Bz, i, split) ; |
|
4723 + ASSIGN (X [i], Bx, Bz, i, Bsplit) ; |
|
4724 } |
|
4725 flops = UMF_lsolve (Numeric, X, Pattern) ; |
|
4726 status = UMFPACK_OK ; |
|
4727 @@ -1043,7 +1046,7 @@ |
|
4728 for (i = 0 ; i < n ; i++) |
|
4729 { |
|
4730 /* W [i] = B [i] ; */ |
|
4731 - ASSIGN (W [i], Bx, Bz, i, split) ; |
|
4732 + ASSIGN (W [i], Bx, Bz, i, Bsplit) ; |
|
4733 } |
|
4734 flops = UMF_lhsolve (Numeric, W, Pattern) ; |
|
4735 for (i = 0 ; i < n ; i++) |
|
4736 @@ -1063,7 +1066,7 @@ |
|
4737 for (i = 0 ; i < n ; i++) |
|
4738 { |
|
4739 /* W [i] = B [i] ; */ |
|
4740 - ASSIGN (W [i], Bx, Bz, i, split) ; |
|
4741 + ASSIGN (W [i], Bx, Bz, i, Bsplit) ; |
|
4742 } |
|
4743 flops = UMF_ltsolve (Numeric, W, Pattern) ; |
|
4744 for (i = 0 ; i < n ; i++) |
|
4745 @@ -1083,7 +1086,7 @@ |
|
4746 for (i = 0 ; i < n ; i++) |
|
4747 { |
|
4748 /* X [i] = B [i] ; */ |
|
4749 - ASSIGN (X [i], Bx, Bz, i, split) ; |
|
4750 + ASSIGN (X [i], Bx, Bz, i, Bsplit) ; |
|
4751 } |
|
4752 flops = UMF_lhsolve (Numeric, X, Pattern) ; |
|
4753 status = UMFPACK_OK ; |
|
4754 @@ -1099,7 +1102,7 @@ |
|
4755 for (i = 0 ; i < n ; i++) |
|
4756 { |
|
4757 /* X [i] = B [i] ; */ |
|
4758 - ASSIGN (X [i], Bx, Bz, i, split) ; |
|
4759 + ASSIGN (X [i], Bx, Bz, i, Bsplit) ; |
|
4760 } |
|
4761 flops = UMF_ltsolve (Numeric, X, Pattern) ; |
|
4762 status = UMFPACK_OK ; |
|
4763 @@ -1115,7 +1118,7 @@ |
|
4764 for (i = 0 ; i < n ; i++) |
|
4765 { |
|
4766 /* W [i] = B [i] ; */ |
|
4767 - ASSIGN (W [i], Bx, Bz, i, split) ; |
|
4768 + ASSIGN (W [i], Bx, Bz, i, Bsplit) ; |
|
4769 } |
|
4770 flops = UMF_usolve (Numeric, W, Pattern) ; |
|
4771 for (i = 0 ; i < n ; i++) |
|
4772 @@ -1134,7 +1137,7 @@ |
|
4773 for (i = 0 ; i < n ; i++) |
|
4774 { |
|
4775 /* X [i] = B [i] ; */ |
|
4776 - ASSIGN (X [i], Bx, Bz, i, split) ; |
|
4777 + ASSIGN (X [i], Bx, Bz, i, Bsplit) ; |
|
4778 } |
|
4779 flops = UMF_usolve (Numeric, X, Pattern) ; |
|
4780 |
|
4781 @@ -1149,7 +1152,7 @@ |
|
4782 for (i = 0 ; i < n ; i++) |
|
4783 { |
|
4784 /* X [i] = B [Cperm [i]] ; */ |
|
4785 - ASSIGN (X [i], Bx, Bz, Cperm [i], split) ; |
|
4786 + ASSIGN (X [i], Bx, Bz, Cperm [i], Bsplit) ; |
|
4787 } |
|
4788 flops = UMF_uhsolve (Numeric, X, Pattern) ; |
|
4789 |
|
4790 @@ -1164,7 +1167,7 @@ |
|
4791 for (i = 0 ; i < n ; i++) |
|
4792 { |
|
4793 /* X [i] = B [Cperm [i]] ; */ |
|
4794 - ASSIGN (X [i], Bx, Bz, Cperm [i], split) ; |
|
4795 + ASSIGN (X [i], Bx, Bz, Cperm [i], Bsplit) ; |
|
4796 } |
|
4797 flops = UMF_utsolve (Numeric, X, Pattern) ; |
|
4798 |
|
4799 @@ -1179,7 +1182,7 @@ |
|
4800 for (i = 0 ; i < n ; i++) |
|
4801 { |
|
4802 /* X [i] = B [i] ; */ |
|
4803 - ASSIGN (X [i], Bx, Bz, i, split) ; |
|
4804 + ASSIGN (X [i], Bx, Bz, i, Bsplit) ; |
|
4805 } |
|
4806 flops = UMF_uhsolve (Numeric, X, Pattern) ; |
|
4807 |
|
4808 @@ -1194,7 +1197,7 @@ |
|
4809 for (i = 0 ; i < n ; i++) |
|
4810 { |
|
4811 /* X [i] = B [i] ; */ |
|
4812 - ASSIGN (X [i], Bx, Bz, i, split) ; |
|
4813 + ASSIGN (X [i], Bx, Bz, i, Bsplit) ; |
|
4814 } |
|
4815 flops = UMF_utsolve (Numeric, X, Pattern) ; |
|
4816 |
|
4817 @@ -1206,7 +1209,7 @@ |
|
4818 |
|
4819 #ifdef COMPLEX |
|
4820 /* copy the solution back, from Entry X [ ] to double Xx [ ] and Xz [ ] */ |
|
4821 - if (split) |
|
4822 + if (AXsplit) |
|
4823 { |
|
4824 for (i = 0 ; i < n ; i++) |
|
4825 { |