Mercurial > octave-nkf
view liboctave/UMFPACK.patch @ 5164:57077d0ddc8e
[project @ 2005-02-25 19:55:24 by jwe]
| author | jwe |
|---|---|
| date | Fri, 25 Feb 2005 19:55:28 +0000 |
| parents | |
| children |
line wrap: on
line source
--- UMFPACKv4.4.orig/AMD/Makefile 2004-01-29 20:40:44.000000000 +0100 +++ UMFPACK/AMD/Makefile 2005-02-01 22:00:38.917972334 +0100 @@ -28,11 +28,18 @@ ( cd Demo ; make cross ) - cat Doc/License +# compile a Octave version +# (not compiled by "make all") +octave: + ( cd OCTAVE ; make ) + - cat Doc/License + # remove object files, but keep the compiled programs and library archives clean: ( cd Source ; make clean ) ( cd Demo ; make clean ) ( cd MATLAB ; make clean ) + ( cd OCTAVE ; make clean ) ( cd Doc ; make clean ) # clean, and then remove compiled programs and library archives @@ -40,6 +47,7 @@ ( cd Source ; make purge ) ( cd Demo ; make purge ) ( cd MATLAB ; make purge ) + ( cd OCTAVE ; make purge ) ( cd Doc ; make purge ) # create PDF documents for the original distribution --- UMFPACKv4.4.orig/UMFPACK/Makefile 2004-01-29 20:40:48.000000000 +0100 +++ UMFPACK/UMFPACK/Makefile 2005-02-01 22:00:38.916972398 +0100 @@ -32,12 +32,19 @@ hb: ( cd Demo ; make hb ) +# compile a Octave version +# (not compiled by "make all") +octave: + ( cd OCTAVE ; make ) + - cat Doc/License + # remove object files, but keep the compiled programs and library archives clean: ( cd ../AMD ; make clean ) ( cd Source ; make clean ) ( cd Demo ; make clean ) ( cd MATLAB ; make clean ) + ( cd OCTAVE ; make clean ) ( cd Doc ; make clean ) # clean, and then remove compiled programs and library archives @@ -46,6 +53,7 @@ ( cd Source ; make purge ) ( cd Demo ; make purge ) ( cd MATLAB ; make purge ) + ( cd OCTAVE ; make purge ) ( cd Doc ; make purge ) # create PDF documents for the original distribution diff -uNr UMFPACKv4.4.orig/UMFPACK/OCTAVE/Contents.m UMFPACK/UMFPACK/OCTAVE/Contents.m --- UMFPACKv4.4.orig/UMFPACK/OCTAVE/Contents.m 1970-01-01 01:00:00.000000000 +0100 +++ UMFPACK/UMFPACK/OCTAVE/Contents.m 2004-12-30 01:58:46.000000000 +0100 @@ -0,1 +1,22 @@ +%Contents of the UMFPACK sparse matrix toolbox: +% +% umfpack computes x=A\b, x=A/b, or lu (A) for a sparse matrix A +% umfpack_make to compile umfpack for use in MATLAB +% umfpack_report prints optional control settings and statistics +% umfpack_demo a long demo +% umfpack_simple a simple demo +% umfpack_btf factorize A using a block triangular form +% umfpack_solve x = A\b or x = b/A +% lu_normest estimates norm (L*U-A, 1) without forming L*U-A +% luflop given L and U, computes # of flops required to compute them +% +% See also: +% amd symmetric minimum degree ordering +% colamd unsymmetric column approx minimum degree ordering +% symamd symmetric approx minimum degree ordering, based on colamd +% +% UMFPACK Version 4.3 (Jan. 16, 2004), Copyright (c) 2004 by Timothy A. +% Davis. All Rights Reserved. Type umfpack_details for License. + +help Contents --- UMFPACKv4.4.orig/UMFPACK/OCTAVE/GNUmakefile 1970-01-01 01:00:00.000000000 +0100 +++ UMFPACK/UMFPACK/OCTAVE/GNUmakefile 2004-12-30 01:58:46.000000000 +0100 @@ -0,0 +1,248 @@ +#------------------------------------------------------------------------------- +# UMFPACK GNUmakefile for the UMFPACK OCTAVE oct-file (GNU "make" only) +#------------------------------------------------------------------------------- + +.PRECIOUS: %.o + +# UMFPACK Version 4.3 (Jan. 16, 2004), Copyright (c) 2004 by Timothy A. +# Davis. All Rights Reserved. See ../README for License. + +all: umfpack luflop + +include ../Make/Make.include + +MKOCT = mkoctfile $(CONFIG) -DNRECIPROCAL -I/usr/include/atlas -I../Include -I../Source -I../../AMD/Include -I../../AMD/Source + +OCT_SPARSE_INC = -I../../../ + + +#------------------------------------------------------------------------------- +# source files +#------------------------------------------------------------------------------- + +# non-user-callable umf_*.[ch] files: +UMFCH = umf_assemble umf_blas3_update \ + umf_build_tuples umf_create_element \ + umf_dump umf_extend_front umf_garbage_collection \ + umf_get_memory umf_init_front umf_kernel \ + umf_kernel_init umf_kernel_wrapup \ + umf_local_search umf_lsolve umf_ltsolve \ + umf_mem_alloc_element umf_mem_alloc_head_block \ + umf_mem_alloc_tail_block umf_mem_free_tail_block \ + umf_mem_init_memoryspace \ + umf_report_vector umf_row_search umf_scale_column \ + umf_set_stats umf_solve umf_symbolic_usage umf_transpose \ + umf_tuple_lengths umf_usolve umf_utsolve umf_valid_numeric \ + umf_valid_symbolic umf_grow_front umf_start_front umf_2by2 \ + umf_store_lu umf_scale + +# non-user-callable umf_*.[ch] files, int/long versions only (no real/complex): +UMFINT = umf_analyze umf_apply_order umf_colamd umf_free umf_fsize \ + umf_is_permutation umf_malloc umf_realloc umf_report_perm \ + umf_singletons + +# non-user-callable and user-callable amd_*.[ch] files (int/long versions only): +AMD = amd_aat amd_1 amd_2 amd_dump amd_postorder amd_post_tree amd_defaults \ + amd_order amd_control amd_info amd_valid + +# non-user-callable, created from umf_ltsolve.c, umf_utsolve.c, +# umf_triplet.c, and umf_assemble.c , with int/long and real/complex versions: +UMF_CREATED = umf_lhsolve umf_uhsolve umf_triplet_map_nox \ + umf_triplet_nomap_x umf_triplet_nomap_nox umf_triplet_map_x \ + umf_assemble_fixq umf_store_lu_drop + +# non-user-callable, int/long and real/complex versions: +UMF = $(UMF_CREATED) $(UMFCH) + +# user-callable umfpack_*.[ch] files (int/long and real/complex): +UMFPACK = umfpack_col_to_triplet umfpack_defaults umfpack_free_numeric \ + umfpack_free_symbolic umfpack_get_numeric umfpack_get_lunz \ + umfpack_get_symbolic umfpack_get_determinant umfpack_numeric \ + umfpack_qsymbolic umfpack_report_control umfpack_report_info \ + umfpack_report_matrix umfpack_report_numeric umfpack_report_perm \ + umfpack_report_status umfpack_report_symbolic umfpack_report_triplet \ + umfpack_report_vector umfpack_solve umfpack_symbolic \ + umfpack_transpose umfpack_triplet_to_col umfpack_scale \ + umfpack_load_numeric umfpack_save_numeric \ + umfpack_load_symbolic umfpack_save_symbolic + +# user-callable, created from umfpack_solve.c (umfpack_wsolve.h exists, though): +# with int/long and real/complex versions: +UMFPACKW = umfpack_wsolve + +USER = $(UMFPACKW) $(UMFPACK) + +# user-callable, only one version for int/long, real/complex, *.[ch] files: +GENERIC = umfpack_timer umfpack_tictoc + +#------------------------------------------------------------------------------- +# include files: +#------------------------------------------------------------------------------- + +AMDH = ../../AMD/Source/amd_internal.h ../../AMD/Include/amd.h + +INC1 = umf_config.h umf_version.h umf_internal.h umf_triplet.h + +INC = ../Include/umfpack.h \ + $(addprefix ../Source/, $(INC1)) \ + $(addprefix ../Source/, $(addsuffix .h,$(UMFCH))) \ + $(addprefix ../Source/, $(addsuffix .h,$(UMFINT))) \ + $(addprefix ../Include/, $(addsuffix .h,$(USER))) \ + $(addprefix ../Include/, $(addsuffix .h,$(GENERIC))) \ + $(AMDH) + +#------------------------------------------------------------------------------- +# Create the umfpack and amd oct-file for OCTAVE (int versions only) +#------------------------------------------------------------------------------- + +OCTI = $(addsuffix .o, $(subst umf_,umf_o_,$(UMFINT))) +OCTDI = $(addsuffix .o, $(subst umf_,umf_od_,$(UMF)) $(subst umfpack_,umfpack_od_,$(USER))) +OCTZI = $(addsuffix .o, $(subst umf_,umf_oz_,$(UMF)) $(subst umfpack_,umfpack_oz_,$(USER)) ) +OCTAMD = $(addsuffix .o, $(subst amd_,amd_o_,$(AMD))) +OCTGN = $(addsuffix .o, $(subst umfpack_,umfpack_o_,$(GENERIC))) + +OCTUMFPACK = $(OCTI) $(OCTDI) $(OCTZI) $(OCTGN) + +OCTUMFPACK_LIB = umfpack_octave.o + +# Note that mkoctfile has an "-o" option, but it doesn't work in conjunction +# with the "-c" option, thus the need for $(MV) commands. +# If it did, then the rules would be much simpler: +# $(MKOCT) -DDINT -c $< -o $@ + +#---------------------------------------- +# integer-only routines (no real/complex): +#---------------------------------------- + +amd_o_%.o: ../../AMD/Source/amd_%.c $(AMDH) + $(MKOCT) -DDINT -c $< + - $(MV) ../../AMD/Source/amd_$*.o $@ + + +umf_o_%.o: ../Source/umf_%.c $(INC) + $(MKOCT) -DDINT -c $< + - $(MV) ../Source/umf_$*.o $@ + +#---------------------------------------- +# Double precision, int version, for OCTAVE +#---------------------------------------- + +umf_od_%.o: ../Source/umf_%.c $(INC) + $(MKOCT) -DDINT -c $< + - $(MV) ../Source/umf_$*.o $@ + +umf_od_%hsolve.o: ../Source/umf_%tsolve.c $(INC) + $(MKOCT) -DDINT -DCONJUGATE_SOLVE -c $< + - $(MV) ../Source/umf_$*tsolve.o $@ + +umf_od_triplet_map_x.o: ../Source/umf_triplet.c $(INC) + $(MKOCT) -DDINT -DDO_MAP -DDO_VALUES -c $< + - $(MV) ../Source/umf_triplet.o $@ + +umf_od_triplet_map_nox.o: ../Source/umf_triplet.c $(INC) + $(MKOCT) -DDINT -DDO_MAP -c $< + - $(MV) ../Source/umf_triplet.o $@ + +umf_od_triplet_nomap_x.o: ../Source/umf_triplet.c $(INC) + $(MKOCT) -DDINT -DDO_VALUES -c $< + - $(MV) ../Source/umf_triplet.o $@ + +umf_od_triplet_nomap_nox.o: ../Source/umf_triplet.c $(INC) + $(MKOCT) -DDINT -c $< + - $(MV) ../Source/umf_triplet.o $@ + +umf_od_assemble_fixq.o: ../Source/umf_assemble.c $(INC) + $(MKOCT) -DDINT -DFIXQ -c $< + - $(MV) ../Source/umf_assemble.o $@ + +umf_od_store_lu_drop.o: ../Source/umf_store_lu.c $(INC) + $(MKOCT) -DDINT -DDROP -c $< + - $(MV) ../Source/umf_store_lu.o $@ + +umfpack_od_wsolve.o: ../Source/umfpack_solve.c $(INC) + $(MKOCT) -DDINT -DWSOLVE -c $< + - $(MV) ../Source/umfpack_solve.o $@ + +umfpack_od_%.o: ../Source/umfpack_%.c $(INC) + $(MKOCT) -DDINT -c $< + - $(MV) ../Source/umfpack_$*.o $@ + +#---------------------------------------- +# Complex double precision, int version, for OCTAVE +#---------------------------------------- + +umf_oz_%.o: ../Source/umf_%.c $(INC) + $(MKOCT) -DZINT -c $< + - $(MV) ../Source/umf_$*.o $@ + +umf_oz_%hsolve.o: ../Source/umf_%tsolve.c $(INC) + $(MKOCT) -DZINT -DCONJUGATE_SOLVE -c $< + - $(MV) ../Source/umf_$*tsolve.o $@ + +umf_oz_triplet_map_x.o: ../Source/umf_triplet.c $(INC) + $(MKOCT) -DZINT -DDO_MAP -DDO_VALUES -c $< + - $(MV) ../Source/umf_triplet.o $@ + +umf_oz_triplet_map_nox.o: ../Source/umf_triplet.c $(INC) + $(MKOCT) -DZINT -DDO_MAP -c $< + - $(MV) ../Source/umf_triplet.o $@ + +umf_oz_triplet_nomap_x.o: ../Source/umf_triplet.c $(INC) + $(MKOCT) -DZINT -DDO_VALUES -c $< + - $(MV) ../Source/umf_triplet.o $@ + +umf_oz_triplet_nomap_nox.o: ../Source/umf_triplet.c $(INC) + $(MKOCT) -DZINT -c $< + - $(MV) ../Source/umf_triplet.o $@ + +umf_oz_assemble_fixq.o: ../Source/umf_assemble.c $(INC) + $(MKOCT) -DZINT -DFIXQ -c $< + - $(MV) ../Source/umf_assemble.o $@ + +umf_oz_store_lu_drop.o: ../Source/umf_store_lu.c $(INC) + $(MKOCT) -DZINT -DDROP -c $< + - $(MV) ../Source/umf_store_lu.o $@ + +umfpack_oz_wsolve.o: ../Source/umfpack_solve.c $(INC) + $(MKOCT) -DZINT -DWSOLVE -c $< + - $(MV) ../Source/umfpack_solve.o $@ + +umfpack_oz_%.o: ../Source/umfpack_%.c $(INC) + $(MKOCT) -DZINT -c $< + - $(MV) ../Source/umfpack_$*.o $@ + +#---------------------------------------- +# Generic routines for OCTAVE +#---------------------------------------- + +umfpack_o_timer.o: ../Source/umfpack_timer.c $(INC) + $(MKOCT) -c $< + - $(MV) ../Source/umfpack_timer.o $@ + +umfpack_o_tictoc.o: ../Source/umfpack_tictoc.c $(INC) + $(MKOCT) -c $< + - $(MV) ../Source/umfpack_tictoc.o $@ + +#---------------------------------------- +# umfpack oct-files +#---------------------------------------- + +umfpack: umfpack.cc $(OCTUMFPACK) $(OCTAMD) + $(MKOCT) $(OCT_SPARSE_INC) umfpack.cc $(OCTUMFPACK) $(OCTAMD) -o umfpack.oct + +luflop: luflop.cc + $(MKOCT) $(OCT_SPARSE_INC) luflop.cc -o luflop.oct + +#---------------------------------------- +# umfpack library to link with octave +#---------------------------------------- + +octave: $(OCTUMFPACK) $(OCTAMD) + ld -r $(OCTUMFPACK) $(OCTAMD) -o ../../../$(OCTUMFPACK_LIB) + +#------------------------------------------------------------------------------- +# Remove all but the files in the original distribution +#------------------------------------------------------------------------------- + +purge: clean + - $(RM) *.oct *.so diff -uNr UMFPACKv4.4.orig/UMFPACK/OCTAVE/luflop.cc UMFPACK/UMFPACK/OCTAVE/luflop.cc --- UMFPACKv4.4.orig/UMFPACK/OCTAVE/luflop.cc 1970-01-01 01:00:00.000000000 +0100 +++ UMFPACK/UMFPACK/OCTAVE/luflop.cc 2004-12-30 01:58:46.000000000 +0100 @@ -0,0 +1,223 @@ +/* + +Copyright (C) 2004 David Bateman + +This program is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +This program is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with this program; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +In addition to the terms of the GPL, you are permitted to link +this program with any Open Source program, as defined by the +Open Source Initiative (www.opensource.org) + +*/ + +/* + +This is the Octave interface to the UMFPACK code, which bore the following +copyright + + UMFPACK Version 4.3 (Jan. 16, 2004), Copyright (c) 2004 by Timothy A. + Davis. All Rights Reserved. See ../README for License. + email: davis@cise.ufl.edu CISE Department, Univ. of Florida. + web: http://www.cise.ufl.edu/research/sparse/umfpack + +*/ + + +#include <cstdlib> +#include <string> + +#include <octave/config.h> +#include <octave/ov.h> +#include <octave/defun-dld.h> +#include <octave/pager.h> +#include <octave/ov-re-mat.h> + +#include "ov-re-sparse.h" +#include "ov-cx-sparse.h" + +DEFUN_DLD (luflop, args, nargout, + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{f} =} luflup (@var{l}, @var{u})\n\ +\n\ +Given an LU factorization, compute how many flops took to compute it. This\n\ +is the same as (assuming U has a zero-free diagonal):\n\ +\n\ +@example\n\ +@group\n\ + Lnz = full (sum (spones (L))) - 1 ;\n\ + Unz = full (sum (spones (U')))' - 1 ;\n\ + f = 2*Lnz*Unz + sum (Lnz) ;\n\ +@end group\n\ +@end example\n\ +\n\ +except that no extra workspace is allocated for spones (L) and spones (U).\n\ +L and U must be sparse.\n\ +\n\ +Note: the above expression has a subtle undercount when exact numerical\n\ +cancelation occurs. Try [L,U,P] = lu (sparse (ones (10))) and then\n\ +luflop (L,U).\n\ +@end deftypefn") +{ + int nargin = args.length (); + octave_value retval; + + // These are here only so that the C++ destructors don't prematurally + // remove the underlying data we are interested in + SparseMatrix Lmatrix, Umatrix; + SparseComplexMatrix CLmatrix, CUmatrix; + int *Lp, *Li, *Up, *Ui; + int Lm, Ln, Um, Un; + + if (nargin != 2) + { + usage ("f = luflop (L, U)"); + return retval; + } + + if (args(0).class_name () == "sparse") + { + if (args(0).is_complex_type ()) + { + CLmatrix = + (((const octave_sparse_complex_matrix&) args(0).get_rep ()) + .sparse_complex_matrix_value ()); + Lp = CLmatrix.cidx (); + Li = CLmatrix.ridx (); + Lm = CLmatrix.rows (); + Ln = CLmatrix.cols (); + } + else + { + Lmatrix = (((const octave_sparse_matrix&) args(0).get_rep ()) + .sparse_matrix_value ()); + Lp = Lmatrix.cidx (); + Li = Lmatrix.ridx (); + Lm = Lmatrix.rows (); + Ln = Lmatrix.cols (); + } + } + else + { + if (args(0).is_complex_type ()) + { + CLmatrix = SparseComplexMatrix (args(0).complex_matrix_value ()); + Lp = CLmatrix.cidx (); + Li = CLmatrix.ridx (); + Lm = CLmatrix.rows (); + Ln = CLmatrix.cols (); + } + else + { + Lmatrix = SparseMatrix (args(0).matrix_value ()); + Lp = Lmatrix.cidx (); + Li = Lmatrix.ridx (); + Lm = Lmatrix.rows (); + Ln = Lmatrix.cols (); + } + } + + + if (args(0).class_name () == "sparse") + { + if (args(1).is_complex_type ()) + { + CUmatrix = + (((const octave_sparse_complex_matrix&) args(1).get_rep ()) + .sparse_complex_matrix_value ()); + Up = CUmatrix.cidx (); + Ui = CUmatrix.ridx (); + Um = CUmatrix.rows (); + Un = CUmatrix.cols (); + } + else + { + Umatrix = + (((const octave_sparse_matrix&) args(1).get_rep ()) + .sparse_matrix_value ()); + Up = Umatrix.cidx (); + Ui = Umatrix.ridx (); + Um = Umatrix.rows (); + Un = Umatrix.cols (); + } + } + else + { + if (args(1).is_complex_type ()) + { + CUmatrix = SparseComplexMatrix (args(1).complex_matrix_value ()); + Up = CUmatrix.cidx (); + Ui = CUmatrix.ridx (); + Um = CUmatrix.rows (); + Un = CUmatrix.cols (); + } + else + { + Umatrix = SparseMatrix (args(1).matrix_value ()); + Up = Umatrix.cidx (); + Ui = Umatrix.ridx (); + Um = Umatrix.rows (); + Un = Umatrix.cols (); + } + } + + + if (error_state) + return retval; + + + int n = Lm; + + if (n != Ln || n != Um || n != Un) + error ("luflop: L and U must be square"); + else + { + + OCTAVE_LOCAL_BUFFER (int, Unz, n); + + // count the nonzeros in each row of U + for (int row = 0 ; row < n ; row++) + Unz [row] = 0 ; + + for (int col = 0 ; col < n ; col++) + { + for (int p = Up [col] ; p < Up [col+1] ; p++) + { + int row = Ui [p] ; + Unz [row]++ ; + } + } + + // count the flops + double flop_count = 0.0 ; + for (int k = 0 ; k < n ; k++) + { + /* off-diagonal nonzeros in column k of L: */ + int Lnz_k = Lp [k+1] - Lp [k] - 1 ; + int Unz_k = Unz [k] - 1 ; + flop_count += (2 * Lnz_k * Unz_k) + Lnz_k ; + } + + // return the result + retval = flop_count; + } + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ diff -uNr UMFPACKv4.4.orig/UMFPACK/OCTAVE/lu_normest.m UMFPACK/UMFPACK/OCTAVE/lu_normest.m --- UMFPACKv4.4.orig/UMFPACK/OCTAVE/lu_normest.m 1970-01-01 01:00:00.000000000 +0100 +++ UMFPACK/UMFPACK/OCTAVE/lu_normest.m 2004-12-30 01:58:46.000000000 +0100 @@ -0,0 +1,106 @@ +function rho = lu_normest (A, L, U) +% LU_NORMEST: estimate the 1-norm of A-L*U without computing L*U +% +% Usage: +% +% rho = lu_normest (A, L, U) +% +% which estimates the computation of the 1-norm: +% +% rho = norm (A-L*U, 1) +% +% Authors: William W. Hager, Math Dept., Univ. of Florida +% Timothy A. Davis, CISE Dept., Univ. of Florida +% Gainesville, FL, 32611, USA. +% based on normest1, contributed on November, 1997 +% +% This code can be quite easily adapted to estimate the 1-norm of any +% matrix E, where E itself is dense or not explicitly represented, but the +% computation of E (and E') times a vector is easy. In this case, our matrix +% of interest is: +% +% E = A-L*U +% +% That is, L*U is the LU factorization of A, where A, L and U +% are sparse. This code works for dense matrices A and L too, +% but it would not be needed in that case, since E is easy to compute +% explicitly. For sparse A, L, and U, computing E explicitly would be quite +% expensive, and thus normest (A-L*U) would be prohibitive. +% +% For a detailed description, see Davis, T. A. and Hager, W. W., +% Modifying a sparse Cholesky factorization, SIAM J. Matrix Analysis and +% Applications, 1999, vol. 20, no. 3, 606-627. + +% The three places that the matrix-vector multiply E*x is used are highlighted. +% Note that E is never formed explicity. + +[m n] = size (A) ; + +if (m ~= n) + % pad A, L, and U with zeros so that they are all square + if (m < n) + U = [ U ; (sparse (n-m,n)) ] ; + L = [ L , (sparse (m,n-m)) ; (sparse (n-m,n)) ] ; + A = [ A ; (sparse (n-m,n)) ] ; + else + U = [ U , (sparse (n,m-n)) ; (sparse (m-n,m)) ] ; + L = [ L , (sparse (m,m-n)) ] ; + A = [ A , (sparse (m,m-n)) ] ; + end +end + +[m n] = size (A) ; + +notvisited = ones (m, 1) ; % nonvisited(j) is zero if j is visited, 1 otherwise +rho = 0 ; % the global rho + +At = A' ; +Lt = L' ; + +for trial = 1:3 % { + + x = notvisited ./ sum (notvisited) ; + rho1 = 0 ; % the current rho for this trial + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %%% COMPUTE Ex1 = E*x EFFICIENTLY: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + Ex1 = (A*x) - L*(U*x) ; + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + rho2 = norm (Ex1, 1) ; + + while rho2 > rho1 % { + + rho1 = rho2 ; + y = 2*(Ex1 >= 0) - 1 ; + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %%% COMPUTE z = E'*y EFFICIENTLY: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + z = (A'*y) - U'*(L'*y) ; + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + [zj, j] = max (abs (z .* notvisited)) ; + j = j (1) ; + if (abs (z (j)) > z'*x) % { + x = zeros (m, 1) ; + x (j) = 1 ; + notvisited (j) = 0 ; + else % } { + break ; + end % } + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %%% COMPUTE Ex1 = E*x EFFICIENTLY: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + Ex1 = (A*x) - L*(U*x) ; + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + rho2 = norm (Ex1, 1) ; + + end % } + + rho = max (rho, rho1) ; + +end % } diff -uNr UMFPACKv4.4.orig/UMFPACK/OCTAVE/Makefile UMFPACK/UMFPACK/OCTAVE/Makefile --- UMFPACKv4.4.orig/UMFPACK/OCTAVE/Makefile 1970-01-01 01:00:00.000000000 +0100 +++ UMFPACK/UMFPACK/OCTAVE/Makefile 2004-12-30 01:58:46.000000000 +0100 @@ -0,0 +1,517 @@ +#------------------------------------------------------------------------------- +# UMFPACK Makefile for the UMFPACK MATLAB mexFunction (old "make" only) +#------------------------------------------------------------------------------- + +# UMFPACK Version 4.3 (Jan. 16, 2004), Copyright (c) 2004 by Timothy A. +# Davis. All Rights Reserved. See ../README for License. + +# This is a very ugly Makefile, and is only provided for those who do not +# have GNU make. Note that it is not used if you have GNU make. It ignores +# dependency checking and just compiles everything. It was created +# automatically, via make -n using the GNUmakefile. That way, I don't have +# maintain two Makefiles. + +all: umfpack luflop + +include ../Make/Make.include + +MKOCT = mkoctfile $(CONFIG) -DNRECIPROCAL -I/usr/include/atlas -I../Include -I../Source -I../../AMD/Include -I../../AMD/Source + +OCT_SPARSE_INC = -I../../../ +OCTUMFPACK_LIB = umfpack_octave.o + +umfpack: + $(MKOCT) -DDINT -c ../Source/umf_analyze.c + $(MV) -f ../Source/umf_analyze.o umf_m_analyze.o + $(MKOCT) -DDINT -c ../Source/umf_apply_order.c + $(MV) -f ../Source/umf_apply_order.o umf_m_apply_order.o + $(MKOCT) -DDINT -c ../Source/umf_colamd.c + $(MV) -f ../Source/umf_colamd.o umf_m_colamd.o + $(MKOCT) -DDINT -c ../Source/umf_free.c + $(MV) -f ../Source/umf_free.o umf_m_free.o + $(MKOCT) -DDINT -c ../Source/umf_fsize.c + $(MV) -f ../Source/umf_fsize.o umf_m_fsize.o + $(MKOCT) -DDINT -c ../Source/umf_is_permutation.c + $(MV) -f ../Source/umf_is_permutation.o umf_m_is_permutation.o + $(MKOCT) -DDINT -c ../Source/umf_malloc.c + $(MV) -f ../Source/umf_malloc.o umf_m_malloc.o + $(MKOCT) -DDINT -c ../Source/umf_realloc.c + $(MV) -f ../Source/umf_realloc.o umf_m_realloc.o + $(MKOCT) -DDINT -c ../Source/umf_report_perm.c + $(MV) -f ../Source/umf_report_perm.o umf_m_report_perm.o + $(MKOCT) -DDINT -c ../Source/umf_singletons.c + $(MV) -f ../Source/umf_singletons.o umf_m_singletons.o + $(MKOCT) -DDINT -DCONJUGATE_SOLVE -c ../Source/umf_ltsolve.c + $(MV) -f ../Source/umf_ltsolve.o umf_od_lhsolve.o + $(MKOCT) -DDINT -DCONJUGATE_SOLVE -c ../Source/umf_utsolve.c + $(MV) -f ../Source/umf_utsolve.o umf_od_uhsolve.o + $(MKOCT) -DDINT -DDO_MAP -c ../Source/umf_triplet.c + $(MV) -f ../Source/umf_triplet.o umf_od_triplet_map_nox.o + $(MKOCT) -DDINT -DDO_VALUES -c ../Source/umf_triplet.c + $(MV) -f ../Source/umf_triplet.o umf_od_triplet_nomap_x.o + $(MKOCT) -DDINT -c ../Source/umf_triplet.c + $(MV) -f ../Source/umf_triplet.o umf_od_triplet_nomap_nox.o + $(MKOCT) -DDINT -DDO_MAP -DDO_VALUES -c ../Source/umf_triplet.c + $(MV) -f ../Source/umf_triplet.o umf_od_triplet_map_x.o + $(MKOCT) -DDINT -DFIXQ -c ../Source/umf_assemble.c + $(MV) -f ../Source/umf_assemble.o umf_od_assemble_fixq.o + $(MKOCT) -DDINT -DDROP -c ../Source/umf_store_lu.c + $(MV) -f ../Source/umf_store_lu.o umf_od_store_lu_drop.o + $(MKOCT) -DDINT -c ../Source/umf_assemble.c + $(MV) -f ../Source/umf_assemble.o umf_od_assemble.o + $(MKOCT) -DDINT -c ../Source/umf_blas3_update.c + $(MV) -f ../Source/umf_blas3_update.o umf_od_blas3_update.o + $(MKOCT) -DDINT -c ../Source/umf_build_tuples.c + $(MV) -f ../Source/umf_build_tuples.o umf_od_build_tuples.o + $(MKOCT) -DDINT -c ../Source/umf_create_element.c + $(MV) -f ../Source/umf_create_element.o umf_od_create_element.o + $(MKOCT) -DDINT -c ../Source/umf_dump.c + $(MV) -f ../Source/umf_dump.o umf_od_dump.o + $(MKOCT) -DDINT -c ../Source/umf_extend_front.c + $(MV) -f ../Source/umf_extend_front.o umf_od_extend_front.o + $(MKOCT) -DDINT -c ../Source/umf_garbage_collection.c + $(MV) -f ../Source/umf_garbage_collection.o umf_od_garbage_collection.o + $(MKOCT) -DDINT -c ../Source/umf_get_memory.c + $(MV) -f ../Source/umf_get_memory.o umf_od_get_memory.o + $(MKOCT) -DDINT -c ../Source/umf_init_front.c + $(MV) -f ../Source/umf_init_front.o umf_od_init_front.o + $(MKOCT) -DDINT -c ../Source/umf_kernel.c + $(MV) -f ../Source/umf_kernel.o umf_od_kernel.o + $(MKOCT) -DDINT -c ../Source/umf_kernel_init.c + $(MV) -f ../Source/umf_kernel_init.o umf_od_kernel_init.o + $(MKOCT) -DDINT -c ../Source/umf_kernel_wrapup.c + $(MV) -f ../Source/umf_kernel_wrapup.o umf_od_kernel_wrapup.o + $(MKOCT) -DDINT -c ../Source/umf_local_search.c + $(MV) -f ../Source/umf_local_search.o umf_od_local_search.o + $(MKOCT) -DDINT -c ../Source/umf_lsolve.c + $(MV) -f ../Source/umf_lsolve.o umf_od_lsolve.o + $(MKOCT) -DDINT -c ../Source/umf_ltsolve.c + $(MV) -f ../Source/umf_ltsolve.o umf_od_ltsolve.o + $(MKOCT) -DDINT -c ../Source/umf_mem_alloc_element.c + $(MV) -f ../Source/umf_mem_alloc_element.o umf_od_mem_alloc_element.o + $(MKOCT) -DDINT -c ../Source/umf_mem_alloc_head_block.c + $(MV) -f ../Source/umf_mem_alloc_head_block.o umf_od_mem_alloc_head_block.o + $(MKOCT) -DDINT -c ../Source/umf_mem_alloc_tail_block.c + $(MV) -f ../Source/umf_mem_alloc_tail_block.o umf_od_mem_alloc_tail_block.o + $(MKOCT) -DDINT -c ../Source/umf_mem_free_tail_block.c + $(MV) -f ../Source/umf_mem_free_tail_block.o umf_od_mem_free_tail_block.o + $(MKOCT) -DDINT -c ../Source/umf_mem_init_memoryspace.c + $(MV) -f ../Source/umf_mem_init_memoryspace.o umf_od_mem_init_memoryspace.o + $(MKOCT) -DDINT -c ../Source/umf_report_vector.c + $(MV) -f ../Source/umf_report_vector.o umf_od_report_vector.o + $(MKOCT) -DDINT -c ../Source/umf_row_search.c + $(MV) -f ../Source/umf_row_search.o umf_od_row_search.o + $(MKOCT) -DDINT -c ../Source/umf_scale_column.c + $(MV) -f ../Source/umf_scale_column.o umf_od_scale_column.o + $(MKOCT) -DDINT -c ../Source/umf_set_stats.c + $(MV) -f ../Source/umf_set_stats.o umf_od_set_stats.o + $(MKOCT) -DDINT -c ../Source/umf_solve.c + $(MV) -f ../Source/umf_solve.o umf_od_solve.o + $(MKOCT) -DDINT -c ../Source/umf_symbolic_usage.c + $(MV) -f ../Source/umf_symbolic_usage.o umf_od_symbolic_usage.o + $(MKOCT) -DDINT -c ../Source/umf_transpose.c + $(MV) -f ../Source/umf_transpose.o umf_od_transpose.o + $(MKOCT) -DDINT -c ../Source/umf_tuple_lengths.c + $(MV) -f ../Source/umf_tuple_lengths.o umf_od_tuple_lengths.o + $(MKOCT) -DDINT -c ../Source/umf_usolve.c + $(MV) -f ../Source/umf_usolve.o umf_od_usolve.o + $(MKOCT) -DDINT -c ../Source/umf_utsolve.c + $(MV) -f ../Source/umf_utsolve.o umf_od_utsolve.o + $(MKOCT) -DDINT -c ../Source/umf_valid_numeric.c + $(MV) -f ../Source/umf_valid_numeric.o umf_od_valid_numeric.o + $(MKOCT) -DDINT -c ../Source/umf_valid_symbolic.c + $(MV) -f ../Source/umf_valid_symbolic.o umf_od_valid_symbolic.o + $(MKOCT) -DDINT -c ../Source/umf_grow_front.c + $(MV) -f ../Source/umf_grow_front.o umf_od_grow_front.o + $(MKOCT) -DDINT -c ../Source/umf_start_front.c + $(MV) -f ../Source/umf_start_front.o umf_od_start_front.o + $(MKOCT) -DDINT -c ../Source/umf_2by2.c + $(MV) -f ../Source/umf_2by2.o umf_od_2by2.o + $(MKOCT) -DDINT -c ../Source/umf_store_lu.c + $(MV) -f ../Source/umf_store_lu.o umf_od_store_lu.o + $(MKOCT) -DDINT -c ../Source/umf_scale.c + $(MV) -f ../Source/umf_scale.o umf_od_scale.o + $(MKOCT) -DDINT -DWSOLVE -c ../Source/umfpack_solve.c + $(MV) -f ../Source/umfpack_solve.o umfpack_od_wsolve.o + $(MKOCT) -DDINT -c ../Source/umfpack_col_to_triplet.c + $(MV) -f ../Source/umfpack_col_to_triplet.o umfpack_od_col_to_triplet.o + $(MKOCT) -DDINT -c ../Source/umfpack_defaults.c + $(MV) -f ../Source/umfpack_defaults.o umfpack_od_defaults.o + $(MKOCT) -DDINT -c ../Source/umfpack_free_numeric.c + $(MV) -f ../Source/umfpack_free_numeric.o umfpack_od_free_numeric.o + $(MKOCT) -DDINT -c ../Source/umfpack_free_symbolic.c + $(MV) -f ../Source/umfpack_free_symbolic.o umfpack_od_free_symbolic.o + $(MKOCT) -DDINT -c ../Source/umfpack_get_numeric.c + $(MV) -f ../Source/umfpack_get_numeric.o umfpack_od_get_numeric.o + $(MKOCT) -DDINT -c ../Source/umfpack_get_lunz.c + $(MV) -f ../Source/umfpack_get_lunz.o umfpack_od_get_lunz.o + $(MKOCT) -DDINT -c ../Source/umfpack_get_symbolic.c + $(MV) -f ../Source/umfpack_get_symbolic.o umfpack_od_get_symbolic.o + $(MKOCT) -DDINT -c ../Source/umfpack_get_determinant.c + $(MV) -f ../Source/umfpack_get_determinant.o umfpack_od_get_determinant.o + $(MKOCT) -DDINT -c ../Source/umfpack_numeric.c + $(MV) -f ../Source/umfpack_numeric.o umfpack_od_numeric.o + $(MKOCT) -DDINT -c ../Source/umfpack_qsymbolic.c + $(MV) -f ../Source/umfpack_qsymbolic.o umfpack_od_qsymbolic.o + $(MKOCT) -DDINT -c ../Source/umfpack_report_control.c + $(MV) -f ../Source/umfpack_report_control.o umfpack_od_report_control.o + $(MKOCT) -DDINT -c ../Source/umfpack_report_info.c + $(MV) -f ../Source/umfpack_report_info.o umfpack_od_report_info.o + $(MKOCT) -DDINT -c ../Source/umfpack_report_matrix.c + $(MV) -f ../Source/umfpack_report_matrix.o umfpack_od_report_matrix.o + $(MKOCT) -DDINT -c ../Source/umfpack_report_numeric.c + $(MV) -f ../Source/umfpack_report_numeric.o umfpack_od_report_numeric.o + $(MKOCT) -DDINT -c ../Source/umfpack_report_perm.c + $(MV) -f ../Source/umfpack_report_perm.o umfpack_od_report_perm.o + $(MKOCT) -DDINT -c ../Source/umfpack_report_status.c + $(MV) -f ../Source/umfpack_report_status.o umfpack_od_report_status.o + $(MKOCT) -DDINT -c ../Source/umfpack_report_symbolic.c + $(MV) -f ../Source/umfpack_report_symbolic.o umfpack_od_report_symbolic.o + $(MKOCT) -DDINT -c ../Source/umfpack_report_triplet.c + $(MV) -f ../Source/umfpack_report_triplet.o umfpack_od_report_triplet.o + $(MKOCT) -DDINT -c ../Source/umfpack_report_vector.c + $(MV) -f ../Source/umfpack_report_vector.o umfpack_od_report_vector.o + $(MKOCT) -DDINT -c ../Source/umfpack_solve.c + $(MV) -f ../Source/umfpack_solve.o umfpack_od_solve.o + $(MKOCT) -DDINT -c ../Source/umfpack_symbolic.c + $(MV) -f ../Source/umfpack_symbolic.o umfpack_od_symbolic.o + $(MKOCT) -DDINT -c ../Source/umfpack_transpose.c + $(MV) -f ../Source/umfpack_transpose.o umfpack_od_transpose.o + $(MKOCT) -DDINT -c ../Source/umfpack_triplet_to_col.c + $(MV) -f ../Source/umfpack_triplet_to_col.o umfpack_od_triplet_to_col.o + $(MKOCT) -DDINT -c ../Source/umfpack_scale.c + $(MV) -f ../Source/umfpack_scale.o umfpack_od_scale.o + $(MKOCT) -DDINT -c ../Source/umfpack_load_numeric.c + $(MV) -f ../Source/umfpack_load_numeric.o umfpack_od_load_numeric.o + $(MKOCT) -DDINT -c ../Source/umfpack_save_numeric.c + $(MV) -f ../Source/umfpack_save_numeric.o umfpack_od_save_numeric.o + $(MKOCT) -DDINT -c ../Source/umfpack_load_symbolic.c + $(MV) -f ../Source/umfpack_load_symbolic.o umfpack_od_load_symbolic.o + $(MKOCT) -DDINT -c ../Source/umfpack_save_symbolic.c + $(MV) -f ../Source/umfpack_save_symbolic.o umfpack_od_save_symbolic.o + $(MKOCT) -DZINT -DCONJUGATE_SOLVE -c ../Source/umf_ltsolve.c + $(MV) -f ../Source/umf_ltsolve.o umf_oz_lhsolve.o + $(MKOCT) -DZINT -DCONJUGATE_SOLVE -c ../Source/umf_utsolve.c + $(MV) -f ../Source/umf_utsolve.o umf_oz_uhsolve.o + $(MKOCT) -DZINT -DDO_MAP -c ../Source/umf_triplet.c + $(MV) -f ../Source/umf_triplet.o umf_oz_triplet_map_nox.o + $(MKOCT) -DZINT -DDO_VALUES -c ../Source/umf_triplet.c + $(MV) -f ../Source/umf_triplet.o umf_oz_triplet_nomap_x.o + $(MKOCT) -DZINT -c ../Source/umf_triplet.c + $(MV) -f ../Source/umf_triplet.o umf_oz_triplet_nomap_nox.o + $(MKOCT) -DZINT -DDO_MAP -DDO_VALUES -c ../Source/umf_triplet.c + $(MV) -f ../Source/umf_triplet.o umf_oz_triplet_map_x.o + $(MKOCT) -DZINT -DFIXQ -c ../Source/umf_assemble.c + $(MV) -f ../Source/umf_assemble.o umf_oz_assemble_fixq.o + $(MKOCT) -DZINT -DDROP -c ../Source/umf_store_lu.c + $(MV) -f ../Source/umf_store_lu.o umf_oz_store_lu_drop.o + $(MKOCT) -DZINT -c ../Source/umf_assemble.c + $(MV) -f ../Source/umf_assemble.o umf_oz_assemble.o + $(MKOCT) -DZINT -c ../Source/umf_blas3_update.c + $(MV) -f ../Source/umf_blas3_update.o umf_oz_blas3_update.o + $(MKOCT) -DZINT -c ../Source/umf_build_tuples.c + $(MV) -f ../Source/umf_build_tuples.o umf_oz_build_tuples.o + $(MKOCT) -DZINT -c ../Source/umf_create_element.c + $(MV) -f ../Source/umf_create_element.o umf_oz_create_element.o + $(MKOCT) -DZINT -c ../Source/umf_dump.c + $(MV) -f ../Source/umf_dump.o umf_oz_dump.o + $(MKOCT) -DZINT -c ../Source/umf_extend_front.c + $(MV) -f ../Source/umf_extend_front.o umf_oz_extend_front.o + $(MKOCT) -DZINT -c ../Source/umf_garbage_collection.c + $(MV) -f ../Source/umf_garbage_collection.o umf_oz_garbage_collection.o + $(MKOCT) -DZINT -c ../Source/umf_get_memory.c + $(MV) -f ../Source/umf_get_memory.o umf_oz_get_memory.o + $(MKOCT) -DZINT -c ../Source/umf_init_front.c + $(MV) -f ../Source/umf_init_front.o umf_oz_init_front.o + $(MKOCT) -DZINT -c ../Source/umf_kernel.c + $(MV) -f ../Source/umf_kernel.o umf_oz_kernel.o + $(MKOCT) -DZINT -c ../Source/umf_kernel_init.c + $(MV) -f ../Source/umf_kernel_init.o umf_oz_kernel_init.o + $(MKOCT) -DZINT -c ../Source/umf_kernel_wrapup.c + $(MV) -f ../Source/umf_kernel_wrapup.o umf_oz_kernel_wrapup.o + $(MKOCT) -DZINT -c ../Source/umf_local_search.c + $(MV) -f ../Source/umf_local_search.o umf_oz_local_search.o + $(MKOCT) -DZINT -c ../Source/umf_lsolve.c + $(MV) -f ../Source/umf_lsolve.o umf_oz_lsolve.o + $(MKOCT) -DZINT -c ../Source/umf_ltsolve.c + $(MV) -f ../Source/umf_ltsolve.o umf_oz_ltsolve.o + $(MKOCT) -DZINT -c ../Source/umf_mem_alloc_element.c + $(MV) -f ../Source/umf_mem_alloc_element.o umf_oz_mem_alloc_element.o + $(MKOCT) -DZINT -c ../Source/umf_mem_alloc_head_block.c + $(MV) -f ../Source/umf_mem_alloc_head_block.o umf_oz_mem_alloc_head_block.o + $(MKOCT) -DZINT -c ../Source/umf_mem_alloc_tail_block.c + $(MV) -f ../Source/umf_mem_alloc_tail_block.o umf_oz_mem_alloc_tail_block.o + $(MKOCT) -DZINT -c ../Source/umf_mem_free_tail_block.c + $(MV) -f ../Source/umf_mem_free_tail_block.o umf_oz_mem_free_tail_block.o + $(MKOCT) -DZINT -c ../Source/umf_mem_init_memoryspace.c + $(MV) -f ../Source/umf_mem_init_memoryspace.o umf_oz_mem_init_memoryspace.o + $(MKOCT) -DZINT -c ../Source/umf_report_vector.c + $(MV) -f ../Source/umf_report_vector.o umf_oz_report_vector.o + $(MKOCT) -DZINT -c ../Source/umf_row_search.c + $(MV) -f ../Source/umf_row_search.o umf_oz_row_search.o + $(MKOCT) -DZINT -c ../Source/umf_scale_column.c + $(MV) -f ../Source/umf_scale_column.o umf_oz_scale_column.o + $(MKOCT) -DZINT -c ../Source/umf_set_stats.c + $(MV) -f ../Source/umf_set_stats.o umf_oz_set_stats.o + $(MKOCT) -DZINT -c ../Source/umf_solve.c + $(MV) -f ../Source/umf_solve.o umf_oz_solve.o + $(MKOCT) -DZINT -c ../Source/umf_symbolic_usage.c + $(MV) -f ../Source/umf_symbolic_usage.o umf_oz_symbolic_usage.o + $(MKOCT) -DZINT -c ../Source/umf_transpose.c + $(MV) -f ../Source/umf_transpose.o umf_oz_transpose.o + $(MKOCT) -DZINT -c ../Source/umf_tuple_lengths.c + $(MV) -f ../Source/umf_tuple_lengths.o umf_oz_tuple_lengths.o + $(MKOCT) -DZINT -c ../Source/umf_usolve.c + $(MV) -f ../Source/umf_usolve.o umf_oz_usolve.o + $(MKOCT) -DZINT -c ../Source/umf_utsolve.c + $(MV) -f ../Source/umf_utsolve.o umf_oz_utsolve.o + $(MKOCT) -DZINT -c ../Source/umf_valid_numeric.c + $(MV) -f ../Source/umf_valid_numeric.o umf_oz_valid_numeric.o + $(MKOCT) -DZINT -c ../Source/umf_valid_symbolic.c + $(MV) -f ../Source/umf_valid_symbolic.o umf_oz_valid_symbolic.o + $(MKOCT) -DZINT -c ../Source/umf_grow_front.c + $(MV) -f ../Source/umf_grow_front.o umf_oz_grow_front.o + $(MKOCT) -DZINT -c ../Source/umf_start_front.c + $(MV) -f ../Source/umf_start_front.o umf_oz_start_front.o + $(MKOCT) -DZINT -c ../Source/umf_2by2.c + $(MV) -f ../Source/umf_2by2.o umf_oz_2by2.o + $(MKOCT) -DZINT -c ../Source/umf_store_lu.c + $(MV) -f ../Source/umf_store_lu.o umf_oz_store_lu.o + $(MKOCT) -DZINT -c ../Source/umf_scale.c + $(MV) -f ../Source/umf_scale.o umf_oz_scale.o + $(MKOCT) -DZINT -DWSOLVE -c ../Source/umfpack_solve.c + $(MV) -f ../Source/umfpack_solve.o umfpack_oz_wsolve.o + $(MKOCT) -DZINT -c ../Source/umfpack_col_to_triplet.c + $(MV) -f ../Source/umfpack_col_to_triplet.o umfpack_oz_col_to_triplet.o + $(MKOCT) -DZINT -c ../Source/umfpack_defaults.c + $(MV) -f ../Source/umfpack_defaults.o umfpack_oz_defaults.o + $(MKOCT) -DZINT -c ../Source/umfpack_free_numeric.c + $(MV) -f ../Source/umfpack_free_numeric.o umfpack_oz_free_numeric.o + $(MKOCT) -DZINT -c ../Source/umfpack_free_symbolic.c + $(MV) -f ../Source/umfpack_free_symbolic.o umfpack_oz_free_symbolic.o + $(MKOCT) -DZINT -c ../Source/umfpack_get_numeric.c + $(MV) -f ../Source/umfpack_get_numeric.o umfpack_oz_get_numeric.o + $(MKOCT) -DZINT -c ../Source/umfpack_get_lunz.c + $(MV) -f ../Source/umfpack_get_lunz.o umfpack_oz_get_lunz.o + $(MKOCT) -DZINT -c ../Source/umfpack_get_symbolic.c + $(MV) -f ../Source/umfpack_get_symbolic.o umfpack_oz_get_symbolic.o + $(MKOCT) -DZINT -c ../Source/umfpack_get_determinant.c + $(MV) -f ../Source/umfpack_get_determinant.o umfpack_oz_get_determinant.o + $(MKOCT) -DZINT -c ../Source/umfpack_numeric.c + $(MV) -f ../Source/umfpack_numeric.o umfpack_oz_numeric.o + $(MKOCT) -DZINT -c ../Source/umfpack_qsymbolic.c + $(MV) -f ../Source/umfpack_qsymbolic.o umfpack_oz_qsymbolic.o + $(MKOCT) -DZINT -c ../Source/umfpack_report_control.c + $(MV) -f ../Source/umfpack_report_control.o umfpack_oz_report_control.o + $(MKOCT) -DZINT -c ../Source/umfpack_report_info.c + $(MV) -f ../Source/umfpack_report_info.o umfpack_oz_report_info.o + $(MKOCT) -DZINT -c ../Source/umfpack_report_matrix.c + $(MV) -f ../Source/umfpack_report_matrix.o umfpack_oz_report_matrix.o + $(MKOCT) -DZINT -c ../Source/umfpack_report_numeric.c + $(MV) -f ../Source/umfpack_report_numeric.o umfpack_oz_report_numeric.o + $(MKOCT) -DZINT -c ../Source/umfpack_report_perm.c + $(MV) -f ../Source/umfpack_report_perm.o umfpack_oz_report_perm.o + $(MKOCT) -DZINT -c ../Source/umfpack_report_status.c + $(MV) -f ../Source/umfpack_report_status.o umfpack_oz_report_status.o + $(MKOCT) -DZINT -c ../Source/umfpack_report_symbolic.c + $(MV) -f ../Source/umfpack_report_symbolic.o umfpack_oz_report_symbolic.o + $(MKOCT) -DZINT -c ../Source/umfpack_report_triplet.c + $(MV) -f ../Source/umfpack_report_triplet.o umfpack_oz_report_triplet.o + $(MKOCT) -DZINT -c ../Source/umfpack_report_vector.c + $(MV) -f ../Source/umfpack_report_vector.o umfpack_oz_report_vector.o + $(MKOCT) -DZINT -c ../Source/umfpack_solve.c + $(MV) -f ../Source/umfpack_solve.o umfpack_oz_solve.o + $(MKOCT) -DZINT -c ../Source/umfpack_symbolic.c + $(MV) -f ../Source/umfpack_symbolic.o umfpack_oz_symbolic.o + $(MKOCT) -DZINT -c ../Source/umfpack_transpose.c + $(MV) -f ../Source/umfpack_transpose.o umfpack_oz_transpose.o + $(MKOCT) -DZINT -c ../Source/umfpack_triplet_to_col.c + $(MV) -f ../Source/umfpack_triplet_to_col.o umfpack_oz_triplet_to_col.o + $(MKOCT) -DZINT -c ../Source/umfpack_scale.c + $(MV) -f ../Source/umfpack_scale.o umfpack_oz_scale.o + $(MKOCT) -DZINT -c ../Source/umfpack_load_numeric.c + $(MV) -f ../Source/umfpack_load_numeric.o umfpack_oz_load_numeric.o + $(MKOCT) -DZINT -c ../Source/umfpack_save_numeric.c + $(MV) -f ../Source/umfpack_save_numeric.o umfpack_oz_save_numeric.o + $(MKOCT) -DZINT -c ../Source/umfpack_load_symbolic.c + $(MV) -f ../Source/umfpack_load_symbolic.o umfpack_oz_load_symbolic.o + $(MKOCT) -DZINT -c ../Source/umfpack_save_symbolic.c + $(MV) -f ../Source/umfpack_save_symbolic.o umfpack_oz_save_symbolic.o + $(MKOCT) -c ../Source/umfpack_timer.c + $(MV) -f ../Source/umfpack_timer.o umfpack_m_timer.o + $(MKOCT) -c ../Source/umfpack_tictoc.c + $(MV) -f ../Source/umfpack_tictoc.o umfpack_m_tictoc.o + $(MKOCT) -DDINT -c ../../AMD/Source/amd_aat.c + $(MV) -f ../../AMD/Source/amd_aat.o amd_m_aat.o + $(MKOCT) -DDINT -c ../../AMD/Source/amd_1.c + $(MV) -f ../../AMD/Source/amd_1.o amd_m_1.o + $(MKOCT) -DDINT -c ../../AMD/Source/amd_2.c + $(MV) -f ../../AMD/Source/amd_2.o amd_m_2.o + $(MKOCT) -DDINT -c ../../AMD/Source/amd_dump.c + $(MV) -f ../../AMD/Source/amd_dump.o amd_m_dump.o + $(MKOCT) -DDINT -c ../../AMD/Source/amd_postorder.c + $(MV) -f ../../AMD/Source/amd_postorder.o amd_m_postorder.o + $(MKOCT) -DDINT -c ../../AMD/Source/amd_post_tree.c + $(MV) -f ../../AMD/Source/amd_post_tree.o amd_m_post_tree.o + $(MKOCT) -DDINT -c ../../AMD/Source/amd_defaults.c + $(MV) -f ../../AMD/Source/amd_defaults.o amd_m_defaults.o + $(MKOCT) -DDINT -c ../../AMD/Source/amd_order.c + $(MV) -f ../../AMD/Source/amd_order.o amd_m_order.o + $(MKOCT) -DDINT -c ../../AMD/Source/amd_control.c + $(MV) -f ../../AMD/Source/amd_control.o amd_m_control.o + $(MKOCT) -DDINT -c ../../AMD/Source/amd_info.c + $(MV) -f ../../AMD/Source/amd_info.o amd_m_info.o + $(MKOCT) -DDINT -c ../../AMD/Source/amd_valid.c + $(MV) -f ../../AMD/Source/amd_valid.o amd_m_valid.o + $(MKOCT) -o umfpack.oct $(OCT_SPARSE_INC) umfpack.cc \ + umf_o_analyze.o umf_o_apply_order.o umf_o_colamd.o umf_o_free.o \ + umf_o_fsize.o umf_o_is_permutation.o umf_o_malloc.o \ + umf_o_realloc.o umf_o_report_perm.o umf_o_singletons.o \ + umf_od_lhsolve.o umf_od_uhsolve.o umf_od_triplet_map_nox.o \ + umf_od_triplet_nomap_x.o umf_od_triplet_nomap_nox.o \ + umf_od_triplet_map_x.o umf_od_assemble_fixq.o \ + umf_od_store_lu_drop.o umf_od_assemble.o umf_od_blas3_update.o \ + umf_od_build_tuples.o umf_od_create_element.o umf_od_dump.o \ + umf_od_extend_front.o umf_od_garbage_collection.o \ + umf_od_get_memory.o umf_od_init_front.o umf_od_kernel.o \ + umf_od_kernel_init.o umf_od_kernel_wrapup.o umf_od_local_search.o \ + umf_od_lsolve.o umf_od_ltsolve.o umf_od_mem_alloc_element.o \ + umf_od_mem_alloc_head_block.o umf_od_mem_alloc_tail_block.o \ + umf_od_mem_free_tail_block.o umf_od_mem_init_memoryspace.o \ + umf_od_report_vector.o umf_od_row_search.o umf_od_scale_column.o \ + umf_od_set_stats.o umf_od_solve.o umf_od_symbolic_usage.o \ + umf_od_transpose.o umf_od_tuple_lengths.o umf_od_usolve.o \ + umf_od_utsolve.o umf_od_valid_numeric.o umf_od_valid_symbolic.o \ + umf_od_grow_front.o umf_od_start_front.o umf_od_2by2.o \ + umf_od_store_lu.o umf_od_scale.o umfpack_od_wsolve.o \ + umfpack_od_col_to_triplet.o umfpack_od_defaults.o \ + umfpack_od_free_numeric.o umfpack_od_free_symbolic.o \ + umfpack_od_get_numeric.o umfpack_od_get_lunz.o \ + umfpack_od_get_symbolic.o umfpack_od_numeric.o \ + umfpack_od_qsymbolic.o umfpack_od_report_control.o \ + umfpack_od_report_info.o umfpack_od_report_matrix.o \ + umfpack_od_report_numeric.o umfpack_od_report_perm.o \ + umfpack_od_report_status.o umfpack_od_report_symbolic.o \ + umfpack_od_report_triplet.o umfpack_od_report_vector.o \ + umfpack_od_solve.o umfpack_od_symbolic.o umfpack_od_transpose.o \ + umfpack_od_triplet_to_col.o umfpack_od_scale.o \ + umfpack_od_load_numeric.o umfpack_od_save_numeric.o \ + umfpack_od_load_symbolic.o umfpack_od_save_symbolic.o \ + umf_oz_lhsolve.o umf_oz_uhsolve.o umf_oz_triplet_map_nox.o \ + umf_oz_triplet_nomap_x.o umf_oz_triplet_nomap_nox.o \ + umf_oz_triplet_map_x.o umf_oz_assemble_fixq.o \ + umf_oz_store_lu_drop.o umf_oz_assemble.o umf_oz_blas3_update.o \ + umf_oz_build_tuples.o umf_oz_create_element.o umf_oz_dump.o \ + umf_oz_extend_front.o umf_oz_garbage_collection.o \ + umf_oz_get_memory.o umf_oz_init_front.o umf_oz_kernel.o \ + umf_oz_kernel_init.o umf_oz_kernel_wrapup.o umf_oz_local_search.o \ + umf_oz_lsolve.o umf_oz_ltsolve.o umf_oz_mem_alloc_element.o \ + umf_oz_mem_alloc_head_block.o umf_oz_mem_alloc_tail_block.o \ + umf_oz_mem_free_tail_block.o umf_oz_mem_init_memoryspace.o \ + umf_oz_report_vector.o umf_oz_row_search.o umf_oz_scale_column.o \ + umf_oz_set_stats.o umf_oz_solve.o umf_oz_symbolic_usage.o \ + umf_oz_transpose.o umf_oz_tuple_lengths.o umf_oz_usolve.o \ + umf_oz_utsolve.o umf_oz_valid_numeric.o umf_oz_valid_symbolic.o \ + umf_oz_grow_front.o umf_oz_start_front.o umf_oz_2by2.o \ + umf_oz_store_lu.o umf_oz_scale.o umfpack_oz_wsolve.o \ + umfpack_oz_col_to_triplet.o umfpack_oz_defaults.o \ + umfpack_oz_free_numeric.o umfpack_oz_free_symbolic.o \ + umfpack_oz_get_numeric.o umfpack_oz_get_lunz.o \ + umfpack_oz_get_symbolic.o umfpack_oz_numeric.o \ + umfpack_oz_qsymbolic.o umfpack_oz_report_control.o \ + umfpack_oz_report_info.o umfpack_oz_report_matrix.o \ + umfpack_oz_report_numeric.o umfpack_oz_report_perm.o \ + umfpack_oz_report_status.o umfpack_oz_report_symbolic.o \ + umfpack_oz_report_triplet.o umfpack_oz_report_vector.o \ + umfpack_oz_solve.o umfpack_oz_symbolic.o umfpack_oz_transpose.o \ + umfpack_oz_triplet_to_col.o umfpack_oz_scale.o \ + umfpack_oz_load_numeric.o umfpack_oz_save_numeric.o \ + umfpack_oz_load_symbolic.o umfpack_oz_save_symbolic.o \ + umfpack_o_timer.o umfpack_o_tictoc.o \ + amd_o_aat.o amd_o_1.o amd_o_2.o amd_o_dump.o \ + amd_o_postorder.o amd_o_post_tree.o amd_o_defaults.o amd_o_order.o \ + amd_o_control.o amd_o_info.o amd_o_valid.o + +luflop: luflop.cc + $(MKOCT) luflop.cc -I$(OCT_SPARSE_INC) -o luflop.oct + +#---------------------------------------- +# umfpack library to link with octave +#---------------------------------------- + +octave: umfpack + ld -r -o ../../../$(OCTUMFPACK_LIB) \ + umf_o_analyze.o umf_o_apply_order.o umf_o_colamd.o umf_o_free.o \ + umf_o_fsize.o umf_o_is_permutation.o umf_o_malloc.o \ + umf_o_realloc.o umf_o_report_perm.o umf_o_singletons.o \ + umf_od_lhsolve.o umf_od_uhsolve.o umf_od_triplet_map_nox.o \ + umf_od_triplet_nomap_x.o umf_od_triplet_nomap_nox.o \ + umf_od_triplet_map_x.o umf_od_assemble_fixq.o \ + umf_od_store_lu_drop.o umf_od_assemble.o umf_od_blas3_update.o \ + umf_od_build_tuples.o umf_od_create_element.o umf_od_dump.o \ + umf_od_extend_front.o umf_od_garbage_collection.o \ + umf_od_get_memory.o umf_od_init_front.o umf_od_kernel.o \ + umf_od_kernel_init.o umf_od_kernel_wrapup.o umf_od_local_search.o \ + umf_od_lsolve.o umf_od_ltsolve.o umf_od_mem_alloc_element.o \ + umf_od_mem_alloc_head_block.o umf_od_mem_alloc_tail_block.o \ + umf_od_mem_free_tail_block.o umf_od_mem_init_memoryspace.o \ + umf_od_report_vector.o umf_od_row_search.o umf_od_scale_column.o \ + umf_od_set_stats.o umf_od_solve.o umf_od_symbolic_usage.o \ + umf_od_transpose.o umf_od_tuple_lengths.o umf_od_usolve.o \ + umf_od_utsolve.o umf_od_valid_numeric.o umf_od_valid_symbolic.o \ + umf_od_grow_front.o umf_od_start_front.o umf_od_2by2.o \ + umf_od_store_lu.o umf_od_scale.o umfpack_od_wsolve.o \ + umfpack_od_col_to_triplet.o umfpack_od_defaults.o \ + umfpack_od_free_numeric.o umfpack_od_free_symbolic.o \ + umfpack_od_get_numeric.o umfpack_od_get_lunz.o \ + umfpack_od_get_symbolic.o umfpack_od_numeric.o \ + umfpack_od_qsymbolic.o umfpack_od_report_control.o \ + umfpack_od_report_info.o umfpack_od_report_matrix.o \ + umfpack_od_report_numeric.o umfpack_od_report_perm.o \ + umfpack_od_report_status.o umfpack_od_report_symbolic.o \ + umfpack_od_report_triplet.o umfpack_od_report_vector.o \ + umfpack_od_solve.o umfpack_od_symbolic.o umfpack_od_transpose.o \ + umfpack_od_triplet_to_col.o umfpack_od_scale.o \ + umfpack_od_load_numeric.o umfpack_od_save_numeric.o \ + umfpack_od_load_symbolic.o umfpack_od_save_symbolic.o \ + umf_oz_lhsolve.o umf_oz_uhsolve.o umf_oz_triplet_map_nox.o \ + umf_oz_triplet_nomap_x.o umf_oz_triplet_nomap_nox.o \ + umf_oz_triplet_map_x.o umf_oz_assemble_fixq.o \ + umf_oz_store_lu_drop.o umf_oz_assemble.o umf_oz_blas3_update.o \ + umf_oz_build_tuples.o umf_oz_create_element.o umf_oz_dump.o \ + umf_oz_extend_front.o umf_oz_garbage_collection.o \ + umf_oz_get_memory.o umf_oz_init_front.o umf_oz_kernel.o \ + umf_oz_kernel_init.o umf_oz_kernel_wrapup.o umf_oz_local_search.o \ + umf_oz_lsolve.o umf_oz_ltsolve.o umf_oz_mem_alloc_element.o \ + umf_oz_mem_alloc_head_block.o umf_oz_mem_alloc_tail_block.o \ + umf_oz_mem_free_tail_block.o umf_oz_mem_init_memoryspace.o \ + umf_oz_report_vector.o umf_oz_row_search.o umf_oz_scale_column.o \ + umf_oz_set_stats.o umf_oz_solve.o umf_oz_symbolic_usage.o \ + umf_oz_transpose.o umf_oz_tuple_lengths.o umf_oz_usolve.o \ + umf_oz_utsolve.o umf_oz_valid_numeric.o umf_oz_valid_symbolic.o \ + umf_oz_grow_front.o umf_oz_start_front.o umf_oz_2by2.o \ + umf_oz_store_lu.o umf_oz_scale.o umfpack_oz_wsolve.o \ + umfpack_oz_col_to_triplet.o umfpack_oz_defaults.o \ + umfpack_oz_free_numeric.o umfpack_oz_free_symbolic.o \ + umfpack_oz_get_numeric.o umfpack_oz_get_lunz.o \ + umfpack_oz_get_symbolic.o umfpack_oz_numeric.o \ + umfpack_oz_qsymbolic.o umfpack_oz_report_control.o \ + umfpack_oz_report_info.o umfpack_oz_report_matrix.o \ + umfpack_oz_report_numeric.o umfpack_oz_report_perm.o \ + umfpack_oz_report_status.o umfpack_oz_report_symbolic.o \ + umfpack_oz_report_triplet.o umfpack_oz_report_vector.o \ + umfpack_oz_solve.o umfpack_oz_symbolic.o umfpack_oz_transpose.o \ + umfpack_oz_triplet_to_col.o umfpack_oz_scale.o \ + umfpack_oz_load_numeric.o umfpack_oz_save_numeric.o \ + umfpack_oz_load_symbolic.o umfpack_oz_save_symbolic.o \ + umfpack_o_timer.o umfpack_o_tictoc.o \ + amd_o_aat.o amd_o_1.o amd_o_2.o amd_o_dump.o \ + amd_o_postorder.o amd_o_post_tree.o amd_o_defaults.o amd_o_order.o \ + amd_o_control.o amd_o_info.o amd_o_valid.o + +#------------------------------------------------------------------------------- +# Remove all but the files in the original distribution +#------------------------------------------------------------------------------- + +purge: clean + - $(RM) *.oct* *.dll diff -uNr UMFPACKv4.4.orig/UMFPACK/OCTAVE/umfpack_btf.m UMFPACK/UMFPACK/OCTAVE/umfpack_btf.m --- UMFPACKv4.4.orig/UMFPACK/OCTAVE/umfpack_btf.m 1970-01-01 01:00:00.000000000 +0100 +++ UMFPACK/UMFPACK/OCTAVE/umfpack_btf.m 2004-12-30 01:58:46.000000000 +0100 @@ -0,0 +1,129 @@ +function x = umfpack_btf (A, b, Control) +% UMFPACK_BTF +% +% x = umfpack_btf (A, b, Control) +% +% solve Ax=b by first permuting the matrix A to block triangular form via dmperm +% and then using UMFPACK to factorize each diagonal block. Adjacent 1-by-1 +% blocks are merged into a single upper triangular block, and solved via +% MATLAB's \ operator. The Control parameter is optional (Type umfpack_details +% and umfpack_report for details on its use). A must be square. +% +% See also: umfpack, umfpack_factorize, umfpack_details, dmperm + +% UMFPACK Version 4.3 (Jan. 16, 2004), Copyright (c) 2004 by Timothy A. +% Davis. All Rights Reserved. Type umfpack_details for License. + +if (nargin < 2) + help umfpack_btf + error ('Usage: x = umfpack_btf (A, b, Control)') ; +end + +[m n] = size (A) ; +if (m ~= n) + help umfpack_btf + error ('umfpack_btf: A must be square') ; +end +[m1 n1] = size (b) ; +if (m1 ~= n) + help umfpack_btf + error ('umfpack_btf: b has the wrong dimensions') ; +end + +if (nargin < 3) + Control = umfpack ; +end + +%------------------------------------------------------------------------------- +% find the block triangular form +%------------------------------------------------------------------------------- + +[p,q,r] = dmperm (A) ; +nblocks = length (r) - 1 ; + +%------------------------------------------------------------------------------- +% solve the system +%------------------------------------------------------------------------------- + +if (nblocks == 1 | sprank (A) < n) + + %--------------------------------------------------------------------------- + % matrix is irreducible or structurally singular + %--------------------------------------------------------------------------- + + x = umfpack_solve (A, '\\', b, Control) ; + +else + + %--------------------------------------------------------------------------- + % A (p,q) is in block triangular form + %--------------------------------------------------------------------------- + + b = b (p,:) ; + A = A (p,q) ; + x = zeros (size (b)) ; + + %--------------------------------------------------------------------------- + % merge adjacent singletons into a single upper triangular block + %--------------------------------------------------------------------------- + + [r, nblocks, is_triangular] = merge_singletons (r) ; + + %--------------------------------------------------------------------------- + % solve the system: x (q) = A\b + %--------------------------------------------------------------------------- + + for k = nblocks:-1:1 + + % get the kth block + k1 = r (k) ; + k2 = r (k+1) - 1 ; + + % solve the system + x (k1:k2,:) = solver (A (k1:k2, k1:k2), b (k1:k2,:), ... + is_triangular (k), Control) ; + + % off-diagonal block back substitution + b (1:k1-1,:) = b (1:k1-1,:) - A (1:k1-1, k1:k2) * x (k1:k2,:) ; + + end + + x (q,:) = x ; + +end + +%------------------------------------------------------------------------------- +% merge_singletons +%------------------------------------------------------------------------------- + +function [r, nblocks, is_triangular] = merge_singletons (r) +% +% Given r from [p,q,r] = dmperm (A), where A is square, return a modified r that +% reflects the merger of adjacent singletons into a single upper triangular +% block. is_triangular (k) is 1 if the kth block is upper triangular. nblocks +% is the number of new blocks. + +nblocks = length (r) - 1 ; +bsize = r (2:nblocks+1) - r (1:nblocks) ; +t = [0 (bsize == 1)] ; +z = (t (1:nblocks) == 0 & t (2:nblocks+1) == 1) | t (2:nblocks+1) == 0 ; +y = [(find (z)) nblocks+1] ; +r = r (y) ; +nblocks = length (y) - 1 ; +is_triangular = y (2:nblocks+1) - y (1:nblocks) > 1 ; + +%------------------------------------------------------------------------------- +% solve Ax=b, but check for small and/or triangular systems +%------------------------------------------------------------------------------- + +function x = solver (A, b, is_triangular, Control) +if (is_triangular) + % back substitution only + x = A \ b ; +elseif (size (A,1) < 4) + % a very small matrix, solve it as a dense linear system + x = full (A) \ b ; +else + % solve it as a sparse linear system + x = umfpack_solve (A, '\\', b, Control) ; +end diff -uNr UMFPACKv4.4.orig/UMFPACK/OCTAVE/umfpack.cc UMFPACK/UMFPACK/OCTAVE/umfpack.cc --- UMFPACKv4.4.orig/UMFPACK/OCTAVE/umfpack.cc 1970-01-01 01:00:00.000000000 +0100 +++ UMFPACK/UMFPACK/OCTAVE/umfpack.cc 2005-02-01 21:58:15.885225363 +0100 @@ -0,0 +1,1399 @@ +/* + +Copyright (C) 2004 David Bateman + +This program is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +This program is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with this program; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +In addition to the terms of the GPL, you are permitted to link +this program with any Open Source program, as defined by the +Open Source Initiative (www.opensource.org) + +*/ + +/* + +This is the Octave interface to the UMFPACK code, which bore the following +copyright + + UMFPACK Version 4.3 (Jan. 16, 2004), Copyright (c) 2004 by Timothy A. + Davis. All Rights Reserved. See ../README for License. + email: davis@cise.ufl.edu CISE Department, Univ. of Florida. + web: http://www.cise.ufl.edu/research/sparse/umfpack + +*/ + +#include <cfloat> +#include <cstdlib> +#include <string> + +#include <octave/config.h> +#include <octave/ov.h> +#include <octave/defun-dld.h> +#include <octave/pager.h> +#include <octave/ov-re-mat.h> + +#include "ov-re-sparse.h" +#include "ov-cx-sparse.h" + +// External UMFPACK functions in C +extern "C" { +#include "umfpack.h" +} + +#define MIN(a,b) (((a) < (b)) ? (a) : (b)) +#define MAX(a,b) (((a) > (b)) ? (a) : (b)) + +// Return an error message +static +void umfpack_error (const char *s, int A_is_complex, int nargout, + octave_value_list retval, NDArray Control, + NDArray Info, int status, int do_info) +{ + if (A_is_complex) + { + umfpack_zi_report_status (Control.fortran_vec (), status) ; + umfpack_zi_report_info (Control.fortran_vec (), Info.fortran_vec ()) ; + } + else + { + umfpack_di_report_status (Control.fortran_vec (), status) ; + umfpack_di_report_info (Control.fortran_vec (), Info.fortran_vec ()) ; + } + if (do_info > 0) + // return Info + retval (do_info) = octave_value (Info); + + error (s); +} + +DEFUN_DLD (umfpack, args, nargout, + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {[@var{x}, @var{Info}] =} umfpack (@var{a}, '\\', @var{b})\n\ +@deftypefnx {Loadable Function} {[@var{x}, @var{Info}] =} umfpack (@var{a}, '\\', @var{b}, @var{Control})\n\ +@deftypefnx {Loadable Function} {[@var{x}, @var{Info}] =} umfpack (@var{a}, @var{Qinit}, '\\', @var{b}, @var{Control})\n\ +@deftypefnx {Loadable Function} {[@var{x}, @var{Info}] =} umfpack (@var{a}, @var{Qinit}, '\\', b)\n\ +@deftypefnx {Loadable Function} {[@var{x}, @var{Info}] =} umfpack (@var{b}, '/', A) ;\n\ +@deftypefnx {Loadable Function} {[@var{x}, @var{Info}] =} umfpack (@var{b}, '/', @var{a}, @var{Control}) ;\n\ +@deftypefnx {Loadable Function} {[@var{x}, @var{Info}] =} umfpack (@var{b}, '/', @var{a}, @var{Qinit}) ;\n\ +@deftypefnx {Loadable Function} {[@var{x}, @var{Info}] =} umfpack (@var{b}, '/', @var{a}, @var{Qinit}, @var{Control}) ;\n\ +@deftypefnx {Loadable Function} {[@var{l}, @var{u}, @var{p}, @var{q}, @var{r}, @var{Info}] =} umfpack (@var{a}) ;\n\ +@deftypefnx {Loadable Function} {[@var{l}, @var{u}, @var{p}, @var{q}, @var{r}, @var{Info}] =} umfpack (@var{a}, @var{Control}) ;\n\ +@deftypefnx {Loadable Function} {[@var{l}, @var{u}, @var{p}, @var{q}, @var{r}, @var{Info}] =} umfpack (@var{a}, @var{Qinit}) ;\n\ +@deftypefnx {Loadable Function} {[@var{l}, @var{u}, @var{p}, @var{q}, @var{r}, @var{Info}] =} umfpack (@var{a}, @var{Qinit}, @var{Control}) ;\n\ +@deftypefnx {Loadable Function} {[@var{l}, @var{u}, @var{p}, @var{q}] =} umfpack (@var{a}) ;\n\ +@deftypefnx {Loadable Function} {[@var{l}, @var{u}, @var{p}, @var{q}] =} umfpack (@var{a}, @var{Control}) ;\n\ +@deftypefnx {Loadable Function} {[@var{l}, @var{u}, @var{p}, @var{q}] =} umfpack (@var{a}, @var{Qinit}) ;\n\ +@deftypefnx {Loadable Function} {[@var{l}, @var{u}, @var{p}, @var{q}] =} umfpack (@var{a}, @var{Qinit}, @var{Control}) ;\n\ +@deftypefnx {Loadable Function} {[P1, Q1, Fr, Ch, Info] =} umfpack (@var{a}, 'symbolic') ;\n\ +@deftypefnx {Loadable Function} {[P1, Q1, Fr, Ch, Info] =} umfpack (@var{a}, 'symbolic', @var{Control}) ;\n\ +@deftypefnx {Loadable Function} {[P1, Q1, Fr, Ch, Info] =} umfpack (@var{a}, @var{Qinit}, 'symbolic') ;\n\ +@deftypefnx {Loadable Function} {[P1, Q1, Fr, Ch, Info] =} umfpack (@var{a}, @var{Qinit}, 'symbolic', @var{Control});\n\ +@deftypefnx {Loadable Function} {@var{Control} =} umfpack ;\n\ +\n\ +UMFPACK v4.3 is a Octave oct-file for solving sparse linear systems.\n\ +\n\ +@iftex\n\ +@tex\n\ +\\vskip 2ex\n\ +\\hfil\\vbox{\n\ +\\offinterlineskip\n\ +\\tabskip=0pt\n\ +\\halign{\n\ +\\vrule height1.75ex depth1.25ex width 0.6pt #\\tabskip=1em &\n\ +\\hfil #\\hfil &\\vrule # & \n\ +\\hfil #\\hfil &\\vrule # width 0.6pt \\tabskip=0pt\\cr\n\ +\\noalign{\\hrule height 0.6pt}\n\ +& UMFPACK v4.3 && OCTAVE approximate equivalent &\\cr\n\ +\\noalign{\\hrule} \n\ +& x = umfpack (A, '\\', b) ; && x = A \\ b &\\cr\n\ +& && &\\cr\n\ +&x = umfpack (b, '/', A) ; && x = b / A &\\cr\n\ +& && &\\cr\n\ +&[L,U,P,Q] = umfpack (A) ; && [m,n] = size (A) ; &\\cr\n\ +& && I = speye (n) ; &\\cr\n\ +& && Q = I (:, colamd (A)) ; &\\cr\n\ +& && [L,U,P] = lu (A*Q) ; &\\cr\n\ +& && &\\cr\n\ +&[L,U,P,Q,R] = umfpack (A) ; && [m,n] = size (A) ; &\\cr\n\ +& && I = speye (n) ; &\\cr\n\ +& && Q = I (:, colamd (A)) ; &\\cr\n\ +& && r = full (sum (abs (A), 2)) ; &\\cr\n\ +& && r (find (r == 0)) = 1 ; &\\cr\n\ +& && R = spdiags (r, 0, m, m) ; &\\cr\n\ +& && [L,U,P] = lu ((R\\A)*Q) ; &\\cr\n\ +& && &\\cr\n\ +&[P,Q,F,C] = umfpack (A, 'symbolic')&& [m,n] = size (A) ; &\\cr\n\ +& && I = speye (n) ; &\\cr\n\ +& && Q = I (:, colamd (A)) ; &\\cr\n\ +& && [count,h,parent,post] = ... &\\cr\n\ +& && symbfact (A*Q, 'col') ; &\\cr\n\ +\\noalign{\\hrule height 0.6pt}\n\ +}}\\hfil\n\ +\\vskip 1ex\n\ +@end tex\n\ +@end iftex\n\ +@ifinfo\n\ +@multitable @columnfractions 0.43 .02 .43\n\ +@item UMFPACK v4.3: @tab | \n\ +@tab OCTAVE approx. equivalent\n\ +@item ------------------------------- @tab | \n\ +@tab --------------------------------\n\ +@item x = umfpack (A, '\\', b) ; @tab | \n\ +@tab x = A \\ b\n\ +@item @tab | \n\ +@tab\n\ +@item x = umfpack (b, '/', A) ; @tab | \n\ +@tab x = b / A\n\ +@item @tab | \n\ +@tab\n\ +@item [L,U,P,Q] = umfpack (A) ; @tab | \n\ +@tab [m,n] = size (A) ;\n\ +@item @tab | \n\ +@tab I = speye (n) ;\n\ +@item @tab | \n\ +@tab Q = I (:, colamd (A)) ;\n\ +@item @tab | \n\ +@tab [L,U,P] = lu (A*Q) ;\n\ +@item @tab | \n\ +@tab\n\ +@item [L,U,P,Q,R] = umfpack (A) ; @tab | \n\ +@tab [m,n] = size (A) ;\n\ +@item @tab | \n\ +@tab I = speye (n) ;\n\ +@item @tab | \n\ +@tab Q = I (:, colamd (A)) ;\n\ +@item @tab | \n\ +@tab r = full (sum (abs (A), 2)) ;\n\ +@item @tab | \n\ +@tab r (find (r == 0)) = 1 ;\n\ +@item @tab | \n\ +@tab R = spdiags (r, 0, m, m) ;\n\ +@item @tab | \n\ +@tab [L,U,P] = lu ((R\\A)*Q) ;\n\ +@item @tab | \n\ +@tab\n\ +@item [P,Q,F,C] = umfpack (A, 'symbolic') @tab | \n\ +@tab [m,n] = size (A) ; \n\ +@item @tab | \n\ +@tab I = speye (n) ;\n\ +@item @tab | \n\ +@tab Q = I (:, colamd (A)) ;\n\ +@item @tab | \n\ +@tab [count,h,parent,post] = ...\n\ +@item @tab | \n\ +@tab symbfact (A*Q, 'col') ;\n\ +@end multitable\n\ +@end ifinfo\n\ +\n\ +A must be sparse. It can be complex, singular, and/or rectangular.\n\ +A must be square for '/' or '\\'. b must be a full real or complex\n\ +vector. For @code{[@var{l}, @var{u}, @var{p}, @var{q}, @var{r}] =\n\ +umfpack (@var{a})}, the factorization is @code{@var{l} * @var{u} =\n\ +@var{p} * (@var{r} \\ @var{a}) * @var{q}}. If @var{a} has a mostly\n\ +symmetric nonzero pattern, then replace @dfn{colamd} with @dfn{amd}\n\ +in the OCTAVE-equivalent column in the table above.\n\ +\n\ +Factor or solve a sparse linear system, returning either the solution\n\ +@var{x} to @code{@var{A} * @var{x} = @var{b}} or @code{@var{A}' * @var{x}'\n\ += @var{b}'}, the factorization LU=PAQ, or LU=P(R\\A)Q. A must be sparse.\n\ +For the solve, A must be square and b must be a dense n-by-1 vector. For LU\n\ +factorization, A can be rectangular. In both cases, A and/or b can be real\n\ +or complex.\n\ +\n\ +UMFPACK analyzes the matrix and selects one of three strategies to factorize\n\ +the matrix. It first finds a set of k initial pivot entries of zero\n\ +Markowitz cost. This forms the first k rows and columns of L and U. The\n\ +remaining submatrix S is then analyzed, based on the symmetry of the nonzero\n\ +pattern of the submatrix and the values on the diagaonal. The strategies\n\ +include:\n\ +\n\ +@table @asis\n\ +@item unsymmetric\n\ +Use a COLAMD pre-ordering, a column elimination tree\n\ +post-ordering, refine the column ordering during factorization,\n\ +and make no effort at selecting pivots on the diagonal.\n\ +@item 2-by-2\n\ +Like the symmetric strategy (see below), except that local\n\ +row permutations are first made to attempt to place large entries\n\ +on the diagonal.\n\ +@item symmetric\n\ +Use an AMD pre-ordering on the matrix @code{@var{s} + @var{s}'}, an\n\ +elimination tree post-ordering, do not refine the column ordering during\n\ +factorization, and attempt to select pivots on the diagonal.\n\ +@end table\n\ +\n\ +Each of the following uses of umfpack (except for 'Control = umfpack') is\n\ +stand-alone. That is, no call to umfpack is required for any subsequent\n\ +call. In each usage, the Info output argument is optional.\n\ +\n\ +Usage:\n\ +\n\ +[x, Info] = umfpack (A, '\\', b) ;\n\ +[x, Info] = umfpack (A, '\\', b, Control) ;\n\ +[x, Info] = umfpack (A, Qinit, '\\', b, Control) ;\n\ +[x, Info] = umfpack (A, Qinit, '\\', b) ;\n\ +\n\ + Solves Ax=b (similar to x = A\\b in OCTAVE).\n\ +\n\ +[x, Info] = umfpack (b, '/', A) ;\n\ +[x, Info] = umfpack (b, '/', A, Control) ;\n\ +[x, Info] = umfpack (b, '/', A, Qinit) ;\n\ +[x, Info] = umfpack (b, '/', A, Qinit, Control) ;\n\ +\n\ + Solves A'x'=b' (similar to x = b/A in OCTAVE).\n\ +\n\ +[L, U, P, Q, R, Info] = umfpack (A) ;\n\ +[L, U, P, Q, R, Info] = umfpack (A, Control) ;\n\ +[L, U, P, Q, R, Info] = umfpack (A, Qinit) ;\n\ +[L, U, P, Q, R, Info] = umfpack (A, Qinit, Control) ;\n\ +\n\ + Returns the LU factorization of A. P and Q are returned as permutation\n\ + matrices. R is a diagonal sparse matrix of scale factors for the rows\n\ + of A, L is lower triangular, and U is upper triangular. The\n\ + factorization is L*U = P*(R\\A)*Q. You can turn off scaling by setting\n\ + Control (17) to zero (in which case R = speye (m)), or by using the\n\ + following syntaxes (in which case Control (17) is ignored):\n\ +\n\ +[L, U, P, Q] = umfpack (A) ;\n\ +[L, U, P, Q] = umfpack (A, Control) ;\n\ +[L, U, P, Q] = umfpack (A, Qinit) ;\n\ +[L, U, P, Q] = umfpack (A, Qinit, Control) ;\n\ +\n\ + Same as above, except that no row scaling is performed. The Info array\n\ + is not returned, either.\n\ +\n\ +[P1, Q1, Fr, Ch, Info] = umfpack (A, 'symbolic') ;\n\ +[P1, Q1, Fr, Ch, Info] = umfpack (A, 'symbolic', Control) ;\n\ +[P1, Q1, Fr, Ch, Info] = umfpack (A, Qinit, 'symbolic') ;\n\ +[P1, Q1, Fr, Ch, Info] = umfpack (A, Qinit, 'symbolic', Control);\n\ +\n\ + Performs only the fill-reducing column pre-ordering (including the\n\ + elimination tree post-ordering) and symbolic factorization. Q1 is the\n\ + initial column permutation (either from colamd, amd, or the input\n\ + ordering Qinit), possibly followed by a column elimination tree post-\n\ + ordering or a symmetric elimination tree post-ordering, depending on\n\ + the strategy used.\n\ +\n\ + For the unsymmetric strategy, P1 is the row ordering induced by Q1\n\ + (row-merge order). For the 2-by-2 strategy, P1 is the row ordering that\n\ + places large entries on the diagonal of P1*A*Q1. For the symmetric\n\ + strategy, P1 = Q1.\n\ +\n\ + Fr is a (nfr+1)-by-4 array containing information about each frontal\n\ + matrix, where nfr <= n is the number of frontal matrices. Fr (:,1) is\n\ + the number of pivot columns in each front, and Fr (:,2) is the parent\n\ + of each front in the supercolumn elimination tree. Fr (k,2) is zero if\n\ + k is a root. The first Fr (1,1) columns of P1*A*Q1 are the pivot\n\ + columns for the first front, the next Fr (2,1) columns of P1*A*Q1\n\ + are the pivot columns for the second front, and so on.\n\ +\n\ + For the unsymmetric strategy, Fr (:,3) is the row index of the first\n\ + row in P1*A*Q1 whose leftmost nonzero entry is in a pivot column for\n\ + the kth front. Fr (:,4) is the leftmost descendent of the kth front.\n\ + Rows in the range Fr (Fr (k,4),3) to Fr (k+1,3)-1 form the entire set\n\ + of candidate pivot rows for the kth front (some of these will typically\n\ + have been selected as pivot rows of fronts Fr (k,3) to k-1, before the\n\ + factorization reaches the kth front. If front k is a leaf node, then\n\ + Fr (k,4) is k.\n\ +\n\ + Ch is a (nchains+1)-by-3 array containing information about each\n\ + 'chain' (unifrontal sequence) of frontal matrices, and where\n\ + nchains <= nfr is the number of chains. The ith chain consists of\n\ + frontal matrices. Chain (i,1) to Chain (i+1,1)-1, and the largest\n\ + front in chain i is Chain (i,2)-by-Chain (i,3).\n\ +\n\ + This use of umfpack is not required to factor or solve a linear system\n\ + in OCTAVE. It analyzes the matrix A and provides information only.\n\ + The OCTAVE statement @code{treeplot (Fr (:,2)')} plots the column\n\ + elimination tree.\n\ +\n\ +Control = umfpack ;\n\ +\n\ + Returns a 20-by-1 vector of default parameter settings for umfpack.\n\ +\n\ +umfpack_report (Control, Info) ;\n\ +\n\ + Prints the current Control settings, and Info\n\ +\n\ +If present, Qinit is a user-supplied 1-by-n permutation vector. It is an\n\ +initial fill-reducing column pre-ordering for A; if not present, then colamd\n\ +or amd are used instead. If present, Control is a user-supplied 20-by-1\n\ +array. Control and Info are optional; if Control is not present, defaults\n\ +are used. If a Control entry is NaN, then the default is used for that entry.\n\ +\n\ +UMFPACK Version 4.3 (Jan. 16, 2004), Copyright @copyright{} 2004 by\n\ +Timothy A. Davis. All Rights Reserved.\n\ +\n\ +UMFPACK License:\n\ +\n\ +@example\n\ +Your use or distribution of UMFPACK or any modified version of\n\ +UMFPACK implies that you agree to this License.\n\ +\n\ +THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY\n\ +EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK.\n\ +\n\ +Permission is hereby granted to use or copy this program, provided\n\ +that the Copyright, this License, and the Availability of the original\n\ +version is retained on all copies. User documentation of any code that\n\ +uses UMFPACK or any modified version of UMFPACK code must cite the\n\ +Copyright, this License, the Availability note, and 'Used by permission.'\n\ +Permission to modify the code and to distribute modified code is granted,\n\ +provided the Copyright, this License, and the Availability note are\n\ +retained, and a notice that the code was modified is included. This\n\ +software was developed with support from the National Science Foundation,\n\ +and is provided to you free of charge.\n\ +@end example\n\ +\n\ +Availability: http://www.cise.ufl.edu/research/sparse/umfpack\n\ +@end deftypefn\n\ +@seealso{lu_normtest, colamd, amd, umfpack_solve}") +{ + octave_value_list retval; + int nargin = args.length (); + int op = 0; + std::string operation; + bool do_solve = false; + int do_info = 0; + + ColumnVector User_Qinit; + SparseComplexMatrix CAmatrix; + ComplexMatrix CBmatrix; + SparseMatrix Amatrix; + Matrix Bmatrix; + NDArray User_Control_matrix; + + bool A_is_complex = false; + bool B_is_complex = false; + bool X_is_complex = false; + bool transpose = false; + bool have_User_Qinit = false; + bool have_User_Control_matrix = false; + bool do_numeric = true; + bool no_scale = false; + + // find the operator + for (int i = 0 ; i < nargin ; i++) + { + if (args(i).is_string ()) + { + op = i; + break; + } + } + + if (op > 0) + { + std::string op_type = args (op).string_value (); + + if (op_type == "\\") + { + + // matrix left divide, x = A\b + + // [x, Info] = umfpack (A, '\', b) ; + // [x, Info] = umfpack (A, '\', b, Control) ; + // [x, Info] = umfpack (A, Qinit, '\', b, Control) ; + // [x, Info] = umfpack (A, Qinit, '\', b) ; + + do_solve = true; + operation = "x = A\\b"; + + if (args(0).class_name () != "sparse") + { + error ("umfpack: input matrix A must be sparse"); + return retval; + } + + if (args(0).is_complex_type ()) + { + CAmatrix = args(0).sparse_complex_matrix_value (); + A_is_complex = true; + } + else + Amatrix = args(0).sparse_matrix_value (); + + + if (args(op+1).is_complex_type ()) + { + CBmatrix = args(op+1).complex_matrix_value (); + B_is_complex = true; + } + else + Bmatrix = args(op+1).matrix_value (); + + if (nargout == 2) + do_info = 1; + + if (op == 2) + { + User_Qinit = args(1).column_vector_value (); + have_User_Qinit = true; + } + + if ((op == 1 && nargin == 4) || (op == 2 && nargin == 5)) + { + User_Control_matrix = args(nargin-1).array_value (); + have_User_Control_matrix = true; + } + + if (error_state) + { + error ("umfpack: incorrect argument type"); + return retval; + } + + if (nargin < 3 || nargin > 5 || nargout > 2) + { + error ("umfpack: wrong number of arguments"); + return retval; + } + + } + else if (op_type == "/") + { + // matrix right divide, x = b/A + + // [x, Info] = umfpack (b, '/', A) ; + // [x, Info] = umfpack (b, '/', A, Control) ; + // [x, Info] = umfpack (b, '/', A, Qinit) ; + // [x, Info] = umfpack (b, '/', A, Qinit, Control) ; + + do_solve = true; + operation = "x = b/A" ; + + transpose = true; + + if (args(2).class_name () != "sparse") + { + error ("umfpack: input matrix A must be sparse"); + return retval; + } + + if (args(2).is_complex_type ()) + { + CAmatrix = args(2).sparse_complex_matrix_value (); + A_is_complex = true; + } + else + Amatrix = args(2).sparse_matrix_value (); + + if (args(0).is_complex_type ()) + { + CBmatrix = args(0).complex_matrix_value (); + B_is_complex = true; + } + else + Bmatrix = args(0).matrix_value (); + + if (nargout == 2) + do_info = 1; + + if (nargin == 5) + { + User_Qinit = args(3).column_vector_value (); + User_Control_matrix = args(4).array_value (); + have_User_Qinit = true; + have_User_Control_matrix = true; + } + else if (nargin == 4) + { + User_Control_matrix = args(3).array_value (); + + if (User_Control_matrix.rows () == 1) + { + User_Qinit = args(3).column_vector_value (); + have_User_Qinit = true; + } + else + have_User_Control_matrix = true; + } + else if (nargin < 3 || nargin > 5 || nargout > 2) + { + error ("umfpack: wrong number of arguments"); + return retval; + } + + if (error_state) + { + error ("umfpack: incorrect argument type"); + return retval; + } + } + else if (op_type == "symbolic") + { + // symbolic factorization only + + // [P Q Fr Ch Info] = umfpack (A, 'symbolic') ; + // [P Q Fr Ch Info] = umfpack (A, 'symbolic', Control) ; + // [P Q Fr Ch Info] = umfpack (A, Qinit, 'symbolic') ; + // [P Q Fr Ch Info] = umfpack (A, Qinit, 'symbolic', Control) ; + + operation = "symbolic factorization"; + do_numeric = false; + + if (args(0).class_name () != "sparse") + { + error ("umfpack: input matrix A must be sparse"); + return retval; + } + + if (args(0).is_complex_type ()) + { + CAmatrix = args(0).sparse_complex_matrix_value (); + A_is_complex = true; + } + else + Amatrix = args(0).sparse_matrix_value (); + + if (nargout == 5) + do_info = 4 ; + + if (op == 2) + { + User_Qinit = args(1).column_vector_value (); + have_User_Qinit = true; + } + if ((op == 1 && nargin == 3) || (op == 2 && nargin == 4)) + { + User_Control_matrix = args(nargin-1).array_value (); + have_User_Control_matrix = true; + } + + if (error_state) + { + error ("umfpack: incorrect argument type"); + return retval; + } + + if (nargin < 2 || nargin > 4 || nargout > 5 || nargout < 4) + { + error ("umfpack: wrong number of arguments") ; + return retval; + } + } + else + { + error ("operator must be '/', '\\', or 'symbolic'") ; + return retval; + } + } + else if (nargin > 0) + { + // LU factorization + + // with scaling: + // [L, U, P, Q, R, Info] = umfpack (A) ; + // [L, U, P, Q, R, Info] = umfpack (A, Qinit) ; + // + // scaling determined by Control settings: + // [L, U, P, Q, R, Info] = umfpack (A, Control) ; + // [L, U, P, Q, R, Info] = umfpack (A, Qinit, Control) ; + // + // with no scaling: + // [L, U, P, Q] = umfpack (A) ; + // [L, U, P, Q] = umfpack (A, Control) ; + // [L, U, P, Q] = umfpack (A, Qinit) ; + // [L, U, P, Q] = umfpack (A, Qinit, Control) ; + + operation = "numeric factorization" ; + + if (args(0).is_complex_type ()) + { + CAmatrix = args(0).sparse_complex_matrix_value (); + A_is_complex = true; + } + else + Amatrix = args(0).sparse_matrix_value (); + + no_scale = nargout <= 4 ; + + if (nargout == 6) + do_info = 5 ; + + if (nargin == 3) + { + User_Qinit = args(1).column_vector_value (); + User_Control_matrix = args(2).array_value (); + have_User_Qinit = true; + have_User_Control_matrix = true; + } + else if (nargin == 2) + { + User_Control_matrix = args(1).array_value (); + + if (User_Control_matrix.rows () == 1) + { + User_Qinit = args(1).column_vector_value (); + have_User_Qinit = true; + } + else + have_User_Control_matrix = true; + } + else if (nargin > 3 || nargout > 6 || nargout < 4) + { + error ("umfpack: wrong number of arguments") ; + return retval; + } + } + else + { + // return default control settings + + // Control = umfpack ; + // umfpack ; + + if (nargout > 1) + { + error ("umfpack: wrong number of arguments") ; + return retval; + } + + NDArray user_control (dim_vector (UMFPACK_CONTROL, 1)); + double *user_control_ptr = user_control.fortran_vec (); + umfpack_di_defaults (user_control_ptr); + retval(0) = user_control; + return retval; + } + + // check inputs + + int n_row = Amatrix.rows (); + int n_col = Amatrix.cols (); + int nn = MAX (n_row, n_col) ; + int n_inner = MIN (n_row, n_col) ; + if (do_solve && n_row != n_col) + { + error ("umfpack: input matrix A must square for '\\' or '/'") ; + return retval; + } + if (n_row == 0 || n_col == 0) + { + error ("umfpack: input matrix A cannot have zero rows or zero columns") ; + return retval; + } + + /* The real/complex status of A determines which version to use, */ + /* (umfpack_di_* or umfpack_zi_*). */ + const int *Ap; + const int *Ai; + const double *Ax; + const double *Bx; + + if (A_is_complex) + { + Ap = CAmatrix.cidx (); + Ai = CAmatrix.ridx (); + Ax = X_CAST (const double *, CAmatrix.data ()); + } + else + { + Ap = Amatrix.cidx (); + Ai = Amatrix.ridx (); + Ax = Amatrix.data (); + } + + if (B_is_complex) + Bx = X_CAST (const double *, CBmatrix.fortran_vec ()); + else + Bx = Bmatrix.fortran_vec (); + + if (do_solve) + { + int b_row = Bmatrix.rows (); + int b_col = Bmatrix.cols (); + + if (!((transpose && b_row == 1 && b_col == nn) || + (!transpose && b_row == nn && b_col == 1))) + { + error ("umfpack: b has the wrong dimensions") ; + return retval; + } + + X_is_complex = A_is_complex || B_is_complex ; + } + + // set the Control parameters + NDArray Control (dim_vector (UMFPACK_CONTROL, 1)); + double *Control_ptr = Control.fortran_vec (); + if (A_is_complex) + umfpack_zi_defaults (Control_ptr) ; + else + umfpack_di_defaults (Control_ptr) ; + + if (have_User_Control_matrix) + { + int size = MIN (UMFPACK_CONTROL, User_Control_matrix.length ()); + for (int i = 0 ; i < size ; i++) + Control (i) = User_Control_matrix (i) ; + } + + if (no_scale) + { + // turn off scaling for [L, U, P, Q] = umfpack (A) ; + // ignoring the input value of Control (24) for the usage + // [L, U, P, Q] = umfpack (A, Control) ; + Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE ; + } + + int print_level; + if (xisnan (Control (UMFPACK_PRL))) + print_level = UMFPACK_DEFAULT_PRL ; + else + print_level = int (Control (UMFPACK_PRL)) ; + + Control (UMFPACK_PRL) = print_level ; + + // check Qinit, if present + int *Qinit = NULL; + if (have_User_Qinit) + { + if(User_Qinit.rows () != 1 || User_Qinit.cols () != n_col) + { + error ("umfpack: Qinit must be 1-by-n_col") ; + return retval; + } + + Qinit = new int [n_col]; + for (int i = 0; i < n_col; i++) + Qinit[i] = static_cast<int> (User_Qinit (i)); + } + + // report the inputs A and Qinit + + if (print_level >= 2) + // print the operation + octave_stdout << "\numfpack: " << operation; + + if (A_is_complex) + { + umfpack_zi_report_control (Control_ptr) ; + + if (print_level >= 3) + octave_stdout << "\nA: " ; + + umfpack_zi_report_matrix (n_row, n_col, Ap, Ai, Ax, NULL, + 1, Control_ptr) ; + if (have_User_Qinit) + { + if (print_level >= 3) + octave_stdout << "\nQinit: " ; + + umfpack_zi_report_perm (n_col, Qinit, Control_ptr) ; + } + } + else + { + umfpack_di_report_control (Control_ptr) ; + + if (print_level >= 3) + octave_stdout << "\nA: " ; + + umfpack_di_report_matrix (n_row, n_col, Ap, Ai, Ax, + 1, Control_ptr) ; + if (have_User_Qinit) + { + if (print_level >= 3) + octave_stdout << "\nQinit: " ; + + umfpack_di_report_perm (n_col, Qinit, Control_ptr) ; + } + } + +#ifndef NO_TRANSPOSE_FORWARD_SLASH + // create the array transpose for x = b/A + if (transpose) + if (A_is_complex) + { + CAmatrix = CAmatrix.transpose (); + Ap = Amatrix.cidx (); + Ai = Amatrix.ridx (); + Ax = X_CAST (const double *, CAmatrix.data ()); + } + else + { + Amatrix = Amatrix.transpose (); + Ap = Amatrix.cidx (); + Ai = Amatrix.ridx (); + Ax = Amatrix.data (); + } +#endif + + // perform the symbolic factorization + + NDArray InfoOut (dim_vector (1, UMFPACK_INFO)); + double * Info = InfoOut.fortran_vec (); + void *Symbolic; + int status, status2; + if (A_is_complex) + status = umfpack_zi_qsymbolic (n_row, n_col, Ap, Ai, Ax, NULL, + Qinit, &Symbolic, Control_ptr, + Info); + else + status = umfpack_di_qsymbolic (n_row, n_col, Ap, Ai, Ax, + Qinit, &Symbolic, Control_ptr, + Info); + + if (status < 0) + { + umfpack_error ("symbolic factorization failed", A_is_complex, + nargout, retval, Control, InfoOut, status, do_info) ; + return retval; + } + + if (have_User_Qinit) + delete [] Qinit; + + // report the Symbolic object + + if (A_is_complex) + umfpack_zi_report_symbolic (Symbolic, Control_ptr) ; + else + umfpack_di_report_symbolic (Symbolic, Control_ptr) ; + + // perform numeric factorization, or just return symbolic factorization + + if (do_numeric) + { + // perform the numeric factorization + void *Numeric; + + if (A_is_complex) + status = umfpack_zi_numeric (Ap, Ai, Ax, NULL, Symbolic, &Numeric, + Control_ptr, Info) ; + else + status = umfpack_di_numeric (Ap, Ai, Ax, Symbolic, &Numeric, + Control_ptr, Info) ; + + // free the symbolic factorization + if (A_is_complex) + umfpack_zi_free_symbolic (&Symbolic) ; + else + umfpack_di_free_symbolic (&Symbolic) ; + + // report the Numeric object + if (status < 0) + { + umfpack_error ("numeric factorization failed", A_is_complex, + nargout, retval, Control, InfoOut, status, do_info); + return retval; + } + + if (A_is_complex) + (void) umfpack_zi_report_numeric (Numeric, Control_ptr) ; + else + (void) umfpack_di_report_numeric (Numeric, Control_ptr) ; + + // return the solution or the factorization + + if (do_solve) + { + int sys; + ComplexNDArray Xcmplx; + NDArray Xreal; + + // solve Ax=b or A'x'=b', and return just the solution x + +#ifndef NO_TRANSPOSE_FORWARD_SLASH + if (transpose) + { + // A.'x.'=b.' gives the same x=b/A as solving A'x'=b' + // since C=A.' was factorized, solve with sys = UMFPACK_A + // since x and b are vectors, x.' and b.' are implicit + if (X_is_complex) + Xcmplx.resize (dim_vector (1, nn)); + else + Xreal.resize (dim_vector (1, nn)); + } + else + { + if (X_is_complex) + Xcmplx.resize (dim_vector (nn, 1)); + else + Xreal.resize (dim_vector (nn, 1)); + } + + sys = UMFPACK_A ; +#else + if (transpose) + { + // If A is real, A'x=b is the same as A.'x=b. + // x and b are vectors, so x and b are the same as x' and b'. + // If A is complex, then A.'x.'=b.' gives the same solution x + // as the complex conjugate transpose. If we used the A'x=b + // option in umfpack_*_solve, we would have to form b' on + // input and x' on output (negating the imaginary part). + // We can save this work by just using the A.'x=b option in + // umfpack_*_solve. Then, forming x.' and b.' is implicit, + // since x and b are just vectors anyway. + // In both cases, the system to solve is A.'x=b + if (X_is_complex) + Xcmplx.resize (dim_vector (1, nn)); + else + Xreal.resize (dim_vector (1, nn)); + + sys = UMFPACK_Aat ; + } + else + { + if (X_is_complex) + Xcmplx.resize (dim_vector (nn, 1)); + else + Xreal.resize (dim_vector (nn, 1)); + sys = UMFPACK_A ; + } +#endif + + // print the right-hand-side, B + if (print_level >= 3) + octave_stdout << "\nright-hand side, b: "; + + if (B_is_complex) + (void) umfpack_zi_report_vector (nn, Bx, NULL, Control_ptr) ; + else + (void) umfpack_di_report_vector (nn, Bx, Control_ptr) ; + + // solve the system + double * Xx; + if (X_is_complex) + Xx = X_CAST (double *, Xcmplx.fortran_vec ()); + else + Xx = Xreal.fortran_vec (); + status2 = UMFPACK_OK ; + + if (A_is_complex) + { + if (!B_is_complex) + { + OCTAVE_LOCAL_BUFFER (double, Bz, nn); + for (int i = 0; i < nn; i++) + Bz[i] = 0.; + + status = umfpack_zi_solve (sys, Ap, Ai, Ax, NULL, Xx, NULL, + Bx, Bz, Numeric, Control_ptr, + Info); + } + else + status = umfpack_zi_solve (sys, Ap, Ai, Ax, NULL, Xx, NULL, + Bx, NULL, Numeric, Control_ptr, + Info); + } + else + { + if (B_is_complex) + { + // Ax=b when b is complex and A is sparse can be split + // into two systems, A*xr=br and A*xi=bi, where r denotes + // the real part and i the imaginary part of x and b. + OCTAVE_LOCAL_BUFFER (double, Tx, nn); + OCTAVE_LOCAL_BUFFER (double, Tz, nn); + + status = umfpack_di_solve (sys, Ap, Ai, Ax, Tx, Bx, + Numeric, Control_ptr, Info); + status2 = umfpack_di_solve (sys, Ap, Ai, Ax, Tz, Bx, + Numeric, Control_ptr, Info) ; + + for (int i = 0; i < nn; i++) + Xcmplx (i) = Complex (Tx[i], Tz[i]); + } + else + status = umfpack_di_solve (sys, Ap, Ai, Ax, Xx, Bx, + Numeric, Control_ptr, Info); + } + + // free the Numeric object + if (A_is_complex) + umfpack_zi_free_numeric (&Numeric) ; + else + umfpack_di_free_numeric (&Numeric) ; + + // check error status + if (status < 0 || status2 < 0) + { + umfpack_error ("solve failed", A_is_complex, nargout, + retval, Control, InfoOut, status, do_info) ; + return retval; + } + + // print the solution, X + if (print_level >= 3) + octave_stdout << "\nsolution, x: "; + + if (X_is_complex) + (void) umfpack_zi_report_vector (nn, Xx, NULL, Control_ptr); + else + (void) umfpack_di_report_vector (nn, Xx, Control_ptr); + + // warn about singular or near-singular matrices + // no warning is given if Control (1) is zero + if (Control (UMFPACK_PRL) >= 1) + { + if (status == UMFPACK_WARNING_singular_matrix) + { + warning ("matrix is singular"); + warning ("Try increasing Control (%d) and Control (%d).", + 1+UMFPACK_PIVOT_TOLERANCE, + 1+UMFPACK_SYM_PIVOT_TOLERANCE); + warning ("(Suppress this warning with Control (%d) = 0.)", + 1+UMFPACK_PRL); + } + else if (InfoOut (UMFPACK_RCOND) < DBL_EPSILON) + { + warning ("matrix is nearly singular, rcond = %g", + InfoOut (UMFPACK_RCOND)); + warning ("Try increasing Control (%d) and Control (%d).", + 1+UMFPACK_PIVOT_TOLERANCE, + 1+UMFPACK_SYM_PIVOT_TOLERANCE); + warning ("(Suppress this warning with Control (%d) = 0.)", + 1+UMFPACK_PRL); + } + } + + // Setup the return value + if (X_is_complex) + retval (0) = octave_value (Xcmplx); + else + retval (0) = octave_value (Xreal); + } + else + { + // get L, U, P, Q, and r + int lnz, unz, ignore1, ignore2, ignore3; + + if (A_is_complex) + status = umfpack_zi_get_lunz (&lnz, &unz, &ignore1, &ignore2, + &ignore3, Numeric) ; + else + status = umfpack_di_get_lunz (&lnz, &unz, &ignore1, &ignore2, + &ignore3, Numeric) ; + + if (status < 0) + { + if (A_is_complex) + umfpack_zi_free_numeric (&Numeric) ; + else + umfpack_di_free_numeric (&Numeric) ; + + umfpack_error ("extracting LU factors failed", A_is_complex, + nargout, retval, Control, InfoOut, status, + do_info); + return retval; + } + + // avoid malloc of zero-sized arrays + lnz = MAX (lnz, 1) ; + unz = MAX (unz, 1) ; + + // get space for the *** ROW *** form of L + SparseMatrix Lreal; + SparseComplexMatrix Limag; + int *Ltp, *Ltj; + double *Ltx; + if (A_is_complex) + { + Limag = SparseComplexMatrix (n_inner, n_row, lnz); + Ltp = Limag.cidx (); + Ltj = Limag.ridx (); + Ltx = X_CAST (double *, Limag.data ()); + } + else + { + Lreal = SparseMatrix (n_inner, n_row, lnz); + Ltp = Lreal.cidx (); + Ltj = Lreal.ridx (); + Ltx = Lreal.data (); + } + + // create permanent copy of the output matrix U + int *Up, *Ui; + double *Ux; + SparseMatrix Ureal; + SparseComplexMatrix Uimag; + + if (A_is_complex) + { + Uimag = SparseComplexMatrix (n_inner, n_col, unz); + Up = Uimag.cidx (); + Ui = Uimag.ridx (); + Ux = X_CAST (double *, Uimag.data ()); + } + else + { + Ureal = SparseMatrix (n_inner, n_col, unz); + Up = Ureal.cidx (); + Ui = Ureal.ridx (); + Ux = Ureal.data (); + } + + // temporary space for the integer permutation vectors + OCTAVE_LOCAL_BUFFER (int, P, n_row); + OCTAVE_LOCAL_BUFFER (int, Q, n_col); + + // get scale factors, if requested + status2 = UMFPACK_OK ; + SparseMatrix Rsout; + double * Rs; + if (!no_scale) + { + // create a diagonal sparse matrix for the scale factors + Rsout = SparseMatrix (n_row, n_row, n_row); + for (int i = 0 ; i < n_row ; i++) + { + Rsout.cidx (i) = i; + Rsout.ridx (i) = i; + } + Rsout.cidx (n_row) = n_row; + Rs = Rsout.data (); + } + else + Rs = (double *) NULL ; + + // get Lt, U, P, Q, and Rs from the numeric object + int do_recip; + if (A_is_complex) + { + status = umfpack_zi_get_numeric (Ltp, Ltj, Ltx, NULL, Up, Ui, + Ux, NULL, P, Q, + (double *) NULL, + (double *) NULL, + &do_recip, Rs, Numeric) ; + umfpack_zi_free_numeric (&Numeric) ; + } + else + { + status = umfpack_di_get_numeric (Ltp, Ltj, Ltx, Up, Ui, + Ux, P, Q, (double *) NULL, + &do_recip, Rs, Numeric) ; + umfpack_di_free_numeric (&Numeric) ; + } + + if (!no_scale) + retval (4) = octave_vale (Rsout); + + // for the oct-file, -DNRECIPROCAL must be set, + // so do_recip must be FALSE + + if (status < 0 || status2 < 0 || do_recip) + { + umfpack_error ("extracting LU factors failed", A_is_complex, + nargout, retval, Control, InfoOut, status, + do_info); + return retval; + } + + if (A_is_complex) + retval (1) = octave_value (Uimag); + else + retval (1) = octave_valye (Ureal); + + // create sparse permutation matrix for P + SparseMatrix Pout (n_row, n_row, n_row); + for (int k = 0 ; k < n_row ; k++) + { + Pout.cidx (k) = k ; + Pout.ridx (P [k]) = k; + Pout.data (k) = 1; + } + Pout.cidx (n_row) = n_row; + retval (2) = octave_value (Pout); + + // create sparse permutation matrix for Q + SparseMatrix Qout (n_col, n_col, n_col); + for (int k = 0 ; k < n_col ; k++) + { + Qout.cidx (k) = k ; + Qout.ridx (k) = Q[k]; + Qout.data (k) = 1; + } + Qout.cidx (n_col) = n_col; + retval (3) = octave_value (Qout); + + // permanent copy of L + if (A_is_complex) + retval (0) = octave_value (Limag.transpose()); + else + retval (0) = octave_value (Lreal.transpose()); + + if (status < 0) + { + umfpack_error ("constructing L failed", A_is_complex, + nargout, retval, Control, InfoOut, status, + do_info) ; + return octave_value (); + } + + // print L, U, P, and Q + if (A_is_complex) + { + if (print_level >= 3) + { + octave_stdout << "\nL: "; + int *Lp = Limag.cidx (); + int *Li = Limag.ridx (); + double *Lx = X_CAST (double *, Limag.data ()); + + (void) umfpack_zi_report_matrix (n_row, n_inner, Lp, Li, + Lx, NULL, 1, Control_ptr) ; + } + + if (print_level >= 3) + octave_stdout << "\nU: "; + (void) umfpack_zi_report_matrix (n_inner, n_col, Up, Ui, + Ux, NULL, 1, Control_ptr) ; + if (print_level >= 3) + octave_stdout << "\nP: "; + (void) umfpack_zi_report_perm (n_row, P, Control_ptr); + if (print_level >= 3) + octave_stdout << "\nQ: "; + (void) umfpack_zi_report_perm (n_col, Q, Control_ptr); + } + else + { + if (print_level >= 3) + { + int *Lp = Lreal.cidx (); + int *Li = Lreal.ridx (); + double *Lx = Lreal.data (); + octave_stdout << "\nL: "; + (void) umfpack_di_report_matrix (n_row, n_inner, Lp, Li, + Lx, 1, Control_ptr); + } + + if (print_level >= 3) + octave_stdout << "\nU: "; + (void) umfpack_di_report_matrix (n_inner, n_col, Up, Ui, + Ux, 1, Control_ptr); + if (print_level >= 3) + octave_stdout << "\nP: "; + (void) umfpack_di_report_perm (n_row, P, Control_ptr); + if (print_level >= 3) + octave_stdout << "\nQ: "; + (void) umfpack_di_report_perm (n_col, Q, Control_ptr); + } + } + } + else + { + // return the symbolic factorization + int ignore1, ignore2, ignore3; + OCTAVE_LOCAL_BUFFER (int, Q, n_col); + OCTAVE_LOCAL_BUFFER (int, P, n_row); + OCTAVE_LOCAL_BUFFER (int, Front_npivcol, nn + 1); + OCTAVE_LOCAL_BUFFER (int, Front_parent, nn + 1); + OCTAVE_LOCAL_BUFFER (int, Front_1strow, nn + 1); + OCTAVE_LOCAL_BUFFER (int, Front_leftmostdesc, nn + 1); + OCTAVE_LOCAL_BUFFER (int, Chain_start, nn + 1); + OCTAVE_LOCAL_BUFFER (int, Chain_maxrows, nn + 1); + OCTAVE_LOCAL_BUFFER (int, Chain_maxcols, nn + 1); + + int nz, nfronts, nchains; + + if (A_is_complex) + { + status = umfpack_zi_get_symbolic (&ignore1, &ignore2, &ignore3, + &nz, &nfronts, &nchains, P, Q, + Front_npivcol, Front_parent, + Front_1strow, + Front_leftmostdesc, + Chain_start, Chain_maxrows, + Chain_maxcols, Symbolic) ; + umfpack_zi_free_symbolic (&Symbolic) ; + } + else + { + status = umfpack_di_get_symbolic (&ignore1, &ignore2, &ignore3, + &nz, &nfronts, &nchains, P, Q, + Front_npivcol, Front_parent, + Front_1strow, + Front_leftmostdesc, + Chain_start, Chain_maxrows, + Chain_maxcols, Symbolic) ; + umfpack_di_free_symbolic (&Symbolic) ; + } + + if (status < 0) + { + umfpack_error ("extracting symbolic factors failed", + A_is_complex, nargout, retval, Control, + InfoOut, status, do_info) ; + return retval; + } + + // create sparse permutation matrix for P + SparseMatrix Pout (n_row, n_row, n_row); + for (int k = 0 ; k < n_row ; k++) + { + Pout.cidx (k) = k; + Pout.ridx (P [k]) = k; + Pout.data (k) = 1; + } + Pout.cidx (n_row) = n_row; + retval (0) = octave_value (Pout); + + // create sparse permutation matrix for Q + SparseMatrix Qout (n_col, n_col, n_col); + for (int k = 0 ; k < n_col ; k++) + { + Qout.cidx (k) = k; + Qout.ridx (k) = Q[k]; + Qout.data (k) = 1; + } + Qout.cidx (n_col) = n_col; + retval (1) = octave_value (Qout); + + // create Fr + Matrix Frout (nfronts + 1, 4); + for (int i = 0 ; i <= nfronts ; i++) + { + // convert parent, 1strow, and leftmostdesc to 1-based + Frout (i, 0) = (double) (Front_npivcol [i]) ; + Frout (i, 1) = (double) (Front_parent [i] + 1) ; + Frout (i, 2) = (double) (Front_1strow [i] + 1) ; + Frout (i, 3) = (double) (Front_leftmostdesc [i] + 1) ; + } + retval (2) = octave_value (Frout); + + // create Ch + Matrix Chout (nchains + 1, 3); + for (int i = 0 ; i <= nchains ; i++) + { + // convert to 1-based + Chout (i, 0) = (double) (Chain_start [i] + 1) ; + Chout (i, 1) = (double) (Chain_maxrows [i]) ; + Chout (i, 2) = (double) (Chain_maxcols [i]) ; + } + Chout (0, nchains) = (double) Chain_start [nchains] + 1 ; + Chout (1, nchains) = 0.; + Chout (2, nchains) = 0.; + retval (3) = octave_value (Chout); + } + + // report Info + if (A_is_complex) + umfpack_zi_report_info (Control_ptr, Info); + else + umfpack_di_report_info (Control_ptr, Info); + + if (do_info > 0) + retval (do_info) = InfoOut; + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ diff -uNr UMFPACKv4.4.orig/UMFPACK/OCTAVE/umfpack_demo.m UMFPACK/UMFPACK/OCTAVE/umfpack_demo.m --- UMFPACKv4.4.orig/UMFPACK/OCTAVE/umfpack_demo.m 1970-01-01 01:00:00.000000000 +0100 +++ UMFPACK/UMFPACK/OCTAVE/umfpack_demo.m 2004-12-30 01:58:46.000000000 +0100 @@ -0,0 +1,191 @@ +function umfpack_demo +% UMFPACK DEMO +% +% A demo of UMFPACK for OCTAVE. +% +% See also umfpack, umfpack_make, umfpack_details, umfpack_report, +% and umfpack_simple. + +% UMFPACK Version 4.3 (Jan. 16, 2004), Copyright (c) 2004 by Timothy A. +% Davis. All Rights Reserved. Type umfpack_details for License. + +%------------------------------------------------------------------------------- +% get default control parameters +%------------------------------------------------------------------------------- + +control = umfpack ; +fprintf ('\nEnter the printing level for UMFPACK''s output statistics:\n') ; +fprintf ('0: none, 1: errors only, 2: statistics, 4: print some of outputs\n') ; +c = input ('5: print all output [default is 1]: ') ; +if (isempty (c)) + c = 1 ; +end +control (1) = c ; + +%------------------------------------------------------------------------------- +% solve a simple system +%------------------------------------------------------------------------------- + +fprintf ('\n--------------------------------------------------------------\n') ; +fprintf ('Factor and solve a small system, Ax=b, using default parameters\n') ; +if (control (1) > 1) + fprintf ('(except for verbose printing enabled)\n') ; +end + +load west0067 ; +A = Problem.A ; +n = size (A, 1) ; +b = rand (n, 1) ; + +fprintf ('Solving Ax=b via UMFPACK:\n') ; +[xu, info] = umfpack (A, '\\', b, control) ; +x = xu ; + +fprintf ('Solving Ax=b via OCTAVE:\n') ; +xm = A\b ; +x = xm ; + +fprintf ('Difference between UMFPACK and OCTAVE solution: %g\n', ... + norm (xu - xm, Inf)) ; + +%------------------------------------------------------------------------------- +% spy the results +%------------------------------------------------------------------------------- + +figure (1) ; +clf + +subplot (2,3,1) +title ('The matrix A') ; +spy (A) + +subplot (2,3,2) +[P1, Q1, Fr, Ch, Info] = umfpack (A, 'symbolic') ; +title ('Supernodal column elimination tree') ; +%% Disable for now !! +%% treeplot (Fr (1:end-1,2)') ; + +subplot (2,3,3) +title ('A, with initial row and column order') ; +spy (P1 * A * Q1) + +subplot (2,3,4) +fprintf ('\n--------------------------------------------------------------\n') ; +fprintf ('\nFactorizing [L, U, P, Q, R] = umfpack (A)\n') ; +[L, U, P, Q, R] = umfpack (A) ; +title ('A, with final row/column order') ; +spy (P*A*Q) + +fprintf ('\nP * (R\\A) * Q - L*U should be zero:\n') ; +fprintf ('norm (P*(R\\A)*Q - L*U, 1) = %g (exact) %g (estimated)\n', ... + norm (P * (R\A) * Q - L*U, 1), lu_normest (P * (R\A) * Q, L, U)) ; + +fprintf ('\nSolution to Ax=b via UMFPACK factorization:\n') ; +fprintf ('x = Q * (U \\ (L \\ (P * (R \\ b))))\n') ; +xu = Q * (U \ (L \ (P * (R \ b)))) ; +x = xu ; + +fprintf ('\nUMFPACK flop count: %d\n', luflop (L, U)) ; + +subplot (2,3,5) +title ('UMFPACK LU factors') ; +spy (spones (L) + spones (U)) + +subplot (2,3,6) +fprintf ('\nFactorizing [L, U, P] = lu (A (:, q))\n') ; +try + q = colamd (A) ; +catch + fprintf ('\n *** colamd not found, using colmmd instead *** \n') ; + q = colmmd (A) ; +end +[L, U, P] = lu (A (:,q)) ; +title ('OCTAVE LU factors') ; +spy (spones (L) + spones (U)) + +fprintf ('\nSolution to Ax=b via OCTAVE factorization:\n') ; +fprintf ('x = U \\ (L \\ (P * b)) ; x (q) = x ;\n') ; +xm = U \ (L \ (P * b)) ; +xm (q) = xm ; + +fprintf ('Difference between UMFPACK and OCTAVE solution: %g\n', ... + norm (xu - xm, Inf)) ; + +fprintf ('\nOCTAVE LU flop count: %d\n', luflop (L, U)) ; + +%------------------------------------------------------------------------------- +% solve A'x=b +%------------------------------------------------------------------------------- + +fprintf ('\n--------------------------------------------------------------\n') ; +fprintf ('Solve A''x=b:\n') ; + +fprintf ('Solving A''x=b via UMFPACK:\n') ; +[xu, info] = umfpack (b', '/', A, control) ; +xu = xu' ; + +fprintf ('Solving A''x=b via OCTAVE:\n') ; +xm = (b'/A)' ; +x = xm ; + +fprintf ('Difference between UMFPACK and OCTAVE solution: %g\n', ... + norm (xu - xm, Inf)) ; + +%------------------------------------------------------------------------------- +% factor A' and then solve Ax=b using the factors of A' +%------------------------------------------------------------------------------- + +fprintf ('\n--------------------------------------------------------------\n') ; +fprintf ('Compute C = A'', and compute the LU factorization of C.\n') ; +fprintf ('Factorizing A'' can sometimes be better than factorizing A itself\n'); +fprintf ('(less work and memory usage). Solve C''x=b; the solution is the\n') ; +fprintf ('same as the solution to Ax=b for the original A.\n'); + +C = A' ; + +% factorize C (P,Q) = L*U +[L, U, P, Q, R, info] = umfpack (C, control) ; + +fprintf ('\nP * (R\\C) * Q - L*U should be zero:\n') ; +fprintf ('norm (P*(R\\C)*Q - L*U, 1) = %g (exact) %g (estimated)\n', ... + norm (P * (R\C) * Q - L*U, 1), lu_normest (P * (R\C) * Q, L, U)) ; + +fprintf ('\nSolution to Ax=b via UMFPACK, using the factors of C:\n') ; +fprintf ('x = R \\ (P'' * (L'' \\ (U'' \\ (Q'' * b)))) ;\n') ; +xu = R \ (P' * (L' \ (U' \ (Q' * b)))) ; +x = xu ; + +fprintf ('Solution to Ax=b via OCTAVE:\n') ; +xm = A\b ; +x = xm ; + +fprintf ('Difference between UMFPACK and OCTAVE solution: %g\n', ... + norm (xu - xm, Inf)) ; + +%------------------------------------------------------------------------------- +% solve Ax=B +%------------------------------------------------------------------------------- + +fprintf ('\n--------------------------------------------------------------\n') ; +fprintf ('\nSolve AX=B, where B is n-by-10, and sparse\n') ; +B = sprandn (n, 10, 0.05) ; +XU = umfpack_solve (A, '\\', B, control) ; +XM = A\B ; + +fprintf ('Difference between UMFPACK and OCTAVE solution: %g\n', ... + norm (XU - XM, Inf)) ; + +fprintf ('\n--------------------------------------------------------------\n') ; +fprintf ('\nSolve AX=B, where B is n-by-10, and sparse, using umfpack_btf\n') ; +XU = umfpack_btf (A, B, control) ; + +fprintf ('Difference between UMFPACK and OCTAVE solution: %g\n', ... + norm (XU - XM, Inf)) ; + +fprintf ('\n--------------------------------------------------------------\n') ; +fprintf ('\nSolve A''X=B, where B is n-by-10, and sparse\n') ; +XU = umfpack_solve (B', '/', A, control) ; +XM = B'/A ; + +fprintf ('Difference between UMFPACK and OCTAVE solution: %g\n', ... + norm (XU - XM, Inf)) ; diff -uNr UMFPACKv4.4.orig/UMFPACK/OCTAVE/umfpack_demo.m.out UMFPACK/UMFPACK/OCTAVE/umfpack_demo.m.out --- UMFPACKv4.4.orig/UMFPACK/OCTAVE/umfpack_demo.m.out 1970-01-01 01:00:00.000000000 +0100 +++ UMFPACK/UMFPACK/OCTAVE/umfpack_demo.m.out 2004-12-30 01:58:46.000000000 +0100 @@ -0,0 +1,72 @@ +>> umfpack_demo + +Enter the printing level for UMFPACK's output statistics: +0: none, 1: errors only, 2: statistics, 4: print some of outputs +5: print all output [default is 1]: + +-------------------------------------------------------------- +Factor and solve a small system, Ax=b, using default parameters +Solving Ax=b via UMFPACK: +Solving Ax=b via MATLAB: +Difference between UMFPACK and MATLAB solution: 1.24345e-14 + +-------------------------------------------------------------- + +Factorizing [L, U, P, Q, R] = umfpack (A) + +P * (R\A) * Q - L*U should be zero: +norm (P*(R\A)*Q - L*U, 1) = 4.2068e-16 (exact) 3.74627e-16 (estimated) + +Solution to Ax=b via UMFPACK factorization: +x = Q * (U \ (L \ (P * (R \ b)))) + +UMFPACK flop count: 2362 + +Factorizing [L, U, P] = lu (A (:, q)) +If you are using a version of MATLAB prior to V6.0, then the +following statement (q = colamd (A)) may fail. Either download +colamd from http://www.cise.ufl.edu/research/sparse, upgrade to +MATLAB V6.0 or later, or replace the statement with +q = colmmd (A) ; + +Solution to Ax=b via MATLAB factorization: +x = U \ (L \ (P * b)) ; x (q) = x ; +Difference between UMFPACK and MATLAB solution: 1.37668e-14 + +MATLAB LU flop count: 3164 + +-------------------------------------------------------------- +Solve A'x=b: +Solving A'x=b via UMFPACK: +Solving A'x=b via MATLAB: +Difference between UMFPACK and MATLAB solution: 3.10862e-15 + +-------------------------------------------------------------- +Compute C = A', and compute the LU factorization of C. +Factorizing A' can sometimes be better than factorizing A itself +(less work and memory usage). Solve C'x=b; the solution is the +same as the solution to Ax=b for the original A. + +P * (R\C) * Q - L*U should be zero: +norm (P*(R\C)*Q - L*U, 1) = 1.31839e-16 (exact) 6.41848e-17 (estimated) + +Solution to Ax=b via UMFPACK, using the factors of C: +x = R \ (P' * (L' \ (U' \ (Q' * b)))) ; +Solution to Ax=b via MATLAB: +Difference between UMFPACK and MATLAB solution: 1.77636e-14 + +-------------------------------------------------------------- + +Solve AX=B, where B is n-by-10, and sparse +Difference between UMFPACK and MATLAB solution: 2.88198e-14 + +-------------------------------------------------------------- + +Solve AX=B, where B is n-by-10, and sparse, using umfpack_btf +Difference between UMFPACK and MATLAB solution: 9.79736e-14 + +-------------------------------------------------------------- + +Solve A'X=B, where B is n-by-10, and sparse +Difference between UMFPACK and MATLAB solution: 1.05244e-13 +>> diary off diff -uNr UMFPACKv4.4.orig/UMFPACK/OCTAVE/umfpack_make.m UMFPACK/UMFPACK/OCTAVE/umfpack_make.m --- UMFPACKv4.4.orig/UMFPACK/OCTAVE/umfpack_make.m 1970-01-01 01:00:00.000000000 +0100 +++ UMFPACK/UMFPACK/OCTAVE/umfpack_make.m 2004-12-30 01:58:46.000000000 +0100 @@ -0,0 +1,356 @@ +function umfpack_make +% UMFPACK_MAKE +% +% Compiles the UMFPACK mexFunction and then runs a simple demo. +% +% UMFPACK Version 4.3 (Jan. 16, 2004), Copyright (c) 2004 by Timothy A. +% Davis. All Rights Reserved. Type umfpack_details for License. +% +% See also: umfpack, umfpack_details, umfpack_report, umfpack_demo, and +% umfpack_simple. + +help umfpack_make + +fprintf ('\n--------------------------------------------------------------\n') ; +fprintf ('Now compiling the UMFPACK and AMD mexFunctions.\n') ; +fprintf ('--------------------------------------------------------------\n') ; + +try + % ispc does not appear in MATLAB 5.3 + pc = ispc ; +catch + % if ispc fails, assume we aren't on a Windows PC. + pc = 0 ; +end + +obj = 'o' ; +blas_lib = '' ; +if (pc) + obj = 'obj' ; +end + +%------------------------------------------------------------------------------- +% BLAS option +%------------------------------------------------------------------------------- + +msg = [ ... + '\nUsing the BLAS is faster, but might not compile correctly.\n', ... + 'If you get an error stating that dgemm, dgemv, dger, zgemm,\n', ... + 'zgemv, and/or zger are not defined, then recompile without the\n', ... + 'BLAS. You can ignore warnings that these routines are implicitly\n', ... + 'declared.\n\nPlease select one of the following options: \n', ... + ' 1: attempt to compile with the BLAS (default)\n', ... + ' 2: do not use the BLAS\n'] ; +fprintf (msg) ; +blas = input (': ') ; +if (isempty (blas)) + blas = 1 ; +end +if (blas == 1) + % try to link to MATLAB's built-in BLAS + blas = '' ; + if (pc) + % the default lcc compiler needs this library to access the BLAS + blas_lib = ' libmwlapack.lib' ; + msg = [ ... + '\nCheck to see if you have a file called libmwlapack.lib in the\n', ... + '<matlab>\\extern\\lib\\win32\\lcc\\ directory, where <matlab> is ', ... + 'the\ndirectory where MATLAB is installed. If a file of that ', ... + 'name is already\nthere, then you don''t have to do anything. ', ... + 'Otherwise, you must first\ncopy the libmwlapack.lib file from ', ... + 'the umfpack\\lcc_lib\\ directory to the\n', ... + '<matlab>\\extern\\lib\\win32\\lcc\\ directory. Next, type\n\n', ... + ' mex -setup\n\n', ... + 'at the MATLAB prompt, and ask MATLAB to select the lcc compiler. ',... + 'You can skip\nall of this if you have already done it, or have ', ... + 'configured mex to use\na different compiler. If you are using ', ... + 'Norton anti-virus software on Windows\n98SE, then you need to ', ... + 'exit MATLAB, turn off virus checking, and restart MATLAB\n', ... + 'before you can use the mex command or compile UMFPACK.\n', ... + 'You may also need to turn off virus checking in other cases.\n', ... + '\nHit enter to continue, or type control-C if you do not wish to '] ; + fprintf (msg) ; + input ('proceed: ') ; + end + fprintf ('\nUsing the BLAS (recommended).\n') ; +else + % No BLAS + fprintf ('\nNot using the BLAS. UMFPACK will be slow.\n') ; + blas = ' -DNBLAS' ; +end + +%------------------------------------------------------------------------------- +% -DNUTIL option (using utMalloc or mxMalloc) +%------------------------------------------------------------------------------- + +utils = '' ; + +if (~pc) + msg = [ ... + '--------------------------------------------------------------\n', ... + '\nUMFPACK uses MATLAB''s memory allocation routines. The internal', ... + '\nutMalloc, utFree, and utRealloc allow for better use of memory,', ... + '\nbut they are internal utility routines that are not documented.\n', ... + 'Thus, they might not always work. Using mxMalloc, mxFree, and\n', ... + 'mxRealloc works, but UMFPACK might run out of memory when solving\n', ... + 'problems that it could otherwise solve. Try using the default.\n', ... + 'If you get an error stating that utMalloc, utFree, and/or\n', ... + 'utRealloc are not defined, then recompile with the mx* routines.\n', ... + '\nPlease select one of the following options:\n', ... + ' 1: attempt to use the ut* routines (default)\n', ... + ' 2: use the standard mx* routines\n'] ; + fprintf (msg) ; + utils = input (': ') ; + if (isempty (utils)) + utils = 1 ; + end + if (utils == 2) + fprintf ('\nNot using utMalloc, utFree, or utRealloc\n') ; + utils = ' -DNUTIL' ; + else + fprintf ('\nUsing utMalloc, utFree, and utRealloc\n') ; + utils = '' ; + end +end + +%------------------------------------------------------------------------------- +% -DNPOSIX option (for sysconf and times timer routines) +%------------------------------------------------------------------------------- + +posix = '' ; + +if (~pc) + msg = [ ... + '--------------------------------------------------------------\n', ... + '\nUMFPACK can use the POSIX routines sysconf () and times ()\n', ... + 'to provide CPU time and wallclock time statistics. If you do not\n', ... + 'have a POSIX-compliant operating system, then UMFPACK won''t\n', ... + 'compile. If you don''t know which option to pick, try the\n', ... + 'default. If you get an error saying that sysconf and/or times\n', ... + 'are not defined, then recompile with the non-POSIX option.\n', ... + '\nPlease select one of the following options:\n', ... + ' 1: use POSIX sysconf and times routines (default)\n', ... + ' 2: do not use POSIX routines\n'] ; + fprintf (msg) ; + posix = input (': ') ; + if (isempty (posix)) + posix = 1 ; + end + if (posix == 2) + fprintf ('\nNot using POSIX sysconf and times routines.\n') ; + posix = ' -DNPOSIX' ; + else + fprintf ('\nUsing POSIX sysconf and times routines.\n') ; + posix = '' ; + end +end + +%------------------------------------------------------------------------------- +% mex command +%------------------------------------------------------------------------------- + +umfdir = sprintf ('..%sSource%s', filesep, filesep) ; +amddir = sprintf ('..%s..%sAMD%sSource%s', filesep, filesep, filesep, filesep) ; +incdir = sprintf ( ... +' -I..%sInclude -I..%sSource -I..%s..%sAMD%sInclude -I..%s..%sAMD%sSource', ... +filesep,filesep, filesep, filesep, filesep, filesep, filesep, filesep) ; + +mx = sprintf ('mex -inline -O%s%s%s%s', blas, utils, posix, incdir) ; +msg = [ ... + '--------------------------------------------------------------\n', ... + '\nCompile options:\n%s\nNow compiling. Please wait.\n'] ; +fprintf (msg, mx) ; + +% The following is adapted from GNUmakefile + +%------------------------------------------------------------------------------- +% source files +%------------------------------------------------------------------------------- + +% non-user-callable umf_*.[ch] files: +umfch = { 'assemble', 'blas3_update', ... + 'build_tuples', 'create_element', ... + 'dump', 'extend_front', 'garbage_collection', ... + 'get_memory', 'init_front', 'kernel', ... + 'kernel_init', 'kernel_wrapup', ... + 'local_search', 'lsolve', 'ltsolve', ... + 'mem_alloc_element', 'mem_alloc_head_block', ... + 'mem_alloc_tail_block', 'mem_free_tail_block', ... + 'mem_init_memoryspace', ... + 'report_vector', 'row_search', 'scale_column', ... + 'set_stats', 'solve', 'symbolic_usage', 'transpose', ... + 'tuple_lengths', 'usolve', 'utsolve', 'valid_numeric', ... + 'valid_symbolic', 'grow_front', 'start_front', '2by2', ... + 'store_lu', 'scale' } ; + +% non-user-callable umf_*.[ch] files, int versions only (no real/complex): +umfint = { 'analyze', 'apply_order', 'colamd', 'free', 'fsize', ... + 'is_permutation', 'malloc', 'realloc', 'report_perm', ... + 'singletons' } ; + +% non-user-callable and user-callable amd_*.[ch] files (int versions only): +amd = { 'aat', '1', '2', 'dump', 'postorder', 'post_tree', 'defaults', ... + 'order', 'control', 'info', 'valid' } ; + +% user-callable umfpack_*.[ch] files (real/complex): +user = { 'col_to_triplet', 'defaults', 'free_numeric', ... + 'free_symbolic', 'get_numeric', 'get_lunz', ... + 'get_symbolic', 'numeric', 'qsymbolic', ... + 'report_control', 'report_info', 'report_matrix', ... + 'report_numeric', 'report_perm', 'report_status', ... + 'report_symbolic', 'report_triplet', ... + 'report_vector', 'solve', 'symbolic', ... + 'transpose', 'triplet_to_col', 'scale' ... + 'load_numeric', 'save_numeric', 'load_symbolic', 'save_symbolic' } ; + +% user-callable umfpack_*.[ch], only one version +generic = { 'timer', 'tictoc' } ; + +M = cell (0) ; + +%------------------------------------------------------------------------------- +% Create the umfpack and amd mexFunctions for MATLAB (int versions only) +%------------------------------------------------------------------------------- + +for k = 1:length(umfint) + M = make (M, '%s -DDINT -c %sumf_%s.c', 'umf_%s.%s', 'umf_%s_%s.%s', ... + mx, umfint {k}, umfint {k}, 'm', obj, umfdir) ; +end + +rules = { [mx ' -DDINT'] , [mx ' -DZINT'] } ; +kinds = { 'md', 'mz' } ; + +for what = 1:2 + + rule = rules {what} ; + kind = kinds {what} ; + + M = make (M, '%s -DCONJUGATE_SOLVE -c %sumf_%s.c', 'umf_%s.%s', ... + 'umf_%s_%s.%s', rule, 'ltsolve', 'lhsolve', kind, obj, umfdir) ; + + M = make (M, '%s -DCONJUGATE_SOLVE -c %sumf_%s.c', 'umf_%s.%s', ... + 'umf_%s_%s.%s', rule, 'utsolve', 'uhsolve', kind, obj, umfdir) ; + + M = make (M, '%s -DDO_MAP -c %sumf_%s.c', 'umf_%s.%s', ... + 'umf_%s_%s_map_nox.%s', rule, 'triplet', 'triplet', kind, obj, umfdir) ; + + M = make (M, '%s -DDO_VALUES -c %sumf_%s.c', 'umf_%s.%s', ... + 'umf_%s_%s_nomap_x.%s', rule, 'triplet', 'triplet', kind, obj, umfdir) ; + + M = make (M, '%s -c %sumf_%s.c', 'umf_%s.%s', ... + 'umf_%s_%s_nomap_nox.%s', rule, 'triplet', 'triplet', kind, obj, ... + umfdir) ; + + M = make (M, '%s -DDO_MAP -DDO_VALUES -c %sumf_%s.c', 'umf_%s.%s', ... + 'umf_%s_%s_map_x.%s', rule, 'triplet', 'triplet', kind, obj, umfdir) ; + + M = make (M, '%s -DFIXQ -c %sumf_%s.c', 'umf_%s.%s', ... + 'umf_%s_%s_fixq.%s', rule, 'assemble', 'assemble', kind, obj, umfdir) ; + + M = make (M, '%s -DDROP -c %sumf_%s.c', 'umf_%s.%s', ... + 'umf_%s_%s_drop.%s', rule, 'store_lu', 'store_lu', kind, obj, umfdir) ; + + for k = 1:length(umfch) + M = make (M, '%s -c %sumf_%s.c', 'umf_%s.%s', 'umf_%s_%s.%s', ... + rule, umfch {k}, umfch {k}, kind, obj, umfdir) ; + end + + M = make (M, '%s -DWSOLVE -c %sumfpack_%s.c', 'umfpack_%s.%s', ... + 'umfpack_%s_w%s.%s', rule, 'solve', 'solve', kind, obj, umfdir) ; + + for k = 1:length(user) + M = make (M, '%s -c %sumfpack_%s.c', 'umfpack_%s.%s', ... + 'umfpack_%s_%s.%s', rule, user {k}, user {k}, kind, obj, umfdir) ; + end +end + +for k = 1:length(generic) + M = make (M, '%s -c %sumfpack_%s.c', 'umfpack_%s.%s', ... + 'umfpack_%s_%s.%s', mx, generic {k}, generic {k}, 'm', obj, umfdir) ; +end + +%---------------------------------------- +% AMD routines (int only) +%---------------------------------------- + +for k = 1:length(amd) + M = make (M, '%s -DDINT -c %samd_%s.c', 'amd_%s.%s', 'amd_%s_%s.%s', ... + mx, amd {k}, amd {k}, 'm', obj, amddir) ; +end + +%---------------------------------------- +% compile the umfpack mexFunction +%---------------------------------------- + +C = sprintf ('%s -output umfpack umfpackmex.c', mx) ; +for i = 1:length (M) + C = [C ' ' (M {i})] ; +end +C = [C ' ' blas_lib] ; +cmd (C) ; + +%---------------------------------------- +% delete the object files +%---------------------------------------- + +for i = 1:length (M) + rmfile (M {i}) ; +end + +%---------------------------------------- +% compile the luflop mexFunction +%---------------------------------------- + +cmd (sprintf ('%s -output luflop luflopmex.c', mx)) ; + +fprintf ('\n\nCompilation has completed. Now trying the umfpack_simple demo.\n'); +umfpack_simple + +%------------------------------------------------------------------------------- +% rmfile: delete a file, but only if it exists +%------------------------------------------------------------------------------- + +function rmfile (file) +if (length (dir (file)) > 0) + delete (file) ; +end + +%------------------------------------------------------------------------------- +% cpfile: copy the src file to the filename dst, overwriting dst if it exists +%------------------------------------------------------------------------------- + +function cpfile (src, dst) +rmfile (dst) +if (length (dir (src)) == 0) + help umfpack_make + error (sprintf ('File does not exist: %s\n', src)) ; +end +copyfile (src, dst) ; + +%------------------------------------------------------------------------------- +% mvfile: move the src file to the filename dst, overwriting dst if it exists +%------------------------------------------------------------------------------- + +function mvfile (src, dst) +cpfile (src, dst) ; +rmfile (src) ; + +%------------------------------------------------------------------------------- +% cmd: display and execute a command +%------------------------------------------------------------------------------- + +function cmd (s) +fprintf ('.') ; +eval (s) ; + +%------------------------------------------------------------------------------- +% make: execute a "make" command for a source file +%------------------------------------------------------------------------------- + +function M = make (M, s, src, dst, rule, file1, file2, kind, obj, srcdir) +cmd (sprintf (s, rule, srcdir, file1)) ; +src = sprintf (src, file1, obj) ; +dst = sprintf (dst, kind, file2, obj) ; +mvfile (src, dst) ; +M {end + 1} = dst ; + diff -uNr UMFPACKv4.4.orig/UMFPACK/OCTAVE/umfpack_report.m UMFPACK/UMFPACK/OCTAVE/umfpack_report.m --- UMFPACKv4.4.orig/UMFPACK/OCTAVE/umfpack_report.m 1970-01-01 01:00:00.000000000 +0100 +++ UMFPACK/UMFPACK/OCTAVE/umfpack_report.m 2004-12-30 01:58:46.000000000 +0100 @@ -0,0 +1,346 @@ +function umfpack_report (Control, Info) +% UMFPACK_REPORT +% +% umfpack_report (Control, Info) ; +% +% Prints the current Control settings for umfpack, and the statistical +% information returned by umfpack in the Info array. If Control is +% an empty matrix, then the default control settings are printed. +% +% Control is 20-by-1, and Info is 90-by-1. Not all entries are used. +% +% Alternative usages: +% +% umfpack_report ([ ], Info) ; print the default control parameters +% and the Info array. +% umfpack_report (Control) ; print the control parameters only. +% umfpack_report ; print the default control parameters +% and an empty Info array. +% +% See also umfpack, umfpack_make, umfpack_details, +% umfpack_demo, and umfpack_simple. + +% UMFPACK Version 4.3 (Jan. 16, 2004), Copyright (c) 2004 by Timothy A. +% Davis. All Rights Reserved. See ../README for License. + +%------------------------------------------------------------------------------- +% get inputs, use defaults if input arguments not present +%------------------------------------------------------------------------------- + +% The contents of Control and Info are defined in umfpack.h +if (nargin < 1) + Control = [] ; +end +if (nargin < 2) + Info = [] ; +end +if (isempty (Control)) + Control = umfpack ; +end +if (isempty (Info)) + Info = [ 0 (-ones (1, 89)) ] ; +end + +%------------------------------------------------------------------------------- +% control settings +%------------------------------------------------------------------------------- + +fprintf ('\nUMFPACK Version 4.3: Control settings:\n\n') ; +fprintf (' Control (1): print level: %d\n', Control (1)) ; +fprintf (' Control (2): dense row parameter: %g\n', Control (2)) ; +fprintf (' "dense" rows have > max (16, (%g)*16*sqrt(n_col)) entries\n', Control (2)) ; +fprintf (' Control (3): dense column parameter: %g\n', Control (3)) ; +fprintf (' "dense" columns have > max (16, (%g)*16*sqrt(n_row)) entries\n', Control (3)) ; +fprintf (' Control (4): pivot tolerance: %g\n', Control (4)) ; +fprintf (' Control (5): max block size for dense matrix kernels: %d\n', Control (5)) ; +prstrat (' Control (6): strategy: %g ', Control (6)) ; +fprintf (' Control (7): initial allocation ratio: %g\n', Control (7)) ; +fprintf (' Control (8): max iterative refinement steps: %d\n', Control (8)) ; +fprintf (' Control (13): 2-by-2 pivot tolerance: %g\n', Control (13)) ; +fprintf (' Control (14): Q fixed during numeric factorization: %g ', Control (14)) ; +if (Control (14) > 0) + fprintf ('(yes)\n') ; +elseif (Control (14) < 0) + fprintf ('(no)\n') ; +else + fprintf ('(auto)\n') ; +end +fprintf (' Control (15): AMD dense row/column parameter: %g\n', Control (15)) ; +fprintf (' "dense" rows/columns in A+A'' have > max (16, (%g)*sqrt(n)) entries.\n', Control (15)) ; +fprintf (' Only used if the AMD ordering is used.\n') ; +fprintf (' Control (16): diagonal pivot tolerance: %g\n', Control (16)) ; +fprintf (' Only used if diagonal pivoting is attempted.\n') ; + +fprintf (' Control (17): scaling option: %g ', Control (17)) ; +if (Control (17) == 0) + fprintf ('(none)\n') ; +elseif (Control (17) == 2) + fprintf ('(scale the matrix by\n') ; + fprintf (' dividing each row by max. abs. value in each row)\n') ; +else + fprintf ('(scale the matrix by\n') ; + fprintf (' dividing each row by sum of abs. values in each row)\n') ; +end + +fprintf (' Control (18): frontal matrix allocation ratio: %g\n', Control (18)) ; +fprintf (' Control (19): drop tolerance: %g\n', Control (19)) ; +fprintf (' Control (20): AMD and COLAMD aggressive absorption: %g ', Control (20)) ; +yes_no (Control (20)) ; + +% compile-time options: + +fprintf ('\n The following options can only be changed at compile-time:\n') ; + +if (Control (9) == 1) + fprintf (' Control (9): compiled to use the BLAS\n') ; +else + fprintf (' Control (9): compiled without the BLAS\n') ; + fprintf (' (you will not get the best possible performance)\n') ; +end + +if (Control (10) == 1) + fprintf (' Control (10): compiled for MATLAB\n') ; +elseif (Control (10) == 2) + fprintf (' Control (10): compiled for MATLAB\n') ; + fprintf (' Uses internal utMalloc, utFree, utRealloc, utPrintf\n') ; + fprintf (' utDivideComplex, and utFdlibm_hypot routines.\n') ; +else + fprintf (' Control (10): not compiled for MATLAB\n') ; + fprintf (' Uses ANSI C malloc, free, realloc, and printf\n') ; + fprintf (' instead of mxMalloc, mxFree, mxRealloc, and mexPrintf.\n') ; + fprintf (' Printing will be in terms of 0-based matrix indexing,\n') ; + fprintf (' not 1-based as is expected in MATLAB. Diary output may\n') ; + fprintf (' not be properly recorded.\n') ; +end + +if (Control (11) == 2) + fprintf (' Control (11): uses POSIX times ( ) to get CPU time and wallclock time.\n') ; +elseif (Control (11) == 1) + fprintf (' Control (11): uses getrusage to get CPU time.\n') ; +else + fprintf (' Control (11): uses ANSI C clock to get CPU time.\n') ; + fprintf (' The CPU time may wrap around, type "help cputime".\n') ; +end + +if (Control (12) == 1) + fprintf (' Control (12): compiled with debugging enabled\n') ; + fprintf (' ###########################################\n') ; + fprintf (' ### This will be exceedingly slow! ########\n') ; + fprintf (' ###########################################\n') ; + if (Control (10) == 1) + fprintf (' Uses mxAssert.\n') ; + elseif (Control (10) == 2) + fprintf (' Uses utAssert.\n') ; + else + fprintf (' Uses ANSI C assert instead of mxAssert.\n') ; + end +else + fprintf (' Control (12): compiled for normal operation (no debugging)\n') ; +end + +%------------------------------------------------------------------------------- +% Info: +%------------------------------------------------------------------------------- + +if (nargin == 1) + return +end + +status = Info (1) ; +fprintf ('\nUMFPACK status: Info (1): %d, ', status) ; + +if (status == 0) + fprintf ('OK\n') ; +elseif (status == 1) + fprintf ('WARNING matrix is singular\n') ; +elseif (status == -1) + fprintf ('ERROR out of memory\n') ; +elseif (status == -3) + fprintf ('ERROR numeric LU factorization is invalid\n') ; +elseif (status == -4) + fprintf ('ERROR symbolic LU factorization is invalid\n') ; +elseif (status == -5) + fprintf ('ERROR required argument is missing\n') ; +elseif (status == -6) + fprintf ('ERROR n <= 0\n') ; +elseif (status <= -7 & status >= -12 | status == -14) + fprintf ('ERROR matrix A is corrupted\n') ; +elseif (status == -13) + fprintf ('ERROR invalid system\n') ; +elseif (status == -15) + fprintf ('ERROR invalid permutation\n') ; +elseif (status == -911) + fprintf ('ERROR internal error!\n') ; + fprintf ('Please report this error to Tim Davis (davis@cise.ufl.edu)\n') ; +else + fprintf ('ERROR unrecognized error. Info array corrupted\n') ; +end + +fprintf (' (a -1 means the entry has not been computed):\n') ; + +fprintf ('\n Basic statistics:\n') ; +fprintf (' Info (2): %d, # of rows of A\n', Info (2)) ; +fprintf (' Info (17): %d, # of columns of A\n', Info (17)) ; +fprintf (' Info (3): %d, nnz (A)\n', Info (3)) ; +fprintf (' Info (4): %d, Unit size, in bytes, for memory usage reported below\n', Info (4)) ; +fprintf (' Info (5): %d, size of int (in bytes)\n', Info (5)) ; +fprintf (' Info (6): %d, size of long (in bytes)\n', Info (6)) ; +fprintf (' Info (7): %d, size of pointer (in bytes)\n', Info (7)) ; +fprintf (' Info (8): %d, size of numerical entry (in bytes)\n', Info (8)) ; + +fprintf ('\n Pivots with zero Markowitz cost removed to obtain submatrix S:\n') ; +fprintf (' Info (57): %d, # of pivots with one entry in pivot column\n', Info (57)) ; +fprintf (' Info (58): %d, # of pivots with one entry in pivot row\n', Info (58)) ; +fprintf (' Info (59): %d, # of rows/columns in submatrix S (if square)\n', Info (59)) ; +fprintf (' Info (60): %d ') ; +if (Info (60) > 0) + fprintf ('submatrix S square and diagonal preserved\n') ; +elseif (Info (60) == 0) + fprintf ('submatrix S not square or diagonal not preserved\n') ; +else + fprintf ('\n') ; +end +fprintf (' Info (9): %d, # of "dense" rows in S\n', Info (9)) ; +fprintf (' Info (10): %d, # of empty rows in S\n', Info (10)) ; +fprintf (' Info (11): %d, # of "dense" columns in S\n', Info (11)) ; +fprintf (' Info (12): %d, # of empty columns in S\n', Info (12)) ; +fprintf (' Info (34): %g, symmetry of pattern of S\n', Info (34)) ; +fprintf (' Info (35): %d, # of off-diagonal nonzeros in S+S''\n', Info (35)) ; +fprintf (' Info (36): %d, nnz (diag (S))\n', Info (36)) ; + +fprintf ('\n 2-by-2 pivoting to place large entries on diagonal:\n') ; +fprintf (' Info (52): %d, # of small diagonal entries of S\n', Info (52)) ; +fprintf (' Info (53): %d, # of unmatched small diagonal entries\n', Info (53)) ; +fprintf (' Info (54): %g, symmetry of P2*S\n', Info (54)) ; +fprintf (' Info (55): %d, # of off-diagonal entries in (P2*S)+(P2*S)''\n', Info (55)) ; +fprintf (' Info (56): %d, nnz (diag (P2*S))\n', Info (56)) ; + +fprintf ('\n AMD results, for strict diagonal pivoting:\n') ; +fprintf (' Info (37): %d, est. nz in L and U\n', Info (37)) ; +fprintf (' Info (38): %g, est. flop count\n', Info (38)) ; +fprintf (' Info (39): %g, # of "dense" rows in S+S''\n', Info (39)) ; +fprintf (' Info (40): %g, est. max. nz in any column of L\n', Info (40)) ; + +fprintf ('\n Final strategy selection, based on the analysis above:\n') ; +prstrat (' Info (19): %d, strategy used ', Info (19)) ; +fprintf (' Info (20): %d, ordering used ', Info (20)) ; +if (Info (20) == 0) + fprintf ('(COLAMD on A)\n') ; +elseif (Info (20) == 1) + fprintf ('(AMD on A+A'')\n') ; +elseif (Info (20) == 2) + fprintf ('(provided by user)\n') ; +else + fprintf ('(undefined ordering option)\n') ; +end +fprintf (' Info (32): %d, Q fixed during numeric factorization: ', Info (32)) ; +yes_no (Info (32)) ; +fprintf (' Info (33): %d, prefer diagonal pivoting: ', Info (33)) ; +yes_no (Info (33)) ; + +fprintf ('\n symbolic analysis time and memory usage:\n') ; +fprintf (' Info (13): %d, defragmentations during symbolic analysis\n', Info (13)) ; +fprintf (' Info (14): %d, memory used during symbolic analysis (Units)\n', Info (14)) ; +fprintf (' Info (15): %d, final size of symbolic factors (Units)\n', Info (15)) ; +fprintf (' Info (16): %.2f, symbolic analysis CPU time (seconds)\n', Info (16)) ; +fprintf (' Info (18): %.2f, symbolic analysis wall clock time (seconds)\n', Info (18)) ; + +fprintf ('\n Estimates computed in the symbolic analysis:\n') ; +fprintf (' Info (21): %d, est. size of LU factors (Units)\n', Info (21)) ; +fprintf (' Info (22): %d, est. total peak memory usage (Units)\n', Info (22)) ; +fprintf (' Info (23): %d, est. factorization flop count\n', Info (23)) ; +fprintf (' Info (24): %d, est. nnz (L)\n', Info (24)) ; +fprintf (' Info (25): %d, est. nnz (U)\n', Info (25)) ; +fprintf (' Info (26): %d, est. initial size, variable-part of LU (Units)\n', Info (26)) ; +fprintf (' Info (27): %d, est. peak size, of variable-part of LU (Units)\n', Info (27)) ; +fprintf (' Info (28): %d, est. final size, of variable-part of LU (Units)\n', Info (28)) ; +fprintf (' Info (29): %d, est. max frontal matrix size (# of entries)\n', Info (29)) ; +fprintf (' Info (30): %d, est. max # of rows in frontal matrix\n', Info (30)) ; +fprintf (' Info (31): %d, est. max # of columns in frontal matrix\n', Info (31)) ; + +fprintf ('\n Computed in the numeric factorization (estimates shown above):\n') ; +fprintf (' Info (41): %d, size of LU factors (Units)\n', Info (41)) ; +fprintf (' Info (42): %d, total peak memory usage (Units)\n', Info (42)) ; +fprintf (' Info (43): %d, factorization flop count\n', Info (43)) ; +fprintf (' Info (44): %d, nnz (L)\n', Info (44)) ; +fprintf (' Info (45): %d, nnz (U)\n', Info (45)) ; +fprintf (' Info (46): %d, initial size of variable-part of LU (Units)\n', Info (46)) ; +fprintf (' Info (47): %d, peak size of variable-part of LU (Units)\n', Info (47)) ; +fprintf (' Info (48): %d, final size of variable-part of LU (Units)\n', Info (48)) ; +fprintf (' Info (49): %d, max frontal matrix size (# of numerical entries)\n', Info (49)) ; +fprintf (' Info (50): %d, max # of rows in frontal matrix\n', Info (50)) ; +fprintf (' Info (51): %d, max # of columns in frontal matrix\n', Info (51)) ; + +fprintf ('\n Computed in the numeric factorization (no estimates computed a priori):\n') ; +fprintf (' Info (61): %d, defragmentations during numeric factorization\n', Info (61)) ; +fprintf (' Info (62): %d, reallocations during numeric factorization\n', Info (62)) ; +fprintf (' Info (63): %d, costly reallocations during numeric factorization\n', Info (63)) ; +fprintf (' Info (64): %d, integer indices in compressed pattern of L and U\n', Info (64)) ; +fprintf (' Info (65): %d, numerical values stored in L and U\n', Info (65)) ; +fprintf (' Info (66): %.2f, numeric factorization CPU time (seconds)\n', Info (66)) ; +fprintf (' Info (76): %.2f, numeric factorization wall clock time (seconds)\n', Info (76)) ; +if (Info (66) > 0.05 & Info (43) > 0) +fprintf (' mflops in numeric factorization phase: %.2f\n', 1e-6 * Info (43) / Info (66)) ; +end +fprintf (' Info (67): %d, nnz (diag (U))\n', Info (67)) ; +fprintf (' Info (68): %g, reciprocal condition number estimate\n', Info (68)) ; +fprintf (' Info (69): %g, matrix was ', Info (69)) ; +if (Info (69) == 0) + fprintf ('not scaled\n') ; +elseif (Info (69) == 2) + fprintf ('scaled (row max)\n') ; +else + fprintf ('scaled (row sum)\n') ; +end +fprintf (' Info (70): %g, min. scale factor of rows of A\n', Info (70)) ; +fprintf (' Info (71): %g, max. scale factor of rows of A\n', Info (71)) ; +fprintf (' Info (72): %g, min. abs. on diagonal of U\n', Info (72)) ; +fprintf (' Info (73): %g, max. abs. on diagonal of U\n', Info (73)) ; +fprintf (' Info (74): %g, initial allocation parameter used\n', Info (74)) ; +fprintf (' Info (75): %g, # of forced updates due to frontal growth\n', Info (75)) ; +fprintf (' Info (77): %d, # of off-diaogonal pivots\n', Info (77)) ; +fprintf (' Info (78): %d, nnz (L), if no small entries dropped\n', Info (78)) ; +fprintf (' Info (79): %d, nnz (U), if no small entries dropped\n', Info (79)) ; +fprintf (' Info (80): %d, # of small entries dropped\n', Info (80)) ; + +fprintf ('\n Computed in the solve step:\n') ; +fprintf (' Info (81): %d, iterative refinement steps taken\n', Info (81)) ; +fprintf (' Info (82): %d, iterative refinement steps attempted\n', Info (82)) ; +fprintf (' Info (83): %g, omega(1), sparse-backward error estimate\n', Info (83)) ; +fprintf (' Info (84): %g, omega(2), sparse-backward error estimate\n', Info (84)) ; +fprintf (' Info (85): %d, solve flop count\n', Info (85)) ; +fprintf (' Info (86): %.2f, solve CPU time (seconds)\n', Info (86)) ; +fprintf (' Info (87): %.2f, solve wall clock time (seconds)\n', Info (87)) ; + +fprintf ('\n Info (88:90): unused\n\n') ; + +%------------------------------------------------------------------------------- + +function prstrat (fmt, strategy) +fprintf (fmt, strategy) ; +if (strategy == 1) + fprintf ('(unsymmetric)\n') ; + fprintf (' Q = COLAMD (A), Q refined during numerical\n') ; + fprintf (' factorization, and no attempt at diagonal pivoting.\n') ; +elseif (strategy == 2) + fprintf ('(symmetric, with 2-by-2 pivoting)\n') ; + fprintf (' P2 = row permutation to place large values on the diagonal\n') ; + fprintf (' Q = AMD (P2*A+(P2*A)''), Q not refined during numeric factorization,\n') ; + fprintf (' and diagonal pivoting attempted.\n') ; +elseif (strategy == 3) + fprintf ('(symmetric)\n') ; + fprintf (' Q = AMD (A+A''), Q not refined during numeric factorization,\n') ; + fprintf (' and diagonal pivoting (P=Q'') attempted.\n') ; +else + strategy = 0 ; + fprintf ('(auto)\n') ; +end + +%------------------------------------------------------------------------------- + +function yes_no (s) +if (s == 0) + fprintf ('(no)\n') ; +else + fprintf ('(yes)\n') ; +end diff -uNr UMFPACKv4.4.orig/UMFPACK/OCTAVE/umfpack_simple.m UMFPACK/UMFPACK/OCTAVE/umfpack_simple.m --- UMFPACKv4.4.orig/UMFPACK/OCTAVE/umfpack_simple.m 1970-01-01 01:00:00.000000000 +0100 +++ UMFPACK/UMFPACK/OCTAVE/umfpack_simple.m 2004-12-30 01:58:46.000000000 +0100 @@ -0,0 +1,61 @@ +% umfpack_simple: a simple demo of UMFPACK +% +% UMFPACK Version 4.3 (Jan. 16, 2004), Copyright (c) 2004 by Timothy A. +% Davis. All Rights Reserved. +% +% UMFPACK License: +% +% Your use or distribution of UMFPACK or any modified version of +% UMFPACK implies that you agree to this License. +% +% THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY +% EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK. +% +% Permission is hereby granted to use or copy this program, provided +% that the Copyright, this License, and the Availability of the original +% version is retained on all copies. User documentation of any code that +% uses UMFPACK or any modified version of UMFPACK code must cite the +% Copyright, this License, the Availability note, and "Used by permission." +% Permission to modify the code and to distribute modified code is granted, +% provided the Copyright, this License, and the Availability note are +% retained, and a notice that the code was modified is included. This +% software was developed with support from the National Science Foundation, +% and is provided to you free of charge. +% +% Availability: http://www.cise.ufl.edu/research/sparse/umfpack +% +% See also: umfpack, umfpack_details + +help umfpack_simple +i = input ('Hit enter to agree to the above License: ', 's') ; +if (~isempty (i)) + error ('terminating') ; +end + +format short + +A = [ + 2 3 0 0 0 + 3 0 4 0 6 + 0 -1 -3 2 0 + 0 0 1 0 0 + 0 4 2 0 1 +] + +A = sparse (A) ; + +b = [8 45 -3 3 19]' + +fprintf ('Solution to Ax=b via UMFPACK:\n') ; +fprintf ('x1 = umfpack (A, ''\\'', b)\n') ; + +x1 = umfpack (A, '\\', b) + +fprintf ('Solution to Ax=b via OCTAVE:\n') ; +fprintf ('x2 = A\\b\n') ; + +x2 = A\b + +fprintf ('norm (x1-x2) should be small: %g\n', norm (x1-x2)) ; + +fprintf ('Type ''umfpack_demo'' for a full demo of UMFPACK\n') ; diff -uNr UMFPACKv4.4.orig/UMFPACK/OCTAVE/umfpack_simple.m.out UMFPACK/UMFPACK/OCTAVE/umfpack_simple.m.out --- UMFPACKv4.4.orig/UMFPACK/OCTAVE/umfpack_simple.m.out 1970-01-01 01:00:00.000000000 +0100 +++ UMFPACK/UMFPACK/OCTAVE/umfpack_simple.m.out 2004-12-30 01:58:46.000000000 +0100 @@ -0,0 +1,79 @@ +octave:4> umfpack_simple +umfpack_simple is the file: /home/dbateman/octave/devel/octave-sparse/UMFPACKv4.3/UMFPACK/OCTAVE/umfpack_simple.m + +umfpack_simple: a simple demo of UMFPACK + +UMFPACK Version 4.3 (Jan. 16, 2004), Copyright (c) 2004 by Timothy A. +Davis. All Rights Reserved. + +UMFPACK License: + + Your use or distribution of UMFPACK or any modified version of + UMFPACK implies that you agree to this License. + + THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY + EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK. + + Permission is hereby granted to use or copy this program, provided + that the Copyright, this License, and the Availability of the original + version is retained on all copies. User documentation of any code that + uses UMFPACK or any modified version of UMFPACK code must cite the + Copyright, this License, the Availability note, and "Used by permission." + Permission to modify the code and to distribute modified code is granted, + provided the Copyright, this License, and the Availability note are + retained, and a notice that the code was modified is included. This + software was developed with support from the National Science Foundation, + and is provided to you free of charge. + +Availability: http://www.cise.ufl.edu/research/sparse/umfpack + +See also: umfpack, umfpack_details + + +Additional help for built-in functions, operators, and variables +is available in the on-line version of the manual. Use the command +`help -i <topic>' to search the manual index. + +Help and information about Octave is also available on the WWW +at http://www.octave.org and via the help@octave.org +mailing list. +Hit enter to agree to the above License: +A = + + 2 3 0 0 0 + 3 0 4 0 6 + 0 -1 -3 2 0 + 0 0 1 0 0 + 0 4 2 0 1 + +b = + + 8 + 45 + -3 + 3 + 19 + +Solution to Ax=b via UMFPACK: +x1 = umfpack (A, '\', b) +x1 = + + 1.00000 + 2.00000 + 3.00000 + 4.00000 + 5.00000 + +Solution to Ax=b via OCTAVE: +x2 = A\b +x2 = + + 1.00000 + 2.00000 + 3.00000 + 4.00000 + 5.00000 + +norm (x1-x2) should be small: 0 +Type 'umfpack_demo' for a full demo of UMFPACK +octave:5> diary off diff -uNr UMFPACKv4.4.orig/UMFPACK/OCTAVE/umfpack_solve.m UMFPACK/UMFPACK/OCTAVE/umfpack_solve.m --- UMFPACKv4.4.orig/UMFPACK/OCTAVE/umfpack_solve.m 1970-01-01 01:00:00.000000000 +0100 +++ UMFPACK/UMFPACK/OCTAVE/umfpack_solve.m 2004-12-30 01:58:46.000000000 +0100 @@ -0,0 +1,97 @@ +function x = umfpack_solve (arg1, op, arg2, Control) +% UMFPACK_SOLVE +% +% x = umfpack_solve (A, '\', b, Control) +% x = umfpack_solve (b, '/', A, Control) +% +% Computes x = A\b, or b/A, where A is square. Uses UMFPACK if A is sparse. +% The Control argument is optional. +% +% See also umfpack, umfpack_make, umfpack_details, umfpack_report, +% and umfpack_simple. + +% UMFPACK Version 4.3 (Jan. 16, 2004), Copyright (c) 2004 by Timothy A. +% Davis. All Rights Reserved. Type umfpack_details for License. + +%------------------------------------------------------------------------------- +% check inputs and get default control parameters +%------------------------------------------------------------------------------- + +if (op == '\\') + A = arg1 ; + b = arg2 ; +elseif (op == '/') + A = arg2 ; + b = arg1 ; +else + help umfack_solve + error ('umfpack_solve: unrecognized operator') ; +end + +[m n] = size (A) ; +if (m ~= n) + help umfpack_solve + error ('umfpack_solve: A must be square') ; +end + +[m1 n1] = size (b) ; +if ((op == '\\' & n ~= m1) | (op == '/' & n1 ~= m)) + help umfpack_solve + error ('umfpack_solve: b has the wrong dimensions') ; +end + +if (nargin < 4) + Control = umfpack ; +end + +%------------------------------------------------------------------------------- +% solve the system +%------------------------------------------------------------------------------- + +if (op == '\\') + + if (~issparse (A)) + + % A is not sparse, so just use MATLAB + x = A\b ; + + elseif (n1 == 1 & ~issparse (b)) + + % the UMFPACK '\' requires b to be a dense column vector + x = umfpack (A, '\\', b, Control) ; + + else + + % factorize with UMFPACK and do the forward/back solves in MATLAB + [L, U, P, Q, R] = umfpack (A, Control) ; + keyboard + x = Q * (U \ (L \ (P * (R \ b)))) ; + + end + +else + + if (~issparse (A)) + + % A is not sparse, so just use MATLAB + x = b/A ; + + elseif (m1 == 1 & ~issparse (b)) + + % the UMFPACK '/' requires b to be a dense column vector + x = umfpack (b, '/', A, Control) ; + + else + + % factorize with UMFPACK and do the forward/back solves in MATLAB + % this mimics the behavior of x = b/A, except for the row scaling + [L, U, P, Q, R] = umfpack (A.', Control) ; + x = (Q * (U \ (L \ (P * (R \ (b.')))))).' ; + + % an alternative method: + % [L, U, P, Q, r] = umfpack (A, Control) ; + % x = (R \ (P' * (L.' \ (U.' \ (Q' * b.'))))).' ; + + end + +end diff -uNr UMFPACKv4.4.orig/AMD/OCTAVE/amd.cc UMFPACK/AMD/OCTAVE/amd.cc --- UMFPACKv4.4.orig/AMD/OCTAVE/amd.cc 1970-01-01 01:00:00.000000000 +0100 +++ UMFPACK/AMD/OCTAVE/amd.cc 2005-02-01 21:51:23.944874478 +0100 @@ -0,0 +1,299 @@ +/* + +Copyright (C) 2004 David Bateman + +This program is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +This program is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with this program; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +In addition to the terms of the GPL, you are permitted to link +this program with any Open Source program, as defined by the +Open Source Initiative (www.opensource.org) + +*/ + +/* + +This is the Octave interface to the UMFPACK AMD code, which bore the following +copyright + + AMD Version 1.1 (Jan. 21, 2004), Copyright (c) 2004 by Timothy A. Davis, + Patrick R. Amestoy, and Iain S. Duff. See ../README for License. + email: davis@cise.ufl.edu CISE Department, Univ. of Florida. + web: http://www.cise.ufl.edu/research/sparse/amd + -------------------------------------------------------------------------- + + Acknowledgements: This work was supported by the National Science + Foundation, under grants ASC-9111263, DMS-9223088, and CCR-0203270. + +*/ + +#include <cstdlib> +#include <string> + +#include <octave/config.h> +#include <octave/ov.h> +#include <octave/defun-dld.h> +#include <octave/pager.h> +#include <octave/ov-re-mat.h> + +#include "ov-re-sparse.h" +#include "ov-cx-sparse.h" + +// External AMD functions in C +extern "C" { +#include "amd.h" +} + +DEFUN_DLD (amd, args, nargout, + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{p} =} amd (@var{s})\n\ +@deftypefnx {Loadable Function} {@var{Control} =} amd ()\n\ +@deftypefnx {Loadable Function} {[@var{p}, @var{info}] =} amd (@var{s})\n\ +\n\ +AMD Approximate minimum degree permutation. Returns the approximate\n\ +minimum degree permutation vector for the sparse matrix\n\ +@code{@var{c} = @var{S} + @var{S}'}. The Cholesky factorization of\n\ +@code{@var{c} (@var{p}, @var{p})}, or @code{@var{s} (@var{p}, @var{p})},\n\ +tends to be sparser than that of @var{c} or @var{s}.\n\ +@var{s} must be square. If @var{s} is full, @code{amd (@var{S})} is\n\ +equivalent to @code{amd (sparse (@var{s}))}.\n\ +\n\ +@table @asis\n\ +@item @var{Control} (1)\n\ +If S is n-by-n, then rows/columns with more than\n\ +@code{@dfn{max} (16, (@var{Control} (1)) * @dfn{sqrt} (@var{n}))} entries\n\ +in @code{@var{s} + @var{S}'} are considered @emph{dense}, and ignored during\n\ +ordering. They are placed last in the output permutation. The default is\n\ +10.0 if @var{Control} is not present.\n\ +@item @var{Control} (2)\n\ +If nonzero, then aggressive absorption is performed. This is the default if\n\ +@var{Control} is not present.\n\ +@item @var{Control} (3)\n\ +If nonzero, print statistics about the ordering.\n\ +@end table\n\ +\n\ +@table @asis\n\ +@item @var{Info} (1)\n\ +status (0: ok, -1: out of memory, -2: matrix invalid)\n\ +@item @var{Info} (2)\n\ +@code{@var{n} = size (@var{a}, 1)}\n\ +@item @var{Info} (3)\n\ +@code{nnz (A)}\n\ +@item @var{Info} (4)\n\ +The symmetry of the matrix @var{s} (0.0 means purely unsymmetric, 1.0 means\n\ +purely symmetric). Computed as: @code{@var{b} = tril (@var{s}, -1) +\n\ +triu (@var{s}, 1); @var{symmetry} = nnz (@var{b} & @var{b}')\n\ +/ nnz (@var{b});}\n\ +@item @var{Info} (5)\n\ +@code{nnz (diag (@var{s}))}\n\ +@item @var{Info} (6)\n\ +@dfn{nnz} in @code{@var{s} + @var{s}'}, excluding the diagonal\n\ +(= @code{nnz (@var{b} + @var{b}')})\n\ +@item @var{Info} (7)\n\ +Number of @emph{dense} rows/columns in @code{@var{s} + @var{s}'}\n\ +@item @var{Info} (8)\n\ +The amount of memory used by AMD, in bytes\n\ +@item @var{Info} (9)\n\ +The number of memory compactions performed by AMD\n\ +@end table\n\ +\n\ +The following statistics are slight upper bounds because of the\n\ +approximate degree in AMD. The bounds are looser if @emph{dense}\n\ +rows/columns are ignored during ordering @code{(@var{Info} (7) > 0)}.\n\ +The statistics are for a subsequent factorization of the matrix\n\ +@code{@var{c} (@var{p},@var{p})}. The LU factorization statistics assume\n \ +no pivoting.\n\ +\n\ +@table @asis\n\ +@item @var{Info} (10)\n\ +The number of nonzeros in L, excluding the diagonal\n\ +@item @var{Info} (11)\n\ +The number of divide operations for LL', LDL', or LU\n\ +@item @var{Info (12)}\n\ +The number of multiply-subtract pairs for LL' or LDL'\n\ +@item @var{Info} (13)\n\ +The number of multiply-subtract pairs for LU\n\ +@item @var{Info} (14)\n\ +The max number of nonzeros in any column of L (incl. diagonal)\n\ +@item @var{Info} (15:20)\n\ +unused, reserved for future use\n\ +@end table\n\ +\n\ +An assembly tree post-ordering is performed, which is typically the same\n\ +as an elimination tree post-ordering. It is not always identical because\n\ +of the approximate degree update used, and because @emph{dense} rows/columns\n\ +do not take part in the post-order. It well-suited for a subsequent\n\ +@dfn{chol}, however. If you require a precise elimination tree\n\ +post-ordering, then do:\n\ +\n\ +@group\n\ + @var{p} = @dfn{amd} (@var{s});\n\ + % skip this if S already symmetric\n\ + @var{c} = spones (@var{s}) + spones (@var{s}');\n\ + [@var{ignore}, @var{q}] = sparsfun ('symetree', @var{c} (@var{p}, @var{p}));\n\ + @var{p} = @var{p} (@var{q});\n\ +@end group\n\ +\n\ +AMD Version 1.1 (Jan. 21, 2004), Copyright @copyright{} 2004 by\n\ +Timothy A. Davis, Patrick R. Amestoy, and Iain S. Duff.\n\ +\n\ +email: davis@@cise.ufl.edu (CISE Department, Univ. of Florida).\n\ +\n\ +web: http://www.cise.ufl.edu/research/sparse/amd\n\ +\n\ +Acknowledgements: This work was supported by the National Science\n\ +Foundation, under grants ASC-9111263, DMS-9223088, and CCR-0203270.\n\ +@end deftypefn") +{ + int nargin = args.length (); + octave_value_list retval; + int spumoni = 0; + + if (nargin > 2 || nargout > 2) + usage ("p = amd (A) or [p, Info] = amd (A, Control)"); + else if (nargin == 0) + { + // Get the default control parameter, and return + NDArray control (dim_vector (AMD_CONTROL, 1)); + double *control_ptr = control.fortran_vec (); + amd_defaults (control_ptr); + if (nargout == 0) + amd_control (control_ptr); + else + retval(0) = control; + } + else + { + NDArray control (dim_vector (AMD_CONTROL, 1)); + double *control_ptr = control.fortran_vec (); + amd_defaults (control_ptr); + + if (nargin > 1) + { + NDArray control_in = args(1).array_value(); + + if (error_state) + { + error ("amd: could not read control vector"); + return retval; + } + + dim_vector dv = control_in.dims (); + if (dv.length() > 2 || (dv(0) != 1 && dv(1) != 1)) + { + error ("amd: control vector isn't a vector"); + return retval; + } + + int nc = dv.numel (); + control (AMD_DENSE) = (nc > 0 ? control_in (AMD_DENSE) : + AMD_DEFAULT_DENSE); + control (AMD_AGGRESSIVE) = (nc > 1 ? control_in (AMD_AGGRESSIVE) : + AMD_DEFAULT_AGGRESSIVE); + spumoni = (nc > 2 ? (control_in (2) != 0) : 0); + } + + if (spumoni > 0) + amd_control (control_ptr); + + int *Ap, *Ai; + int n, m, nz; + + // These are here only so that the C++ destructors don't prematurally + // remove the underlying data we are interested in + SparseMatrix sm; + SparseComplexMatrix scm; + Matrix mm; + ComplexMatrix cm; + + if (args(0).class_name () != "sparse" && spumoni > 0) + octave_stdout << " input matrix A is full (sparse copy" + << " of A will be created)" << std::endl; + + if (args(0).is_complex_type ()) + { + scm = args(0).sparse_complex_matrix_value (); + Ai = scm.ridx (); + Ap = scm.cidx (); + m = scm.rows (); + n = scm.cols (); + nz = scm.nnz (); + } + else + { + sm = args(0).sparse_matrix_value (); + Ai = sm.ridx (); + Ap = sm.cidx (); + m = sm.rows (); + n = sm.cols (); + nz = sm.nnz (); + } + + if (spumoni > 0) + octave_stdout << " input matrix A is " << m << "-by-" << n + << std::endl; + + if (m != n) + { + error ("amd: A must be square"); + return retval; + } + + if (spumoni > 0) + octave_stdout << " input matrix A has " << nz << + " nonzero entries" << std::endl; + + // allocate workspace for output permutation + Array<int> P(n+1); + int *P_ptr = P.fortran_vec (); + NDArray info (dim_vector (AMD_INFO, 1)); + double *info_ptr = info.fortran_vec (); + int result; + + // order the matrix + result = amd_order (n, Ap, Ai, P_ptr, control_ptr, info_ptr); + + // print results (including return value) + if (spumoni > 0) + amd_info (info_ptr); + + // check error conditions + if (result == AMD_OUT_OF_MEMORY) + error ("amd: out of memory"); + else if (result == AMD_INVALID) + error ("amd: input matrix A is corrupted"); + else + { + // copy the outputs to Octave + + // output permutation, P + NDArray perm (dim_vector (1, n)); + for (int i = 0; i < n; i++) + perm (i) = double (P(i) + 1); // 1-based indexing for Octave + + retval (0) = perm; + + // Info + if (nargout > 1) + retval (1) = info; + } + } + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ diff -uNr UMFPACKv4.4.orig/AMD/OCTAVE/amd_demo.m UMFPACK/AMD/OCTAVE/amd_demo.m --- UMFPACKv4.4.orig/AMD/OCTAVE/amd_demo.m 1970-01-01 01:00:00.000000000 +0100 +++ UMFPACK/AMD/OCTAVE/amd_demo.m 2004-12-30 01:58:45.000000000 +0100 @@ -0,0 +1,61 @@ +function amd_demo +% AMD DEMO +% +% A demo of AMD for OCTAVE. +% +% -------------------------------------------------------------------------- +% AMD Version 1.1 (Jan. 21, 2004), Copyright (c) 2004 by Timothy A. Davis, +% Patrick R. Amestoy, and Iain S. Duff. See ../README for License. +% email: davis@cise.ufl.edu CISE Department, Univ. of Florida. +% web: http://www.cise.ufl.edu/research/sparse/amd +% -------------------------------------------------------------------------- +% +% See also: amd + +% Get the Harwell/Boeing can_24 matrix. This is an example matrix from the +% MATLAB-accessible UF sparse matrix collection, and can be loaded into +% MATLAB with the statment "Problem = UFget ('HB/can_24')", after obtaining +% the UFget function and its supporting routines at +% http://www.cise.ufl.edu/sparse/mat . + +load can_24.mat +A = Problem.A ; +n = size (A,1) ; + +figure (1) +clf +hold off +subplot (2,1,1) ; +% remove the "_" from the name before printing it in the plot title +title (sprintf ('%s', strrep (Problem.name, '_', '-'))) ; +fprintf ('Matrix name: %s\n', Problem.name) ; +fprintf ('Matrix title: %s\n', Problem.title) ; +spy (A) + +% print the details during AMD ordering and SYMBFACT +control = amd (); +control (3) = 1; + +% order the matrix. Note that the Info argument is optional. +fprintf ('\nIf the next step fails, then you have\n') ; +fprintf ('not yet compiled the AMD mexFunction.\n') ; +[p, Info] = amd (A, control) ; + +% order again, but this time print some statistics +[p, Info] = amd (A, [10 1 1]) ; + +fprintf ('Permutation vector:\n') ; +fprintf (' %2d', p) ; +fprintf ('\n\n') ; + +subplot (2,1,2) ; +title ('Permuted matrix') ; +spy (A (p,p)) + +% approximations from amd: +lnz2 = n + Info (10) ; +fl2 = n + Info (11) + 2 * Info (12) ; +fprintf ('\nResults from AMD''s approximate analysis:\n') ; +fprintf ('number of nonzeros in L (including diagonal): %d\n', lnz2) ; +fprintf ('floating point operation count for chol (A (p,p)): %d\n\n', fl2) ; + diff -uNr UMFPACKv4.4.orig/AMD/OCTAVE/amd_demo.m.out UMFPACK/AMD/OCTAVE/amd_demo.m.out --- UMFPACKv4.4.orig/AMD/OCTAVE/amd_demo.m.out 1970-01-01 01:00:00.000000000 +0100 +++ UMFPACK/AMD/OCTAVE/amd_demo.m.out 2004-12-30 01:58:45.000000000 +0100 @@ -0,1 +1,21 @@ +octave:3> amd_demo +ans = 1 +Matrix name: HB/can_24 +Matrix title: 1SYMMETRIC PATTERN FROM CANNES,LUCIEN MARRO,JUNE 1981. + +If the next step fails, then you have +not yet compiled the AMD mexFunction. + input matrix A is 24-by-24 + input matrix A has 160 nonzero entries + input matrix A is 24-by-24 + input matrix A has 160 nonzero entries +Permutation vector: + 23 21 11 24 13 6 17 9 15 5 16 8 2 10 14 18 1 3 4 7 12 19 22 20 + + +Results from AMD's approximate analysis: +number of nonzeros in L (including diagonal): 121 +floating point operation count for chol (A (p,p)): 671 + +octave:4> quit --- UMFPACKv4.4.orig/AMD/OCTAVE/GNUmakefile 1970-01-01 01:00:00.000000000 +0100 +++ UMFPACK/AMD/OCTAVE/GNUmakefile 2004-12-30 01:58:45.000000000 +0100 @@ -0,0 +1,33 @@ +#------------------------------------------------------------------------------ +# GNUmakefile for the AMD MATLAB mexFunction +#------------------------------------------------------------------------------ + +all: amd + +include ../Make/Make.include + +MKOCT = mkoctfile -I../Include + +OCT_SPARSE_INC = -I../../../ + +AMD = amd_aat amd_1 amd_2 amd_dump amd_postorder amd_post_tree amd_defaults \ + amd_order amd_control amd_info amd_valid + +INC = ../Include/amd.h ../Source/amd_internal.h + +OCTAMD = $(addsuffix .o, $(subst amd_,amd_o_,$(AMD))) + +amd_o_%.o: ../Source/amd_%.c $(INC) + $(MKOCT) -DDINT -c $< -o $@ + - $(MV) ../Source/amd_$*.o + +# Note temporary addition of octave sparse path +amd: amd.cc $(OCTAMD) $(INC) + $(MKOCT) amd.cc $(OCTAMD) $(OCT_SPARSE_INC) -o amd.oct + +#------------------------------------------------------------------------------ +# Remove all but the files in the original distribution +#------------------------------------------------------------------------------ + +purge: clean + - $(RM) amd.oct diff -uNr UMFPACKv4.4.orig/AMD/OCTAVE/Makefile UMFPACK/AMD/OCTAVE/Makefile --- UMFPACKv4.4.orig/AMD/OCTAVE/Makefile 1970-01-01 01:00:00.000000000 +0100 +++ UMFPACK/AMD/OCTAVE/Makefile 2004-12-30 01:58:45.000000000 +0100 @@ -0,0 +1,41 @@ +#------------------------------------------------------------------------------ +# compile the AMD mexFunction for MATLAB (original make only) +#------------------------------------------------------------------------------ + +# This is a very ugly Makefile, and is only provided for those who do not +# have GNU make. Note that it is not used if you have GNU make. It ignores +# dependency checking and just compiles everything. It was created +# automatically, via make -n using the GNUmakefile. That way, I don't have +# maintain two Makefiles. + +all: amd + +include ../Make/Make.include + +MKOCT = mkoctfile -I../Include + +OCT_SPARSE_INC = ../../../ + +amd: + $(MKOCT) -DDINT -o amd_o_aat.o -c ../Source/amd_aat.c + $(MKOCT) -DDINT -o amd_o_1.o -c ../Source/amd_1.c + $(MKOCT) -DDINT -o amd_o_2.o -c ../Source/amd_2.c + $(MKOCT) -DDINT -o amd_o_dump.o -c ../Source/amd_dump.c + $(MKOCT) -DDINT -o amd_o_postorder.o -c ../Source/amd_postorder.c + $(MKOCT) -DDINT -o amd_o_post_tree.o -c ../Source/amd_post_tree.c + $(MKOCT) -DDINT -o amd_o_defaults.o -c ../Source/amd_defaults.c + $(MKOCT) -DDINT -o amd_o_order.o -c ../Source/amd_order.c + $(MKOCT) -DDINT -o amd_o_control.o -c ../Source/amd_control.c + $(MKOCT) -DDINT -o amd_o_info.o -c ../Source/amd_info.c + $(MKOCT) -DDINT -o amd_o_valid.o -c ../Source/amd_valid.c + $(MKOCT) -output amd.oct amd_mex.c amd_o_aat.o \ + amd_o_1.o amd_o_2.o amd_o_dump.o amd_o_postorder.o \ + amd_o_post_tree.o amd_o_defaults.o amd_o_order.o amd_o_control.o \ + amd_o_info.o amd_o_valid.o $(OCT_SPARSE_INC) -o amd.oct + +#------------------------------------------------------------------------------ +# Remove all but the files in the original distribution +#------------------------------------------------------------------------------ + +purge: clean + - $(RM) amd.oct --- UMFPACKv4.4.orig/UMFPACK/Source/umf_solve.c 2005-01-17 17:21:29.000000000 +0100 +++ UMFPACK/UMFPACK/Source/umf_solve.c 2005-02-02 00:01:22.904346651 +0100 @@ -79,7 +79,8 @@ double *Z2, *Y, *B2, *Rs ; Int *Rperm, *Cperm, i, n, p, step, j, nz, status, p2, do_scale ; #ifdef COMPLEX - Int split ; + Int AXsplit ; + Int Bsplit ; #endif #ifndef NRECIPROCAL Int do_recip = Numeric->do_recip ; @@ -141,7 +142,7 @@ return (UMFPACK_ERROR_argument_missing) ; } /* A, B, and X in split format if Az, Bz, and Xz present */ - split = SPLIT (Az) && SPLIT (Bz) && SPLIT (Xz) ; + AXsplit = SPLIT (Az) || SPLIT(Xz); Z = (Entry *) (SolveWork + 4*n) ; /* Entry Z [0..n-1] */ S = (Entry *) (SolveWork + 6*n) ; /* Entry S [0..n-1] */ Y = (double *) (SolveWork + 8*n) ; /* double Y [0..n-1] */ @@ -150,10 +151,12 @@ } else { - /* A is ignored, only look at X and B for split/packed cases */ - split = SPLIT (Bz) && SPLIT (Xz) ; + /* A is ignored, only look at X for split/packed cases */ + AXsplit = SPLIT(Xz); } - if (split) + Bsplit = SPLIT (Bz); + + if (AXsplit) { X = (Entry *) (SolveWork + 2*n) ; /* Entry X [0..n-1] */ } @@ -209,7 +212,7 @@ for (p = 0 ; p < p2 ; p++) { /* Y [Ai [p]] += ABS (Ax [p]) ; */ - ASSIGN (aij, Ax, Az, p, split) ; + ASSIGN (aij, Ax, Az, p, AXsplit) ; ABS (d, aij) ; Y [Ai [p]] += d ; } @@ -219,7 +222,7 @@ for (i = 0 ; i < n ; i++) { /* B2 [i] = ABS (B [i]) ; */ - ASSIGN (bi, Bx, Bz, i, split) ; + ASSIGN (bi, Bx, Bz, i, Bsplit) ; ABS (B2 [i], bi) ; } @@ -276,7 +279,7 @@ /* multiply by the scale factors */ for (i = 0 ; i < n ; i++) { - ASSIGN (X [i], Bx, Bz, i, split) ; + ASSIGN (X [i], Bx, Bz, i, Bsplit) ; SCALE (X [i], Rs [i]) ; } } @@ -286,7 +289,7 @@ /* divide by the scale factors */ for (i = 0 ; i < n ; i++) { - ASSIGN (X [i], Bx, Bz, i, split) ; + ASSIGN (X [i], Bx, Bz, i, Bsplit) ; SCALE_DIV (X [i], Rs [i]) ; } } @@ -302,7 +305,7 @@ for (i = 0 ; i < n ; i++) { /* W [i] = B [Rperm [i]] ; */ - ASSIGN (W [i], Bx, Bz, Rperm [i], split) ; + ASSIGN (W [i], Bx, Bz, Rperm [i], Bsplit) ; } } } @@ -311,7 +314,7 @@ for (i = 0 ; i < n ; i++) { /* Z [i] = B [i] ; */ - ASSIGN (Z [i], Bx, Bz, i, split) ; + ASSIGN (Z [i], Bx, Bz, i, Bsplit) ; } flops += MULTSUB_FLOPS * nz ; for (i = 0 ; i < n ; i++) @@ -321,7 +324,7 @@ for (p = Ap [i] ; p < p2 ; p++) { /* Z [Ai [p]] -= Ax [p] * xi ; */ - ASSIGN (aij, Ax, Az, p, split) ; + ASSIGN (aij, Ax, Az, p, AXsplit) ; MULT_SUB (Z [Ai [p]], aij, xi) ; } } @@ -390,7 +393,7 @@ for (i = 0 ; i < n ; i++) { /* W [i] = B [i] ; */ - ASSIGN (W [i], Bx, Bz, i, split) ; + ASSIGN (W [i], Bx, Bz, i, Bsplit) ; Z2 [i] = 0. ; } flops += (MULT_FLOPS + DECREMENT_FLOPS + ABS_FLOPS + 1) * nz ; @@ -403,7 +406,7 @@ i = Ai [p] ; /* axx = Ax [p] * xj ; */ - ASSIGN (aij, Ax, Az, p, split) ; + ASSIGN (aij, Ax, Az, p, AXsplit) ; MULT (axx, aij, xj) ; /* W [i] -= axx ; */ @@ -493,7 +496,7 @@ /* yi += ABS (Ax [p]) * Rs [Ai [p]] ; */ /* note that abs (aij) is the same as * abs (conj (aij)) */ - ASSIGN (aij, Ax, Az, p, split) ; + ASSIGN (aij, Ax, Az, p, AXsplit) ; ABS (d, aij) ; yi += (d * Rs [Ai [p]]) ; } @@ -513,7 +516,7 @@ /* yi += ABS (Ax [p]) / Rs [Ai [p]] ; */ /* note that abs (aij) is the same as * abs (conj (aij)) */ - ASSIGN (aij, Ax, Az, p, split) ; + ASSIGN (aij, Ax, Az, p, AXsplit) ; ABS (d, aij) ; yi += (d / Rs [Ai [p]]) ; } @@ -534,7 +537,7 @@ /* yi += ABS (Ax [p]) ; */ /* note that abs (aij) is the same as * abs (conj (aij)) */ - ASSIGN (aij, Ax, Az, p, split) ; + ASSIGN (aij, Ax, Az, p, AXsplit) ; ABS (d, aij) ; yi += d ; } @@ -546,7 +549,7 @@ for (i = 0 ; i < n ; i++) { /* B2 [i] = ABS (B [i]) ; */ - ASSIGN (bi, Bx, Bz, i, split) ; + ASSIGN (bi, Bx, Bz, i, Bsplit) ; ABS (B2 [i], bi) ; } @@ -568,7 +571,7 @@ for (i = 0 ; i < n ; i++) { /* W [i] = B [Cperm [i]] ; */ - ASSIGN (W [i], Bx, Bz, Cperm [i], split) ; + ASSIGN (W [i], Bx, Bz, Cperm [i], Bsplit) ; } } else @@ -577,7 +580,7 @@ for (i = 0 ; i < n ; i++) { /* Z [i] = B [i] ; */ - ASSIGN (Z [i], Bx, Bz, i, split) ; + ASSIGN (Z [i], Bx, Bz, i, Bsplit) ; } flops += MULTSUB_FLOPS * nz ; for (i = 0 ; i < n ; i++) @@ -587,7 +590,7 @@ for (p = Ap [i] ; p < p2 ; p++) { /* zi -= conjugate (Ax [p]) * X [Ai [p]] ; */ - ASSIGN (aij, Ax, Az, p, split) ; + ASSIGN (aij, Ax, Az, p, Bsplit) ; MULT_SUB_CONJ (zi, X [Ai [p]], aij) ; } Z [i] = zi ; @@ -696,13 +699,13 @@ for (i = 0 ; i < n ; i++) { /* wi = B [i] ; */ - ASSIGN (wi, Bx, Bz, i, split) ; + ASSIGN (wi, Bx, Bz, i, Bsplit) ; z2i = 0. ; p2 = Ap [i+1] ; for (p = Ap [i] ; p < p2 ; p++) { /* axx = conjugate (Ax [p]) * X [Ai [p]] ; */ - ASSIGN (aij, Ax, Az, p, split) ; + ASSIGN (aij, Ax, Az, p, AXsplit) ; MULT_CONJ (axx, X [Ai [p]], aij) ; /* wi -= axx ; */ @@ -766,7 +769,7 @@ /* yi += ABS (Ax [p]) * Rs [Ai [p]] ; */ /* note that A.' is the array transpose, * so no conjugate */ - ASSIGN (aij, Ax, Az, p, split) ; + ASSIGN (aij, Ax, Az, p, AXsplit) ; ABS (d, aij) ; yi += (d * Rs [Ai [p]]) ; } @@ -786,7 +789,7 @@ /* yi += ABS (Ax [p]) / Rs [Ai [p]] ; */ /* note that A.' is the array transpose, * so no conjugate */ - ASSIGN (aij, Ax, Az, p, split) ; + ASSIGN (aij, Ax, Az, p, AXsplit) ; ABS (d, aij) ; yi += (d / Rs [Ai [p]]) ; } @@ -807,7 +810,7 @@ /* yi += ABS (Ax [p]) */ /* note that A.' is the array transpose, * so no conjugate */ - ASSIGN (aij, Ax, Az, p, split) ; + ASSIGN (aij, Ax, Az, p, AXsplit) ; ABS (d, aij) ; yi += d ; } @@ -819,7 +822,7 @@ for (i = 0 ; i < n ; i++) { /* B2 [i] = ABS (B [i]) ; */ - ASSIGN (bi, Bx, Bz, i, split) ; + ASSIGN (bi, Bx, Bz, i, Bsplit) ; ABS (B2 [i], bi) ; } @@ -841,7 +844,7 @@ for (i = 0 ; i < n ; i++) { /* W [i] = B [Cperm [i]] ; */ - ASSIGN (W [i], Bx, Bz, Cperm [i], split) ; + ASSIGN (W [i], Bx, Bz, Cperm [i], Bsplit) ; } } else @@ -850,7 +853,7 @@ for (i = 0 ; i < n ; i++) { /* Z [i] = B [i] ; */ - ASSIGN (Z [i], Bx, Bz, i, split) ; + ASSIGN (Z [i], Bx, Bz, i, Bsplit) ; } flops += MULTSUB_FLOPS * nz ; for (i = 0 ; i < n ; i++) @@ -860,7 +863,7 @@ for (p = Ap [i] ; p < p2 ; p++) { /* zi -= Ax [p] * X [Ai [p]] ; */ - ASSIGN (aij, Ax, Az, p, split) ; + ASSIGN (aij, Ax, Az, p, AXsplit) ; MULT_SUB (zi, aij, X [Ai [p]]) ; } Z [i] = zi ; @@ -969,13 +972,13 @@ for (i = 0 ; i < n ; i++) { /* wi = B [i] ; */ - ASSIGN (wi, Bx, Bz, i, split) ; + ASSIGN (wi, Bx, Bz, i, Bsplit) ; z2i = 0. ; p2 = Ap [i+1] ; for (p = Ap [i] ; p < p2 ; p++) { /* axx = Ax [p] * X [Ai [p]] ; */ - ASSIGN (aij, Ax, Az, p, split) ; + ASSIGN (aij, Ax, Az, p, AXsplit) ; MULT (axx, aij, X [Ai [p]]) ; /* wi -= axx ; */ @@ -1011,7 +1014,7 @@ for (i = 0 ; i < n ; i++) { /* X [i] = B [Rperm [i]] ; */ - ASSIGN (X [i], Bx, Bz, Rperm [i], split) ; + ASSIGN (X [i], Bx, Bz, Rperm [i], Bsplit) ; } flops = UMF_lsolve (Numeric, X, Pattern) ; status = UMFPACK_OK ; @@ -1027,7 +1030,7 @@ for (i = 0 ; i < n ; i++) { /* X [i] = B [i] ; */ - ASSIGN (X [i], Bx, Bz, i, split) ; + ASSIGN (X [i], Bx, Bz, i, Bsplit) ; } flops = UMF_lsolve (Numeric, X, Pattern) ; status = UMFPACK_OK ; @@ -1043,7 +1046,7 @@ for (i = 0 ; i < n ; i++) { /* W [i] = B [i] ; */ - ASSIGN (W [i], Bx, Bz, i, split) ; + ASSIGN (W [i], Bx, Bz, i, Bsplit) ; } flops = UMF_lhsolve (Numeric, W, Pattern) ; for (i = 0 ; i < n ; i++) @@ -1063,7 +1066,7 @@ for (i = 0 ; i < n ; i++) { /* W [i] = B [i] ; */ - ASSIGN (W [i], Bx, Bz, i, split) ; + ASSIGN (W [i], Bx, Bz, i, Bsplit) ; } flops = UMF_ltsolve (Numeric, W, Pattern) ; for (i = 0 ; i < n ; i++) @@ -1083,7 +1086,7 @@ for (i = 0 ; i < n ; i++) { /* X [i] = B [i] ; */ - ASSIGN (X [i], Bx, Bz, i, split) ; + ASSIGN (X [i], Bx, Bz, i, Bsplit) ; } flops = UMF_lhsolve (Numeric, X, Pattern) ; status = UMFPACK_OK ; @@ -1099,7 +1102,7 @@ for (i = 0 ; i < n ; i++) { /* X [i] = B [i] ; */ - ASSIGN (X [i], Bx, Bz, i, split) ; + ASSIGN (X [i], Bx, Bz, i, Bsplit) ; } flops = UMF_ltsolve (Numeric, X, Pattern) ; status = UMFPACK_OK ; @@ -1115,7 +1118,7 @@ for (i = 0 ; i < n ; i++) { /* W [i] = B [i] ; */ - ASSIGN (W [i], Bx, Bz, i, split) ; + ASSIGN (W [i], Bx, Bz, i, Bsplit) ; } flops = UMF_usolve (Numeric, W, Pattern) ; for (i = 0 ; i < n ; i++) @@ -1134,7 +1137,7 @@ for (i = 0 ; i < n ; i++) { /* X [i] = B [i] ; */ - ASSIGN (X [i], Bx, Bz, i, split) ; + ASSIGN (X [i], Bx, Bz, i, Bsplit) ; } flops = UMF_usolve (Numeric, X, Pattern) ; @@ -1149,7 +1152,7 @@ for (i = 0 ; i < n ; i++) { /* X [i] = B [Cperm [i]] ; */ - ASSIGN (X [i], Bx, Bz, Cperm [i], split) ; + ASSIGN (X [i], Bx, Bz, Cperm [i], Bsplit) ; } flops = UMF_uhsolve (Numeric, X, Pattern) ; @@ -1164,7 +1167,7 @@ for (i = 0 ; i < n ; i++) { /* X [i] = B [Cperm [i]] ; */ - ASSIGN (X [i], Bx, Bz, Cperm [i], split) ; + ASSIGN (X [i], Bx, Bz, Cperm [i], Bsplit) ; } flops = UMF_utsolve (Numeric, X, Pattern) ; @@ -1179,7 +1182,7 @@ for (i = 0 ; i < n ; i++) { /* X [i] = B [i] ; */ - ASSIGN (X [i], Bx, Bz, i, split) ; + ASSIGN (X [i], Bx, Bz, i, Bsplit) ; } flops = UMF_uhsolve (Numeric, X, Pattern) ; @@ -1194,7 +1197,7 @@ for (i = 0 ; i < n ; i++) { /* X [i] = B [i] ; */ - ASSIGN (X [i], Bx, Bz, i, split) ; + ASSIGN (X [i], Bx, Bz, i, Bsplit) ; } flops = UMF_utsolve (Numeric, X, Pattern) ; @@ -1206,7 +1209,7 @@ #ifdef COMPLEX /* copy the solution back, from Entry X [ ] to double Xx [ ] and Xz [ ] */ - if (split) + if (AXsplit) { for (i = 0 ; i < n ; i++) {