add pivoting Cholesky from LAPACK 3.2

git-svn-id: https://svn.r-project.org/R/trunk@60504 00db46b3-68df-0310-9c12-caf00c1e9a41

doc/NEWS.Rd

@@ -194,6 +194,9 @@
       
       \item More matrix algebra functions (e.g. \code{chol()} and
       \code{solve()}) accept logical matrices (and coerce to numeric).
+      
+      \item \code{chol(pivot = TRUE, LINPACK = FALSE)} is now available
+      using LAPACK 3.2 subroutine \code{DPSTRF}.
     }
   }
 
@@ -457,6 +460,8 @@
       implementations can be slow or flawed.  Simply undefine
       \code{HAVE_LONG_DOUBLE} in \file{src/include/config.h} after
       configuration.
+
+      \item If an external LAPACK is used, it must be version 3.2 or later.
     }
   }
 

doc/manual/R-admin.texi

@@ -3217,16 +3217,18 @@ run-time library path or @command{ld.so} cache).
 Provision is made for using an external LAPACK library, principally to
 cope with @acronym{BLAS} libraries which contain a copy of LAPACK (such
 as @code{sunperf} on Solaris, @code{vecLib} on OS X and ACML on
-@cputype{ix86}/@cputype{x86_64} Linux).  However, the likely performance
-gains are thought to be small (and may be negative), and the default is
-not to search for a suitable LAPACK library, and this is definitely
-@strong{not} recommended.  You can specify a specific LAPACK library or
-a search for a generic library by the configuration option
-@option{--with-lapack}.  The default for @option{--with-lapack} is to
-check the @acronym{BLAS} library and then look for an external library
-@samp{-llapack}.  Sites searching for the fastest possible linear
-algebra may want to build a LAPACK library using the ATLAS-optimized
-subset of LAPACK.  To do so specify something like
+@cputype{ix86}/@cputype{x86_64} Linux).  At least LAPACK version 3.2 is
+required. 
+
+However, the likely performance gains are thought to be small (and may
+be negative), and the default is not to search for a suitable LAPACK
+library, and this is definitely @strong{not} recommended.  You can
+specify a specific LAPACK library or a search for a generic library by
+the configuration option @option{--with-lapack}.  The default for
+@option{--with-lapack} is to check the @acronym{BLAS} library and then
+look for an external library @samp{-llapack}.  Sites searching for the
+fastest possible linear algebra may want to build a LAPACK library using
+the ATLAS-optimized subset of LAPACK.  To do so specify something like
 
 @example
 --with-lapack="-L/path/to/ATLAS/libs -llapack -lcblas"
@@ -3245,7 +3247,7 @@ Since ACML contains a full LAPACK, if selected as the @acronym{BLAS} it
 can be used as the LAPACK @emph{via} @option{--with-lapack}.
 
 If you do use @option{--with-lapack}, be aware of potential problems
-with bugs in the LAPACK 3.0 sources (or in the posted corrections to those
+with bugs in the LAPACK sources (or in the posted corrections to those
 sources).  In particular, bugs in @code{DGEEV} and @code{DGESDD} have
 resulted in error messages such as
 
@@ -3255,9 +3257,7 @@ DGEBRD gave error code -10
 
 @noindent
 .  Other potential problems are incomplete versions of the libraries,
-seen several times over the years.  For problems compiling LAPACK using
-recent versions of @code{gcc} on @cputype{ix86} Linux, see @ref{New
-platforms}.
+seen several times in Linux distributions over the years.
 
 Please @strong{do} bear in mind that using @option{--with-lapack} is
 `definitely @strong{not} recommended': it is provided @strong{only}

src/include/R_ext/Lapack.h

@@ -2767,6 +2767,11 @@ La_extern void
 F77_NAME(dtzrzf)(int *m, int *n, double *a, int *
 	lda, double *tau, double *work, int *lwork, int *info);
 
+La_extern void
+F77_NAME(dpstrf)(const char* uplo, const int* n,
+		 double* a, const int* lda, int* piv, int* rank,
+		 double* tol, double *work, int* info);
+
 
 La_extern int
 F77_NAME(lsame)(char *ca, char *cb);

src/library/base/R/chol.R

@@ -18,7 +18,7 @@
 
 chol <- function(x, ...) UseMethod("chol")
 
-chol.default <- function(x, pivot = FALSE, LINPACK = pivot, ...)
+chol.default <- function(x, pivot = FALSE, LINPACK = pivot, tol = -1, ...)
 {
     if (is.complex(x))
         stop("complex matrices not permitted at present")
@@ -30,7 +30,10 @@ chol.default <- function(x, pivot = FALSE, LINPACK = pivot, ...)
 	if(length(x) != 1L) stop("non-matrix argument to 'chol'")
 	n <- 1L
     }
+
+    ## FIXME add pivoting version using DPSTRF (added in LAPACK 3.2)
     if(!pivot && !LINPACK) return(.Internal(La_chol(as.matrix(x))))
+    if(!LINPACK) return(.Internal(La_chol_piv(as.matrix(x), tol)))
 
     ## sanity checks
     n <- as.integer(n)

src/library/base/man/chol.Rd

@@ -14,15 +14,17 @@
 \usage{
 chol(x, \dots)
 
-\method{chol}{default}(x, pivot = FALSE,  LINPACK = pivot, \dots)
+\method{chol}{default}(x, pivot = FALSE,  LINPACK = pivot, tol = -1, \dots)
 }
 \arguments{
   \item{x}{an object for which a method exists.  The default method
     applies to numeric (or logical) symmetric, positive-definite matrices.}
   \item{\dots}{arguments to be based to or from methods.}
   \item{pivot}{Should pivoting be used?}
-  \item{LINPACK}{logical.  Should LINPACK be used in the non-pivoting
-    case (for compatibility with \R < 1.7.0)?}
+  \item{LINPACK}{logical.  Should LINPACK be used? (For compatibility
+    with \R < 1.7.0 in the non-pivoting case.)}
+  \item{tol}{A numeric tolerance for use with \code{pivot = TRUE,
+      LINPACK = FALSE)}.}
 }
 \value{
   The upper triangular factor of the Choleski decomposition, i.e., the
@@ -35,15 +37,15 @@ chol(x, \dots)
   \code{chol} is generic: the description here applies to the default
   method.
   
-  This is an interface to the LAPACK routine DPOTRF and the
+  This is an interface to the LAPACK routines DPOTRF and DPSTRF and the
   LINPACK routines DPOFA and DCHDC.
 
   Note that only the upper triangular part of \code{x} is used, so
   that \eqn{R'R = x} when \code{x} is symmetric.
 
   If \code{pivot = FALSE} and \code{x} is not non-negative definite an
   error occurs.  If \code{x} is positive semi-definite (i.e., some zero
-  eigenvalues) an error will also occur, as a numerical tolerance is used.
+  eigenvalues) an error will also occur as a numerical tolerance is used.
 
   If \code{pivot = TRUE}, then the Choleski decomposition of a positive
   semi-definite \code{x} can be computed.  The rank of \code{x} is
@@ -55,6 +57,12 @@ chol(x, \dots)
   or, alternatively, \code{t(Q) \%*\% Q} equals \code{x[pivot,
     pivot]}.  See the examples.
 
+  Pivoting with LAPACK requires LAPACK >= 3.2 and was added in \R
+  2.16.0.  The value of \code{tol} is passed to LAPACK, with negative
+  values selecting the default tolerance of (usually) \code{nrow(x) *
+  .Machine$double.neg.eps * max(diag(x)}.  The algorithm terminates once
+  the pivot is less than \code{tol}.
+
   The LINPACK interface is restricted to matrices \code{x} with less
   than \eqn{2^{31}}{2^31} elements.
 }
@@ -105,6 +113,7 @@ crossprod(cm)  #-- = 'm'
 # now for something positive semi-definite
 x <- matrix(c(1:5, (1:5)^2), 5, 2)
 x <- cbind(x, x[, 1] + 3*x[, 2])
+colnames(x) <- letters[20:22]
 m <- crossprod(x)
 qr(m)$rank # is 2, as it should be
 
@@ -117,6 +126,12 @@ try(chol(m))
 pivot <- attr(Q, "pivot")
 crossprod(Q[, order(pivot)]) # recover m
 
+(Q <- chol(m, TRUE, FALSE)) # NB rank is correct
+## we can use this by
+pivot <- attr(Q, "pivot")
+crossprod(Q[, order(pivot)]) # recover m
+
+
 ## now for a non-positive-definite matrix
 ( m <- matrix(c(5,-5,-5,3),2,2) )
 try(chol(m))  # fails

src/modules/lapack/Lapack.c

@@ -945,6 +945,51 @@ static SEXP La_chol(SEXP A)
     return ans;
 }
 
+static SEXP La_chol_piv(SEXP A, SEXP stol)
+{
+    if (!isMatrix(A)) error(_("'a' must be a numeric matrix"));
+    double tol = asReal(stol);
+
+    SEXP ans = PROTECT((TYPEOF(A) == REALSXP) ? duplicate(A):
+		       coerceVector(A, REALSXP));
+    SEXP adims = getAttrib(A, R_DimSymbol);
+    if (TYPEOF(adims) != INTSXP) error("non-integer dims");
+    int m = INTEGER(adims)[0], n = INTEGER(adims)[1];
+
+    if (m != n) error(_("'a' must be a square matrix"));
+    if (m <= 0) error(_("'a' must have dims > 0"));
+    size_t N = n;
+    for (int j = 0; j < n; j++) 	/* zero the lower triangle */
+	for (int i = j+1; i < n; i++) REAL(ans)[i + N * j] = 0.;
+
+    SEXP piv = PROTECT(allocVector(INTSXP, m));
+    int *ip = INTEGER(piv);
+    double *work = (double *) R_alloc(2 * (size_t)m, sizeof(double));
+    int rank, info;
+    F77_CALL(dpstrf)("U", &m, REAL(ans), &m, ip, &rank, &tol, work, &info);
+    if (info != 0) {
+	if (info > 0)
+	    warning(_("the matrix is either rank-deficient or indefinite"));
+	else
+	    error(_("argument %d of Lapack routine %s had invalid value"),
+		  -info, "dpotrf");
+    }
+    setAttrib(ans, install("pivot"), piv);
+    setAttrib(ans, install("rank"), ScalarInteger(rank));
+    SEXP cn, dn = getAttrib(ans, R_DimNamesSymbol);
+    if (!isNull(dn) && !isNull(cn = VECTOR_ELT(dn, 1))) {
+	// need to pivot the colnames
+	SEXP dn2 = PROTECT(duplicate(dn));
+	SEXP cn2 = VECTOR_ELT(dn2, 1);
+	for(int i = 0; i < m; i++) 
+	    SET_STRING_ELT(cn2, i, STRING_ELT(cn, ip[i] - 1)); // base 1
+	setAttrib(ans, R_DimNamesSymbol, dn2);
+	UNPROTECT(1);
+    }
+    UNPROTECT(2);
+    return ans;
+}
+
 static SEXP La_chol2inv(SEXP A, SEXP size)
 {
     int sz = asInteger(size);
@@ -1289,6 +1334,7 @@ static SEXP mod_do_lapack(SEXP call, SEXP op, SEXP args, SEXP env)
 
     case 200: ans = La_chol(CAR(args)); break;
     case 201: ans = La_chol2inv(CAR(args), CADR(args)); break;
+    case 202: ans = La_chol_piv(CAR(args), CADR(args)); break;
 
     case 300: ans = qr_coef_real(CAR(args), CADR(args)); break;
     case 301: ans = qr_qy_real(CAR(args), CADR(args), CADDR(args)); break;