Mercurial > hg > octave-max
changeset 5506:b4cfbb0ec8c4
[project @ 2005-10-23 19:09:32 by dbateman]
line wrap: on
line diff
--- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,9 @@ +2005-10-23 David Bateman <dbateman@free.fr> + + * configure.in (OCTAVE_UMFPACK_SEPERATE_SPLIT): Check for metis + separately. + * PROJECTS: Remove completed sparse matrix tasks. + 2005-10-17 Paul Kienzle <pkienzle@users.sf.net> * octave.test/system/system.exp: rmdir no longer prints a
--- a/PROJECTS +++ b/PROJECTS @@ -107,17 +107,12 @@ Sparse Matrices: --------------- - * Sparse Cholesky factorization for Fchol function and also for the - in built polymorphic solvers. - * QR factorization functions, also for use in lssolve functions. Write svds function based on this. Write sprank function based on svds. * Once dmperm is implemented, perhaps use the technique to detect permuted triangular matrices with permutations in both rows and cols. - * Sparse inverse function, based on Andy's code from octave-forge. - * Implement fourth argument to the sprand and sprandn that the leading brand implements. @@ -163,17 +158,12 @@ * Port the sparse testing code into the DejaGNU testing code. Add tests for spkron, matrix_type, luinc.. - * Treat the dispatching of the splu, spdet, functions, etc within octave - itself. This either means that classes need implementing or at the - minimum the octave-forge dispatch function is implemented within octave - * Other missing Functions - eigs - - dmperm Tim Davis is apparently working on something + - dmperm Use Tim Davis' BTF code when available - symmmd Superseded by symamd - colmmd Superseded by colamd - sprandsym - - symbfact - etreeplot - treeplot - gplot
--- a/configure.in +++ b/configure.in @@ -29,7 +29,7 @@ EXTERN_CXXFLAGS="$CXXFLAGS" AC_INIT -AC_REVISION($Revision: 1.486 $) +AC_REVISION($Revision: 1.487 $) AC_PREREQ(2.57) AC_CONFIG_SRCDIR([src/octave.cc]) AC_CONFIG_HEADER(config.h) @@ -823,6 +823,8 @@ AC_SUBST(CHOLMOD_LIBS) CHOLMOD_INCLUDE=cholmod.h AC_SUBST(CHOLMOD_INCLUDE) +METIS_INCLUDE=metis.h +AC_SUBST(METIS_INCLUDE) AC_ARG_WITH(cholmod, [ --without-cholmod don't use CHOLMOD, disable some sparse functionality], @@ -832,18 +834,35 @@ test "$with_ccolamd" = "yes" && test "$with_amd" = "yes"; then with_cholmod=no ACX_CHECK_HEADER_IN_DIRS(cholmod.h, [umfpack ufsparse],[ - AC_CHECK_LIB(metis, METIS_NodeND, [ + if test x"$acx_include_dir" != x; then + CHOLMOD_INCLUDE=$acx_include_dir/cholmod.h + fi + + ACX_CHECK_HEADER_IN_DIRS(metis.h, [metis umfpack ufsparse],[ + AC_CHECK_LIB(metis, METIS_NodeND, with_metis=yes, with_metis=no)], + with_metis=no) + + if test "$with_metis" = yes; then + if test x"$acx_include_dir" != x; then + METIS_INCLUDE=$acx_include_dir/metis.h + fi + AC_DEFINE(HAVE_METIS, 1, [Define if the METIS library is used.]) AC_CHECK_LIB(cholmod, cholmod_start, [CHOLMOD_LIBS="-lcholmod -lmetis"; with_cholmod=yes], [ - AC_CHECK_LIB(cholmod_start, cholmod, [CHOLMOD_LIBS="-lcholmod -cblas -lmetis"; + AC_CHECK_LIB(cholmod_start, cholmod, + [CHOLMOD_LIBS="-lcholmod -cblas -lmetis"; with_cholmod=yes], [], + AMD_LIBS $COLAMD_LIBS $CCOLAMD_LIBS $BLAS_LIBS $FLIBS -lmetis)], + $AMD_LIBS $COLAMD_LIBS $CCOLAMD_LIBS $BLAS_LIBS $FLIBS -lmetis) + else + AC_CHECK_LIB(cholmod, cholmod_start, [CHOLMOD_LIBS="-lcholmod"; + with_cholmod=yes], [ + AC_CHECK_LIB(cholmod_start, cholmod, [CHOLMOD_LIBS="-lcholmod -cblas"; with_cholmod=yes], [], - AMD_LIBS $COLAMD_LIBS $CCOLAMD_LIBS $BLAS_LIBS $FLIBS -lmetis)], - $AMD_LIBS $COLAMD_LIBS $CCOLAMD_LIBS $BLAS_LIBS $FLIBS -lmetis)]) + AMD_LIBS $COLAMD_LIBS $CCOLAMD_LIBS $BLAS_LIBS $FLIBS)], + $AMD_LIBS $COLAMD_LIBS $CCOLAMD_LIBS $BLAS_LIBS $FLIBS) + fi if test "$with_cholmod" = yes; then - if test x"$acx_include_dir" != x; then - CHOLMOD_INCLUDE=$acx_include_dir/cholmod.h - fi AC_DEFINE(HAVE_CHOLMOD, 1, [Define if the CHOLMOD library is used.]) else warn_cholmod="CHOLMOD not found. This will result in some lack of functionality for sparse matrices."
--- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -1,3 +1,8 @@ +2005-10-23 David Bateman <dbateman@free.fr> + + * sparse.txi: Updates for new ufsparse licensing, new functions and + various typos. + 2005-09-19 Rafael Laboissiere <rafael@debian.org> * interpreter/octave-config.1: Use bold instead of italics to
--- a/doc/interpreter/octave.texi +++ b/doc/interpreter/octave.texi @@ -388,7 +388,6 @@ * Sparse Linear Algebra:: * Iterative Techniques:: * Oct-Files:: -* License:: * Function Reference:: Quadrature
--- a/doc/interpreter/sparse.txi +++ b/doc/interpreter/sparse.txi @@ -11,7 +11,6 @@ * Sparse Linear Algebra:: Linear Algebra on Sparse Matrices * Iterative Techniques:: Iterative Techniques applied to Sparse Matrices * Oct-Files:: Using Sparse Matrices in Oct-files -* License:: Licensing of Third Party Software * Function Reference:: Documentation from the Specific Sparse Functions @end menu @@ -23,11 +22,11 @@ the speed of the computer and its available memory place limitation on the problem size. -There are many classes mathematical problems which give rise to +There are many classes of mathematical problems which give rise to matrices, where a large number of the elements are zero. In this case it makes sense to have a special matrix type to handle this class of problems where only the non-zero elements of the matrix are -stored. Not only done this reduce the amount of memory to store the +stored. Not only does this reduce the amount of memory to store the matrix, but it also means that operations on this type of matrix can take advantage of the a-priori knowledge of the positions of the non-zero elements to accelerate their calculations. @@ -69,7 +68,7 @@ of the non-zero elements of the matrix must equally be stored. An obvious way to do this is by storing the elements of the matrix as -triplets, with two elements being the of the position in the array +triplets, with two elements being their position in the array (rows and column) and the third being the data itself. This is conceptually easy to grasp, but requires more storage than is strictly needed. @@ -122,14 +121,14 @@ @example @group - @var{cidx} = [0, 2, 2, 4] + @var{cidx} = [0, 1, 2, 2, 4] @var{ridx} = [0, 0, 1, 2] @var{data} = [1, 2, 3, 4] @end group @end example Note that this is the representation of these elements with the first row -and column assumed to start as zero. Thus the number of elements in the +and column assumed to start at zero. Thus the number of elements in the @var{i}-th column is given by @code{@var{cidx} (@var{i} + 1) - @var{cidx} (@var{i})}. @@ -141,7 +140,7 @@ A further constraint on the sparse matrix storage used by Octave is that all elements in the rows are stored in increasing order of their row -index. This makes certain operations, later be faster. However, it imposes +index, which makes certain operations faster. However, it imposes the need to sort the elements on the creation of sparse matrices. Having dis-ordered elements is potentially an advantage in that it makes operations such as concatenating two sparse matrices together easier and faster, however @@ -158,7 +157,7 @@ @dfn{speye}, @dfn{sprand}, @dfn{spdiag}, etc. @item Constructed from matrices or vectors The function @dfn{sparse} allows a sparse matrix to be constructed from -three vectors representing the row, column and data. ALternatively, the +three vectors representing the row, column and data. Alternatively, the function @dfn{spconvert} uses a three column matrix format to allow easy importation of data from elsewhere. @item Created and then filled @@ -177,20 +176,18 @@ Another typical sparse matrix that is often needed is a random distribution of random elements. The functions @dfn{sprand} and @dfn{sprandn} perform -this for uniform and random distributions of elements. They have exactly -the same calling convention as @code{sprand (@var{r}, @var{c}, @var{d})}, -which creates an @var{r}-by-@var{c} sparse matrix a density of filled +this for uniform and normal random distributions of elements. They have exactly +the same calling convention, where @code{sprand (@var{r}, @var{c}, @var{d})}, +creates an @var{r}-by-@var{c} sparse matrix with a density of filled elements of @var{d}. -Other functions of interest that directly creates a sparse matrix, are +Other functions of interest that directly creates a sparse matrices, are @dfn{spdiag} or its generalization @dfn{spdiags}, that can take the definition of the diagonals of the matrix and create the sparse matrix that corresponds to this. For example -@c FIXME, when spdiag is overloaded to diag the below won't work. - @example -s = spdiag (randn(1,n), -1); +s = diag (sparse(randn(1,n)), -1); @end example creates a sparse (@var{n}+1)-by-(@var{n}+1) sparse matrix with a single @@ -222,8 +219,8 @@ creates an @var{r}-by-@var{c} sparse matrix with a random distribution of 2 elements per row. The elements of the vectors do not need to be sorted in -any particular order as Octave will store them prior to storing the -data. However, per sorting the data will make teh creation of the sparse +any particular order as Octave will sort them prior to storing the +data. However, pre-sorting the data will make the creation of the sparse matrix faster. The function @dfn{spconvert} takes a three or four column real matrix. @@ -255,11 +252,11 @@ It should be noted, that due to the way that the Octave assignment functions are written that the assignment will reallocate -the memory used by the sparse matrix at each iteration. Therefore the -@dfn{spalloc} function ignores the @var{nz} argument and does not -preassign the memory for the matrix. Therefore, it is vitally +the memory used by the sparse matrix at each iteration of the above loop. +Therefore the @dfn{spalloc} function ignores the @var{nz} argument and +does not preassign the memory for the matrix. Therefore, it is vitally important that code using to above structure should be as vectorized -as possible to minimize the number of assignments and reduce the +much as possible to minimize the number of assignments and reduce the number of memory allocations. The above problem can be avoided in oct-files. However, the @@ -285,10 +282,10 @@ @node ReturnType, MathConsiderations, Functions, Operators and Functions @subsubsection The Return Types of Operators and Functions -The two basic reason to use sparse matrices are to reduce the memory +The two basic reasons to use sparse matrices are to reduce the memory usage and to not have to do calculations on zero elements. The two are closely related in that the computation time on a sparse matrix operator -or function is roughly linear with the numberof non-zero elements. +or function is roughly linear with the number of non-zero elements. Therefore, there is a certain density of non-zero elements of a matrix where it no longer makes sense to store it as a sparse matrix, but rather @@ -373,7 +370,7 @@ A particular problem of sparse matrices comes about due to the fact that as the zeros are not stored, the sign-bit of these zeros is equally not -stored... In certain cases the sign-bit of zero is important. For example +stored. In certain cases the sign-bit of zero is important. For example @example a = 0 ./ [-1, 1; 1, -1]; @@ -391,6 +388,12 @@ efficient, and so the user is warned that calculations where the sign-bit of zero is important must not be done using sparse matrices. +In general any function or operator used on a sparse matrix will result +in a sparse matrix with the same or a larger number of non-zero elements +than the original matrix. This is particularly true for the important +case of sparse matrix factorizations. + + Also discuss issues of fill-in. Discuss symamd etc, and mention randperm that is included elsewhere in the docs... @@ -415,7 +418,7 @@ Octave includes a poly-morphic solver for sparse matrices, where the exact solver used to factorize the matrix, depends on the properties -of the sparse matrix, itself. The cost of determining the matrix type +of the sparse matrix itself. The cost of determining the matrix type is small relative to the cost of factorizing the matrix itself, but in any case the matrix type is cached once it is calculated, so that it is not re-determined each time it is used in a linear equation. @@ -462,10 +465,7 @@ or backward subsitution, and goto 9 @item If the matrix is hermitian with a real positive diagonal, attempt -sparse Cholesky factorization. - -FIXME: Detection of positive definite matrices written and tested, but - Cholesky factorization isn't yet written +sparse Cholesky factorization using CHOLMOD. @item If the sparse Cholesky factorization failed or the matrix is not hermitian with a real positive diagonal, factorize using UMFPACK. @@ -497,14 +497,14 @@ function. This overcomes the cost of discovering the type of the matrix. However, it should be noted incorrectly identifying the type of the matrix will lead to unpredictable results, and so @code{matrix_type} should be -use dwith care. +used with care. @node Iterative Techniques, Oct-Files, Sparse Linear Algebra, Sparse Matrices @section Iterative Techniques applied to sparse matrices WRITE ME -@node Oct-Files, License, Iterative Techniques, Sparse Matrices +@node Oct-Files, Function Reference, Iterative Techniques, Sparse Matrices @section Using Sparse Matrices in Oct-files An oct-file is a means of writing an Octave function in a compilable @@ -762,7 +762,7 @@ @} sm.cidx(j+1) = ii; @} - sm.maybe_mutate (); // If don't know a-priori the final no of nz. + sm.maybe_compress (); // If don't know a-priori the final no of nz. @end example which is probably the most efficient means of creating the sparse matrix. @@ -837,72 +837,7 @@ The conversion to an octave-value is handled by the sparse @code{octave_value} constructors, and so no special care is needed. -@node License, Function Reference, Oct-Files, Sparse Matrices -@section Licensing of Third Party Software - -There are several third party software packages used by the Octave -sparse matrix. - -@table @asis -@item COLAMD -is used for the @dfn{colamd} and @dfn{symamd} functions. - -@table @asis -@item Authors -The authors of the code itself are Stefan I. Larimore and Timothy A. -Davis (davis@@cise.ufl.edu), University of Florida. The algorithm was -developed in collaboration with John Gilbert, Xerox PARC, and Esmond -Ng, Oak Ridge National Laboratory. - -@item License -Copyright @copyright{} 1998-2003 by the University of Florida. -All Rights Reserved. - -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, copy, modify, and/or distribute -this program, provided that the Copyright, this License, and the -Availability of the original version is retained on all copies and made -accessible to the end-user of any code or package that includes COLAMD -or any modified version of COLAMD. - -@item Availability -@url{http://www.cise.ufl.edu/research/sparse/colamd/} -@end table - -@item UMFPACK -is used in various places with the sparse types, including the -LU decomposition and solvers. - -@table @asis -@item License -UMFPACK Version 4.4 (Jan. 28, 2005), Copyright @copyright{} 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. - -@item Availability -@url{http://www.cise.ufl.edu/research/sparse/umfpack/} -@end table -@end table - -@node Function Reference, , License, Sparse Matrices +@node Function Reference, , Oct-Files, Sparse Matrices @section Function Reference @iftex @@ -965,10 +900,14 @@ @end table @subsubsection Sparse matrix reordering @table @asis +@item ccolamd +Constrained column approximate minimum degree permutation. @item colamd Column approximate minimum degree permutation. @item colperm Returns the column permutations such that the columns of `S (:, P)' are ordered in terms of increase number of non-zero elements. +@item csymamd +For a symmetric positive definite matrix S, returns the permutation vector p such that `S (P, P)' tends to have a sparser Cholesky factor than S. @item dmperm Perform a Deulmage-Mendelsohn permutation on the sparse matrix S. @item symamd @@ -988,12 +927,22 @@ Identify the matrix type or mark a matrix as a particular type. @item normest @emph{Not implemented} +@item spchol +Compute the Cholesky factor, R, of the symmetric positive definite. +@item spcholinv +Use the Cholesky factorization to compute the inverse of the +sparse symmetric positive definite matrix A. +@item spchol2inv +Invert a sparse symmetric, positive definite square matrix from its +Cholesky decomposition, U. @item spdet Compute the determinant of sparse matrix A using UMFPACK. @item spinv Compute the inverse of the square matrix A. @item spkron Form the kronecker product of two sparse matrices. +@item splchol +Compute the Cholesky factor, L, of the symmetric positive definite. @item splu Compute the LU decomposition of the sparse matrix A, using subroutines from UMFPACK. @item sprank @@ -1054,18 +1003,20 @@ Compute atan (Y / X) for corresponding sparse matrix elements of Y and X. @item spdiag Return a diagonal matrix with the sparse vector V on diagonal K. -@item spreshape -Return a sparse matrix with M rows and N columns whose elements are taken from the sparse matrix A. @end table @subsection Functions Alphabetically @end iftex @menu +* ccolamd:: Constrained column approximate minimum degree permutation. * colamd:: Column approximate minimum degree permutation. * colperm:: Returns the column permutations such that the columns of `S (:, P)' are ordered in terms of increase number of non-zero elements. +* csymamd:: For a symmetric positive definite matrix S, returns the + permutation vector p such that `S (P, P)' tends to have a + sparser Cholesky factor than S. * dmperm:: Perfrom a Deulmage-Mendelsohn permutation on the sparse matrix S. * etree:: Returns the elimination tree for the matrix S. @@ -1086,6 +1037,12 @@ * sparse:: SPARSE: create a sparse matrix * spatan2:: Compute atan (Y / X) for corresponding sparse matrix elements of Y and X. +* spchol:: Compute the Cholesky factor, R, of the symmetric + positive definite. +* spcholinv:: Use the Cholesky factorization to compute the inverse of the + sparse symmetric positive definite matrix A. +* spchol2inv:: Invert a sparse symmetric, positive definite square matrix + from its Cholesky decomposition, U. * spconvert:: This function converts for a simple sparse matrix format easily produced by other programs into Octave's internal sparse format. @@ -1101,6 +1058,8 @@ a sparse matrix with the same structure as X. * spinv:: Compute the inverse of the square matrix A. * spkron:: Form the kronecker product of two sparse matrices. +* splchol:: Compute the Cholesky factor, L, of the symmetric positive + definite. * splu:: Compute the LU decomposition of the sparse matrix A, using subroutines from UMFPACK. * spmax:: For a vector argument, return the maximum value. @@ -1111,8 +1070,6 @@ * spprod:: Product of elements along dimension DIM. * sprand:: Generate a random sparse matrix. * sprandn:: Generate a random sparse matrix. -* spreshape:: Return a sparse matrix with M rows and N columns whose - elements are taken from the sparse matrix A. * spstats:: Return the stats for the non-zero elements of the sparse matrix S COUNT is the number of non-zeros in each column, MEAN is the mean of the non-zeros in each column, and VAR @@ -1129,17 +1086,27 @@ matrix S. @end menu -@node colamd, colperm, , Function Reference +@node colamd, ccolamd, , Function Reference @subsubsection colamd @DOCSTRING(colamd) -@node colperm, dmperm, colamd, Function Reference +@node ccolamd, colperm, colamd, Function Reference +@subsubsection ccolamd + +@DOCSTRING(ccolamd) + +@node colperm, csymamd, ccolamd, Function Reference @subsubsection colperm @DOCSTRING(colperm) -@node dmperm, etree, colperm, Function Reference +@node csymamd, dmperm, colperm, Function Reference +@subsubsection csymamd + +@DOCSTRING(csymamd) + +@node dmperm, etree, csymamd, Function Reference @subsubsection dmperm @DOCSTRING(dmperm) @@ -1194,12 +1161,27 @@ @DOCSTRING(sparse) -@node spatan2, spconvert, sparse, Function Reference +@node spatan2, spchol, sparse, Function Reference @subsubsection spatan2 @DOCSTRING(spatan2) -@node spconvert, spcumprod, spatan2, Function Reference +@node spchol, spcholinv, spatan2, Function Reference +@subsubsection spchol + +@DOCSTRING(spchol) + +@node spcholinv, spchol2inv, spchol, Function Reference +@subsubsection spcholinv + +@DOCSTRING(spcholinv) + +@node spchol2inv, spconvert, spcholinv, Function Reference +@subsubsection spchol2inv + +@DOCSTRING(spchol2inv) + +@node spconvert, spcumprod, spchol2inv, Function Reference @subsubsection spconvert @DOCSTRING(spconvert) @@ -1249,12 +1231,17 @@ @DOCSTRING(spinv) -@node spkron, splu, spinv, Function Reference +@node spkron, splchol, spinv, Function Reference @subsubsection spkron @DOCSTRING(spkron) -@node splu, spmax, spkron, Function Reference +@node splchol, splu, spkron, Function Reference +@subsubsection splchol + +@DOCSTRING(splchol) + +@node splu, spmax, splchol, Function Reference @subsubsection splu @DOCSTRING(splu) @@ -1289,17 +1276,12 @@ @DOCSTRING(sprand) -@node sprandn, spreshape, sprand, Function Reference +@node sprandn, spstats, sprand, Function Reference @subsubsection sprandn @DOCSTRING(sprandn) -@node spreshape, spstats, sprandn, Function Reference -@subsubsection spreshape - -@DOCSTRING(spreshape) - -@node spstats, spsum, spreshape, Function Reference +@node spstats, spsum, sprandn, Function Reference @subsubsection spstats @DOCSTRING(spstats)
--- a/liboctave/CSparse.cc +++ b/liboctave/CSparse.cc @@ -41,6 +41,8 @@ #include "oct-spparms.h" #include "SparseCmplxLU.h" #include "oct-sparse.h" +#include "sparse-util.h" +#include "SparseCmplxCHOL.h" // Fortran functions we call. extern "C" @@ -585,24 +587,378 @@ { octave_idx_type info; double rcond; - return inverse (info, rcond, 0, 0); + SparseType mattype (*this); + return inverse (mattype, info, rcond, 0, 0); +} + +SparseComplexMatrix +SparseComplexMatrix::inverse (SparseType& mattype) const +{ + octave_idx_type info; + double rcond; + return inverse (mattype, info, rcond, 0, 0); } SparseComplexMatrix -SparseComplexMatrix::inverse (octave_idx_type& info) const +SparseComplexMatrix::inverse (SparseType& mattype, octave_idx_type& info) const { double rcond; - return inverse (info, rcond, 0, 0); + return inverse (mattype, info, rcond, 0, 0); +} + +SparseComplexMatrix +SparseComplexMatrix::dinverse (SparseType &mattyp, octave_idx_type& info, + double& rcond, const bool force, + const bool calccond) const +{ + SparseComplexMatrix retval; + + octave_idx_type nr = rows (); + octave_idx_type nc = cols (); + info = 0; + + if (nr == 0 || nc == 0 || nr != nc) + (*current_liboctave_error_handler) ("inverse requires square matrix"); + else + { + // Print spparms("spumoni") info if requested + int typ = mattyp.type (); + mattyp.info (); + + if (typ == SparseType::Diagonal || + typ == SparseType::Permuted_Diagonal) + { + if (typ == SparseType::Permuted_Diagonal) + retval = transpose(); + else + retval = *this; + + // Force make_unique to be called + Complex *v = retval.data(); + + if (calccond) + { + double dmax = 0., dmin = octave_Inf; + for (octave_idx_type i = 0; i < nr; i++) + { + double tmp = std::abs(v[i]); + if (tmp > dmax) + dmax = tmp; + if (tmp < dmin) + dmin = tmp; + } + rcond = dmin / dmax; + } + + for (octave_idx_type i = 0; i < nr; i++) + v[i] = 1.0 / v[i]; + } + else + (*current_liboctave_error_handler) ("incorrect matrix type"); + } + + return retval; +} + +SparseComplexMatrix +SparseComplexMatrix::tinverse (SparseType &mattyp, octave_idx_type& info, + double& rcond, const bool force, + const bool calccond) const +{ + SparseComplexMatrix retval; + + octave_idx_type nr = rows (); + octave_idx_type nc = cols (); + info = 0; + + if (nr == 0 || nc == 0 || nr != nc) + (*current_liboctave_error_handler) ("inverse requires square matrix"); + else + { + // Print spparms("spumoni") info if requested + int typ = mattyp.type (); + mattyp.info (); + + if (typ == SparseType::Upper || typ == SparseType::Permuted_Upper || + typ == SparseType::Lower || typ == SparseType::Permuted_Lower) + { + double anorm = 0.; + double ainvnorm = 0.; + + if (calccond) + { + // Calculate the 1-norm of matrix for rcond calculation + for (octave_idx_type j = 0; j < nr; j++) + { + double atmp = 0.; + for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) + atmp += std::abs(data(i)); + if (atmp > anorm) + anorm = atmp; + } + } + + if (typ == SparseType::Upper || typ == SparseType::Lower) + { + octave_idx_type nz = nnz(); + octave_idx_type cx = 0; + octave_idx_type nz2 = nz; + retval = SparseComplexMatrix (nr, nc, nz2); + + for (octave_idx_type i = 0; i < nr; i++) + { + OCTAVE_QUIT; + // place the 1 in the identity position + octave_idx_type cx_colstart = cx; + + if (cx == nz2) + { + nz2 *= 2; + retval.change_capacity (nz2); + } + + retval.xcidx(i) = cx; + retval.xridx(cx) = i; + retval.xdata(cx) = 1.0; + cx++; + + // iterate accross columns of input matrix + for (octave_idx_type j = i+1; j < nr; j++) + { + Complex v = 0.; + // iterate to calculate sum + octave_idx_type colXp = retval.xcidx(i); + octave_idx_type colUp = cidx(j); + octave_idx_type rpX, rpU; + do + { + OCTAVE_QUIT; + rpX = retval.xridx(colXp); + rpU = ridx(colUp); + + if (rpX < rpU) + colXp++; + else if (rpX > rpU) + colUp++; + else + { + v -= retval.xdata(colXp) * data(colUp); + colXp++; + colUp++; + } + } while ((rpX<j) && (rpU<j) && + (colXp<cx) && (colUp<nz)); + + // get A(m,m) + colUp = cidx(j+1) - 1; + Complex pivot = data(colUp); + if (pivot == 0.) + (*current_liboctave_error_handler) + ("division by zero"); + + if (v != 0.) + { + if (cx == nz2) + { + nz2 *= 2; + retval.change_capacity (nz2); + } + + retval.xridx(cx) = j; + retval.xdata(cx) = v / pivot; + cx++; + } + } + + // get A(m,m) + octave_idx_type colUp = cidx(i+1) - 1; + Complex pivot = data(colUp); + if (pivot == 0.) + (*current_liboctave_error_handler) ("division by zero"); + + if (pivot != 1.0) + for (octave_idx_type j = cx_colstart; j < cx; j++) + retval.xdata(j) /= pivot; + } + retval.xcidx(nr) = cx; + retval.maybe_compress (); + } + else + { + octave_idx_type nz = nnz(); + octave_idx_type cx = 0; + octave_idx_type nz2 = nz; + retval = SparseComplexMatrix (nr, nc, nz2); + + OCTAVE_LOCAL_BUFFER (Complex, work, nr); + OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nr); + + octave_idx_type *perm = mattyp.triangular_perm(); + if (typ == SparseType::Permuted_Upper) + { + for (octave_idx_type i = 0; i < nr; i++) + rperm[perm[i]] = i; + } + else + { + for (octave_idx_type i = 0; i < nr; i++) + rperm[i] = perm[i]; + for (octave_idx_type i = 0; i < nr; i++) + perm[rperm[i]] = i; + } + + for (octave_idx_type i = 0; i < nr; i++) + { + OCTAVE_QUIT; + octave_idx_type iidx = rperm[i]; + + for (octave_idx_type j = 0; j < nr; j++) + work[j] = 0.; + + // place the 1 in the identity position + work[iidx] = 1.0; + + // iterate accross columns of input matrix + for (octave_idx_type j = iidx+1; j < nr; j++) + { + Complex v = 0.; + octave_idx_type jidx = perm[j]; + // iterate to calculate sum + for (octave_idx_type k = cidx(jidx); + k < cidx(jidx+1); k++) + { + OCTAVE_QUIT; + v -= work[ridx(k)] * data(k); + } + + // get A(m,m) + Complex pivot = data(cidx(jidx+1) - 1); + if (pivot == 0.) + (*current_liboctave_error_handler) + ("division by zero"); + + work[j] = v / pivot; + } + + // get A(m,m) + octave_idx_type colUp = cidx(perm[iidx]+1) - 1; + Complex pivot = data(colUp); + if (pivot == 0.) + (*current_liboctave_error_handler) + ("division by zero"); + + octave_idx_type new_cx = cx; + for (octave_idx_type j = iidx; j < nr; j++) + if (work[j] != 0.0) + { + new_cx++; + if (pivot != 1.0) + work[j] /= pivot; + } + + if (cx < new_cx) + { + nz2 = (2*nz2 < new_cx ? new_cx : 2*nz2); + retval.change_capacity (nz2); + } + + retval.xcidx(i) = cx; + for (octave_idx_type j = iidx; j < nr; j++) + if (work[j] != 0.) + { + retval.xridx(cx) = j; + retval.xdata(cx++) = work[j]; + } + } + + retval.xcidx(nr) = cx; + retval.maybe_compress (); + } + + if (calccond) + { + // Calculate the 1-norm of inverse matrix for rcond calculation + for (octave_idx_type j = 0; j < nr; j++) + { + double atmp = 0.; + for (octave_idx_type i = retval.cidx(j); + i < retval.cidx(j+1); i++) + atmp += std::abs(retval.data(i)); + if (atmp > ainvnorm) + ainvnorm = atmp; + } + + rcond = 1. / ainvnorm / anorm; + } + } + else + (*current_liboctave_error_handler) ("incorrect matrix type"); + } + + return retval; } SparseComplexMatrix -SparseComplexMatrix::inverse (octave_idx_type& info, double& rcond, int force, - int calc_cond) const -{ - info = -1; - (*current_liboctave_error_handler) - ("SparseComplexMatrix::inverse not implemented yet"); - return SparseComplexMatrix (); +SparseComplexMatrix::inverse (SparseType& mattype, octave_idx_type& info, + double& rcond, int force, int calc_cond) const +{ + int typ = mattype.type (false); + SparseComplexMatrix ret; + + if (typ == SparseType::Unknown) + typ = mattype.type (*this); + + if (typ == SparseType::Diagonal || typ == SparseType::Permuted_Diagonal) + ret = dinverse (mattype, info, rcond, true, calc_cond); + else if (typ == SparseType::Upper || typ == SparseType::Permuted_Upper) + ret = tinverse (mattype, info, rcond, true, calc_cond).transpose(); + else if (typ == SparseType::Lower || typ == SparseType::Permuted_Lower) + ret = transpose().tinverse (mattype, info, rcond, true, calc_cond); + else if (typ != SparseType::Rectangular) + { + if (mattype.is_hermitian()) + { + SparseType tmp_typ (SparseType::Upper); + SparseComplexCHOL fact (*this, info, false); + rcond = fact.rcond(); + if (info == 0) + { + double rcond2; + SparseMatrix Q = fact.Q(); + SparseComplexMatrix InvL = fact.L().transpose(). + tinverse(tmp_typ, info, rcond2, true, false); + ret = Q * InvL.hermitian() * InvL * Q.transpose(); + } + else + { + // Matrix is either singular or not positive definite + mattype.mark_as_unsymmetric (); + typ = SparseType::Full; + } + } + + if (!mattype.is_hermitian()) + { + octave_idx_type n = rows(); + ColumnVector Qinit(n); + for (octave_idx_type i = 0; i < n; i++) + Qinit(i) = i; + + SparseType tmp_typ (SparseType::Upper); + SparseComplexLU fact (*this, Qinit, -1.0, false); + rcond = fact.rcond(); + double rcond2; + SparseComplexMatrix InvL = fact.L().transpose(). + tinverse(tmp_typ, info, rcond2, true, false); + SparseComplexMatrix InvU = fact.U(). + tinverse(tmp_typ, info, rcond2, true, false).transpose(); + ret = fact.Pc().transpose() * InvU * InvL * fact.Pr(); + } + } + else + (*current_liboctave_error_handler) ("inverse requires square matrix"); + + return ret; } ComplexDET @@ -4646,14 +5002,152 @@ if (typ == SparseType::Hermitian) { - // XXX FIXME XXX Write the cholesky solver and only fall - // through if cholesky factorization fails - +#ifdef HAVE_CHOLMOD + cholmod_common Common; + cholmod_common *cm = &Common; + + // Setup initial parameters + CHOLMOD_NAME(start) (cm); + cm->prefer_zomplex = FALSE; + + double spu = Voctave_sparse_controls.get_key ("spumoni"); + if (spu == 0.) + { + cm->print = -1; + cm->print_function = NULL; + } + else + { + cm->print = (int)spu + 2; + cm->print_function =&SparseCholPrint; + } + + cm->error_handler = &SparseCholError; + cm->complex_divide = CHOLMOD_NAME(divcomplex); + cm->hypotenuse = CHOLMOD_NAME(hypot); + +#ifdef HAVE_METIS + // METIS 4.0.1 uses malloc and free, and will terminate MATLAB if + // it runs out of memory. Use CHOLMOD's memory guard for METIS, + // which mxMalloc's a huge block of memory (and then immediately + // mxFree's it) before calling METIS + cm->metis_memory = 2.0; + +#if defined(METIS_VERSION) +#if (METIS_VERSION >= METIS_VER(4,0,2)) + // METIS 4.0.2 uses function pointers for malloc and free + METIS_malloc = cm->malloc_memory; + METIS_free = cm->free_memory; + // Turn off METIS memory guard. It is not needed, because mxMalloc + // will safely terminate the mexFunction and free any workspace + // without killing all of octave. + cm->metis_memory = 0.0; +#endif +#endif +#endif + + cm->final_ll = TRUE; + + cholmod_sparse Astore; + cholmod_sparse *A = &Astore; + double dummy; + A->nrow = nr; + A->ncol = nc; + + A->p = cidx(); + A->i = ridx(); + A->nzmax = nonzero(); + A->packed = TRUE; + A->sorted = TRUE; + A->nz = NULL; +#ifdef IDX_TYPE_LONG + A->itype = CHOLMOD_LONG; +#else + A->itype = CHOLMOD_INT; +#endif + A->dtype = CHOLMOD_DOUBLE; + A->stype = 1; + A->xtype = CHOLMOD_COMPLEX; + + if (nr < 1) + A->x = &dummy; + else + A->x = data(); + + cholmod_dense Bstore; + cholmod_dense *B = &Bstore; + B->nrow = b.rows(); + B->ncol = b.cols(); + B->d = B->nrow; + B->nzmax = B->nrow * B->ncol; + B->dtype = CHOLMOD_DOUBLE; + B->xtype = CHOLMOD_REAL; + if (nc < 1 || b.cols() < 1) + B->x = &dummy; + else + // We won't alter it, honest :-) + B->x = const_cast<double *>(b.fortran_vec()); + + cholmod_factor *L; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + L = CHOLMOD_NAME(analyze) (A, cm); + CHOLMOD_NAME(factorize) (A, L, cm); + rcond = CHOLMOD_NAME(rcond)(L, cm); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + if (rcond == 0.0) + { + // Either its indefinite or singular. Try UMFPACK + mattype.mark_as_unsymmetric (); + typ = SparseType::Full; + } + else + { + volatile double rcond_plus_one = rcond + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcond)) + { + err = -2; + + if (sing_handler) + sing_handler (rcond); + else + (*current_liboctave_error_handler) + ("SparseMatrix::solve matrix singular to machine precision, rcond = %g", + rcond); + +#ifdef HAVE_LSSOLVE + return retval; +#endif + } + + cholmod_dense *X; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + X = CHOLMOD_NAME(solve) (CHOLMOD_A, L, B, cm); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + retval.resize (b.rows (), b.cols()); + for (octave_idx_type j = 0; j < b.cols(); j++) + { + octave_idx_type jr = j * b.rows(); + for (octave_idx_type i = 0; i < b.rows(); i++) + retval.xelem(i,j) = static_cast<Complex *>(X->x)[jr + i]; + } + + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CHOLMOD_NAME(free_dense) (&X, cm); + CHOLMOD_NAME(free_factor) (&L, cm); + CHOLMOD_NAME(finish) (cm); + CHOLMOD_NAME(print_common) (" ", cm); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } +#else (*current_liboctave_warning_handler) - ("SparseMatrix::solve XXX FIXME XXX Cholesky code not done"); + ("CHOLMOD not installed"); mattype.mark_as_unsymmetric (); typ = SparseType::Full; +#endif } if (typ == SparseType::Full) @@ -4772,19 +5266,173 @@ else { // Print spparms("spumoni") info if requested - int typ = mattype.type (); + volatile int typ = mattype.type (); mattype.info (); if (typ == SparseType::Hermitian) { - // XXX FIXME XXX Write the cholesky solver and only fall - // through if cholesky factorization fails - +#ifdef HAVE_CHOLMOD + cholmod_common Common; + cholmod_common *cm = &Common; + + // Setup initial parameters + CHOLMOD_NAME(start) (cm); + cm->prefer_zomplex = FALSE; + + double spu = Voctave_sparse_controls.get_key ("spumoni"); + if (spu == 0.) + { + cm->print = -1; + cm->print_function = NULL; + } + else + { + cm->print = (int)spu + 2; + cm->print_function =&SparseCholPrint; + } + + cm->error_handler = &SparseCholError; + cm->complex_divide = CHOLMOD_NAME(divcomplex); + cm->hypotenuse = CHOLMOD_NAME(hypot); + +#ifdef HAVE_METIS + // METIS 4.0.1 uses malloc and free, and will terminate MATLAB if + // it runs out of memory. Use CHOLMOD's memory guard for METIS, + // which mxMalloc's a huge block of memory (and then immediately + // mxFree's it) before calling METIS + cm->metis_memory = 2.0; + +#if defined(METIS_VERSION) +#if (METIS_VERSION >= METIS_VER(4,0,2)) + // METIS 4.0.2 uses function pointers for malloc and free + METIS_malloc = cm->malloc_memory; + METIS_free = cm->free_memory; + // Turn off METIS memory guard. It is not needed, because mxMalloc + // will safely terminate the mexFunction and free any workspace + // without killing all of octave. + cm->metis_memory = 0.0; +#endif +#endif +#endif + + cm->final_ll = TRUE; + + cholmod_sparse Astore; + cholmod_sparse *A = &Astore; + double dummy; + A->nrow = nr; + A->ncol = nc; + + A->p = cidx(); + A->i = ridx(); + A->nzmax = nonzero(); + A->packed = TRUE; + A->sorted = TRUE; + A->nz = NULL; +#ifdef IDX_TYPE_LONG + A->itype = CHOLMOD_LONG; +#else + A->itype = CHOLMOD_INT; +#endif + A->dtype = CHOLMOD_DOUBLE; + A->stype = 1; + A->xtype = CHOLMOD_COMPLEX; + + if (nr < 1) + A->x = &dummy; + else + A->x = data(); + + cholmod_sparse Bstore; + cholmod_sparse *B = &Bstore; + B->nrow = b.rows(); + B->ncol = b.cols(); + B->p = b.cidx(); + B->i = b.ridx(); + B->nzmax = b.nonzero(); + B->packed = TRUE; + B->sorted = TRUE; + B->nz = NULL; +#ifdef IDX_TYPE_LONG + B->itype = CHOLMOD_LONG; +#else + B->itype = CHOLMOD_INT; +#endif + B->dtype = CHOLMOD_DOUBLE; + B->stype = 0; + B->xtype = CHOLMOD_REAL; + + if (b.rows() < 1 || b.cols() < 1) + B->x = &dummy; + else + B->x = b.data(); + + cholmod_factor *L; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + L = CHOLMOD_NAME(analyze) (A, cm); + CHOLMOD_NAME(factorize) (A, L, cm); + rcond = CHOLMOD_NAME(rcond)(L, cm); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + if (rcond == 0.0) + { + // Either its indefinite or singular. Try UMFPACK + mattype.mark_as_unsymmetric (); + typ = SparseType::Full; + } + else + { + volatile double rcond_plus_one = rcond + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcond)) + { + err = -2; + + if (sing_handler) + sing_handler (rcond); + else + (*current_liboctave_error_handler) + ("SparseMatrix::solve matrix singular to machine precision, rcond = %g", + rcond); + +#ifdef HAVE_LSSOLVE + return retval; +#endif + } + + cholmod_sparse *X; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + X = CHOLMOD_NAME(spsolve) (CHOLMOD_A, L, B, cm); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + retval = SparseComplexMatrix + (static_cast<octave_idx_type>(X->nrow), + static_cast<octave_idx_type>(X->ncol), + static_cast<octave_idx_type>(X->nzmax)); + for (octave_idx_type j = 0; + j <= static_cast<octave_idx_type>(X->ncol); j++) + retval.xcidx(j) = static_cast<octave_idx_type *>(X->p)[j]; + for (octave_idx_type j = 0; + j < static_cast<octave_idx_type>(X->nzmax); j++) + { + retval.xridx(j) = static_cast<octave_idx_type *>(X->i)[j]; + retval.xdata(j) = static_cast<Complex *>(X->x)[j]; + } + + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CHOLMOD_NAME(free_sparse) (&X, cm); + CHOLMOD_NAME(free_factor) (&L, cm); + CHOLMOD_NAME(finish) (cm); + CHOLMOD_NAME(print_common) (" ", cm); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } +#else (*current_liboctave_warning_handler) - ("SparseMatrix::solve XXX FIXME XXX Cholesky code not done"); + ("CHOLMOD not installed"); mattype.mark_as_unsymmetric (); typ = SparseType::Full; +#endif } if (typ == SparseType::Full) @@ -4932,19 +5580,157 @@ else { // Print spparms("spumoni") info if requested - int typ = mattype.type (); + volatile int typ = mattype.type (); mattype.info (); if (typ == SparseType::Hermitian) { - // XXX FIXME XXX Write the cholesky solver and only fall - // through if cholesky factorization fails - +#ifdef HAVE_CHOLMOD + cholmod_common Common; + cholmod_common *cm = &Common; + + // Setup initial parameters + CHOLMOD_NAME(start) (cm); + cm->prefer_zomplex = FALSE; + + double spu = Voctave_sparse_controls.get_key ("spumoni"); + if (spu == 0.) + { + cm->print = -1; + cm->print_function = NULL; + } + else + { + cm->print = (int)spu + 2; + cm->print_function =&SparseCholPrint; + } + + cm->error_handler = &SparseCholError; + cm->complex_divide = CHOLMOD_NAME(divcomplex); + cm->hypotenuse = CHOLMOD_NAME(hypot); + +#ifdef HAVE_METIS + // METIS 4.0.1 uses malloc and free, and will terminate MATLAB if + // it runs out of memory. Use CHOLMOD's memory guard for METIS, + // which mxMalloc's a huge block of memory (and then immediately + // mxFree's it) before calling METIS + cm->metis_memory = 2.0; + +#if defined(METIS_VERSION) +#if (METIS_VERSION >= METIS_VER(4,0,2)) + // METIS 4.0.2 uses function pointers for malloc and free + METIS_malloc = cm->malloc_memory; + METIS_free = cm->free_memory; + // Turn off METIS memory guard. It is not needed, because mxMalloc + // will safely terminate the mexFunction and free any workspace + // without killing all of octave. + cm->metis_memory = 0.0; +#endif +#endif +#endif + + cm->final_ll = TRUE; + + cholmod_sparse Astore; + cholmod_sparse *A = &Astore; + double dummy; + A->nrow = nr; + A->ncol = nc; + + A->p = cidx(); + A->i = ridx(); + A->nzmax = nonzero(); + A->packed = TRUE; + A->sorted = TRUE; + A->nz = NULL; +#ifdef IDX_TYPE_LONG + A->itype = CHOLMOD_LONG; +#else + A->itype = CHOLMOD_INT; +#endif + A->dtype = CHOLMOD_DOUBLE; + A->stype = 1; + A->xtype = CHOLMOD_COMPLEX; + + if (nr < 1) + A->x = &dummy; + else + A->x = data(); + + cholmod_dense Bstore; + cholmod_dense *B = &Bstore; + B->nrow = b.rows(); + B->ncol = b.cols(); + B->d = B->nrow; + B->nzmax = B->nrow * B->ncol; + B->dtype = CHOLMOD_DOUBLE; + B->xtype = CHOLMOD_COMPLEX; + if (nc < 1 || b.cols() < 1) + B->x = &dummy; + else + // We won't alter it, honest :-) + B->x = const_cast<Complex *>(b.fortran_vec()); + + cholmod_factor *L; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + L = CHOLMOD_NAME(analyze) (A, cm); + CHOLMOD_NAME(factorize) (A, L, cm); + rcond = CHOLMOD_NAME(rcond)(L, cm); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + if (rcond == 0.0) + { + // Either its indefinite or singular. Try UMFPACK + mattype.mark_as_unsymmetric (); + typ = SparseType::Full; + } + else + { + volatile double rcond_plus_one = rcond + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcond)) + { + err = -2; + + if (sing_handler) + sing_handler (rcond); + else + (*current_liboctave_error_handler) + ("SparseMatrix::solve matrix singular to machine precision, rcond = %g", + rcond); + +#ifdef HAVE_LSSOLVE + return retval; +#endif + } + + cholmod_dense *X; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + X = CHOLMOD_NAME(solve) (CHOLMOD_A, L, B, cm); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + retval.resize (b.rows (), b.cols()); + for (octave_idx_type j = 0; j < b.cols(); j++) + { + octave_idx_type jr = j * b.rows(); + for (octave_idx_type i = 0; i < b.rows(); i++) + retval.xelem(i,j) = static_cast<Complex *>(X->x)[jr + i]; + } + + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CHOLMOD_NAME(free_dense) (&X, cm); + CHOLMOD_NAME(free_factor) (&L, cm); + CHOLMOD_NAME(finish) (cm); + CHOLMOD_NAME(print_common) (" ", cm); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } +#else (*current_liboctave_warning_handler) - ("SparseMatrix::solve XXX FIXME XXX Cholesky code not done"); + ("CHOLMOD not installed"); mattype.mark_as_unsymmetric (); typ = SparseType::Full; +#endif } if (typ == SparseType::Full) @@ -5041,19 +5827,173 @@ else { // Print spparms("spumoni") info if requested - int typ = mattype.type (); + volatile int typ = mattype.type (); mattype.info (); if (typ == SparseType::Hermitian) { - // XXX FIXME XXX Write the cholesky solver and only fall - // through if cholesky factorization fails - +#ifdef HAVE_CHOLMOD + cholmod_common Common; + cholmod_common *cm = &Common; + + // Setup initial parameters + CHOLMOD_NAME(start) (cm); + cm->prefer_zomplex = FALSE; + + double spu = Voctave_sparse_controls.get_key ("spumoni"); + if (spu == 0.) + { + cm->print = -1; + cm->print_function = NULL; + } + else + { + cm->print = (int)spu + 2; + cm->print_function =&SparseCholPrint; + } + + cm->error_handler = &SparseCholError; + cm->complex_divide = CHOLMOD_NAME(divcomplex); + cm->hypotenuse = CHOLMOD_NAME(hypot); + +#ifdef HAVE_METIS + // METIS 4.0.1 uses malloc and free, and will terminate MATLAB if + // it runs out of memory. Use CHOLMOD's memory guard for METIS, + // which mxMalloc's a huge block of memory (and then immediately + // mxFree's it) before calling METIS + cm->metis_memory = 2.0; + +#if defined(METIS_VERSION) +#if (METIS_VERSION >= METIS_VER(4,0,2)) + // METIS 4.0.2 uses function pointers for malloc and free + METIS_malloc = cm->malloc_memory; + METIS_free = cm->free_memory; + // Turn off METIS memory guard. It is not needed, because mxMalloc + // will safely terminate the mexFunction and free any workspace + // without killing all of octave. + cm->metis_memory = 0.0; +#endif +#endif +#endif + + cm->final_ll = TRUE; + + cholmod_sparse Astore; + cholmod_sparse *A = &Astore; + double dummy; + A->nrow = nr; + A->ncol = nc; + + A->p = cidx(); + A->i = ridx(); + A->nzmax = nonzero(); + A->packed = TRUE; + A->sorted = TRUE; + A->nz = NULL; +#ifdef IDX_TYPE_LONG + A->itype = CHOLMOD_LONG; +#else + A->itype = CHOLMOD_INT; +#endif + A->dtype = CHOLMOD_DOUBLE; + A->stype = 1; + A->xtype = CHOLMOD_COMPLEX; + + if (nr < 1) + A->x = &dummy; + else + A->x = data(); + + cholmod_sparse Bstore; + cholmod_sparse *B = &Bstore; + B->nrow = b.rows(); + B->ncol = b.cols(); + B->p = b.cidx(); + B->i = b.ridx(); + B->nzmax = b.nonzero(); + B->packed = TRUE; + B->sorted = TRUE; + B->nz = NULL; +#ifdef IDX_TYPE_LONG + B->itype = CHOLMOD_LONG; +#else + B->itype = CHOLMOD_INT; +#endif + B->dtype = CHOLMOD_DOUBLE; + B->stype = 0; + B->xtype = CHOLMOD_COMPLEX; + + if (b.rows() < 1 || b.cols() < 1) + B->x = &dummy; + else + B->x = b.data(); + + cholmod_factor *L; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + L = CHOLMOD_NAME(analyze) (A, cm); + CHOLMOD_NAME(factorize) (A, L, cm); + rcond = CHOLMOD_NAME(rcond)(L, cm); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + if (rcond == 0.0) + { + // Either its indefinite or singular. Try UMFPACK + mattype.mark_as_unsymmetric (); + typ = SparseType::Full; + } + else + { + volatile double rcond_plus_one = rcond + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcond)) + { + err = -2; + + if (sing_handler) + sing_handler (rcond); + else + (*current_liboctave_error_handler) + ("SparseMatrix::solve matrix singular to machine precision, rcond = %g", + rcond); + +#ifdef HAVE_LSSOLVE + return retval; +#endif + } + + cholmod_sparse *X; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + X = CHOLMOD_NAME(spsolve) (CHOLMOD_A, L, B, cm); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + retval = SparseComplexMatrix + (static_cast<octave_idx_type>(X->nrow), + static_cast<octave_idx_type>(X->ncol), + static_cast<octave_idx_type>(X->nzmax)); + for (octave_idx_type j = 0; + j <= static_cast<octave_idx_type>(X->ncol); j++) + retval.xcidx(j) = static_cast<octave_idx_type *>(X->p)[j]; + for (octave_idx_type j = 0; + j < static_cast<octave_idx_type>(X->nzmax); j++) + { + retval.xridx(j) = static_cast<octave_idx_type *>(X->i)[j]; + retval.xdata(j) = static_cast<Complex *>(X->x)[j]; + } + + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CHOLMOD_NAME(free_sparse) (&X, cm); + CHOLMOD_NAME(free_factor) (&L, cm); + CHOLMOD_NAME(finish) (cm); + CHOLMOD_NAME(print_common) (" ", cm); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } +#else (*current_liboctave_warning_handler) - ("SparseMatrix::solve XXX FIXME XXX Cholesky code not done"); + ("CHOLMOD not installed"); mattype.mark_as_unsymmetric (); typ = SparseType::Full; +#endif } if (typ == SparseType::Full)
--- a/liboctave/CSparse.h +++ b/liboctave/CSparse.h @@ -120,9 +120,22 @@ friend SparseComplexMatrix conj (const SparseComplexMatrix& a); +private: + SparseComplexMatrix dinverse (SparseType &mattyp, octave_idx_type& info, + double& rcond, const bool force = false, + const bool calccond = true) const; + + SparseComplexMatrix tinverse (SparseType &mattyp, octave_idx_type& info, + double& rcond, const bool force = false, + const bool calccond = true) const; + +public: SparseComplexMatrix inverse (void) const; - SparseComplexMatrix inverse (octave_idx_type& info) const; - SparseComplexMatrix inverse (octave_idx_type& info, double& rcond, int force = 0, + SparseComplexMatrix inverse (SparseType& mattype) const; + SparseComplexMatrix inverse (SparseType& mattype, + octave_idx_type& info) const; + SparseComplexMatrix inverse (SparseType& mattype, octave_idx_type& info, + double& rcond, int force = 0, int calc_cond = 1) const; ComplexDET determinant (void) const;
--- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,3 +1,45 @@ +2005-10-23 David Bateman <dbateman@free.fr> + + * Sparse-op-defs.h (SPARSE_SPARSE_MUL): Check whether trailing zero + elements need to be removed. + + * oct-sparse.h.in: Include metis headers and some macros for long/int + versions of cholmod. + + * CSparse.cc (tinverse): New private function for the inversion of + an upper triangular matrix. + (dinverse): ditto for diagonal matrices. + (inverse): Add SparseType as an argument. Implement matrix inverse + using tinverse and dinverse. + (fsolve): Use cholmod to implement Cholesky solver. + * CSparse.h (tinverse, dinverse): Declarations + (inverse): Alter declaration to include SparseType. + + * dSparse.cc (tinverse, dinverse, inverse, fsolve): ditto. + * dSparse.h (tinverse, dinverse, inverse): ditto. + + * SparseType.cc: Fix complex constructor for hermitian matrices. + + * sparse-util.cc: New file for sparse utility functions. + * sparse-util.h: New file with declarations of sparse utility + functions. + + * sparse-base-chol.cc: New file with sparse cholesky class based + on cholmod. + * sparse-base-chol.h: New file with declaration of sparse cholesky + class based on cholmod. + + * SparseCmplxCHOL.cc: Instantiate sparse cholesky class for Complex. + * SparseCmplxCHOL.h: Declaration of sparse cholesky class. + + * SparsedbleCHOL.cc: ditto. + * SparsedbleCHOL.h: ditto. + + * Makefile.in (MATRIX_INC): Include sparse-base-chol.h. + (INCLUDES): Include sparse-util.h + (TEMPLATE_SRC): Include sparse-base-chol.cc + (MATRIX_SRC): Include SparseCmplxCHOL.cc and SparsedbleCHOL.cc + 2005-10-12 John W. Eaton <jwe@octave.org> * oct-env.cc (octave_env::have_x11_display): New function.
--- a/liboctave/Makefile.in +++ b/liboctave/Makefile.in @@ -35,8 +35,9 @@ dNDArray.h dRowVector.h dbleAEPBAL.h dbleCHOL.h dbleDET.h \ dbleHESS.h dbleLU.h dbleQR.h dbleQRP.h dbleSCHUR.h dbleSVD.h \ boolSparse.h CSparse.h dSparse.h MSparse-defs.h MSparse.h \ - Sparse.h sparse-base-lu.h SparseCmplxLU.h \ - SparsedbleLU.h Sparse-op-defs.h SparseType.h \ + Sparse.h sparse-base-lu.h SparseCmplxLU.h SparsedbleLU.h \ + sparse-base-chol.h SparseCmplxCHOL.h SparsedbleCHOL.h \ + Sparse-op-defs.h SparseType.h \ int8NDArray.h uint8NDArray.h int16NDArray.h uint16NDArray.h \ int32NDArray.h uint32NDArray.h int64NDArray.h uint64NDArray.h \ intNDArray.h @@ -66,8 +67,8 @@ oct-passwd.h oct-rand.h oct-rl-edit.h oct-rl-hist.h \ oct-shlib.h oct-sort.h oct-spparms.h oct-syscalls.h \ oct-sparse.h oct-time.h oct-types.h pathlen.h pathsearch.h \ - prog-args.h so-array.h sparse-sort.h statdefs.h str-vec.h \ - sun-utils.h sysdir.h systime.h syswait.h \ + prog-args.h so-array.h sparse-sort.h statdefs.h str-vec.h \ + sparse-util.h sun-utils.h sysdir.h systime.h syswait.h \ $(OPTS_INC) \ $(MATRIX_INC) \ $(MX_OP_INC) \ @@ -76,7 +77,7 @@ TEMPLATE_SRC := Array.cc ArrayN.cc DiagArray2.cc \ MArray.cc MArray2.cc MArrayN.cc MDiagArray2.cc \ - base-lu.cc oct-sort.cc sparse-base-lu.cc + base-lu.cc oct-sort.cc sparse-base-lu.cc sparse-base-chol.cc TI_SRC := Array-C.cc Array-b.cc Array-ch.cc Array-i.cc Array-d.cc \ Array-s.cc Array-so.cc Array-str.cc Array-idx-vec.cc \ @@ -94,7 +95,7 @@ dbleDET.cc dbleHESS.cc dbleLU.cc dbleQR.cc dbleQRP.cc \ dbleSCHUR.cc dbleSVD.cc boolSparse.cc CSparse.cc dSparse.cc \ MSparse.cc Sparse.cc SparseCmplxLU.cc SparsedbleLU.cc \ - SparseType.cc \ + SparseCmplxCHOL.cc SparsedbleCHOL.cc SparseType.cc \ int8NDArray.cc uint8NDArray.cc int16NDArray.cc uint16NDArray.cc \ int32NDArray.cc uint32NDArray.cc int64NDArray.cc uint64NDArray.cc @@ -113,7 +114,7 @@ lo-utils.cc mach-info.cc oct-alloc.cc oct-env.cc \ oct-fftw.cc oct-group.cc oct-passwd.cc oct-rand.cc oct-shlib.cc \ oct-spparms.cc oct-syscalls.cc oct-time.cc prog-args.cc \ - so-array.cc sparse-sort.cc str-vec.cc \ + so-array.cc sparse-sort.cc sparse-util.cc str-vec.cc \ $(TEMPLATE_SRC) \ $(TI_SRC) \ $(MATRIX_SRC) \
--- a/liboctave/Sparse-op-defs.h +++ b/liboctave/Sparse-op-defs.h @@ -1599,6 +1599,7 @@ } \ retval.cidx(i+1) = ii; \ } \ + retval.maybe_compress (); \ return retval; \ } \ }
new file mode 100644 --- /dev/null +++ b/liboctave/SparseCmplxCHOL.cc @@ -0,0 +1,75 @@ +/* + +Copyright (C) 2005 David Bateman +Copyright (C) 1998-2005 Andy Adler + +Octave 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. + +Octave 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, Inc., 51 Franklin Street, Fifth Floor, +Boston, MA 02110-1301, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "SparseCmplxCHOL.h" + +// Instantiate the base CHOL class for the type we need +#define OCTAVE_CHOLMOD_TYPE CHOLMOD_COMPLEX +#include "sparse-base-chol.h" +#include "sparse-base-chol.cc" +template class sparse_base_chol <SparseComplexMatrix, Complex, SparseMatrix>; + +// Compute the inverse of a matrix using the Cholesky factorization. +SparseComplexMatrix +chol2inv (const SparseComplexMatrix& r) +{ + octave_idx_type r_nr = r.rows (); + octave_idx_type r_nc = r.cols (); + SparseComplexMatrix retval; + + if (r_nr == r_nc) + { + SparseType mattype (r); + int typ = mattype.type (false); + double rcond; + octave_idx_type info; + SparseComplexMatrix rinv; + + if (typ == SparseType::Upper) + { + rinv = r.inverse(mattype, info, rcond, true, false); + retval = rinv.transpose() * rinv; + } + else if (typ == SparseType::Lower) + { + rinv = r.transpose().inverse(mattype, info, rcond, true, false); + retval = rinv.transpose() * rinv; + } + else + (*current_liboctave_error_handler) + ("spchol2inv requires triangular matrix"); + } + else + (*current_liboctave_error_handler) ("spchol2inv requires square matrix"); + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/liboctave/SparseCmplxCHOL.h @@ -0,0 +1,102 @@ +/* + +Copyright (C) 2005 David Bateman +Copyright (C) 1998-2005 Andy Adler + +Octave 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. + +Octave 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, Inc., 51 Franklin Street, Fifth Floor, +Boston, MA 02110-1301, USA. + +*/ + +#if !defined (octave_sparse_complex_CHOL_h) +#define octave_sparse_complex_CHOL_h 1 + +#include "sparse-base-chol.h" +#include "dSparse.h" +#include "CSparse.h" + +class +SparseComplexCHOL : + public sparse_base_chol <SparseComplexMatrix, Complex, SparseMatrix> +{ +public: + + SparseComplexCHOL (void) : + sparse_base_chol<SparseComplexMatrix, Complex, SparseMatrix> () { } + + SparseComplexCHOL (const SparseComplexMatrix& a, bool natural = true) : + sparse_base_chol<SparseComplexMatrix, Complex, SparseMatrix> + (a, natural) { } + + SparseComplexCHOL (const SparseComplexMatrix& a, octave_idx_type& info, + bool natural = true) : + sparse_base_chol<SparseComplexMatrix, Complex, SparseMatrix> + (a, info, natural) { } + + SparseComplexCHOL (const SparseComplexCHOL& a) : + sparse_base_chol<SparseComplexMatrix, Complex, SparseMatrix> (a) { } + + ~SparseComplexCHOL (void) { } + + SparseComplexCHOL& operator = (const SparseComplexCHOL& a) + { + if (this != &a) + sparse_base_chol <SparseComplexMatrix, Complex, SparseMatrix> :: + operator = (a); + + return *this; + } + + SparseComplexMatrix chol_matrix (void) const { return R(); } + + SparseComplexMatrix L (void) const + { return sparse_base_chol<SparseComplexMatrix, Complex, + SparseMatrix>:: L (); } + + SparseComplexMatrix R (void) const + { return sparse_base_chol<SparseComplexMatrix, Complex, + SparseMatrix>:: R (); } + + octave_idx_type P (void) const + { return sparse_base_chol<SparseComplexMatrix, Complex, + SparseMatrix>:: P (); } + + ColumnVector perm (void) const + { return sparse_base_chol<SparseComplexMatrix, Complex, + SparseMatrix>:: perm (); } + + SparseMatrix Q (void) const + { return sparse_base_chol<SparseComplexMatrix, Complex, + SparseMatrix>:: Q (); } + + double rcond (void) const + { return sparse_base_chol<SparseComplexMatrix, Complex, + SparseMatrix>:: rcond (); } + + // Compute the inverse of a matrix using the Cholesky factorization. + SparseComplexMatrix inverse (void) const + { return sparse_base_chol<SparseComplexMatrix, Complex, + SparseMatrix>:: inverse (); } +}; + +SparseComplexMatrix chol2inv (const SparseComplexMatrix& r); + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- a/liboctave/SparseType.cc +++ b/liboctave/SparseType.cc @@ -112,7 +112,7 @@ typ = tmp_typ; } - if (typ == Full) + if (typ == SparseType::Full) { // Search for banded, upper and lower triangular matrices bool singular = false; @@ -179,7 +179,7 @@ maybe_hermitian = true; } - if (typ == Full) + if (typ == SparseType::Full) { // Search for a permuted triangular matrix, and test if // permutation is singular @@ -288,7 +288,7 @@ { if (a.ridx(k) == j) { - if (a.data(i) == conj (a.data(k))) + if (a.data(i) == a.data(k)) found = true; break; } @@ -559,7 +559,7 @@ { if (a.ridx(k) == j) { - if (a.data(i) == a.data(k)) + if (a.data(i) == conj(a.data(k))) found = true; break; }
new file mode 100644 --- /dev/null +++ b/liboctave/SparsedbleCHOL.cc @@ -0,0 +1,75 @@ +/* + +Copyright (C) 2005 David Bateman +Copyright (C) 1998-2005 Andy Adler + +Octave 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. + +Octave 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, Inc., 51 Franklin Street, Fifth Floor, +Boston, MA 02110-1301, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "SparsedbleCHOL.h" + +// Instantiate the base CHOL class for the type we need +#define OCTAVE_CHOLMOD_TYPE CHOLMOD_REAL +#include "sparse-base-chol.h" +#include "sparse-base-chol.cc" +template class sparse_base_chol <SparseMatrix, double, SparseMatrix>; + +// Compute the inverse of a matrix using the Cholesky factorization. +SparseMatrix +chol2inv (const SparseMatrix& r) +{ + octave_idx_type r_nr = r.rows (); + octave_idx_type r_nc = r.cols (); + SparseMatrix retval; + + if (r_nr == r_nc) + { + SparseType mattype (r); + int typ = mattype.type (false); + double rcond; + octave_idx_type info; + SparseMatrix rinv; + + if (typ == SparseType::Upper) + { + rinv = r.inverse(mattype, info, rcond, true, false); + retval = rinv.transpose() * rinv; + } + else if (typ == SparseType::Lower) + { + rinv = r.transpose().inverse(mattype, info, rcond, true, false); + retval = rinv.transpose() * rinv; + } + else + (*current_liboctave_error_handler) + ("spchol2inv requires triangular matrix"); + } + else + (*current_liboctave_error_handler) ("spchol2inv requires square matrix"); + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/liboctave/SparsedbleCHOL.h @@ -0,0 +1,90 @@ +/* + +Copyright (C) 2005 David Bateman +Copyright (C) 1998-2005 Andy Adler + +Octave 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. + +Octave 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, Inc., 51 Franklin Street, Fifth Floor, +Boston, MA 02110-1301, USA. + +*/ + +#if !defined (octave_sparse_CHOL_h) +#define octave_sparse_CHOL_h 1 + +#include "sparse-base-chol.h" +#include "dSparse.h" + +class +SparseCHOL : public sparse_base_chol <SparseMatrix, double, SparseMatrix> +{ +public: + + SparseCHOL (void) : sparse_base_chol<SparseMatrix, double, SparseMatrix> () { } + + SparseCHOL (const SparseMatrix& a, bool natural = true) : + sparse_base_chol<SparseMatrix, double, SparseMatrix> (a, natural) { } + + SparseCHOL (const SparseMatrix& a, octave_idx_type& info, + bool natural = true) : + sparse_base_chol<SparseMatrix, double, SparseMatrix> (a, info, natural) { } + + SparseCHOL (const SparseCHOL& a) : + sparse_base_chol<SparseMatrix, double, SparseMatrix> (a) { } + + ~SparseCHOL (void) { } + + SparseCHOL& operator = (const SparseCHOL& a) + { + if (this != &a) + sparse_base_chol <SparseMatrix, double, SparseMatrix> :: operator = (a); + + return *this; + } + + SparseMatrix chol_matrix (void) const { return R(); } + + SparseMatrix L (void) const + { return sparse_base_chol<SparseMatrix, double, SparseMatrix>:: L (); } + + SparseMatrix R (void) const + { return sparse_base_chol<SparseMatrix, double, SparseMatrix>:: R (); } + + octave_idx_type P (void) const + { return sparse_base_chol<SparseMatrix, double, SparseMatrix>:: P (); } + + ColumnVector perm (void) const + { return sparse_base_chol<SparseMatrix, double, SparseMatrix>:: perm (); } + + SparseMatrix Q (void) const + { return sparse_base_chol<SparseMatrix, double, SparseMatrix>:: Q (); } + + double rcond (void) const + { return sparse_base_chol<SparseMatrix, double, SparseMatrix>:: rcond (); } + + // Compute the inverse of a matrix using the Cholesky factorization. + SparseMatrix inverse (void) const + { return sparse_base_chol<SparseMatrix, double, SparseMatrix>:: + inverse (); } +}; + +SparseMatrix chol2inv (const SparseMatrix& r); + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- a/liboctave/dSparse.cc +++ b/liboctave/dSparse.cc @@ -42,6 +42,8 @@ #include "SparsedbleLU.h" #include "SparseType.h" #include "oct-sparse.h" +#include "sparse-util.h" +#include "SparsedbleCHOL.h" // Fortran functions we call. extern "C" @@ -672,23 +674,378 @@ { octave_idx_type info; double rcond; - return inverse (info, rcond, 0, 0); + SparseType mattype (*this); + return inverse (mattype, info, rcond, 0, 0); +} + +SparseMatrix +SparseMatrix::inverse (SparseType& mattype) const +{ + octave_idx_type info; + double rcond; + return inverse (mattype, info, rcond, 0, 0); } SparseMatrix -SparseMatrix::inverse (octave_idx_type& info) const +SparseMatrix::inverse (SparseType& mattype, octave_idx_type& info) const { double rcond; - return inverse (info, rcond, 0, 0); + return inverse (mattype, info, rcond, 0, 0); +} + +SparseMatrix +SparseMatrix::dinverse (SparseType &mattyp, octave_idx_type& info, + double& rcond, const bool force, + const bool calccond) const +{ + SparseMatrix retval; + + octave_idx_type nr = rows (); + octave_idx_type nc = cols (); + info = 0; + + if (nr == 0 || nc == 0 || nr != nc) + (*current_liboctave_error_handler) ("inverse requires square matrix"); + else + { + // Print spparms("spumoni") info if requested + int typ = mattyp.type (); + mattyp.info (); + + if (typ == SparseType::Diagonal || + typ == SparseType::Permuted_Diagonal) + { + if (typ == SparseType::Permuted_Diagonal) + retval = transpose(); + else + retval = *this; + + // Force make_unique to be called + double *v = retval.data(); + + if (calccond) + { + double dmax = 0., dmin = octave_Inf; + for (octave_idx_type i = 0; i < nr; i++) + { + double tmp = fabs(v[i]); + if (tmp > dmax) + dmax = tmp; + if (tmp < dmin) + dmin = tmp; + } + rcond = dmin / dmax; + } + + for (octave_idx_type i = 0; i < nr; i++) + v[i] = 1.0 / v[i]; + } + else + (*current_liboctave_error_handler) ("incorrect matrix type"); + } + + return retval; +} + +SparseMatrix +SparseMatrix::tinverse (SparseType &mattyp, octave_idx_type& info, + double& rcond, const bool force, + const bool calccond) const +{ + SparseMatrix retval; + + octave_idx_type nr = rows (); + octave_idx_type nc = cols (); + info = 0; + + if (nr == 0 || nc == 0 || nr != nc) + (*current_liboctave_error_handler) ("inverse requires square matrix"); + else + { + // Print spparms("spumoni") info if requested + int typ = mattyp.type (); + mattyp.info (); + + if (typ == SparseType::Upper || typ == SparseType::Permuted_Upper || + typ == SparseType::Lower || typ == SparseType::Permuted_Lower) + { + double anorm = 0.; + double ainvnorm = 0.; + + if (calccond) + { + // Calculate the 1-norm of matrix for rcond calculation + for (octave_idx_type j = 0; j < nr; j++) + { + double atmp = 0.; + for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) + atmp += fabs(data(i)); + if (atmp > anorm) + anorm = atmp; + } + } + + if (typ == SparseType::Upper || typ == SparseType::Lower) + { + octave_idx_type nz = nnz(); + octave_idx_type cx = 0; + octave_idx_type nz2 = nz; + retval = SparseMatrix (nr, nc, nz2); + + for (octave_idx_type i = 0; i < nr; i++) + { + OCTAVE_QUIT; + // place the 1 in the identity position + octave_idx_type cx_colstart = cx; + + if (cx == nz2) + { + nz2 *= 2; + retval.change_capacity (nz2); + } + + retval.xcidx(i) = cx; + retval.xridx(cx) = i; + retval.xdata(cx) = 1.0; + cx++; + + // iterate accross columns of input matrix + for (octave_idx_type j = i+1; j < nr; j++) + { + double v = 0.; + // iterate to calculate sum + octave_idx_type colXp = retval.xcidx(i); + octave_idx_type colUp = cidx(j); + octave_idx_type rpX, rpU; + do + { + OCTAVE_QUIT; + rpX = retval.xridx(colXp); + rpU = ridx(colUp); + + if (rpX < rpU) + colXp++; + else if (rpX > rpU) + colUp++; + else + { + v -= retval.xdata(colXp) * data(colUp); + colXp++; + colUp++; + } + } while ((rpX<j) && (rpU<j) && + (colXp<cx) && (colUp<nz)); + + // get A(m,m) + colUp = cidx(j+1) - 1; + double pivot = data(colUp); + if (pivot == 0.) + (*current_liboctave_error_handler) + ("division by zero"); + + if (v != 0.) + { + if (cx == nz2) + { + nz2 *= 2; + retval.change_capacity (nz2); + } + + retval.xridx(cx) = j; + retval.xdata(cx) = v / pivot; + cx++; + } + } + + // get A(m,m) + octave_idx_type colUp = cidx(i+1) - 1; + double pivot = data(colUp); + if (pivot == 0.) + (*current_liboctave_error_handler) ("division by zero"); + + if (pivot != 1.0) + for (octave_idx_type j = cx_colstart; j < cx; j++) + retval.xdata(j) /= pivot; + } + retval.xcidx(nr) = cx; + retval.maybe_compress (); + } + else + { + octave_idx_type nz = nnz(); + octave_idx_type cx = 0; + octave_idx_type nz2 = nz; + retval = SparseMatrix (nr, nc, nz2); + + OCTAVE_LOCAL_BUFFER (double, work, nr); + OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nr); + + octave_idx_type *perm = mattyp.triangular_perm(); + if (typ == SparseType::Permuted_Upper) + { + for (octave_idx_type i = 0; i < nr; i++) + rperm[perm[i]] = i; + } + else + { + for (octave_idx_type i = 0; i < nr; i++) + rperm[i] = perm[i]; + for (octave_idx_type i = 0; i < nr; i++) + perm[rperm[i]] = i; + } + + for (octave_idx_type i = 0; i < nr; i++) + { + OCTAVE_QUIT; + octave_idx_type iidx = rperm[i]; + + for (octave_idx_type j = 0; j < nr; j++) + work[j] = 0.; + + // place the 1 in the identity position + work[iidx] = 1.0; + + // iterate accross columns of input matrix + for (octave_idx_type j = iidx+1; j < nr; j++) + { + double v = 0.; + octave_idx_type jidx = perm[j]; + // iterate to calculate sum + for (octave_idx_type k = cidx(jidx); + k < cidx(jidx+1); k++) + { + OCTAVE_QUIT; + v -= work[ridx(k)] * data(k); + } + + // get A(m,m) + double pivot = data(cidx(jidx+1) - 1); + if (pivot == 0.) + (*current_liboctave_error_handler) + ("division by zero"); + + work[j] = v / pivot; + } + + // get A(m,m) + octave_idx_type colUp = cidx(perm[iidx]+1) - 1; + double pivot = data(colUp); + if (pivot == 0.) + (*current_liboctave_error_handler) + ("division by zero"); + + octave_idx_type new_cx = cx; + for (octave_idx_type j = iidx; j < nr; j++) + if (work[j] != 0.0) + { + new_cx++; + if (pivot != 1.0) + work[j] /= pivot; + } + + if (cx < new_cx) + { + nz2 = (2*nz2 < new_cx ? new_cx : 2*nz2); + retval.change_capacity (nz2); + } + + retval.xcidx(i) = cx; + for (octave_idx_type j = iidx; j < nr; j++) + if (work[j] != 0.) + { + retval.xridx(cx) = j; + retval.xdata(cx++) = work[j]; + } + } + + retval.xcidx(nr) = cx; + retval.maybe_compress (); + } + + if (calccond) + { + // Calculate the 1-norm of inverse matrix for rcond calculation + for (octave_idx_type j = 0; j < nr; j++) + { + double atmp = 0.; + for (octave_idx_type i = retval.cidx(j); + i < retval.cidx(j+1); i++) + atmp += fabs(retval.data(i)); + if (atmp > ainvnorm) + ainvnorm = atmp; + } + + rcond = 1. / ainvnorm / anorm; + } + } + else + (*current_liboctave_error_handler) ("incorrect matrix type"); + } + + return retval; } SparseMatrix -SparseMatrix::inverse (octave_idx_type& info, double& rcond, int force, int calc_cond) const -{ - info = -1; - (*current_liboctave_error_handler) - ("SparseMatrix::inverse not implemented yet"); - return SparseMatrix (); +SparseMatrix::inverse (SparseType &mattype, octave_idx_type& info, + double& rcond, int force, int calc_cond) const +{ + int typ = mattype.type (false); + SparseMatrix ret; + + if (typ == SparseType::Unknown) + typ = mattype.type (*this); + + if (typ == SparseType::Diagonal || typ == SparseType::Permuted_Diagonal) + ret = dinverse (mattype, info, rcond, true, calc_cond); + else if (typ == SparseType::Upper || typ == SparseType::Permuted_Upper) + ret = tinverse (mattype, info, rcond, true, calc_cond).transpose(); + else if (typ == SparseType::Lower || typ == SparseType::Permuted_Lower) + ret = transpose().tinverse (mattype, info, rcond, true, calc_cond); + else if (typ != SparseType::Rectangular) + { + if (mattype.is_hermitian()) + { + SparseType tmp_typ (SparseType::Upper); + SparseCHOL fact (*this, info, false); + rcond = fact.rcond(); + if (info == 0) + { + double rcond2; + SparseMatrix Q = fact.Q(); + SparseMatrix InvL = fact.L().transpose().tinverse(tmp_typ, + info, rcond2, true, false); + ret = Q * InvL.transpose() * InvL * Q.transpose(); + } + else + { + // Matrix is either singular or not positive definite + mattype.mark_as_unsymmetric (); + typ = SparseType::Full; + } + } + + if (!mattype.is_hermitian()) + { + octave_idx_type n = rows(); + ColumnVector Qinit(n); + for (octave_idx_type i = 0; i < n; i++) + Qinit(i) = i; + + SparseType tmp_typ (SparseType::Upper); + SparseLU fact (*this, Qinit, -1.0, false); + rcond = fact.rcond(); + double rcond2; + SparseMatrix InvL = fact.L().transpose().tinverse(tmp_typ, + info, rcond2, true, false); + SparseMatrix InvU = fact.U().tinverse(tmp_typ, info, rcond2, + true, false).transpose(); + ret = fact.Pc().transpose() * InvU * InvL * fact.Pr(); + } + } + else + (*current_liboctave_error_handler) ("inverse requires square matrix"); + + return ret; } DET @@ -4859,19 +5216,157 @@ else { // Print spparms("spumoni") info if requested - int typ = mattype.type (); + volatile int typ = mattype.type (); mattype.info (); if (typ == SparseType::Hermitian) { - // XXX FIXME XXX Write the cholesky solver and only fall - // through if cholesky factorization fails - +#ifdef HAVE_CHOLMOD + cholmod_common Common; + cholmod_common *cm = &Common; + + // Setup initial parameters + CHOLMOD_NAME(start) (cm); + cm->prefer_zomplex = FALSE; + + double spu = Voctave_sparse_controls.get_key ("spumoni"); + if (spu == 0.) + { + cm->print = -1; + cm->print_function = NULL; + } + else + { + cm->print = (int)spu + 2; + cm->print_function =&SparseCholPrint; + } + + cm->error_handler = &SparseCholError; + cm->complex_divide = CHOLMOD_NAME(divcomplex); + cm->hypotenuse = CHOLMOD_NAME(hypot); + +#ifdef HAVE_METIS + // METIS 4.0.1 uses malloc and free, and will terminate MATLAB if + // it runs out of memory. Use CHOLMOD's memory guard for METIS, + // which mxMalloc's a huge block of memory (and then immediately + // mxFree's it) before calling METIS + cm->metis_memory = 2.0; + +#if defined(METIS_VERSION) +#if (METIS_VERSION >= METIS_VER(4,0,2)) + // METIS 4.0.2 uses function pointers for malloc and free + METIS_malloc = cm->malloc_memory; + METIS_free = cm->free_memory; + // Turn off METIS memory guard. It is not needed, because mxMalloc + // will safely terminate the mexFunction and free any workspace + // without killing all of octave. + cm->metis_memory = 0.0; +#endif +#endif +#endif + + cm->final_ll = TRUE; + + cholmod_sparse Astore; + cholmod_sparse *A = &Astore; + double dummy; + A->nrow = nr; + A->ncol = nc; + + A->p = cidx(); + A->i = ridx(); + A->nzmax = nonzero(); + A->packed = TRUE; + A->sorted = TRUE; + A->nz = NULL; +#ifdef IDX_TYPE_LONG + A->itype = CHOLMOD_LONG; +#else + A->itype = CHOLMOD_INT; +#endif + A->dtype = CHOLMOD_DOUBLE; + A->stype = 1; + A->xtype = CHOLMOD_REAL; + + if (nr < 1) + A->x = &dummy; + else + A->x = data(); + + cholmod_dense Bstore; + cholmod_dense *B = &Bstore; + B->nrow = b.rows(); + B->ncol = b.cols(); + B->d = B->nrow; + B->nzmax = B->nrow * B->ncol; + B->dtype = CHOLMOD_DOUBLE; + B->xtype = CHOLMOD_REAL; + if (nc < 1 || b.cols() < 1) + B->x = &dummy; + else + // We won't alter it, honest :-) + B->x = const_cast<double *>(b.fortran_vec()); + + cholmod_factor *L; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + L = CHOLMOD_NAME(analyze) (A, cm); + CHOLMOD_NAME(factorize) (A, L, cm); + rcond = CHOLMOD_NAME(rcond)(L, cm); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + if (rcond == 0.0) + { + // Either its indefinite or singular. Try UMFPACK + mattype.mark_as_unsymmetric (); + typ = SparseType::Full; + } + else + { + volatile double rcond_plus_one = rcond + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcond)) + { + err = -2; + + if (sing_handler) + sing_handler (rcond); + else + (*current_liboctave_error_handler) + ("SparseMatrix::solve matrix singular to machine precision, rcond = %g", + rcond); + +#ifdef HAVE_LSSOLVE + return retval; +#endif + } + + cholmod_dense *X; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + X = CHOLMOD_NAME(solve) (CHOLMOD_A, L, B, cm); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + retval.resize (b.rows (), b.cols()); + for (octave_idx_type j = 0; j < b.cols(); j++) + { + octave_idx_type jr = j * b.rows(); + for (octave_idx_type i = 0; i < b.rows(); i++) + retval.xelem(i,j) = static_cast<double *>(X->x)[jr + i]; + } + + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CHOLMOD_NAME(free_dense) (&X, cm); + CHOLMOD_NAME(free_factor) (&L, cm); + CHOLMOD_NAME(finish) (cm); + CHOLMOD_NAME(print_common) (" ", cm); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } +#else (*current_liboctave_warning_handler) - ("SparseMatrix::solve XXX FIXME XXX Cholesky code not done"); + ("CHOLMOD not installed"); mattype.mark_as_unsymmetric (); typ = SparseType::Full; +#endif } if (typ == SparseType::Full) @@ -4963,19 +5458,172 @@ else { // Print spparms("spumoni") info if requested - int typ = mattype.type (); + volatile int typ = mattype.type (); mattype.info (); if (typ == SparseType::Hermitian) { - // XXX FIXME XXX Write the cholesky solver and only fall - // through if cholesky factorization fails - +#ifdef HAVE_CHOLMOD + cholmod_common Common; + cholmod_common *cm = &Common; + + // Setup initial parameters + CHOLMOD_NAME(start) (cm); + cm->prefer_zomplex = FALSE; + + double spu = Voctave_sparse_controls.get_key ("spumoni"); + if (spu == 0.) + { + cm->print = -1; + cm->print_function = NULL; + } + else + { + cm->print = (int)spu + 2; + cm->print_function =&SparseCholPrint; + } + + cm->error_handler = &SparseCholError; + cm->complex_divide = CHOLMOD_NAME(divcomplex); + cm->hypotenuse = CHOLMOD_NAME(hypot); + +#ifdef HAVE_METIS + // METIS 4.0.1 uses malloc and free, and will terminate MATLAB if + // it runs out of memory. Use CHOLMOD's memory guard for METIS, + // which mxMalloc's a huge block of memory (and then immediately + // mxFree's it) before calling METIS + cm->metis_memory = 2.0; + +#if defined(METIS_VERSION) +#if (METIS_VERSION >= METIS_VER(4,0,2)) + // METIS 4.0.2 uses function pointers for malloc and free + METIS_malloc = cm->malloc_memory; + METIS_free = cm->free_memory; + // Turn off METIS memory guard. It is not needed, because mxMalloc + // will safely terminate the mexFunction and free any workspace + // without killing all of octave. + cm->metis_memory = 0.0; +#endif +#endif +#endif + + cm->final_ll = TRUE; + + cholmod_sparse Astore; + cholmod_sparse *A = &Astore; + double dummy; + A->nrow = nr; + A->ncol = nc; + + A->p = cidx(); + A->i = ridx(); + A->nzmax = nonzero(); + A->packed = TRUE; + A->sorted = TRUE; + A->nz = NULL; +#ifdef IDX_TYPE_LONG + A->itype = CHOLMOD_LONG; +#else + A->itype = CHOLMOD_INT; +#endif + A->dtype = CHOLMOD_DOUBLE; + A->stype = 1; + A->xtype = CHOLMOD_REAL; + + if (nr < 1) + A->x = &dummy; + else + A->x = data(); + + cholmod_sparse Bstore; + cholmod_sparse *B = &Bstore; + B->nrow = b.rows(); + B->ncol = b.cols(); + B->p = b.cidx(); + B->i = b.ridx(); + B->nzmax = b.nonzero(); + B->packed = TRUE; + B->sorted = TRUE; + B->nz = NULL; +#ifdef IDX_TYPE_LONG + B->itype = CHOLMOD_LONG; +#else + B->itype = CHOLMOD_INT; +#endif + B->dtype = CHOLMOD_DOUBLE; + B->stype = 0; + B->xtype = CHOLMOD_REAL; + + if (b.rows() < 1 || b.cols() < 1) + B->x = &dummy; + else + B->x = b.data(); + + cholmod_factor *L; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + L = CHOLMOD_NAME(analyze) (A, cm); + CHOLMOD_NAME(factorize) (A, L, cm); + rcond = CHOLMOD_NAME(rcond)(L, cm); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + if (rcond == 0.0) + { + // Either its indefinite or singular. Try UMFPACK + mattype.mark_as_unsymmetric (); + typ = SparseType::Full; + } + else + { + volatile double rcond_plus_one = rcond + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcond)) + { + err = -2; + + if (sing_handler) + sing_handler (rcond); + else + (*current_liboctave_error_handler) + ("SparseMatrix::solve matrix singular to machine precision, rcond = %g", + rcond); + +#ifdef HAVE_LSSOLVE + return retval; +#endif + } + + cholmod_sparse *X; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + X = CHOLMOD_NAME(spsolve) (CHOLMOD_A, L, B, cm); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + retval = SparseMatrix (static_cast<octave_idx_type>(X->nrow), + static_cast<octave_idx_type>(X->ncol), + static_cast<octave_idx_type>(X->nzmax)); + for (octave_idx_type j = 0; + j <= static_cast<octave_idx_type>(X->ncol); j++) + retval.xcidx(j) = static_cast<octave_idx_type *>(X->p)[j]; + for (octave_idx_type j = 0; + j < static_cast<octave_idx_type>(X->nzmax); j++) + { + retval.xridx(j) = static_cast<octave_idx_type *>(X->i)[j]; + retval.xdata(j) = static_cast<double *>(X->x)[j]; + } + + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CHOLMOD_NAME(free_sparse) (&X, cm); + CHOLMOD_NAME(free_factor) (&L, cm); + CHOLMOD_NAME(finish) (cm); + CHOLMOD_NAME(print_common) (" ", cm); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } +#else (*current_liboctave_warning_handler) - ("SparseMatrix::solve XXX FIXME XXX Cholesky code not done"); + ("CHOLMOD not installed"); mattype.mark_as_unsymmetric (); typ = SparseType::Full; +#endif } if (typ == SparseType::Full) @@ -5099,19 +5747,157 @@ else { // Print spparms("spumoni") info if requested - int typ = mattype.type (); + volatile int typ = mattype.type (); mattype.info (); if (typ == SparseType::Hermitian) { - // XXX FIXME XXX Write the cholesky solver and only fall - // through if cholesky factorization fails - +#ifdef HAVE_CHOLMOD + cholmod_common Common; + cholmod_common *cm = &Common; + + // Setup initial parameters + CHOLMOD_NAME(start) (cm); + cm->prefer_zomplex = FALSE; + + double spu = Voctave_sparse_controls.get_key ("spumoni"); + if (spu == 0.) + { + cm->print = -1; + cm->print_function = NULL; + } + else + { + cm->print = (int)spu + 2; + cm->print_function =&SparseCholPrint; + } + + cm->error_handler = &SparseCholError; + cm->complex_divide = CHOLMOD_NAME(divcomplex); + cm->hypotenuse = CHOLMOD_NAME(hypot); + +#ifdef HAVE_METIS + // METIS 4.0.1 uses malloc and free, and will terminate MATLAB if + // it runs out of memory. Use CHOLMOD's memory guard for METIS, + // which mxMalloc's a huge block of memory (and then immediately + // mxFree's it) before calling METIS + cm->metis_memory = 2.0; + +#if defined(METIS_VERSION) +#if (METIS_VERSION >= METIS_VER(4,0,2)) + // METIS 4.0.2 uses function pointers for malloc and free + METIS_malloc = cm->malloc_memory; + METIS_free = cm->free_memory; + // Turn off METIS memory guard. It is not needed, because mxMalloc + // will safely terminate the mexFunction and free any workspace + // without killing all of octave. + cm->metis_memory = 0.0; +#endif +#endif +#endif + + cm->final_ll = TRUE; + + cholmod_sparse Astore; + cholmod_sparse *A = &Astore; + double dummy; + A->nrow = nr; + A->ncol = nc; + + A->p = cidx(); + A->i = ridx(); + A->nzmax = nonzero(); + A->packed = TRUE; + A->sorted = TRUE; + A->nz = NULL; +#ifdef IDX_TYPE_LONG + A->itype = CHOLMOD_LONG; +#else + A->itype = CHOLMOD_INT; +#endif + A->dtype = CHOLMOD_DOUBLE; + A->stype = 1; + A->xtype = CHOLMOD_REAL; + + if (nr < 1) + A->x = &dummy; + else + A->x = data(); + + cholmod_dense Bstore; + cholmod_dense *B = &Bstore; + B->nrow = b.rows(); + B->ncol = b.cols(); + B->d = B->nrow; + B->nzmax = B->nrow * B->ncol; + B->dtype = CHOLMOD_DOUBLE; + B->xtype = CHOLMOD_COMPLEX; + if (nc < 1 || b.cols() < 1) + B->x = &dummy; + else + // We won't alter it, honest :-) + B->x = const_cast<Complex *>(b.fortran_vec()); + + cholmod_factor *L; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + L = CHOLMOD_NAME(analyze) (A, cm); + CHOLMOD_NAME(factorize) (A, L, cm); + rcond = CHOLMOD_NAME(rcond)(L, cm); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + if (rcond == 0.0) + { + // Either its indefinite or singular. Try UMFPACK + mattype.mark_as_unsymmetric (); + typ = SparseType::Full; + } + else + { + volatile double rcond_plus_one = rcond + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcond)) + { + err = -2; + + if (sing_handler) + sing_handler (rcond); + else + (*current_liboctave_error_handler) + ("SparseMatrix::solve matrix singular to machine precision, rcond = %g", + rcond); + +#ifdef HAVE_LSSOLVE + return retval; +#endif + } + + cholmod_dense *X; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + X = CHOLMOD_NAME(solve) (CHOLMOD_A, L, B, cm); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + retval.resize (b.rows (), b.cols()); + for (octave_idx_type j = 0; j < b.cols(); j++) + { + octave_idx_type jr = j * b.rows(); + for (octave_idx_type i = 0; i < b.rows(); i++) + retval.xelem(i,j) = static_cast<Complex *>(X->x)[jr + i]; + } + + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CHOLMOD_NAME(free_dense) (&X, cm); + CHOLMOD_NAME(free_factor) (&L, cm); + CHOLMOD_NAME(finish) (cm); + CHOLMOD_NAME(print_common) (" ", cm); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } +#else (*current_liboctave_warning_handler) - ("SparseMatrix::solve XXX FIXME XXX Cholesky code not done"); + ("CHOLMOD not installed"); mattype.mark_as_unsymmetric (); typ = SparseType::Full; +#endif } if (typ == SparseType::Full) @@ -5223,19 +6009,173 @@ else { // Print spparms("spumoni") info if requested - int typ = mattype.type (); + volatile int typ = mattype.type (); mattype.info (); if (typ == SparseType::Hermitian) { - // XXX FIXME XXX Write the cholesky solver and only fall - // through if cholesky factorization fails - +#ifdef HAVE_CHOLMOD + cholmod_common Common; + cholmod_common *cm = &Common; + + // Setup initial parameters + CHOLMOD_NAME(start) (cm); + cm->prefer_zomplex = FALSE; + + double spu = Voctave_sparse_controls.get_key ("spumoni"); + if (spu == 0.) + { + cm->print = -1; + cm->print_function = NULL; + } + else + { + cm->print = (int)spu + 2; + cm->print_function =&SparseCholPrint; + } + + cm->error_handler = &SparseCholError; + cm->complex_divide = CHOLMOD_NAME(divcomplex); + cm->hypotenuse = CHOLMOD_NAME(hypot); + +#ifdef HAVE_METIS + // METIS 4.0.1 uses malloc and free, and will terminate MATLAB if + // it runs out of memory. Use CHOLMOD's memory guard for METIS, + // which mxMalloc's a huge block of memory (and then immediately + // mxFree's it) before calling METIS + cm->metis_memory = 2.0; + +#if defined(METIS_VERSION) +#if (METIS_VERSION >= METIS_VER(4,0,2)) + // METIS 4.0.2 uses function pointers for malloc and free + METIS_malloc = cm->malloc_memory; + METIS_free = cm->free_memory; + // Turn off METIS memory guard. It is not needed, because mxMalloc + // will safely terminate the mexFunction and free any workspace + // without killing all of octave. + cm->metis_memory = 0.0; +#endif +#endif +#endif + + cm->final_ll = TRUE; + + cholmod_sparse Astore; + cholmod_sparse *A = &Astore; + double dummy; + A->nrow = nr; + A->ncol = nc; + + A->p = cidx(); + A->i = ridx(); + A->nzmax = nonzero(); + A->packed = TRUE; + A->sorted = TRUE; + A->nz = NULL; +#ifdef IDX_TYPE_LONG + A->itype = CHOLMOD_LONG; +#else + A->itype = CHOLMOD_INT; +#endif + A->dtype = CHOLMOD_DOUBLE; + A->stype = 1; + A->xtype = CHOLMOD_REAL; + + if (nr < 1) + A->x = &dummy; + else + A->x = data(); + + cholmod_sparse Bstore; + cholmod_sparse *B = &Bstore; + B->nrow = b.rows(); + B->ncol = b.cols(); + B->p = b.cidx(); + B->i = b.ridx(); + B->nzmax = b.nonzero(); + B->packed = TRUE; + B->sorted = TRUE; + B->nz = NULL; +#ifdef IDX_TYPE_LONG + B->itype = CHOLMOD_LONG; +#else + B->itype = CHOLMOD_INT; +#endif + B->dtype = CHOLMOD_DOUBLE; + B->stype = 0; + B->xtype = CHOLMOD_COMPLEX; + + if (b.rows() < 1 || b.cols() < 1) + B->x = &dummy; + else + B->x = b.data(); + + cholmod_factor *L; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + L = CHOLMOD_NAME(analyze) (A, cm); + CHOLMOD_NAME(factorize) (A, L, cm); + rcond = CHOLMOD_NAME(rcond)(L, cm); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + if (rcond == 0.0) + { + // Either its indefinite or singular. Try UMFPACK + mattype.mark_as_unsymmetric (); + typ = SparseType::Full; + } + else + { + volatile double rcond_plus_one = rcond + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcond)) + { + err = -2; + + if (sing_handler) + sing_handler (rcond); + else + (*current_liboctave_error_handler) + ("SparseMatrix::solve matrix singular to machine precision, rcond = %g", + rcond); + +#ifdef HAVE_LSSOLVE + return retval; +#endif + } + + cholmod_sparse *X; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + X = CHOLMOD_NAME(spsolve) (CHOLMOD_A, L, B, cm); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + retval = SparseComplexMatrix + (static_cast<octave_idx_type>(X->nrow), + static_cast<octave_idx_type>(X->ncol), + static_cast<octave_idx_type>(X->nzmax)); + for (octave_idx_type j = 0; + j <= static_cast<octave_idx_type>(X->ncol); j++) + retval.xcidx(j) = static_cast<octave_idx_type *>(X->p)[j]; + for (octave_idx_type j = 0; + j < static_cast<octave_idx_type>(X->nzmax); j++) + { + retval.xridx(j) = static_cast<octave_idx_type *>(X->i)[j]; + retval.xdata(j) = static_cast<Complex *>(X->x)[j]; + } + + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CHOLMOD_NAME(free_sparse) (&X, cm); + CHOLMOD_NAME(free_factor) (&L, cm); + CHOLMOD_NAME(finish) (cm); + CHOLMOD_NAME(print_common) (" ", cm); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } +#else (*current_liboctave_warning_handler) - ("SparseMatrix::solve XXX FIXME XXX Cholesky code not done"); + ("CHOLMOD not installed"); mattype.mark_as_unsymmetric (); typ = SparseType::Full; +#endif } if (typ == SparseType::Full)
--- a/liboctave/dSparse.h +++ b/liboctave/dSparse.h @@ -112,11 +112,23 @@ { return MSparse<double>::transpose (); } + SparseMatrix hermitian (void) const { return transpose (); } +private: + SparseMatrix dinverse (SparseType &mattyp, octave_idx_type& info, + double& rcond, const bool force = false, + const bool calccond = true) const; + + SparseMatrix tinverse (SparseType &mattyp, octave_idx_type& info, + double& rcond, const bool force = false, + const bool calccond = true) const; + +public: SparseMatrix inverse (void) const; - SparseMatrix inverse (octave_idx_type& info) const; - SparseMatrix inverse (octave_idx_type& info, double& rcond, int force = 0, - int calc_cond = 1) const; + SparseMatrix inverse (SparseType& mattype) const; + SparseMatrix inverse (SparseType& mattype, octave_idx_type& info) const; + SparseMatrix inverse (SparseType& mattype, octave_idx_type& info, + double& rcond, int force = 0, int calc_cond = 1) const; DET determinant (void) const; DET determinant (octave_idx_type& info) const;
--- a/liboctave/oct-sparse.h.in +++ b/liboctave/oct-sparse.h.in @@ -24,6 +24,10 @@ #if !defined (oct_sparse_h) #define oct_sparse_h 1 +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + #ifdef __cplusplus extern "C" { #endif @@ -43,9 +47,27 @@ #include <@CCOLAMD_INCLUDE@> #endif +#ifdef HAVE_METIS +/* External METIS functions in C */ +#include <@METIS_INCLUDE@> +#endif + #ifdef HAVE_CHOLMOD /* External CHOLMOD functions in C */ #include <@CHOLMOD_INCLUDE@> + +#ifndef FALSE +#define FALSE 0 +#endif +#ifndef TRUE +#define TRUE 1 +#endif + +#ifdef IDX_TYPE_LONG +#define CHOLMOD_NAME(name) cholmod_l_ ## name +#else +#define CHOLMOD_NAME(name) cholmod_ ## name +#endif #endif #ifdef __cplusplus
new file mode 100644 --- /dev/null +++ b/liboctave/sparse-base-chol.cc @@ -0,0 +1,304 @@ +/* + +Copyright (C) 2005 David Bateman +Copyright (C) 1998-2005 Andy Adler + +Octave 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. + +Octave 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, Inc., 51 Franklin Street, Fifth Floor, +Boston, MA 02110-1301, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "sparse-base-chol.h" +#include "sparse-util.h" +#include "lo-error.h" +#include "oct-sparse.h" +#include "oct-spparms.h" +#include "quit.h" +#include "SparseType.h" + +// Can't use CHOLMOD_NAME(drop)(0.0, S, cm). It doesn't treat complex matrices +template <class chol_type, class chol_elt, class p_type> +void +sparse_base_chol<chol_type, chol_elt, p_type>::sparse_base_chol_rep::drop_zeros + (const cholmod_sparse* S) +{ + chol_elt sik; + octave_idx_type *Sp, *Si; + chol_elt *Sx; + octave_idx_type pdest, k, ncol, p, pend; + + if (S == NULL) + return; + + Sp = static_cast<octave_idx_type *>(S->p); + Si = static_cast<octave_idx_type *>(S->i); + Sx = static_cast<chol_elt *>(S->x); + pdest = 0; + ncol = S->ncol; + + for (k = 0; k < ncol; k++) + { + p = Sp [k]; + pend = Sp [k+1]; + Sp [k] = pdest; + for (; p < pend; p++) + { + sik = Sx [p]; + if (CHOLMOD_IS_NONZERO (sik)) + { + if (p != pdest) + { + Si [pdest] = Si [p]; + Sx [pdest] = sik; + } + pdest++; + } + } + } + Sp [ncol] = pdest; +} + +template <class chol_type, class chol_elt, class p_type> +octave_idx_type +sparse_base_chol<chol_type, chol_elt, p_type>::sparse_base_chol_rep::init + (const chol_type& a, bool natural) +{ + octave_idx_type info = 0; +#ifdef HAVE_CHOLMOD + octave_idx_type a_nr = a.rows (); + octave_idx_type a_nc = a.cols (); + + if (a_nr != a_nc) + { + (*current_liboctave_error_handler) + ("SparseCHOL requires square matrix"); + return -1; + } + + cholmod_common *cm = &Common; + + // Setup initial parameters + CHOLMOD_NAME(start) (cm); + cm->prefer_zomplex = FALSE; + + double spu = Voctave_sparse_controls.get_key ("spumoni"); + if (spu == 0.) + { + cm->print = -1; + cm->print_function = NULL; + } + else + { + cm->print = (int)spu + 2; + cm->print_function =&SparseCholPrint; + } + + cm->error_handler = &SparseCholError; + cm->complex_divide = CHOLMOD_NAME(divcomplex); + cm->hypotenuse = CHOLMOD_NAME(hypot); + +#ifdef HAVE_METIS + // METIS 4.0.1 uses malloc and free, and will terminate MATLAB if it runs + // out of memory. Use CHOLMOD's memory guard for METIS, which mxMalloc's + // a huge block of memory (and then immediately mxFree's it) before calling + // METIS + cm->metis_memory = 2.0; + +#if defined(METIS_VERSION) +#if (METIS_VERSION >= METIS_VER(4,0,2)) + // METIS 4.0.2 uses function pointers for malloc and free + METIS_malloc = cm->malloc_memory; + METIS_free = cm->free_memory; + // Turn off METIS memory guard. It is not needed, because mxMalloc will + // safely terminate the mexFunction and free any workspace without killing + // all of MATLAB. + cm->metis_memory = 0.0; +#endif +#endif +#endif + + cm->final_asis = FALSE; + cm->final_super = FALSE; + cm->final_ll = TRUE; + cm->final_pack = TRUE; + cm->final_monotonic = TRUE; + cm->final_resymbol = FALSE; + + cholmod_sparse A; + cholmod_sparse *ac = &A; + double dummy; + ac->nrow = a_nr; + ac->ncol = a_nc; + + ac->p = a.cidx(); + ac->i = a.ridx(); + ac->nzmax = a.nonzero(); + ac->packed = TRUE; + ac->sorted = TRUE; + ac->nz = NULL; +#ifdef IDX_TYPE_LONG + ac->itype = CHOLMOD_LONG; +#else + ac->itype = CHOLMOD_INT; +#endif + ac->dtype = CHOLMOD_DOUBLE; + ac->stype = 1; +#ifdef OCTAVE_CHOLMOD_TYPE + ac->xtype = OCTAVE_CHOLMOD_TYPE; +#else + ac->xtype = CHOLMOD_REAL; +#endif + + if (a_nr < 1) + ac->x = &dummy; + else + ac->x = a.data(); + + // use natural ordering if no q output parameter + if (natural) + { + cm->nmethods = 1 ; + cm->method [0].ordering = CHOLMOD_NATURAL ; + cm->postorder = FALSE ; + } + + cholmod_factor *Lfactor; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + Lfactor = CHOLMOD_NAME(analyze) (ac, cm); + CHOLMOD_NAME(factorize) (ac, Lfactor, cm); + cond = CHOLMOD_NAME(rcond) (Lfactor, cm); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + minor_p = Lfactor->minor; + is_pd = cm->status == CHOLMOD_OK; + info = (is_pd ? 0 : cm->status); + + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + Lsparse = CHOLMOD_NAME(factor_to_sparse) (Lfactor, cm); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + if (minor_p > 0 && minor_p < a_nr) + { + size_t n1 = a_nr + 1; + Lsparse->p = CHOLMOD_NAME(realloc) (minor_p+1, sizeof(octave_idx_type), + Lsparse->p, &n1, cm); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CHOLMOD_NAME(reallocate_sparse) + (static_cast<octave_idx_type *>(Lsparse->p)[minor_p], Lsparse, cm); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + Lsparse->ncol = minor_p; + } + + drop_zeros (Lsparse); + + if (! natural) + { + perms.resize (a_nr); + for (octave_idx_type i = 0; i < a_nr; i++) + perms(i) = static_cast<octave_idx_type *>(Lfactor->Perm)[i]; + } + + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CHOLMOD_NAME(free_factor) (&Lfactor, cm); + CHOLMOD_NAME(finish) (cm); + CHOLMOD_NAME(print_common) (" ", cm); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; +#else + (*current_liboctave_error_handler) + ("Missing CHOLMOD. Sparse cholesky factorization disabled"); +#endif + return info; +} + +template <class chol_type, class chol_elt, class p_type> +chol_type +sparse_base_chol<chol_type, chol_elt, p_type>::L (void) const +{ +#ifdef HAVE_CHOLMOD + cholmod_sparse *m = rep->L(); + octave_idx_type nc = m->ncol; + octave_idx_type nnz = m->nzmax; + chol_type ret (m->nrow, nc, nnz); + for (octave_idx_type j = 0; j < nc+1; j++) + ret.xcidx(j) = static_cast<octave_idx_type *>(m->p)[j]; + for (octave_idx_type i = 0; i < nnz; i++) + { + ret.xridx(i) = static_cast<octave_idx_type *>(m->i)[i]; + ret.xdata(i) = static_cast<chol_elt *>(m->x)[i]; + } + return ret; +#else + return chol_type(); +#endif +} + +template <class chol_type, class chol_elt, class p_type> +p_type +sparse_base_chol<chol_type, chol_elt, p_type>:: +sparse_base_chol_rep::Q (void) const +{ +#ifdef HAVE_CHOLMOD + octave_idx_type n = Lsparse->nrow; + p_type p (n, n, n); + + for (octave_idx_type i = 0; i < n; i++) + { + p.xcidx(i) = i; + p.xridx(i) = static_cast<int>(perms(i)); + p.xdata(i) = 1; + } + p.xcidx(n) = n; + + return p; +#else + return p_type(); +#endif +} + +template <class chol_type, class chol_elt, class p_type> +chol_type +sparse_base_chol<chol_type, chol_elt, p_type>::inverse (void) const +{ + chol_type retval; +#ifdef HAVE_CHOLMOD + cholmod_sparse *m = rep->L(); + octave_idx_type n = m->ncol; + ColumnVector perms = rep->perm(); + chol_type ret; + double rcond2; + octave_idx_type info; + SparseType mattype (SparseType::Upper); + chol_type linv = L().transpose().inverse(mattype, info, rcond2, 1, 0); + + if (perms.length() == n) + { + p_type Qc = Q(); + retval = Qc * linv.transpose() * linv * Qc.transpose(); + } + else + retval = linv.transpose() * linv; +#endif + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/liboctave/sparse-base-chol.h @@ -0,0 +1,159 @@ +/* + +Copyright (C) 2005 David Bateman +Copyright (C) 1998-2005 Andy Adler + +Octave 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. + +Octave 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, Inc., 51 Franklin Street, Fifth Floor, +Boston, MA 02110-1301, USA. + +*/ + +#if !defined (octave_sparse_base_chol_h) +#define octave_sparse_base_chol_h 1 + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "oct-sparse.h" +#include "dColVector.h" + +template <class chol_type, class chol_elt, class p_type> +class +sparse_base_chol +{ +protected: + class sparse_base_chol_rep + { + public: + sparse_base_chol_rep (void) : count (1), Lsparse (NULL), + is_pd (false), minor_p (0) { } + + sparse_base_chol_rep (const chol_type& a, + const bool natural) : count (1) + { init (a, natural); } + + sparse_base_chol_rep (const chol_type& a, octave_idx_type& info, + const bool natural) : count (1) + { info = init (a, natural); } + +#ifndef HAVE_CHOLMOD + ~sparse_base_chol_rep (void) { } +#else + ~sparse_base_chol_rep (void) + { CHOLMOD_NAME(free_sparse) (&Lsparse, &Common); } +#endif + + cholmod_sparse * L (void) const { return Lsparse; } + + octave_idx_type P (void) const + { return (minor_p == static_cast<octave_idx_type>(Lsparse->ncol) ? + 0 : minor_p + 1); } + + ColumnVector perm (void) const { return perms + 1; } + + p_type Q (void) const; + + bool is_positive_definite (void) const { return is_pd; } + + double rcond (void) const { return cond; } + + int count; + + private: + cholmod_sparse *Lsparse; + + cholmod_common Common; + + bool is_pd; + + octave_idx_type minor_p; + + ColumnVector perms; + + double cond; + + octave_idx_type init (const chol_type& a, bool natural = true); + + void drop_zeros (const cholmod_sparse* S); + + // No assignment + sparse_base_chol_rep& operator = (const sparse_base_chol_rep& a); + }; + + private: + sparse_base_chol_rep *rep; + +public: + + sparse_base_chol (void) : rep (new typename + sparse_base_chol<chol_type, chol_elt, p_type>::sparse_base_chol_rep ()) { } + + sparse_base_chol (const chol_type& a, const bool n) : rep (new typename + sparse_base_chol<chol_type, chol_elt, p_type>:: + sparse_base_chol_rep (a, n)) { } + + sparse_base_chol (const chol_type& a, octave_idx_type& info, const bool n) : + rep (new typename sparse_base_chol<chol_type, chol_elt, p_type>:: + sparse_base_chol_rep (a, info, n)) { } + + sparse_base_chol (const sparse_base_chol<chol_type, chol_elt, p_type>& a) : + rep (a.rep) { rep->count++; } + + ~sparse_base_chol (void) + { + if (--rep->count <= 0) + delete rep; + } + + sparse_base_chol& operator = (const sparse_base_chol& a) + { + if (this != &a) + { + if (--rep->count <= 0) + delete rep; + + rep = a.rep; + rep->count++; + } + + return *this; + } + + chol_type L (void) const; + + chol_type R (void) const { return L().hermitian (); } + + octave_idx_type P (void) const { return rep->P(); } + + ColumnVector perm (void) const { return rep->perm(); } + + p_type Q (void) const { return rep->Q(); } + + bool is_positive_definite (void) const + { return rep->is_positive_definite(); } + + double rcond (void) const { return rep->rcond(); } + + chol_type inverse (void) const; +}; + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/liboctave/sparse-util.cc @@ -0,0 +1,57 @@ +/* + +Copyright (C) 2005 David Bateman +Copyright (C) 1998-2005 Andy Adler + +Octave 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. + +Octave 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, Inc., 51 Franklin Street, Fifth Floor, +Boston, MA 02110-1301, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <stdio.h> +#include <stdarg.h> +#include "lo-error.h" +#include "sparse-util.h" + +void +SparseCholError (int status, char *file, int line, char *message) +{ + (*current_liboctave_warning_handler)("warning %i, at line %i in file %s", + status, line, file); + + (*current_liboctave_warning_handler)(message); +} + +int +SparseCholPrint (const char *fmt, ...) +{ + va_list args; + va_start (args, fmt); + int ret = vfprintf (stderr, fmt, args); + fflush (stderr); + va_end (args); + return ret; +} + + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/liboctave/sparse-util.h @@ -0,0 +1,35 @@ +/* + +Copyright (C) 2005 David Bateman +Copyright (C) 1998-2005 Andy Adler + +Octave 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. + +Octave 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, Inc., 51 Franklin Street, Fifth Floor, +Boston, MA 02110-1301, USA. + +*/ + +#if !defined (octave_sparse_util_h) +#define octave_sparse_util_h 1 + +extern void SparseCholError (int status, char *file, int line, char *message); +extern int SparseCholPrint (const char *fmt, ...); + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- a/src/ChangeLog +++ b/src/ChangeLog @@ -1,3 +1,21 @@ +2005-10-23 David Bateman <dbateman@free.fr> + + * Makefile.in (DLD_XSRC): Add spchol.cc. + * sparse-xpow.cc (xpow): Change call to sparse inverse to include + SparseType. + * DLD-FUNCTIONS/colamd.c (Fsymbfact): Remove. + * DLD-FUNCTIONS/matrix_type.cc (Fmatrix_type): 64-bit fixes and fix + for permutation of upper/lower triangular matrices. + * DLD-FUNCTIONS/splu.cc (Fspinv): Implemtation of sparse inverse + function. + * DLD-FUNCTIONS/spchol.cc (sparse_chol): Static function for core of + the sparse cholesky factorization. + (Fspchol): New function for sparse cholesky factorization R'R. + (Fsplchol): New function for sparse cholesky factorization LL'. + (Fspcholinv): New cholesky inverse function. + (Fspchol2inv): New cholesky inverse function. + (Fsymbfact): Implementation of symbolic factorization using cholmod. + 2005-10-21 John W. Eaton <jwe@octave.org> * DLD-FUNCTIONS/gplot.l (read_until): Special case STRING.
--- a/src/DLD-FUNCTIONS/ccolamd.cc +++ b/src/DLD-FUNCTIONS/ccolamd.cc @@ -71,7 +71,7 @@ \n\ @table @code\n\ @item @var{knobs}(1)\n\ -if nonzero, the ordering is optimized for @code{lu(S(:,p)). It will be a\n\ +if nonzero, the ordering is optimized for @code{lu(S(:,p))}. It will be a\n\ poor ordering for @code{chol(@var{s}(:,@var{p})'*@var{s}(:,@var{p}))}. This\n\ is the most important knob for ccolamd.\n\ \n\ @@ -100,7 +100,7 @@ not present or empty. @code{ccolamd (@var{s},[],1:n)} returns @code{1:n}\n\ \n\ @code{@var{p} = ccolamd(@var{s})} is about the same as @code{@var{p} =\n\ -colamd(@var{s}). @var{knobs} and its default values differ. @code{colamd}\n\ +colamd(@var{s})}. @var{knobs} and its default values differ. @code{colamd}\n\ always does aggressive absorption, and it finds an ordering suitable for\n\ both @code{lu(@var{s}(:,@var{p}))} and @code{chol(@var{S}(:,@var{p})'*\n\ @var{s}(:,@var{p}))}; it cannot optimize its ordering for @code{lu(@var{s}\n\
--- a/src/DLD-FUNCTIONS/colamd.cc +++ b/src/DLD-FUNCTIONS/colamd.cc @@ -760,17 +760,6 @@ return retval; } -DEFUN_DLD (symbfact, args, nargout, - "-*- texinfo -*-\n\ -@deftypefn {Loadable Function} {[@var{count}, @var{h}, @var{parent}, @var{post}, @var{r}]} = symbfact (@var{s}, @var{typ})\n\ -\n\ -Performs a symbolic factorization analysis on the sparse matrix @var{s}.\n\ -@end deftypefn") -{ - error ("symbfact: not implemented yet"); - return octave_value (); -} - /* ;;; Local Variables: *** ;;; mode: C++ ***
--- a/src/DLD-FUNCTIONS/matrix_type.cc +++ b/src/DLD-FUNCTIONS/matrix_type.cc @@ -171,8 +171,8 @@ // XXX FIXME, why do I have to explicitly call the constructor? SparseType mattyp = SparseType (); - int nl = 0; - int nu = 0; + octave_idx_type nl = 0; + octave_idx_type nu = 0; if (error_state) error ("Matrix type must be a string"); @@ -238,7 +238,7 @@ error ("matrix_type: Invalid permutation vector"); else { - int len = perm.length (); + octave_idx_type len = perm.length (); dim_vector dv = args(0).dims (); if (len != dv(0)) @@ -247,8 +247,8 @@ { OCTAVE_LOCAL_BUFFER (octave_idx_type, p, len); - for (int i = 0; i < len; i++) - p[i] = (int) (perm (i)); + for (octave_idx_type i = 0; i < len; i++) + p[i] = (octave_idx_type) (perm (i)) - 1; if (str_typ == "upper") mattyp.mark_as_permuted (len, p);
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/spchol.cc @@ -0,0 +1,672 @@ +/* + +Copyright (C) 2005 David Bateman +Copyright (C) 1998-2005 Andy Adler + +Octave 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. + +Octave 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, Inc., 51 Franklin Street, Fifth Floor, +Boston, MA 02110-1301, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "defun-dld.h" +#include "error.h" +#include "gripes.h" +#include "oct-obj.h" +#include "utils.h" + +#include "SparseCmplxCHOL.h" +#include "SparsedbleCHOL.h" +#include "ov-re-sparse.h" +#include "ov-cx-sparse.h" +#include "oct-spparms.h" +#include "sparse-util.h" + +static octave_value_list +sparse_chol (const octave_value_list& args, const int nargout, + const std::string& name, const bool LLt) +{ + octave_value_list retval; + int nargin = args.length (); + + if (nargin != 1 || nargout > 3) + { + print_usage (name); + return retval; + } + + octave_value arg = args(0); + + octave_idx_type nr = arg.rows (); + octave_idx_type nc = arg.columns (); + bool natural = (nargout != 3); + + int arg_is_empty = empty_arg (name.c_str(), nr, nc); + + if (arg_is_empty < 0) + return retval; + if (arg_is_empty > 0) + return octave_value (Matrix ()); + + if (arg.is_real_type ()) + { + SparseMatrix m = arg.sparse_matrix_value (); + + if (! error_state) + { + octave_idx_type info; + SparseCHOL fact (m, info, natural); + if (nargout == 3) + retval(2) = fact.Q(); + + if (nargout > 1 || info == 0) + { + retval(1) = fact.P(); + if (LLt) + retval(0) = fact.L(); + else + retval(0) = fact.R(); + } + else + error ("%s: matrix not positive definite", name.c_str()); + } + } + else if (arg.is_complex_type ()) + { + SparseComplexMatrix m = arg.sparse_complex_matrix_value (); + + if (! error_state) + { + octave_idx_type info; + SparseComplexCHOL fact (m, info, natural); + + if (nargout == 3) + retval(2) = fact.Q(); + + if (nargout > 1 || info == 0) + { + retval(1) = fact.P(); + if (LLt) + retval(0) = fact.L(); + else + retval(0) = fact.R(); + } + else + error ("%s: matrix not positive definite", name.c_str()); + } + } + else + gripe_wrong_type_arg (name.c_str(), arg); + + return retval; +} + +// PKG_ADD: dispatch ("chol", "spchol", "sparse matrix") +// PKG_ADD: dispatch ("chol", "spchol", "sparse complex matrix") +// PKG_ADD: dispatch ("chol", "spchol", "sparse bool matrix") +DEFUN_DLD (spchol, args, nargout, + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{r} =} spchol (@var{a})\n\ +@deftypefnx {Loadable Function} {[@var{r}, @var{p}] =} spchol (@var{a})\n\ +@deftypefnx {Loadable Function} {[@var{r}, @var{p}, @var{q}] =} spchol (@var{a})\n\ +@cindex Cholesky factorization\n\ +Compute the Cholesky factor, @var{r}, of the symmetric positive definite\n\ +sparse matrix @var{a}, where\n\ +@iftex\n\ +@tex\n\ +$ R^T R = A $.\n\ +@end tex\n\ +@end iftex\n\ +@ifinfo\n\ +\n\ +@example\n\ +r' * r = a.\n\ +@end example\n\ +@end ifinfo\n\ +\n\ +If called with 2 or more outputs @var{p} is the 0 when @var{r} is positive\n\ +definite and @var{p} is a positive integer otherwise.\n\ +\n\ +If called with 3 outputs then a sparsity preserving row/column permutation\n\ +is applied to @var{a} prior to the factorization. That is @var{r}\n\ +is the factorization of @code{@var{a}(@var{q},@var{q})} such that\n\ +@iftex\n\ +@tex\n\ +$ R^T R = Q A Q^T$.\n\ +@end tex\n\ +@end iftex\n\ +@ifinfo\n\ +\n\ +@example\n\ +r' * r = q * a * q'.\n\ +@end example\n\ +@end ifinfo\n\ +\n\ +Note that @code{splchol} factorizations is faster and use less memory.\n\ +@end deftypefn\n\ +@seealso{spcholinv, spchol2inv, splchol}") +{ + return sparse_chol (args, nargout, "spchol", false); +} + +// PKG_ADD: dispatch ("lchol", "splchol", "sparse matrix") +// PKG_ADD: dispatch ("lchol", "splchol", "sparse complex matrix") +// PKG_ADD: dispatch ("lchol", "splchol", "sparse bool matrix") +DEFUN_DLD (splchol, args, nargout, + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{l} =} splchol (@var{a})\n\ +@deftypefnx {Loadable Function} {[@var{l}, @var{p}] =} splchol (@var{a})\n\ +@deftypefnx {Loadable Function} {[@var{l}, @var{p}, @var{q}] =} splchol (@var{a})\n\ +@cindex Cholesky factorization\n\ +Compute the Cholesky factor, @var{l}, of the symmetric positive definite\n\ +sparse matrix @var{a}, where\n\ +@iftex\n\ +@tex\n\ +$ L L^T = A $.\n\ +@end tex\n\ +@end iftex\n\ +@ifinfo\n\ +\n\ +@example\n\ +l * l' = a.\n\ +@end example\n\ +@end ifinfo\n\ +\n\ +If called with 2 or more outputs @var{p} is the 0 when @var{l} is positive\n\ +definite and @var{l} is a positive integer otherwise.\n\ +\n\ +If called with 3 outputs that a sparsity preserving row/column permutation\n\ +is applied to @var{a} prior to the factorization. That is @var{l}\n\ +is the factorization of @code{@var{a}(@var{q},@var{q})} such that\n\ +@iftex\n\ +@tex\n\ +$ L R^T = A (Q, Q)$.\n\ +@end tex\n\ +@end iftex\n\ +@ifinfo\n\ +\n\ +@example\n\ +r * r' = a (q, q).\n\ +@end example\n\ +@end ifinfo\n\ +\n\ +Note that @code{splchol} factorizations is faster and use less memory\n\ +than @code{spchol}. @code{splchol(@var{a})} is equivalent to\n\ +@code{spchol(@var{a})'}.\n\ +@end deftypefn\n\ +@seealso{spcholinv, spchol2inv, splchol}") +{ + return sparse_chol (args, nargout, "splchol", true); +} + +// PKG_ADD: dispatch ("cholinv", "spcholinv", "sparse matrix") +// PKG_ADD: dispatch ("cholinv", "spcholinv", "sparse complex matrix") +// PKG_ADD: dispatch ("cholinv", "spcholinv", "sparse bool matrix") +DEFUN_DLD (spcholinv, args, , + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {} spcholinv (@var{a})\n\ +Use the Cholesky factorization to compute the inverse of the\n\ +sparse symmetric positive definite matrix @var{a}.\n\ +@seealso{spchol, spchol2inv}\n\ +@end deftypefn") +{ + octave_value retval; + + int nargin = args.length (); + + if (nargin == 1) + { + octave_value arg = args(0); + + octave_idx_type nr = arg.rows (); + octave_idx_type nc = arg.columns (); + + if (nr == 0 || nc == 0) + retval = Matrix (); + else + { + if (arg.is_real_type ()) + { + SparseMatrix m = arg.sparse_matrix_value (); + + if (! error_state) + { + octave_idx_type info; + SparseCHOL chol (m, info); + if (info == 0) + retval = chol.inverse (); + else + error ("spcholinv: matrix not positive definite"); + } + } + else if (arg.is_complex_type ()) + { + SparseComplexMatrix m = arg.sparse_complex_matrix_value (); + + if (! error_state) + { + octave_idx_type info; + SparseComplexCHOL chol (m, info); + if (info == 0) + retval = chol.inverse (); + else + error ("spcholinv: matrix not positive definite"); + } + } + else + gripe_wrong_type_arg ("spcholinv", arg); + } + } + else + print_usage ("spcholinv"); + + return retval; +} + +// PKG_ADD: dispatch ("chol2inv", "spchol2inv", "sparse matrix") +// PKG_ADD: dispatch ("chol2inv", "spchol2inv", "sparse complex matrix") +// PKG_ADD: dispatch ("chol2inv", "spchol2inv", "sparse bool matrix") +DEFUN_DLD (spchol2inv, args, , + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {} spchol2inv (@var{u})\n\ +Invert a sparse symmetric, positive definite square matrix from its\n\ +Cholesky decomposition, @var{u}. Note that @var{u} should be an\n\ +upper-triangular matrix with positive diagonal elements.\n\ +@code{chol2inv (@var{u})} provides @code{inv (@var{u}'*@var{u})} but\n\ +it is much faster than using @code{inv}.\n\ +@seealso{spchol, spcholinv}\n\ +@end deftypefn") +{ + octave_value retval; + + int nargin = args.length (); + + if (nargin == 1) + { + octave_value arg = args(0); + + octave_idx_type nr = arg.rows (); + octave_idx_type nc = arg.columns (); + + if (nr == 0 || nc == 0) + retval = Matrix (); + else + { + if (arg.is_real_type ()) + { + SparseMatrix r = arg.sparse_matrix_value (); + + if (! error_state) + retval = chol2inv (r); + } + else if (arg.is_complex_type ()) + { + SparseComplexMatrix r = arg.sparse_complex_matrix_value (); + + if (! error_state) + retval = chol2inv (r); + } + else + gripe_wrong_type_arg ("spchol2inv", arg); + } + } + else + print_usage ("spchol2inv"); + + return retval; +} + +DEFUN_DLD (symbfact, args, nargout, + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {[@var{count}, @var{h}, @var{parent}, @var{post}, @var{r}]} = symbfact (@var{s}, @var{typ}, @var{mode})\n\ +\n\ +Performs a symbolic factorization analysis on the sparse matrix @var{s}.\n\ +Where\n\ +\n\ +@table @asis\n\ +@item @var{s}\n\ +@var{s} is a complex or real sparse matrix.\n\ +\n\ +@item @var{typ}\n\ +Is the type of the factorization and can be one of\n\ +\n\ +@table @code\n\ +@item sym\n\ +Factorize @var{s}. This is the default.\n\ +\n\ +@item col\n\ +Factorize @code{@var{s}' * @var{s}}.\n\ +@item row\n\ +Factorize @code{@var{s} * @var{s}'}.\n\ +@item lo\n\ +Factorize @code{@var{s}'}\n\ +@end table\n\ +\n\ +@item @var{mode}\n\ +The default is to return the Cholesky factorization for @var{r}, and if\n\ +@var{mode} is 'L', the conjugate transpose of the Choleksy factorization\n\ +is returned. The conjugate transpose version is faster and uses less\n\ +memory, but returns the same values for @var{count}, @var{h}, @var{parent}\n\ +and @var{post} outputs.\n\ +@end table\n\ +\n\ +The output variables are\n\ +\n\ +@table @asis\n\ +@item @var{count}\n\ +The row counts of the Cholesky factorization as determined by @var{typ}.\n\ +\n\ +@item @var{h}\n\ +The height of the elimination tree.\n\ +\n\ +@item @var{parent}\n\ +The elimination tree itself.\n\ +\n\ +@item @var{post}\n\ +A sparse boolean matrix whose structure is that of the Cholesky\n\ +factorization as determined by @var{typ}.\n\ +@end table\n\ +@end deftypefn") +{ + octave_value_list retval; + int nargin = args.length (); + + if (nargin < 1 || nargin > 3 || nargout > 5) + { + print_usage ("symbfact"); + return retval; + } + + cholmod_common Common; + cholmod_common *cm = &Common; + CHOLMOD_NAME(start) (cm); + + double spu = Voctave_sparse_controls.get_key ("spumoni"); + if (spu == 0.) + { + cm->print = -1; + cm->print_function = NULL; + } + else + { + cm->print = (int)spu + 2; + cm->print_function =&SparseCholPrint; + } + + cm->error_handler = &SparseCholError; + cm->complex_divide = CHOLMOD_NAME(divcomplex); + cm->hypotenuse = CHOLMOD_NAME(hypot); + +#ifdef HAVE_METIS + // METIS 4.0.1 uses malloc and free, and will terminate MATLAB if it runs + // out of memory. Use CHOLMOD's memory guard for METIS, which mxMalloc's + // a huge block of memory (and then immediately mxFree's it) before calling + // METIS + cm->metis_memory = 2.0; + +#if defined(METIS_VERSION) +#if (METIS_VERSION >= METIS_VER(4,0,2)) + // METIS 4.0.2 uses function pointers for malloc and free + METIS_malloc = cm->malloc_memory; + METIS_free = cm->free_memory; + // Turn off METIS memory guard. It is not needed, because mxMalloc will + // safely terminate the mexFunction and free any workspace without killing + // all of MATLAB. + cm->metis_memory = 0.0; +#endif +#endif +#endif + + double dummy; + cholmod_sparse Astore; + cholmod_sparse *A = &Astore; + A->packed = TRUE; + A->sorted = TRUE; + A->nz = NULL; +#ifdef IDX_TYPE_LONG + A->itype = CHOLMOD_LONG; +#else + A->itype = CHOLMOD_INT; +#endif + A->dtype = CHOLMOD_DOUBLE; + A->stype = 1; + A->x = &dummy; + + if (args(0).is_real_type ()) + { + const SparseMatrix a = args(0).sparse_matrix_value(); + A->nrow = a.rows(); + A->ncol = a.cols(); + A->p = a.cidx(); + A->i = a.ridx(); + A->nzmax = a.nonzero(); + A->xtype = CHOLMOD_REAL; + + if (a.rows() > 0 && a.cols() > 0) + A->x = a.data(); + } + else if (args(0).is_complex_type ()) + { + const SparseComplexMatrix a = args(0).sparse_complex_matrix_value(); + A->nrow = a.rows(); + A->ncol = a.cols(); + A->p = a.cidx(); + A->i = a.ridx(); + A->nzmax = a.nonzero(); + A->xtype = CHOLMOD_COMPLEX; + + if (a.rows() > 0 && a.cols() > 0) + A->x = a.data(); + } + else + gripe_wrong_type_arg ("symbfact", arg(0)); + + octave_idx_type coletree = FALSE; + octave_idx_type n = A->nrow; + + if (nargin > 1) + { + char ch; + std::string str = args(1).string_value(); + ch = tolower (str.c_str()[0]); + if (ch == 'r') + A->stype = 0; + else if (ch == 'c') + { + n = A->ncol; + coletree = TRUE; + A->stype = 0; + } + else if (ch == 's') + A->stype = 1; + else if (ch == 's') + A->stype = -1; + else + error ("Unrecognized typ in symbolic factorization"); + } + + if (A->stype && A->nrow != A->ncol) + error ("Matrix must be square"); + + if (!error_state) + { + OCTAVE_LOCAL_BUFFER (octave_idx_type, Parent, n); + OCTAVE_LOCAL_BUFFER (octave_idx_type, Post, n); + OCTAVE_LOCAL_BUFFER (octave_idx_type, ColCount, n); + OCTAVE_LOCAL_BUFFER (octave_idx_type, First, n); + OCTAVE_LOCAL_BUFFER (octave_idx_type, Level, n); + + cholmod_sparse *F = CHOLMOD_NAME(transpose) (A, 0, cm); + cholmod_sparse *Aup, *Alo; + + if (A->stype == 1 || coletree) + { + Aup = A ; + Alo = F ; + } + else + { + Aup = F ; + Alo = A ; + } + + CHOLMOD_NAME(etree) (Aup, Parent, cm); + + if (cm->status < CHOLMOD_OK) + { + error("matrix corrupted"); + goto symbfact_error; + } + + if (CHOLMOD_NAME(postorder) (Parent, n, NULL, Post, cm) != n) + { + error("postorder failed"); + goto symbfact_error; + } + + CHOLMOD_NAME(rowcolcounts) (Alo, NULL, 0, Parent, Post, NULL, + ColCount, First, Level, cm); + + if (cm->status < CHOLMOD_OK) + { + error("matrix corrupted"); + goto symbfact_error; + } + + if (nargout > 4) + { + cholmod_sparse *A1, *A2; + + if (A->stype == 1) + { + A1 = A; + A2 = NULL; + } + else if (A->stype == -1) + { + A1 = F; + A2 = NULL; + } + else if (coletree) + { + A1 = F; + A2 = A; + } + else + { + A1 = A; + A2 = F; + } + + // count the total number of entries in L + octave_idx_type lnz = 0 ; + for (octave_idx_type j = 0 ; j < n ; j++) + lnz += ColCount [j] ; + + + // allocate the output matrix L (pattern-only) + SparseBoolMatrix L (n, n, lnz); + + // initialize column pointers + lnz = 0; + for (octave_idx_type j = 0 ; j < n ; j++) + { + L.xcidx(j) = lnz; + lnz += ColCount [j]; + } + L.xcidx(n) = lnz; + + + /* create a copy of the column pointers */ + octave_idx_type *W = First; + for (octave_idx_type j = 0 ; j < n ; j++) + W [j] = L.xcidx(j); + + // get workspace for computing one row of L + cholmod_sparse *R = cholmod_allocate_sparse (n, 1, n, FALSE, TRUE, + 0, CHOLMOD_PATTERN, cm); + octave_idx_type *Rp = static_cast<octave_idx_type *>(R->p); + octave_idx_type *Ri = static_cast<octave_idx_type *>(R->i); + + // compute L one row at a time + for (octave_idx_type k = 0 ; k < n ; k++) + { + // get the kth row of L and store in the columns of L + cholmod_row_subtree (A1, A2, k, Parent, R, cm) ; + for (octave_idx_type p = 0 ; p < Rp [1] ; p++) + L.xridx (W [Ri [p]]++) = k ; + + // add the diagonal entry + L.xridx (W [k]++) = k ; + } + + // free workspace + cholmod_free_sparse (&R, cm) ; + + + // transpose L to get R, or leave as is + if (nargin < 3) + L = L.transpose (); + + // fill numerical values of L with one's + for (octave_idx_type p = 0 ; p < lnz ; p++) + L.xdata(p) = true; + + retval(4) = L; + } + + ColumnVector tmp (n); + if (nargout > 3) + { + for (octave_idx_type i = 0; i < n; i++) + tmp(i) = Post[i] + 1; + retval(3) = tmp; + } + + if (nargout > 2) + { + for (octave_idx_type i = 0; i < n; i++) + tmp(i) = Parent[i] + 1; + retval(2) = tmp; + } + + if (nargout > 1) + { + /* compute the elimination tree height */ + octave_idx_type height = 0 ; + for (int i = 0 ; i < n ; i++) + height = (height > Level[i] ? height : Level[i]); + height++ ; + retval(1) = (double)height; + } + + for (octave_idx_type i = 0; i < n; i++) + tmp(i) = ColCount[i]; + retval(0) = tmp; + } + + symbfact_error: + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +
--- a/src/DLD-FUNCTIONS/splu.cc +++ b/src/DLD-FUNCTIONS/splu.cc @@ -398,23 +398,102 @@ DEFUN_DLD (spinv, args, nargout, "-*- texinfo -*-\n\ @deftypefn {Loadable Function} {[@var{x}, @var{rcond}] = } spinv (@var{a}, @var{Q})\n\ -@deftypefnx {Loadable Function} {[@var{x}, @var{rcond}, @var{Q}] = } spinv (@var{a}, @var{Q})\n\ -Compute the inverse of the square matrix @var{a}. Return an estimate\n\ +Compute the inverse of the sparse square matrix @var{a}. Return an estimate\n\ of the reciprocal condition number if requested, otherwise warn of an\n\ ill-conditioned matrix if the reciprocal condition number is small.\n\ +This function takes advantage of the sparsity of the matrix to accelerate\n\ +the calculation of the inverse.\n\ \n\ -An optional second input argument @var{Q} is the optional pre-ordering of\n\ -the matrix, such that @code{@var{x} = inv (@var{a} (:, @var{Q}))}. @var{Q}\n\ -can equally be a matrix, in which case @code{@var{x} = inv (@var{a} *\n\ -@var{Q}))}.\n\ -\n\ -If a third output argument is given then the permuations to achieve a sparse\n\ -inverse are returned. It is not required that the return column permutations\n\ -@var{Q} and the same as the user supplied permutations\n\ +In general @var{x} will be a full matrix, and so if possible forming the\n\ +inverse of a sparse matrix should be avoided. It is significantly more\n\ +accurate and faster to do @code{@var{y} = @var{a} \\ @var{b}}, rather\n\ +than @code{@var{y} = spinv (@var{a}) * @var{b}}.\n\ @end deftypefn") { - error ("spinv: not implemented yet"); - return octave_value (); + octave_value_list retval; + + int nargin = args.length (); + + if (nargin != 1) + { + print_usage ("inv"); + return retval; + } + + octave_value arg = args(0); + const octave_value& rep = arg.get_rep (); + + octave_idx_type nr = arg.rows (); + octave_idx_type nc = arg.columns (); + + int arg_is_empty = empty_arg ("spinverse", nr, nc); + + if (arg_is_empty < 0) + return retval; + else if (arg_is_empty > 0) + return octave_value (Matrix ()); + + if (nr != nc) + { + gripe_square_matrix_required ("spinverse"); + return retval; + } + + if (arg.is_real_type ()) + { + + SparseMatrix m = arg.sparse_matrix_value (); + + if (! error_state) + { + octave_idx_type info; + double rcond = 0.0; + SparseType mattyp = ((octave_sparse_matrix &)rep).sparse_type (); + SparseMatrix result = m.inverse (mattyp, info, rcond, 1); + ((octave_sparse_matrix &)(arg.get_rep())).sparse_type (mattyp); + + if (nargout > 1) + retval(1) = rcond; + + retval(0) = result; + + volatile double xrcond = rcond; + xrcond += 1.0; + if (nargout < 2 && (info == -1 || xrcond == 1.0)) + warning ("spinverse: matrix singular to machine precision,\ + rcond = %g", rcond); + } + } + else if (arg.is_complex_type ()) + { + SparseComplexMatrix m = arg.sparse_complex_matrix_value (); + + if (! error_state) + { + octave_idx_type info; + double rcond = 0.0; + + SparseType mattyp = + ((octave_sparse_complex_matrix &)rep).sparse_type (); + SparseComplexMatrix result = m.inverse (mattyp, info, rcond, 1); + ((octave_sparse_matrix &)rep).sparse_type (mattyp); + + if (nargout > 1) + retval(1) = rcond; + + retval(0) = result; + + volatile double xrcond = rcond; + xrcond += 1.0; + if (nargout < 2 && (info == -1 || xrcond == 1.0)) + warning ("spinverse: matrix singular to machine precision,\ + rcond = %g", rcond); + } + } + else + gripe_wrong_type_arg ("spinverse", arg); + + return retval; } /*
--- a/src/Makefile.in +++ b/src/Makefile.in @@ -49,8 +49,8 @@ getpwent.cc getrusage.cc givens.cc hess.cc inv.cc kron.cc \ lpsolve.cc lsode.cc lu.cc luinc.cc matrix_type.cc minmax.cc \ pinv.cc qr.cc quad.cc qz.cc rand.cc schur.cc sort.cc sparse.cc \ - spdet.cc spkron.cc splu.cc spparms.cc sqrtm.cc svd.cc syl.cc \ - time.cc gplot.l __glpk__.cc __qp__.cc + spchol.cc spdet.cc spkron.cc splu.cc spparms.cc sqrtm.cc svd.cc \ + syl.cc time.cc gplot.l __glpk__.cc __qp__.cc DLD_SRC := $(addprefix DLD-FUNCTIONS/, $(DLD_XSRC))
--- a/src/sparse-xpow.cc +++ b/src/sparse-xpow.cc @@ -90,8 +90,9 @@ octave_idx_type info; double rcond = 0.0; + SparseType mattyp (a); - atmp = a.inverse (info, rcond, 1); + atmp = a.inverse (mattyp, info, rcond, 1); if (info == -1) warning ("inverse: matrix singular to machine\ @@ -162,8 +163,9 @@ octave_idx_type info; double rcond = 0.0; + SparseType mattyp (a); - atmp = a.inverse (info, rcond, 1); + atmp = a.inverse (mattyp, info, rcond, 1); if (info == -1) warning ("inverse: matrix singular to machine\