octave-nkf: liboctave/UMFPACK.patch comparison

comparison 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

comparison

equal deleted inserted replaced

-:9f3299378193
+:57077d0ddc8e
+--- 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++)
+	{

Mercurial > octave-nkf

comparison liboctave/UMFPACK.patch @ 5164:57077d0ddc8e