# HG changeset patch # User jwe # Date 1110848336 0 # Node ID dbeafbc0ff64e0aad4f04636a242d5ef6894c3fe # Parent 49de4e624020bb71ce361167ad2f6f5b6d012009 [project @ 2005-03-15 00:58:55 by jwe] diff -r 49de4e624020 -r dbeafbc0ff64 ChangeLog --- a/ChangeLog Mon Mar 14 20:51:59 2005 +0000 +++ b/ChangeLog Tue Mar 15 00:58:56 2005 +0000 @@ -1,7 +1,19 @@ +2005-03-14 John W. Eaton + + * configure.in: Check for umfpack/umfpack.h instead of just umfpack.h. + +2004-06-22 David Bateman + + * configure.in: Check for UMFPACK library and header files. + +2005-03-14 John W. Eaton + + * configure.in: Also print a warning if HDF5 library is not found. + 2005-03-10 John W. Eaton * mkoctfile.in: Accept -R DIR. - +p 2005-03-09 John W. Eaton * examples/Makefile.in (bin-dist): Delete target. diff -r 49de4e624020 -r dbeafbc0ff64 configure.in --- a/configure.in Mon Mar 14 20:51:59 2005 +0000 +++ b/configure.in Tue Mar 15 00:58:56 2005 +0000 @@ -29,7 +29,7 @@ EXTERN_CXXFLAGS="$CXXFLAGS" AC_INIT -AC_REVISION($Revision: 1.464 $) +AC_REVISION($Revision: 1.465 $) AC_PREREQ(2.57) AC_CONFIG_SRCDIR([src/octave.cc]) AC_CONFIG_HEADER(config.h) @@ -396,6 +396,12 @@ AC_DEFINE(HAVE_H5GGET_NUM_OBJS, 1, [Define if HDF5 has H5Gget_num_objs.])])])])]) fi +if $WITH_HDF5; then + true +else + warn_hdf5="HDF5 library not found. Octave will not be able to save or load HDF5 data files." +fi + # Checks for FFTW header and library. # subdirectories of libcruft to build if they aren't found on the system: @@ -415,7 +421,7 @@ with_fftw3=no AC_CHECK_HEADER(fftw3.h, [have_fftw3_header=yes; break]) if test "$have_fftw3_header" = yes; then - AC_CHECK_LIB(fftw3, fftw_plan_dft_1d, [FFTW_LIBS="-lfftw3"; with_fftw3=yes]) + AC_CHECK_LIB(fftw3, fftw_plan_dft_1d, [FFTW_LIBS="-lfftw3"; with_fftw3=yes]) fi fi @@ -424,6 +430,44 @@ AC_DEFINE(HAVE_FFTW3, 1, [Define if the FFTW3 library is used.]) fi +# Check for UMFPACK library. + +UMFPACK_LIBS= +AC_SUBST(UMFPACK_LIBS) + +AC_ARG_WITH(umfpack, + [ --without-umfpack don't use UMFPACK, disable some sparse functionality], + with_umfpack=$withval, with_umfpack=yes) + +if test "$with_umfpack" = "yes"; then + have_umfpack_header=no + with_umfpack=no + AC_CHECK_HEADER(umfpack/umfpack.h, [have_umfpack_header=yes; break]) + if test "$have_umfpack_header" = yes; then + AC_CHECK_LIB(amd, amd_postorder, [ + LIBS="-lamd $LIBS" + AC_CHECK_LIB(umfpack, umfpack_zi_get_determinant, [ + LIBS="-lumfpack $LIBS" + UMFPACK_LIBS="-lumfpack"; with_umfpack=yes])]) + if test "$with_umfpack" = yes; then + # For now the code needed for this is not in umfpack, will add + # a test later that will probably have to be based on version + # numbers as there is no interface changes that are visible at + # compile time. + with_umfpack_split=no + fi + fi +fi + +if test "$with_umfpack" = yes; then + AC_DEFINE(HAVE_UMFPACK, 1, [Define if the UMFPACK library is used.]) + if test "$with_umfpack_split" = yes; then + AC_DEFINE(UMFPACK_SEPARATE_SPLIT, 1, [Define if the UMFPACK Complex solver allow matrix and RHS to be split independently]) + fi +else + warn_umfpack="UMFPACK not found. This will result in some lack of functionality for sparse matrices." +fi + WITH_MPI=true AC_ARG_WITH(mpi, [ --without-mpi don't use MPI], @@ -1479,6 +1523,7 @@ Fortran libraries: $FLIBS BLAS libraries: $BLAS_LIBS FFTW libraries: $FFTW_LIBS + UMFPACK libraries: $UMFPACK_LIBS HDF5 libraries: $HDF5_LIBS MPI libraries: $MPI_LIBS LIBS: $LIBS @@ -1564,6 +1609,16 @@ warn_msg_printed=true fi +if test -n "$warn_umfpack"; then + AC_MSG_WARN($warn_umfpack) + warn_msg_printed=true +fi + +if test -n "$warn_hdf5"; then + AC_MSG_WARN($warn_hdf5) + warn_msg_printed=true +fi + if test -n "$warn_gnuplot"; then ## If you change this text, be sure to also change the corresponding diff -r 49de4e624020 -r dbeafbc0ff64 liboctave/CSparse.cc --- a/liboctave/CSparse.cc Mon Mar 14 20:51:59 2005 +0000 +++ b/liboctave/CSparse.cc Tue Mar 15 00:58:56 2005 +0000 @@ -40,10 +40,12 @@ #include "oct-spparms.h" #include "SparseCmplxLU.h" +#ifdef HAVE_UMFPACK // External UMFPACK functions in C extern "C" { -#include "umfpack.h" +#include } +#endif // Fortran functions we call. extern "C" @@ -627,6 +629,7 @@ SparseComplexMatrix::determinant (int& err, double& rcond, int calc_cond) const { ComplexDET retval; +#ifdef HAVE_UMFPACK int nr = rows (); int nc = cols (); @@ -742,6 +745,9 @@ } } } +#else + (*current_liboctave_error_handler) ("UMFPACK not installed"); +#endif return retval; } @@ -4446,6 +4452,7 @@ void *Numeric; err = 0; +#ifdef HAVE_UMFPACK // Setup the control parameters Control = Matrix (UMFPACK_CONTROL, 1); double *control = Control.fortran_vec (); @@ -4542,6 +4549,9 @@ if (err != 0) umfpack_zi_free_numeric (&Numeric); +#else + (*current_liboctave_error_handler) ("UMFPACK not installed"); +#endif return Numeric; } @@ -4580,6 +4590,7 @@ if (typ == SparseType::Full) { +#ifdef HAVE_UMFPACK Matrix Control, Info; void *Numeric = factorize (err, rcond, Control, Info, sing_handler); @@ -4594,16 +4605,20 @@ const int *Ap = cidx (); const int *Ai = ridx (); const Complex *Ax = data (); +#ifdef UMFPACK_SEPARATE_SPLIT const double *Bx = b.fortran_vec (); OCTAVE_LOCAL_BUFFER (double, Bz, b_nr); for (int i = 0; i < b_nr; i++) Bz[i] = 0.; - +#else + OCTAVE_LOCAL_BUFFER (Complex, Bz, b_nr); +#endif retval.resize (b_nr, b_nc); Complex *Xx = retval.fortran_vec (); for (int j = 0, iidx = 0; j < b_nc; j++, iidx += b_nr) { +#ifdef UMFPACK_SEPARATE_SPLIT status = umfpack_zi_solve (UMFPACK_A, Ap, Ai, X_CAST (const double *, Ax), NULL, @@ -4611,6 +4626,20 @@ NULL, &Bx[iidx], Bz, Numeric, control, info); +#else + for (int i = 0; i < b_nr; i++) + Bz[i] = b.elem (i, j); + + status = umfpack_zi_solve (UMFPACK_A, Ap, Ai, + X_CAST (const double *, Ax), + NULL, + X_CAST (double *, &Xx[iidx]), + NULL, + X_CAST (const double *, Bz), + NULL, Numeric, + control, info); +#endif + if (status < 0) { (*current_liboctave_error_handler) @@ -4647,6 +4676,9 @@ umfpack_zi_free_numeric (&Numeric); } +#else + (*current_liboctave_error_handler) ("UMFPACK not installed"); +#endif } else if (typ != SparseType::Hermitian) (*current_liboctave_error_handler) ("incorrect matrix type"); @@ -4689,6 +4721,7 @@ if (typ == SparseType::Full) { +#ifdef HAVE_UMFPACK Matrix Control, Info; void *Numeric = factorize (err, rcond, Control, Info, sing_handler); @@ -4703,10 +4736,14 @@ const int *Ai = ridx (); const Complex *Ax = data (); +#ifdef UMFPACK_SEPARATE_SPLIT OCTAVE_LOCAL_BUFFER (double, Bx, b_nr); OCTAVE_LOCAL_BUFFER (double, Bz, b_nr); for (int i = 0; i < b_nr; i++) Bz[i] = 0.; +#else + OCTAVE_LOCAL_BUFFER (Complex, Bz, b_nr); +#endif // Take a first guess that the number of non-zero terms // will be as many as in b @@ -4720,6 +4757,7 @@ for (int j = 0; j < b_nc; j++) { +#ifdef UMFPACK_SEPARATE_SPLIT for (int i = 0; i < b_nr; i++) Bx[i] = b.elem (i, j); @@ -4729,6 +4767,18 @@ X_CAST (double *, Xx), NULL, Bx, Bz, Numeric, control, info); +#else + for (int i = 0; i < b_nr; i++) + Bz[i] = b.elem (i, j); + + status = umfpack_zi_solve (UMFPACK_A, Ap, Ai, + X_CAST (const double *, Ax), + NULL, + X_CAST (double *, Xx), NULL, + X_CAST (double *, Bz), NULL, + Numeric, control, + info); +#endif if (status < 0) { (*current_liboctave_error_handler) @@ -4786,6 +4836,9 @@ umfpack_zi_free_numeric (&Numeric); } +#else + (*current_liboctave_error_handler) ("UMFPACK not installed"); +#endif } else if (typ != SparseType::Hermitian) (*current_liboctave_error_handler) ("incorrect matrix type"); @@ -4828,6 +4881,7 @@ if (typ == SparseType::Full) { +#ifdef HAVE_UMFPACK Matrix Control, Info; void *Numeric = factorize (err, rcond, Control, Info, sing_handler); @@ -4891,6 +4945,9 @@ umfpack_zi_free_numeric (&Numeric); } +#else + (*current_liboctave_error_handler) ("UMFPACK not installed"); +#endif } else if (typ != SparseType::Hermitian) (*current_liboctave_error_handler) ("incorrect matrix type"); @@ -4933,6 +4990,7 @@ if (typ == SparseType::Full) { +#ifdef HAVE_UMFPACK Matrix Control, Info; void *Numeric = factorize (err, rcond, Control, Info, sing_handler); @@ -5026,6 +5084,9 @@ umfpack_zi_free_numeric (&Numeric); } +#else + (*current_liboctave_error_handler) ("UMFPACK not installed"); +#endif } else if (typ != SparseType::Hermitian) (*current_liboctave_error_handler) ("incorrect matrix type"); diff -r 49de4e624020 -r dbeafbc0ff64 liboctave/ChangeLog --- a/liboctave/ChangeLog Mon Mar 14 20:51:59 2005 +0000 +++ b/liboctave/ChangeLog Tue Mar 15 00:58:56 2005 +0000 @@ -1,3 +1,27 @@ +2005-03-14 John W. Eaton + + * Makefile.in (DISTFILES): Don't include $(UMFPACK_EXTRAS). + (DISTDIRS): Don't include UMFPACK. + (LIBOCTAVE_OBJECTS): Don't include $(UMFPACK_OBJ). + (UMFPACK_SPECIAL_1, UMFPACK_SPECIAL): No need for special include + flags for these files. + Don't include include $(srcdir)/UMFPACK.files. + Don't include include $(srcdir)/UMFPACK.rules. + + * UMFPACK.README, UMFPACK.files, UMFPACK.patch, UMFPACK.rules: + Delete files. + * UMFPACK: Delete directory tree. + + * dSparse.cc: Include instead of just "umfpack.h". + * CSparse.cc: Likewise. + * SparsedbleLU.cc: Likewise. + * SparseCmplxLU.cc: Likewise. + +2005-03-14 David Bateman + + * CSParse.cc, SparseCmplxLU.cc, SparsedbleLU.cc, dSparse.cc: + Allow compilation to succeed if UMFPACK is not available. + 2005-03-09 John W. Eaton * Makefile.in (bin-dist): Delete target. diff -r 49de4e624020 -r dbeafbc0ff64 liboctave/Makefile.in --- a/liboctave/Makefile.in Mon Mar 14 20:51:59 2005 +0000 +++ b/liboctave/Makefile.in Tue Mar 15 00:58:56 2005 +0000 @@ -24,8 +24,6 @@ include $(srcdir)/COLAMD.files -include $(srcdir)/UMFPACK.files - MATRIX_INC := Array.h Array2.h Array3.h ArrayN.h DiagArray2.h \ Array-flags.h Array-util.h ArrayN-idx.h MArray-defs.h \ MArray.h MArray2.h MDiagArray2.h Matrix.h MArrayN.h \ @@ -155,10 +153,10 @@ DISTFILES := Makefile.in ChangeLog mk-ops.awk mx-ops vx-ops \ sparse-mk-ops.awk sparse-mx-ops \ $(SOURCES) $(INCLUDES) $(EXTRAS) $(OPTS_INC_DATA) \ - $(COLAMD_EXTRAS) $(UMFPACK_EXTRAS) + $(COLAMD_EXTRAS) # Complete directory trees to distribute. -DISTDIRS := COLAMD UMFPACK +DISTDIRS := COLAMD MAKEDEPS_2 := $(SOURCES) $(COLAMD_SRC) MAKEDEPS_1 := $(patsubst %.cc, %.d, $(MAKEDEPS_2)) @@ -166,7 +164,6 @@ LIBOCTAVE_OBJECTS := \ $(COLAMD_OBJ) \ - $(UMFPACK_OBJ) \ $(LIBOCTAVE_CXX_SOURCES:.cc=.o) \ $(LIBOCTAVE_C_SOURCES:.c=.o) \ $(LIBOCT_READLINE_CXX_SOURCES:.cc=.o) \ @@ -182,20 +179,11 @@ endif endif -UMFPACK_SPECIAL_1 := CSparse.o dSparse.o SparseCmplxLU.o SparsedbleLU.o -UMFPACK_SPECIAL := \ - $(UMFPACK_SPECIAL_1) \ - $(addprefix pic/, $(UMFPACK_SPECIAL_1)) \ - $(UMFPACK_SPECIAL_1:.o=.d) -$(UMFPACK_SPECIAL): INCFLAGS += $(UMFPACK_INCFLAGS) - all: libraries .PHONY: all objects: $(OBJECTS) -include $(srcdir)/UMFPACK.rules - stmp-pic: pic @if [ -f stmp-pic ]; then \ true; \ diff -r 49de4e624020 -r dbeafbc0ff64 liboctave/SparseCmplxLU.cc --- a/liboctave/SparseCmplxLU.cc Mon Mar 14 20:51:59 2005 +0000 +++ b/liboctave/SparseCmplxLU.cc Tue Mar 15 00:58:56 2005 +0000 @@ -37,14 +37,17 @@ template class sparse_base_lu ; +#ifdef HAVE_UMFPACK // Include the UMFPACK functions extern "C" { -#include "umfpack.h" +#include } +#endif SparseComplexLU::SparseComplexLU (const SparseComplexMatrix& a, double piv_thres) { +#ifdef HAVE_UMFPACK int nr = a.rows (); int nc = a.cols (); @@ -210,12 +213,16 @@ } } } +#else + (*current_liboctave_error_handler) ("UMFPACK not installed"); +#endif } SparseComplexLU::SparseComplexLU (const SparseComplexMatrix& a, const ColumnVector& Qinit, double piv_thres, bool FixedQ) { +#ifdef HAVE_UMFPACK int nr = a.rows (); int nc = a.cols (); @@ -395,6 +402,9 @@ } } } +#else + (*current_liboctave_error_handler) ("UMFPACK not installed"); +#endif } /* diff -r 49de4e624020 -r dbeafbc0ff64 liboctave/SparsedbleLU.cc --- a/liboctave/SparsedbleLU.cc Mon Mar 14 20:51:59 2005 +0000 +++ b/liboctave/SparsedbleLU.cc Tue Mar 15 00:58:56 2005 +0000 @@ -37,13 +37,16 @@ template class sparse_base_lu ; +#ifdef HAVE_UMFPACK // Include the UMFPACK functions extern "C" { -#include "umfpack.h" +#include } +#endif SparseLU::SparseLU (const SparseMatrix& a, double piv_thres) { +#ifdef HAVE_UMFPACK int nr = a.rows (); int nc = a.cols (); @@ -204,11 +207,15 @@ } } } +#else + (*current_liboctave_error_handler) ("UMFPACK not installed"); +#endif } SparseLU::SparseLU (const SparseMatrix& a, const ColumnVector& Qinit, double piv_thres, bool FixedQ) { +#ifdef HAVE_UMFPACK int nr = a.rows (); int nc = a.cols (); @@ -382,6 +389,9 @@ } } } +#else + (*current_liboctave_error_handler) ("UMFPACK not installed"); +#endif } /* diff -r 49de4e624020 -r dbeafbc0ff64 liboctave/UMFPACK.README --- a/liboctave/UMFPACK.README Mon Mar 14 20:51:59 2005 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,47 +0,0 @@ -This directory contains a modified copy of UMFPACKv4.4 in the -directory UMFPACK. The complete changes from the original version may -be found in the file UMFPACK.patch. - -UMFPACK was written by Timothy A. Davis (davis@cise.ufl.edu), -University of Florida. - -UMFPACK includes a modified version of COLAMD V2.0, by Stefan -I. Larimore and Timothy A. Davis, University of Florida. The COLAMD -algorithm was developed in collaboration with John Gilbert, Xerox Palo -Alto Research Center, and Esmond Ng, Lawrence Berkeley National -Laboratory. - -UMFPACK also includes AMD, by Timothy A. Davis, Patrick R. Amestoy, -and Iain S. Duff. - -UMFPACK Version 2.2.1 (MA38 in the Harwell Subroutine Library) is -co-authored with Iain S. Duff, Rutherford Appleton Laboratory. - -The copyright and license information is: - - UMFPACK Version 4.4 (Jan. 28, 2005), Copyright (c) 2005 by Timothy A. - Davis. All Rights Reserved. - - 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. - -John W. Eaton -jwe@bevo.che.wisc.edu -University of Wisconsin-Madison -Department of Chemical & Biological Engineering - -Wed Feb 01 22:15:20 2005 diff -r 49de4e624020 -r dbeafbc0ff64 liboctave/UMFPACK.files --- a/liboctave/UMFPACK.files Mon Mar 14 20:51:59 2005 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,79 +0,0 @@ - -# 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 - -UMFPACK_BASE := \ - $(subst umf_, umf_o_, $(UMFINT)) \ - $(subst umf_, umf_od_, $(UMF)) \ - $(subst umfpack_, umfpack_od_, $(USER)) \ - $(subst umf_, umf_oz_, $(UMF)) \ - $(subst umfpack_, umfpack_oz_, $(USER)) \ - $(subst amd_, amd_o_, $(AMD)) \ - $(subst umfpack_, umfpack_o_, $(GENERIC)) - -UMFPACK_OBJ := $(addsuffix .o, $(UMFPACK_BASE)) - -UMFPACK_PICOBJ := $(addprefix pic/, $(UMFPACK_OBJ)) - -UMFPACK_INCFLAGS := \ - -I$(top_srcdir)/liboctave/UMFPACK/UMFPACK/Include \ - -I$(top_srcdir)/liboctave/UMFPACK/UMFPACK/Source \ - -I$(top_srcdir)/liboctave/UMFPACK/AMD/Include \ - -I$(top_srcdir)/liboctave/UMFPACK/AMD/Source - -$(UMFPACK_OBJ) $(UMFPACK_PICOBJ): INCFLAGS += $(UMFPACK_INCFLAGS) - -UMFPACK_EXTRAS := UMFPACK.files UMFPACK.rules UMFPACK.patch UMFPACK.README - diff -r 49de4e624020 -r dbeafbc0ff64 liboctave/UMFPACK.patch --- a/liboctave/UMFPACK.patch Mon Mar 14 20:51:59 2005 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,4825 +0,0 @@ ---- 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 -+#include -+ -+#include -+#include -+#include -+#include -+#include -+ -+#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 -+#include -+#include -+ -+#include -+#include -+#include -+#include -+#include -+ -+#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 (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', ... -+ '\\extern\\lib\\win32\\lcc\\ directory, where 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', ... -+ '\\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 ' 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 -+#include -+ -+#include -+#include -+#include -+#include -+#include -+ -+#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 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++) - { diff -r 49de4e624020 -r dbeafbc0ff64 liboctave/UMFPACK.rules --- a/liboctave/UMFPACK.rules Mon Mar 14 20:51:59 2005 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,183 +0,0 @@ - -#---------------------------------------- -# integer-only routines (no real/complex): -#---------------------------------------- - -amd_o_%.o: $(srcdir)/UMFPACK/AMD/Source/amd_%.c - $(CC) -c $(CPPFLAGS) $(ALL_CFLAGS) -DDINT $< -o $@ - - -umf_o_%.o: $(srcdir)/UMFPACK/UMFPACK/Source/umf_%.c - $(CC) -c $(CPPFLAGS) $(ALL_CFLAGS) -DDINT $< -o $@ - -#---------------------------------------- -# Double precision, int version, for OCTAVE -#---------------------------------------- - -umf_od_%.o: $(srcdir)/UMFPACK/UMFPACK/Source/umf_%.c - $(CC) -c $(CPPFLAGS) $(ALL_CFLAGS) -DDINT $< -o $@ - -umf_od_%hsolve.o: $(srcdir)/UMFPACK/UMFPACK/Source/umf_%tsolve.c - $(CC) -c $(CPPFLAGS) $(ALL_CFLAGS) -DDINT -DCONJUGATE_SOLVE $< -o $@ - -umf_od_triplet_map_x.o: $(srcdir)/UMFPACK/UMFPACK/Source/umf_triplet.c - $(CC) -c $(CPPFLAGS) $(ALL_CFLAGS) -DDINT -DDO_MAP -DDO_VALUES $< -o $@ - -umf_od_triplet_map_nox.o: $(srcdir)/UMFPACK/UMFPACK/Source/umf_triplet.c - $(CC) -c $(CPPFLAGS) $(ALL_CFLAGS) -DDINT -DDO_MAP $< -o $@ - -umf_od_triplet_nomap_x.o: $(srcdir)/UMFPACK/UMFPACK/Source/umf_triplet.c - $(CC) -c $(CPPFLAGS) $(ALL_CFLAGS) -DDINT -DDO_VALUES $< -o $@ - -umf_od_triplet_nomap_nox.o: $(srcdir)/UMFPACK/UMFPACK/Source/umf_triplet.c - $(CC) -c $(CPPFLAGS) $(ALL_CFLAGS) -DDINT $< -o $@ - -umf_od_assemble_fixq.o: $(srcdir)/UMFPACK/UMFPACK/Source/umf_assemble.c - $(CC) -c $(CPPFLAGS) $(ALL_CFLAGS) -DDINT -DFIXQ $< -o $@ - -umf_od_store_lu_drop.o: $(srcdir)/UMFPACK/UMFPACK/Source/umf_store_lu.c - $(CC) -c $(CPPFLAGS) $(ALL_CFLAGS) -DDINT -DDROP $< -o $@ - -umfpack_od_wsolve.o: $(srcdir)/UMFPACK/UMFPACK/Source/umfpack_solve.c - $(CC) -c $(CPPFLAGS) $(ALL_CFLAGS) -DDINT -DWSOLVE $< -o $@ - -umfpack_od_%.o: $(srcdir)/UMFPACK/UMFPACK/Source/umfpack_%.c - $(CC) -c $(CPPFLAGS) $(ALL_CFLAGS) -DDINT $< -o $@ - -#---------------------------------------- -# Complex double precision, int version, for OCTAVE -#---------------------------------------- - -umf_oz_%.o: $(srcdir)/UMFPACK/UMFPACK/Source/umf_%.c - $(CC) -c $(CPPFLAGS) $(ALL_CFLAGS) -DZINT $< -o $@ - -umf_oz_%hsolve.o: $(srcdir)/UMFPACK/UMFPACK/Source/umf_%tsolve.c - $(CC) -c $(CPPFLAGS) $(ALL_CFLAGS) -DZINT -DCONJUGATE_SOLVE $< -o $@ - -umf_oz_triplet_map_x.o: $(srcdir)/UMFPACK/UMFPACK/Source/umf_triplet.c - $(CC) -c $(CPPFLAGS) $(ALL_CFLAGS) -DZINT -DDO_MAP -DDO_VALUES $< -o $@ - -umf_oz_triplet_map_nox.o: $(srcdir)/UMFPACK/UMFPACK/Source/umf_triplet.c - $(CC) -c $(CPPFLAGS) $(ALL_CFLAGS) -DZINT -DDO_MAP $< -o $@ - -umf_oz_triplet_nomap_x.o: $(srcdir)/UMFPACK/UMFPACK/Source/umf_triplet.c - $(CC) -c $(CPPFLAGS) $(ALL_CFLAGS) -DZINT -DDO_VALUES $< -o $@ - -umf_oz_triplet_nomap_nox.o: $(srcdir)/UMFPACK/UMFPACK/Source/umf_triplet.c - $(CC) -c $(CPPFLAGS) $(ALL_CFLAGS) -DZINT $< -o $@ - -umf_oz_assemble_fixq.o: $(srcdir)/UMFPACK/UMFPACK/Source/umf_assemble.c - $(CC) -c $(CPPFLAGS) $(ALL_CFLAGS) -DZINT -DFIXQ $< -o $@ - -umf_oz_store_lu_drop.o: $(srcdir)/UMFPACK/UMFPACK/Source/umf_store_lu.c - $(CC) -c $(CPPFLAGS) $(ALL_CFLAGS) -DZINT -DDROP $< -o $@ - -umfpack_oz_wsolve.o: $(srcdir)/UMFPACK/UMFPACK/Source/umfpack_solve.c - $(CC) -c $(CPPFLAGS) $(ALL_CFLAGS) -DZINT -DWSOLVE $< -o $@ - -umfpack_oz_%.o: $(srcdir)/UMFPACK/UMFPACK/Source/umfpack_%.c - $(CC) -c $(CPPFLAGS) $(ALL_CFLAGS) -DZINT $< -o $@ - -#---------------------------------------- -# Generic routines for OCTAVE -#---------------------------------------- - -umfpack_o_timer.o: $(srcdir)/UMFPACK/UMFPACK/Source/umfpack_timer.c - $(CC) -c $(CPPFLAGS) $(ALL_CFLAGS) $< -o $@ - -umfpack_o_tictoc.o: $(srcdir)/UMFPACK/UMFPACK/Source/umfpack_tictoc.c - $(CC) -c $(CPPFLAGS) $(ALL_CFLAGS) $< -o $@ - -######################################################################## -# PIC rules -######################################################################## - -#---------------------------------------- -# integer-only routines (no real/complex): -#---------------------------------------- - -pic/amd_o_%.o: $(srcdir)/UMFPACK/AMD/Source/amd_%.c - $(CC) -c $(CPPFLAGS) $(CPICFLAG) $(ALL_CFLAGS) -DDINT $< -o $@ - - -pic/umf_o_%.o: $(srcdir)/UMFPACK/UMFPACK/Source/umf_%.c - $(CC) -c $(CPPFLAGS) $(CPICFLAG) $(ALL_CFLAGS) -DDINT $< -o $@ - -#---------------------------------------- -# Double precision, int version, for OCTAVE -#---------------------------------------- - -pic/umf_od_%.o: $(srcdir)/UMFPACK/UMFPACK/Source/umf_%.c - $(CC) -c $(CPPFLAGS) $(CPICFLAG) $(ALL_CFLAGS) -DDINT $< -o $@ - -pic/umf_od_%hsolve.o: $(srcdir)/UMFPACK/UMFPACK/Source/umf_%tsolve.c - $(CC) -c $(CPPFLAGS) $(CPICFLAG) $(ALL_CFLAGS) -DDINT -DCONJUGATE_SOLVE $< -o $@ - -pic/umf_od_triplet_map_x.o: $(srcdir)/UMFPACK/UMFPACK/Source/umf_triplet.c - $(CC) -c $(CPPFLAGS) $(CPICFLAG) $(ALL_CFLAGS) -DDINT -DDO_MAP -DDO_VALUES $< -o $@ - -pic/umf_od_triplet_map_nox.o: $(srcdir)/UMFPACK/UMFPACK/Source/umf_triplet.c - $(CC) -c $(CPPFLAGS) $(CPICFLAG) $(ALL_CFLAGS) -DDINT -DDO_MAP $< -o $@ - -pic/umf_od_triplet_nomap_x.o: $(srcdir)/UMFPACK/UMFPACK/Source/umf_triplet.c - $(CC) -c $(CPPFLAGS) $(CPICFLAG) $(ALL_CFLAGS) -DDINT -DDO_VALUES $< -o $@ - -pic/umf_od_triplet_nomap_nox.o: $(srcdir)/UMFPACK/UMFPACK/Source/umf_triplet.c - $(CC) -c $(CPPFLAGS) $(CPICFLAG) $(ALL_CFLAGS) -DDINT $< -o $@ - -pic/umf_od_assemble_fixq.o: $(srcdir)/UMFPACK/UMFPACK/Source/umf_assemble.c - $(CC) -c $(CPPFLAGS) $(CPICFLAG) $(ALL_CFLAGS) -DDINT -DFIXQ $< -o $@ - -pic/umf_od_store_lu_drop.o: $(srcdir)/UMFPACK/UMFPACK/Source/umf_store_lu.c - $(CC) -c $(CPPFLAGS) $(CPICFLAG) $(ALL_CFLAGS) -DDINT -DDROP $< -o $@ - -pic/umfpack_od_wsolve.o: $(srcdir)/UMFPACK/UMFPACK/Source/umfpack_solve.c - $(CC) -c $(CPPFLAGS) $(CPICFLAG) $(ALL_CFLAGS) -DDINT -DWSOLVE $< -o $@ - -pic/umfpack_od_%.o: $(srcdir)/UMFPACK/UMFPACK/Source/umfpack_%.c - $(CC) -c $(CPPFLAGS) $(CPICFLAG) $(ALL_CFLAGS) -DDINT $< -o $@ - -#---------------------------------------- -# Complex double precision, int version, for OCTAVE -#---------------------------------------- - -pic/umf_oz_%.o: $(srcdir)/UMFPACK/UMFPACK/Source/umf_%.c - $(CC) -c $(CPPFLAGS) $(CPICFLAG) $(ALL_CFLAGS) -DZINT $< -o $@ - -pic/umf_oz_%hsolve.o: $(srcdir)/UMFPACK/UMFPACK/Source/umf_%tsolve.c - $(CC) -c $(CPPFLAGS) $(CPICFLAG) $(ALL_CFLAGS) -DZINT -DCONJUGATE_SOLVE $< -o $@ - -pic/umf_oz_triplet_map_x.o: $(srcdir)/UMFPACK/UMFPACK/Source/umf_triplet.c - $(CC) -c $(CPPFLAGS) $(CPICFLAG) $(ALL_CFLAGS) -DZINT -DDO_MAP -DDO_VALUES $< -o $@ - -pic/umf_oz_triplet_map_nox.o: $(srcdir)/UMFPACK/UMFPACK/Source/umf_triplet.c - $(CC) -c $(CPPFLAGS) $(CPICFLAG) $(ALL_CFLAGS) -DZINT -DDO_MAP $< -o $@ - -pic/umf_oz_triplet_nomap_x.o: $(srcdir)/UMFPACK/UMFPACK/Source/umf_triplet.c - $(CC) -c $(CPPFLAGS) $(CPICFLAG) $(ALL_CFLAGS) -DZINT -DDO_VALUES $< -o $@ - -pic/umf_oz_triplet_nomap_nox.o: $(srcdir)/UMFPACK/UMFPACK/Source/umf_triplet.c - $(CC) -c $(CPPFLAGS) $(CPICFLAG) $(ALL_CFLAGS) -DZINT $< -o $@ - -pic/umf_oz_assemble_fixq.o: $(srcdir)/UMFPACK/UMFPACK/Source/umf_assemble.c - $(CC) -c $(CPPFLAGS) $(CPICFLAG) $(ALL_CFLAGS) -DZINT -DFIXQ $< -o $@ - -pic/umf_oz_store_lu_drop.o: $(srcdir)/UMFPACK/UMFPACK/Source/umf_store_lu.c - $(CC) -c $(CPPFLAGS) $(CPICFLAG) $(ALL_CFLAGS) -DZINT -DDROP $< -o $@ - -pic/umfpack_oz_wsolve.o: $(srcdir)/UMFPACK/UMFPACK/Source/umfpack_solve.c - $(CC) -c $(CPPFLAGS) $(CPICFLAG) $(ALL_CFLAGS) -DZINT -DWSOLVE $< -o $@ - -pic/umfpack_oz_%.o: $(srcdir)/UMFPACK/UMFPACK/Source/umfpack_%.c - $(CC) -c $(CPPFLAGS) $(CPICFLAG) $(ALL_CFLAGS) -DZINT $< -o $@ - -#---------------------------------------- -# Generic routines for OCTAVE -#---------------------------------------- - -pic/umfpack_o_timer.o: $(srcdir)/UMFPACK/UMFPACK/Source/umfpack_timer.c - $(CC) -c $(CPPFLAGS) $(CPICFLAG) $(ALL_CFLAGS) $< -o $@ - -pic/umfpack_o_tictoc.o: $(srcdir)/UMFPACK/UMFPACK/Source/umfpack_tictoc.c - $(CC) -c $(CPPFLAGS) $(CPICFLAG) $(ALL_CFLAGS) $< -o $@ - diff -r 49de4e624020 -r dbeafbc0ff64 liboctave/dSparse.cc --- a/liboctave/dSparse.cc Mon Mar 14 20:51:59 2005 +0000 +++ b/liboctave/dSparse.cc Tue Mar 15 00:58:56 2005 +0000 @@ -41,10 +41,12 @@ #include "SparsedbleLU.h" #include "SparseType.h" +#ifdef HAVE_UMFPACK // External UMFPACK functions in C extern "C" { -#include "umfpack.h" -} +#include +} +#endif // Fortran functions we call. extern "C" @@ -714,6 +716,7 @@ { DET retval; +#ifdef HAVE_UMFPACK int nr = rows (); int nc = cols (); @@ -821,6 +824,9 @@ } } } +#else + (*current_liboctave_error_handler) ("UMFPACK not installed"); +#endif return retval; } @@ -4671,6 +4677,7 @@ void *Numeric; err = 0; +#ifdef HAVE_UMFPACK // Setup the control parameters Control = Matrix (UMFPACK_CONTROL, 1); double *control = Control.fortran_vec (); @@ -4766,6 +4773,10 @@ if (err != 0) umfpack_di_free_numeric (&Numeric); +#else + (*current_liboctave_error_handler) ("UMFPACK not installed"); +#endif + return Numeric; } @@ -4803,6 +4814,7 @@ if (typ == SparseType::Full) { +#ifdef HAVE_UMFPACK Matrix Control, Info; void *Numeric = factorize (err, rcond, Control, Info, sing_handler); @@ -4862,6 +4874,9 @@ umfpack_di_free_numeric (&Numeric); } +#else + (*current_liboctave_error_handler) ("UMFPACK not installed"); +#endif } else if (typ != SparseType::Hermitian) (*current_liboctave_error_handler) ("incorrect matrix type"); @@ -4903,6 +4918,7 @@ if (typ == SparseType::Full) { +#ifdef HAVE_UMFPACK Matrix Control, Info; void *Numeric = factorize (err, rcond, Control, Info, sing_handler); @@ -4994,6 +5010,9 @@ umfpack_di_free_numeric (&Numeric); } +#else + (*current_liboctave_error_handler) ("UMFPACK not installed"); +#endif } else if (typ != SparseType::Hermitian) (*current_liboctave_error_handler) ("incorrect matrix type"); @@ -5035,6 +5054,7 @@ if (typ == SparseType::Full) { +#ifdef HAVE_UMFPACK Matrix Control, Info; void *Numeric = factorize (err, rcond, Control, Info, sing_handler); @@ -5113,6 +5133,9 @@ umfpack_di_free_numeric (&Numeric); } +#else + (*current_liboctave_error_handler) ("UMFPACK not installed"); +#endif } else if (typ != SparseType::Hermitian) (*current_liboctave_error_handler) ("incorrect matrix type"); @@ -5155,6 +5178,7 @@ if (typ == SparseType::Full) { +#ifdef HAVE_UMFPACK Matrix Control, Info; void *Numeric = factorize (err, rcond, Control, Info, sing_handler); @@ -5256,6 +5280,9 @@ umfpack_di_free_numeric (&Numeric); } +#else + (*current_liboctave_error_handler) ("UMFPACK not installed"); +#endif } else if (typ != SparseType::Hermitian) (*current_liboctave_error_handler) ("incorrect matrix type");