[med-svn] [r-cran-nnls] 03/04: Import upstream version 1.4.
Alba Crespi
albac-guest at moszumanska.debian.org
Mon May 18 18:48:46 UTC 2015
This is an automated email from the git hooks/post-receive script.
albac-guest pushed a commit to annotated tag upstream/1.4
in repository r-cran-nnls.
commit 697cbbb63f5e8e6d305ac0c334ad07efa86cf225
Author: Alba Crespi <crespialba+debian at gmail.com>
Date: Mon May 18 19:05:06 2015 +0100
Import upstream version 1.4.
---
DESCRIPTION | 14 ++
MD5 | 10 +
NAMESPACE | 13 ++
R/nnls.R | 94 ++++++++++
R/zzz.R | 6 +
inst/COPYRIGHTS | 16 ++
man/nnls-package.Rd | 31 ++++
man/nnls.Rd | 155 ++++++++++++++++
man/nnnpls.Rd | 161 ++++++++++++++++
src/lawson_hanson_nnls.f | 474 +++++++++++++++++++++++++++++++++++++++++++++++
src/nnnpls.f | 353 +++++++++++++++++++++++++++++++++++
11 files changed, 1327 insertions(+)
diff --git a/DESCRIPTION b/DESCRIPTION
new file mode 100644
index 0000000..a496fed
--- /dev/null
+++ b/DESCRIPTION
@@ -0,0 +1,14 @@
+Package: nnls
+Type: Package
+Title: The Lawson-Hanson algorithm for non-negative least squares
+ (NNLS)
+Version: 1.4
+Author: Katharine M. Mullen and Ivo H. M. van Stokkum
+Maintainer: Katharine Mullen <mullenkate at gmail.com>
+Description: An R interface to the Lawson-Hanson implementation of an
+ algorithm for non-negative least squares (NNLS). Also allows
+ the combination of non-negative and non-positive constraints.
+License: GPL (>= 2)
+Packaged: 2012-03-19 16:06:34 UTC; kmullen
+Repository: CRAN
+Date/Publication: 2012-03-19 16:28:59
diff --git a/MD5 b/MD5
new file mode 100644
index 0000000..d9b4602
--- /dev/null
+++ b/MD5
@@ -0,0 +1,10 @@
+3e268652f4df305e6a43b114007e2d39 *DESCRIPTION
+7928ddd104b66affd60aa4b41b97a59c *NAMESPACE
+d5a96eeca107e0fd9a921fcdbacbdedf *R/nnls.R
+0e1136ce55b874bcc3c56a0b2435f043 *R/zzz.R
+9bdb2d8803a5d7eea4a4744ca983e39a *inst/COPYRIGHTS
+7e59ea0c95fe3a0d1eb424f5de5df8fd *man/nnls-package.Rd
+8324b62a5ee954a869a553c1635e2cb7 *man/nnls.Rd
+a674f8fab264cc0d0dc526bd913be57e *man/nnnpls.Rd
+b2be94c3abd8d61c2c52b355739fb41c *src/lawson_hanson_nnls.f
+32649ce3a5d68d5e49f70193ebffe2fe *src/nnnpls.f
diff --git a/NAMESPACE b/NAMESPACE
new file mode 100644
index 0000000..9c18437
--- /dev/null
+++ b/NAMESPACE
@@ -0,0 +1,13 @@
+export(nnls,nnnpls)
+
+S3method("print", "nnls")
+S3method("residuals", "nnls")
+S3method("fitted", "nnls")
+S3method("coef", "nnls")
+S3method("deviance", "nnls")
+
+S3method("print", "nnnpls")
+S3method("residuals", "nnnpls")
+S3method("fitted", "nnnpls")
+S3method("coef", "nnnpls")
+S3method("deviance", "nnnpls")
diff --git a/R/nnls.R b/R/nnls.R
new file mode 100644
index 0000000..d7b37eb
--- /dev/null
+++ b/R/nnls.R
@@ -0,0 +1,94 @@
+nnls <- function(A, b) {
+ MDA <- M <- nrow(A)
+ N <- ncol(A)
+ RNORM <- MODE <- NSETP <- 0
+ W <- INDEX <- X <- rep(0, N)
+ ZZ <- rep(0, M)
+ sol <- .Fortran("nnls", A = as.numeric(A), MDA = as.integer(MDA), M =
+ as.integer(M), N = as.integer(N), B = as.numeric(b),
+ X = as.numeric(X), RNORM = as.numeric(RNORM), W =
+ as.numeric(W), ZZ = as.numeric(ZZ), INDEX =
+ as.integer(INDEX), MODE = as.integer(MODE),
+ NSETP = as.integer(NSETP), PACKAGE="nnls")
+ fitted <- A %*% sol$X
+ resid <- b - fitted
+ index <- sol$INDEX
+ nsetp <- sol$NSETP
+ if(nsetp > 0)
+ passive <- index[1:nsetp]
+ else passive <- vector()
+ if(nsetp == N)
+ bound <- vector()
+ else
+ bound <- index[(nsetp+1):N]
+ nnls.out <- list(x=sol$X, deviance=sol$RNORM^2,
+ residuals=resid, fitted = fitted,mode=sol$MODE,
+ passive = passive, bound = bound, nsetp = nsetp)
+ class(nnls.out) <- "nnls"
+ nnls.out
+}
+nnnpls <- function(A, b, con) {
+ MDA <- M <- nrow(A)
+ N <- ncol(A)
+ RNORM <- MODE <- NSETP <- 0
+ W <- INDEX <- X <- rep(0, N)
+ ZZ <- rep(0, M)
+ sol <- .Fortran("nnnpls", A = as.numeric(A), MDA = as.integer(MDA), M =
+ as.integer(M), N = as.integer(N),
+ CON = as.numeric(con),
+ B = as.numeric(b),
+ X = as.numeric(X),
+ RNORM = as.numeric(RNORM), W =
+ as.numeric(W), ZZ = as.numeric(ZZ), INDEX =
+ as.integer(INDEX), MODE = as.integer(MODE),
+ NSETP = as.integer(NSETP), PACKAGE="nnls")
+ fitted <- A %*% sol$X
+ resid <- b - fitted
+ index <- sol$INDEX
+ nsetp <- sol$NSETP
+ if(nsetp > 0)
+ passive <- index[1:nsetp]
+ else passive <- vector()
+ if(nsetp == N)
+ bound <- vector()
+ else
+ bound <- index[(nsetp+1):N]
+ nnnpls.out <- list(x=sol$X, deviance=sol$RNORM^2,
+ residuals=resid, fitted = fitted,mode=sol$MODE,
+ passive = passive, bound = bound, nsetp = nsetp)
+ class(nnnpls.out) <- "nnnpls"
+ nnnpls.out
+}
+print.nnnpls <- function(x, digits = max(3, getOption("digits") - 3), ...)
+{
+ cat("Nonnegative-nonpositive least squares model\n")
+
+ cat("x estimates:", x$x, "\n")
+ cat("residual sum-of-squares: ", format(x$deviance, digits = digits),
+ "\n", sep = '')
+ stopmess <- switch(x$mode, "The solution has been computed sucessfully.",
+ "The dimensions of the problem are bad",
+ "Iteration count exceded. More than 3*N iterations.")
+
+ cat("reason terminated: ", stopmess, "\n", sep='')
+ invisible(x)
+}
+print.nnls <- function(x, digits = max(3, getOption("digits") - 3), ...)
+{
+ cat("Nonnegative least squares model\n")
+
+ cat("x estimates:", x$x, "\n")
+ cat("residual sum-of-squares: ", format(x$deviance, digits = digits),
+ "\n", sep = '')
+ stopmess <- switch(x$mode, "The solution has been computed sucessfully.",
+ "The dimensions of the problem are bad",
+ "Iteration count exceded. More than 3*N iterations.")
+
+ cat("reason terminated: ", stopmess, "\n", sep='')
+ invisible(x)
+}
+
+residuals.nnls <- residuals.nnnpls <- function(object, ...) object$residuals
+coef.nnls <- coef.nnnpls <- function(object, ...) object$x
+fitted.nnls <- fitted.nnnpls <- function(object, ...) object$fitted
+deviance.nnls <- deviance.nnnpls <- function(object, ...) object$deviance
diff --git a/R/zzz.R b/R/zzz.R
new file mode 100644
index 0000000..e6480c5
--- /dev/null
+++ b/R/zzz.R
@@ -0,0 +1,6 @@
+".onLoad" <- function (lib, pkg)
+{
+
+ library.dynam(pkg, pkg, lib)
+
+}
diff --git a/inst/COPYRIGHTS b/inst/COPYRIGHTS
new file mode 100644
index 0000000..5f87fdd
--- /dev/null
+++ b/inst/COPYRIGHTS
@@ -0,0 +1,16 @@
+src/nnls.f
+ From the public domain source code distributed with the book
+ Lawson CL, Hanson RJ (1995). Solving Least Squares
+ Problems. Classics in Applied Mathematics. SIAM, Philadelphia, and
+ downloaded from http://www.netlib.org/lawson-hanson/
+src/nnnpls.f
+ From the public domain source code distributed with the book
+ Lawson CL, Hanson RJ (1995). Solving Least Squares
+ Problems. Classics in Applied Mathematics. SIAM, Philadelphia, and
+ downloaded from http://www.netlib.org/lawson-hanson/
+ with trivial modifications to allow for non-positive constraints in
+ combination with non-negative constraints.
+
+
+
+
diff --git a/man/nnls-package.Rd b/man/nnls-package.Rd
new file mode 100644
index 0000000..b43655d
--- /dev/null
+++ b/man/nnls-package.Rd
@@ -0,0 +1,31 @@
+\name{nnls-package}
+\alias{nnls-package}
+\docType{package}
+\title{The Lawson-Hanson NNLS implementation of non-negative least squares}
+\description{
+ An R interface to the Lawson-Hanson
+ NNLS implementation of an algorithm
+ for non-negative linear least squares
+ that solves the least squares problem
+ \eqn{\min{\parallel A x = b \parallel_2}}
+ with the constraint \eqn{x \ge 0} where
+ \eqn{x \in R^n, b \in R^m} and \eqn{A} is an
+ \eqn{m \times n} matrix.
+ Also allows the combination of non-negative and non-positive
+ constraints on \eqn{x}.
+}
+
+\references{
+Lawson CL, Hanson RJ (1974). Solving Least Squares Problems. Prentice
+Hall, Englewood Cliffs, NJ.
+
+Lawson CL, Hanson RJ (1995). Solving Least Squares Problems. Classics
+in Applied Mathematics. SIAM, Philadelphia.
+}
+
+\keyword{ package }
+\seealso{ \link{nnls}, \link{nnnpls},
+ the method \code{"L-BFGS-B"} for \link{optim},
+ \link[quadprog]{solve.QP}, \link[bvls]{bvls}
+}
+
diff --git a/man/nnls.Rd b/man/nnls.Rd
new file mode 100644
index 0000000..6330b0d
--- /dev/null
+++ b/man/nnls.Rd
@@ -0,0 +1,155 @@
+\name{nnls}
+\alias{nnls}
+\title{The Lawson-Hanson NNLS implemention of non-negative least squares}
+\description{
+ An R interface to the Lawson-Hanson
+ NNLS implementation of an algorithm
+ for non-negative linear least squares
+ that solves
+ \eqn{\min{\parallel A x - b \parallel_2}} with the
+ constraint \eqn{x \ge 0}, where
+ \eqn{x \in R^n, b \in R^m} and \eqn{A} is an \eqn{m \times n} matrix.
+}
+\usage{
+nnls(A, b)
+}
+\arguments{
+ \item{A}{numeric matrix with \code{m} rows and \code{n} columns}
+ \item{b}{numeric vector of length \code{m} }
+}
+\value{
+ \code{nnls} returns
+ an object of class \code{"nnls"}.
+
+ The generic accessor functions \code{coefficients},
+ \code{fitted.values}, \code{deviance} and \code{residuals} extract
+ various useful features of the value returned by \code{nnls}.
+
+ An object of class \code{"nnls"} is a list containing the
+ following components:
+
+
+ \item{x}{the parameter estimates.}
+ \item{deviance}{the residual sum-of-squares.}
+ \item{residuals}{the residuals, that is response minus fitted values.}
+ \item{fitted}{the fitted values.}
+ \item{mode}{a character vector containing a message regarding why
+ termination occured.}
+ \item{passive}{vector of the indices of \code{x} that are not bound
+ at zero. }
+ \item{bound}{vector of the indices of \code{x} that are bound
+ at zero.}
+ \item{nsetp}{the number of elements of \code{x} that are not bound
+ at zero. }
+}
+\references{
+Lawson CL, Hanson RJ (1974). Solving Least Squares Problems. Prentice
+Hall, Englewood Cliffs, NJ.
+
+Lawson CL, Hanson RJ (1995). Solving Least Squares Problems. Classics
+in Applied Mathematics. SIAM, Philadelphia.
+}
+
+\source{
+ This is an R interface to the Fortran77 code distributed
+ with the book referenced below by Lawson CL, Hanson RJ (1995),
+ obtained from Netlib (file \file{lawson-hanson/all}),
+ with a trivial modification to return the variable
+ NSETP.
+}
+\seealso{\link{nnnpls}, the method \code{"L-BFGS-B"} for \link{optim},
+ \link[quadprog]{solve.QP}, \link[bvls]{bvls}
+}
+\examples{
+## simulate a matrix A
+## with 3 columns, each containing an exponential decay
+t <- seq(0, 2, by = .04)
+k <- c(.5, .6, 1)
+A <- matrix(nrow = 51, ncol = 3)
+Acolfunc <- function(k, t) exp(-k*t)
+for(i in 1:3) A[,i] <- Acolfunc(k[i],t)
+
+## simulate a matrix X
+## with 3 columns, each containing a Gaussian shape
+## the Gaussian shapes are non-negative
+X <- matrix(nrow = 51, ncol = 3)
+wavenum <- seq(18000,28000, by=200)
+location <- c(25000, 22000, 20000)
+delta <- c(3000,3000,3000)
+Xcolfunc <- function(wavenum, location, delta)
+ exp( - log(2) * (2 * (wavenum - location)/delta)^2)
+for(i in 1:3) X[,i] <- Xcolfunc(wavenum, location[i], delta[i])
+
+## set seed for reproducibility
+set.seed(3300)
+
+## simulated data is the product of A and X with some
+## spherical Gaussian noise added
+matdat <- A \%*\% t(X) + .005 * rnorm(nrow(A) * nrow(X))
+
+## estimate the rows of X using NNLS criteria
+nnls_sol <- function(matdat, A) {
+ X <- matrix(0, nrow = 51, ncol = 3)
+ for(i in 1:ncol(matdat))
+ X[i,] <- coef(nnls(A,matdat[,i]))
+ X
+}
+X_nnls <- nnls_sol(matdat,A)
+
+matplot(X_nnls,type="b",pch=20)
+abline(0,0, col=grey(.6))
+
+\dontrun{
+## can solve the same problem with L-BFGS-B algorithm
+## but need starting values for x
+bfgs_sol <- function(matdat, A) {
+ startval <- rep(0, ncol(A))
+ fn1 <- function(par1, b, A) sum( ( b - A \%*\% par1)^2)
+ X <- matrix(0, nrow = 51, ncol = 3)
+ for(i in 1:ncol(matdat))
+ X[i,] <- optim(startval, fn = fn1, b=matdat[,i], A=A,
+ lower = rep(0,3), method="L-BFGS-B")$par
+ X
+}
+X_bfgs <- bfgs_sol(matdat,A)
+
+## the RMS deviation under NNLS is less than under L-BFGS-B
+sqrt(sum((X - X_nnls)^2)) < sqrt(sum((X - X_bfgs)^2))
+
+## and L-BFGS-B is much slower
+system.time(nnls_sol(matdat,A))
+system.time(bfgs_sol(matdat,A))
+
+## can also solve the same problem by reformulating it as a
+## quadratic program (this requires the quadprog package; if you
+## have quadprog installed, uncomment lines below starting with
+## only 1 "#" )
+
+# library(quadprog)
+
+# quadprog_sol <- function(matdat, A) {
+# X <- matrix(0, nrow = 51, ncol = 3)
+# bvec <- rep(0, ncol(A))
+# Dmat <- crossprod(A,A)
+# Amat <- diag(ncol(A))
+# for(i in 1:ncol(matdat)) {
+# dvec <- crossprod(A,matdat[,i])
+# X[i,] <- solve.QP(dvec = dvec, bvec = bvec, Dmat=Dmat,
+# Amat=Amat)$solution
+# }
+# X
+# }
+# X_quadprog <- quadprog_sol(matdat,A)
+
+## the RMS deviation under NNLS is about the same as under quadprog
+# sqrt(sum((X - X_nnls)^2))
+# sqrt(sum((X - X_quadprog)^2))
+
+## and quadprog requires about the same amount of time
+# system.time(nnls_sol(matdat,A))
+# system.time(quadprog_sol(matdat,A))
+
+}
+
+}
+\keyword{optimize}
diff --git a/man/nnnpls.Rd b/man/nnnpls.Rd
new file mode 100644
index 0000000..5a19789
--- /dev/null
+++ b/man/nnnpls.Rd
@@ -0,0 +1,161 @@
+\name{nnnpls}
+\alias{nnnpls}
+\title{An implementation of least squares with non-negative and non-positive
+ constraints}
+\description{
+ An implementation of an algorithm for linear least squares problems
+ with non-negative and non-positive
+ constraints based on the Lawson-Hanson
+ NNLS algorithm. Solves \eqn{\min{\parallel A x - b \parallel_2}}
+ with the constraint \eqn{x_i \ge 0}
+ if \eqn{con_i \ge 0} and \eqn{x_i \le 0} otherwise, where
+ \eqn{x, con \in R^n, b \in R^m}, and \eqn{A} is an \eqn{m \times n} matrix.
+}
+\usage{
+nnnpls(A, b, con)
+}
+\arguments{
+ \item{A}{numeric matrix with \code{m} rows and \code{n} columns}
+ \item{b}{numeric vector of length \code{m} }
+ \item{con}{numeric vector of length \code{m} where element \code{i}
+ is negative if and only if element \code{i} of the solution vector
+ \code{x} should be constrained to non-positive, as opposed to
+ non-negative, values. }
+}
+\value{
+ \code{nnnpls} returns
+ an object of class \code{"nnnpls"}.
+
+ The generic accessor functions \code{coefficients},
+ \code{fitted.values}, \code{deviance} and \code{residuals} extract
+ various useful features of the value returned by \code{nnnpls}.
+
+ An object of class \code{"nnnpls"} is a list containing the
+ following components:
+
+ \item{x}{the parameter estimates.}
+ \item{deviance}{the residual sum-of-squares.}
+ \item{residuals}{the residuals, that is response minus fitted values.}
+ \item{fitted}{the fitted values.}
+ \item{mode}{a character vector containing a message regarding why
+ termination occured.}
+ \item{passive}{vector of the indices of \code{x} that are not bound
+ at zero. }
+ \item{bound}{vector of the indices of \code{x} that are bound
+ at zero.}
+ \item{nsetp}{the number of elements of \code{x} that are not bound
+ at zero. }
+}
+\references{
+Lawson CL, Hanson RJ (1974). Solving Least Squares Problems. Prentice
+Hall, Englewood Cliffs, NJ.
+
+Lawson CL, Hanson RJ (1995). Solving Least Squares Problems. Classics
+in Applied Mathematics. SIAM, Philadelphia.
+}
+
+
+\source{
+ This is an R interface to Fortran77 code distributed
+ with the book referenced below by Lawson CL, Hanson RJ (1995),
+ obtained from Netlib (file \file{lawson-hanson/all}), with some
+ trivial modifications to allow for the combination of constraints to
+ non-negative and non-positive values, and to return the variable
+ NSETP.
+}
+\seealso{
+\link{nnls}, the method \code{"L-BFGS-B"} for \link{optim},
+\link[quadprog]{solve.QP}, \link[bvls]{bvls}
+
+}
+\examples{
+## simulate a matrix A
+## with 3 columns, each containing an exponential decay
+t <- seq(0, 2, by = .04)
+k <- c(.5, .6, 1)
+A <- matrix(nrow = 51, ncol = 3)
+Acolfunc <- function(k, t) exp(-k*t)
+for(i in 1:3) A[,i] <- Acolfunc(k[i],t)
+
+## simulate a matrix X
+## with 3 columns, each containing a Gaussian shape
+## 2 of the Gaussian shapes are non-negative and 1 is non-positive
+X <- matrix(nrow = 51, ncol = 3)
+wavenum <- seq(18000,28000, by=200)
+location <- c(25000, 22000, 20000)
+delta <- c(3000,3000,3000)
+Xcolfunc <- function(wavenum, location, delta)
+ exp( - log(2) * (2 * (wavenum - location)/delta)^2)
+for(i in 1:3) X[,i] <- Xcolfunc(wavenum, location[i], delta[i])
+X[,2] <- -X[,2]
+
+## set seed for reproducibility
+set.seed(3300)
+
+## simulated data is the product of A and X with some
+## spherical Gaussian noise added
+matdat <- A \%*\% t(X) + .005 * rnorm(nrow(A) * nrow(X))
+
+## estimate the rows of X using NNNPLS criteria
+nnnpls_sol <- function(matdat, A) {
+ X <- matrix(0, nrow = 51, ncol = 3)
+ for(i in 1:ncol(matdat))
+ X[i,] <- coef(nnnpls(A,matdat[,i],con=c(1,-1,1)))
+ X
+}
+X_nnnpls <- nnnpls_sol(matdat,A)
+
+\dontrun{
+## can solve the same problem with L-BFGS-B algorithm
+## but need starting values for x and
+## impose a very low/high bound where none is desired
+bfgs_sol <- function(matdat, A) {
+ startval <- rep(0, ncol(A))
+ fn1 <- function(par1, b, A) sum( ( b - A \%*\% par1)^2)
+ X <- matrix(0, nrow = 51, ncol = 3)
+ for(i in 1:ncol(matdat))
+ X[i,] <- optim(startval, fn = fn1, b=matdat[,i], A=A,
+ lower=rep(0, -1000, 0), upper=c(1000,0,1000),
+ method="L-BFGS-B")$par
+ X
+}
+X_bfgs <- bfgs_sol(matdat,A)
+
+## the RMS deviation under NNNPLS is less than under L-BFGS-B
+sqrt(sum((X - X_nnnpls)^2)) < sqrt(sum((X - X_bfgs)^2))
+
+## and L-BFGS-B is much slower
+system.time(nnnpls_sol(matdat,A))
+system.time(bfgs_sol(matdat,A))
+
+## can also solve the same problem by reformulating it as a
+## quadratic program (this requires the quadprog package; if you
+## have quadprog installed, uncomment lines below starting with
+## only 1 "#" )
+
+# library(quadprog)
+
+# quadprog_sol <- function(matdat, A) {
+# X <- matrix(0, nrow = 51, ncol = 3)
+# bvec <- rep(0, ncol(A))
+# Dmat <- crossprod(A,A)
+# Amat <- diag(c(1,-1,1))
+# for(i in 1:ncol(matdat)) {
+# dvec <- crossprod(A,matdat[,i])
+# X[i,] <- solve.QP(dvec = dvec, bvec = bvec, Dmat=Dmat,
+# Amat=Amat)$solution
+# }
+# X
+# }
+# X_quadprog <- quadprog_sol(matdat,A)
+
+## the RMS deviation under NNNPLS is about the same as under quadprog
+# sqrt(sum((X - X_nnnpls)^2))
+# sqrt(sum((X - X_quadprog)^2))
+
+## and quadprog requires about the same amount of time
+# system.time(nnnpls_sol(matdat,A))
+# system.time(quadprog_sol(matdat,A))
+}
+}
+\keyword{optimize}
diff --git a/src/lawson_hanson_nnls.f b/src/lawson_hanson_nnls.f
new file mode 100644
index 0000000..b2cd997
--- /dev/null
+++ b/src/lawson_hanson_nnls.f
@@ -0,0 +1,474 @@
+c Note that calls to 'write' have been commented out
+c since such call trigger warnings in R -KMM, March 2012
+c Also, made DUMMY double precision DUMMY(*)
+
+ double precision FUNCTION DIFF(X,Y)
+c
+c Function used in tests that depend on machine precision.
+c
+c The original version of this code was developed by
+c Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
+c 1973 JUN 7, and published in the book
+c "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
+c Revised FEB 1995 to accompany reprinting of the book by SIAM.
+C
+ double precision X, Y
+ DIFF=X-Y
+ RETURN
+ END
+ SUBROUTINE G1 (A,B,CTERM,STERM,SIG)
+c
+C COMPUTE ORTHOGONAL ROTATION MATRIX..
+c
+c The original version of this code was developed by
+c Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
+c 1973 JUN 12, and published in the book
+c "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
+c Revised FEB 1995 to accompany reprinting of the book by SIAM.
+C
+C COMPUTE.. MATRIX (C, S) SO THAT (C, S)(A) = (SQRT(A**2+B**2))
+C (-S,C) (-S,C)(B) ( 0 )
+C COMPUTE SIG = SQRT(A**2+B**2)
+C SIG IS COMPUTED LAST TO ALLOW FOR THE POSSIBILITY THAT
+C SIG MAY BE IN THE SAME LOCATION AS A OR B .
+C ------------------------------------------------------------------
+ double precision A, B, CTERM, ONE, SIG, STERM, XR, YR, ZERO
+ parameter(ONE = 1.0d0, ZERO = 0.0d0)
+C ------------------------------------------------------------------
+ if (abs(A) .gt. abs(B)) then
+ XR=B/A
+ YR=sqrt(ONE+XR**2)
+ CTERM=sign(ONE/YR,A)
+ STERM=CTERM*XR
+ SIG=abs(A)*YR
+ RETURN
+ endif
+
+ if (B .ne. ZERO) then
+ XR=A/B
+ YR=sqrt(ONE+XR**2)
+ STERM=sign(ONE/YR,B)
+ CTERM=STERM*XR
+ SIG=abs(B)*YR
+ RETURN
+ endif
+
+ SIG=ZERO
+ CTERM=ZERO
+ STERM=ONE
+ RETURN
+ END
+C SUBROUTINE H12 (MODE,LPIVOT,L1,M,U,IUE,UP,C,ICE,ICV,NCV)
+C
+C CONSTRUCTION AND/OR APPLICATION OF A SINGLE
+C HOUSEHOLDER TRANSFORMATION.. Q = I + U*(U**T)/B
+C
+c The original version of this code was developed by
+c Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
+c 1973 JUN 12, and published in the book
+c "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
+c Revised FEB 1995 to accompany reprinting of the book by SIAM.
+C ------------------------------------------------------------------
+c Subroutine Arguments
+c
+C MODE = 1 OR 2 Selects Algorithm H1 to construct and apply a
+c Householder transformation, or Algorithm H2 to apply a
+c previously constructed transformation.
+C LPIVOT IS THE INDEX OF THE PIVOT ELEMENT.
+C L1,M IF L1 .LE. M THE TRANSFORMATION WILL BE CONSTRUCTED TO
+C ZERO ELEMENTS INDEXED FROM L1 THROUGH M. IF L1 GT. M
+C THE SUBROUTINE DOES AN IDENTITY TRANSFORMATION.
+C U(),IUE,UP On entry with MODE = 1, U() contains the pivot
+c vector. IUE is the storage increment between elements.
+c On exit when MODE = 1, U() and UP contain quantities
+c defining the vector U of the Householder transformation.
+c on entry with MODE = 2, U() and UP should contain
+c quantities previously computed with MODE = 1. These will
+c not be modified during the entry with MODE = 2.
+C C() ON ENTRY with MODE = 1 or 2, C() CONTAINS A MATRIX WHICH
+c WILL BE REGARDED AS A SET OF VECTORS TO WHICH THE
+c HOUSEHOLDER TRANSFORMATION IS TO BE APPLIED.
+c ON EXIT C() CONTAINS THE SET OF TRANSFORMED VECTORS.
+C ICE STORAGE INCREMENT BETWEEN ELEMENTS OF VECTORS IN C().
+C ICV STORAGE INCREMENT BETWEEN VECTORS IN C().
+C NCV NUMBER OF VECTORS IN C() TO BE TRANSFORMED. IF NCV .LE. 0
+C NO OPERATIONS WILL BE DONE ON C().
+C ------------------------------------------------------------------
+ SUBROUTINE H12 (MODE,LPIVOT,L1,M,U,IUE,UP,C,ICE,ICV,NCV)
+C ------------------------------------------------------------------
+ integer I, I2, I3, I4, ICE, ICV, INCR, IUE, J
+ integer L1, LPIVOT, M, MODE, NCV
+ double precision B, C(*), CL, CLINV, ONE, SM
+c double precision U(IUE,M)
+ double precision U(IUE,*)
+ double precision UP
+ parameter(ONE = 1.0d0)
+C ------------------------------------------------------------------
+ IF (0.GE.LPIVOT.OR.LPIVOT.GE.L1.OR.L1.GT.M) RETURN
+ CL=abs(U(1,LPIVOT))
+ IF (MODE.EQ.2) GO TO 60
+C ****** CONSTRUCT THE TRANSFORMATION. ******
+ DO 10 J=L1,M
+ 10 CL=MAX(abs(U(1,J)),CL)
+ IF (CL) 130,130,20
+ 20 CLINV=ONE/CL
+ SM=(U(1,LPIVOT)*CLINV)**2
+ DO 30 J=L1,M
+ 30 SM=SM+(U(1,J)*CLINV)**2
+ CL=CL*SQRT(SM)
+ IF (U(1,LPIVOT)) 50,50,40
+ 40 CL=-CL
+ 50 UP=U(1,LPIVOT)-CL
+ U(1,LPIVOT)=CL
+ GO TO 70
+C ****** APPLY THE TRANSFORMATION I+U*(U**T)/B TO C. ******
+C
+ 60 IF (CL) 130,130,70
+ 70 IF (NCV.LE.0) RETURN
+ B= UP*U(1,LPIVOT)
+C B MUST BE NONPOSITIVE HERE. IF B = 0., RETURN.
+C
+ IF (B) 80,130,130
+ 80 B=ONE/B
+ I2=1-ICV+ICE*(LPIVOT-1)
+ INCR=ICE*(L1-LPIVOT)
+ DO 120 J=1,NCV
+ I2=I2+ICV
+ I3=I2+INCR
+ I4=I3
+ SM=C(I2)*UP
+ DO 90 I=L1,M
+ SM=SM+C(I3)*U(1,I)
+ 90 I3=I3+ICE
+ IF (SM) 100,120,100
+ 100 SM=SM*B
+ C(I2)=C(I2)+SM*UP
+ DO 110 I=L1,M
+ C(I4)=C(I4)+SM*U(1,I)
+ 110 I4=I4+ICE
+ 120 CONTINUE
+ 130 RETURN
+ END
+C SUBROUTINE NNLS (A,MDA,M,N,B,X,RNORM,W,ZZ,INDEX,MODE)
+C
+C Algorithm NNLS: NONNEGATIVE LEAST SQUARES
+C
+c The original version of this code was developed by
+c Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
+c 1973 JUN 15, and published in the book
+c "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
+c Revised FEB 1995 to accompany reprinting of the book by SIAM.
+c
+C GIVEN AN M BY N MATRIX, A, AND AN M-VECTOR, B, COMPUTE AN
+C N-VECTOR, X, THAT SOLVES THE LEAST SQUARES PROBLEM
+C
+C A * X = B SUBJECT TO X .GE. 0
+C ------------------------------------------------------------------
+c Subroutine Arguments
+c
+C A(),MDA,M,N MDA IS THE FIRST DIMENSIONING PARAMETER FOR THE
+C ARRAY, A(). ON ENTRY A() CONTAINS THE M BY N
+C MATRIX, A. ON EXIT A() CONTAINS
+C THE PRODUCT MATRIX, Q*A , WHERE Q IS AN
+C M BY M ORTHOGONAL MATRIX GENERATED IMPLICITLY BY
+C THIS SUBROUTINE.
+C B() ON ENTRY B() CONTAINS THE M-VECTOR, B. ON EXIT B() CON-
+C TAINS Q*B.
+C X() ON ENTRY X() NEED NOT BE INITIALIZED. ON EXIT X() WILL
+C CONTAIN THE SOLUTION VECTOR.
+C RNORM ON EXIT RNORM CONTAINS THE EUCLIDEAN NORM OF THE
+C RESIDUAL VECTOR.
+C W() AN N-ARRAY OF WORKING SPACE. ON EXIT W() WILL CONTAIN
+C THE DUAL SOLUTION VECTOR. W WILL SATISFY W(I) = 0.
+C FOR ALL I IN SET P AND W(I) .LE. 0. FOR ALL I IN SET Z
+C ZZ() AN M-ARRAY OF WORKING SPACE.
+C INDEX() AN INTEGER WORKING ARRAY OF LENGTH AT LEAST N.
+C ON EXIT THE CONTENTS OF THIS ARRAY DEFINE THE SETS
+C P AND Z AS FOLLOWS..
+C
+C INDEX(1) THRU INDEX(NSETP) = SET P.
+C INDEX(IZ1) THRU INDEX(IZ2) = SET Z.
+C IZ1 = NSETP + 1 = NPP1
+C IZ2 = N
+C MODE THIS IS A SUCCESS-FAILURE FLAG WITH THE FOLLOWING
+C MEANINGS.
+C 1 THE SOLUTION HAS BEEN COMPUTED SUCCESSFULLY.
+C 2 THE DIMENSIONS OF THE PROBLEM ARE BAD.
+C EITHER M .LE. 0 OR N .LE. 0.
+C 3 ITERATION COUNT EXCEEDED. MORE THAN 3*N ITERATIONS.
+C
+C ------------------------------------------------------------------
+ SUBROUTINE NNLS (A,MDA,M,N,B,X,RNORM,W,ZZ,INDEX,MODE,NSETP)
+C ------------------------------------------------------------------
+ integer I, II, IP, ITER, ITMAX, IZ, IZ1, IZ2, IZMAX, J, JJ, JZ, L
+ integer M, MDA, MODE,N, NPP1, NSETP, RTNKEY
+c integer INDEX(N)
+c double precision A(MDA,N), B(M), W(N), X(N), ZZ(M)
+ integer INDEX(*)
+ double precision A(MDA,*), B(*), W(*), X(*), ZZ(*), DUMMY(1)
+ double precision ALPHA, ASAVE, CC, DIFF, FACTOR, RNORM
+ double precision SM, SS, T, TEMP, TWO, UNORM, UP, WMAX
+ double precision ZERO, ZTEST
+ parameter(FACTOR = 0.01d0)
+ parameter(TWO = 2.0d0, ZERO = 0.0d0)
+C ------------------------------------------------------------------
+ MODE=1
+ IF (M .le. 0 .or. N .le. 0) then
+ MODE=2
+ RETURN
+ endif
+ ITER=0
+ ITMAX=3*N
+C
+C INITIALIZE THE ARRAYS INDEX() AND X().
+C
+ DO 20 I=1,N
+ X(I)=ZERO
+ 20 INDEX(I)=I
+C
+ IZ2=N
+ IZ1=1
+ NSETP=0
+ NPP1=1
+C ****** MAIN LOOP BEGINS HERE ******
+ 30 CONTINUE
+C QUIT IF ALL COEFFICIENTS ARE ALREADY IN THE SOLUTION.
+C OR IF M COLS OF A HAVE BEEN TRIANGULARIZED.
+C
+ IF (IZ1 .GT.IZ2.OR.NSETP.GE.M) GO TO 350
+C
+C COMPUTE COMPONENTS OF THE DUAL (NEGATIVE GRADIENT) VECTOR W().
+C
+ DO 50 IZ=IZ1,IZ2
+ J=INDEX(IZ)
+ SM=ZERO
+ DO 40 L=NPP1,M
+ 40 SM=SM+A(L,J)*B(L)
+ W(J)=SM
+ 50 continue
+C FIND LARGEST POSITIVE W(J).
+ 60 continue
+ WMAX=ZERO
+ DO 70 IZ=IZ1,IZ2
+ J=INDEX(IZ)
+ IF (W(J) .gt. WMAX) then
+ WMAX=W(J)
+ IZMAX=IZ
+ endif
+ 70 CONTINUE
+C
+C IF WMAX .LE. 0. GO TO TERMINATION.
+C THIS INDICATES SATISFACTION OF THE KUHN-TUCKER CONDITIONS.
+C
+ IF (WMAX .le. ZERO) go to 350
+ IZ=IZMAX
+ J=INDEX(IZ)
+C
+C THE SIGN OF W(J) IS OK FOR J TO BE MOVED TO SET P.
+C BEGIN THE TRANSFORMATION AND CHECK NEW DIAGONAL ELEMENT TO AVOID
+C NEAR LINEAR DEPENDENCE.
+C
+ ASAVE=A(NPP1,J)
+ CALL H12 (1,NPP1,NPP1+1,M,A(1,J),1,UP,DUMMY,1,1,0)
+ UNORM=ZERO
+ IF (NSETP .ne. 0) then
+ DO 90 L=1,NSETP
+ 90 UNORM=UNORM+A(L,J)**2
+ endif
+ UNORM=sqrt(UNORM)
+ IF (DIFF(UNORM+ABS(A(NPP1,J))*FACTOR,UNORM) .gt. ZERO) then
+C
+C COL J IS SUFFICIENTLY INDEPENDENT. COPY B INTO ZZ, UPDATE ZZ
+C AND SOLVE FOR ZTEST ( = PROPOSED NEW VALUE FOR X(J) ).
+C
+ DO 120 L=1,M
+ 120 ZZ(L)=B(L)
+ CALL H12 (2,NPP1,NPP1+1,M,A(1,J),1,UP,ZZ,1,1,1)
+ ZTEST=ZZ(NPP1)/A(NPP1,J)
+C
+C SEE IF ZTEST IS POSITIVE
+C
+ IF (ZTEST .gt. ZERO) go to 140
+ endif
+C
+C REJECT J AS A CANDIDATE TO BE MOVED FROM SET Z TO SET P.
+C RESTORE A(NPP1,J), SET W(J)=0., AND LOOP BACK TO TEST DUAL
+C COEFFS AGAIN.
+C
+ A(NPP1,J)=ASAVE
+ W(J)=ZERO
+ GO TO 60
+C
+C THE INDEX J=INDEX(IZ) HAS BEEN SELECTED TO BE MOVED FROM
+C SET Z TO SET P. UPDATE B, UPDATE INDICES, APPLY HOUSEHOLDER
+C TRANSFORMATIONS TO COLS IN NEW SET Z, ZERO SUBDIAGONAL ELTS IN
+C COL J, SET W(J)=0.
+C
+ 140 continue
+ DO 150 L=1,M
+ 150 B(L)=ZZ(L)
+C
+ INDEX(IZ)=INDEX(IZ1)
+ INDEX(IZ1)=J
+ IZ1=IZ1+1
+ NSETP=NPP1
+ NPP1=NPP1+1
+C
+ IF (IZ1 .le. IZ2) then
+ DO 160 JZ=IZ1,IZ2
+ JJ=INDEX(JZ)
+ CALL H12 (2,NSETP,NPP1,M,A(1,J),1,UP,A(1,JJ),1,MDA,1)
+ 160 continue
+ endif
+C
+ IF (NSETP .ne. M) then
+ DO 180 L=NPP1,M
+ 180 A(L,J)=ZERO
+ endif
+C
+ W(J)=ZERO
+C SOLVE THE TRIANGULAR SYSTEM.
+C STORE THE SOLUTION TEMPORARILY IN ZZ().
+ RTNKEY = 1
+ GO TO 400
+ 200 CONTINUE
+C
+C ****** SECONDARY LOOP BEGINS HERE ******
+C
+C ITERATION COUNTER.
+C
+ 210 continue
+ ITER=ITER+1
+ IF (ITER .gt. ITMAX) then
+ MODE=3
+c write (*,'(/a)') ' NNLS quitting on iteration count.'
+ GO TO 350
+ endif
+C
+C SEE IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE.
+C IF NOT COMPUTE ALPHA.
+C
+ ALPHA=TWO
+ DO 240 IP=1,NSETP
+ L=INDEX(IP)
+ IF (ZZ(IP) .le. ZERO) then
+ T=-X(L)/(ZZ(IP)-X(L))
+ IF (ALPHA .gt. T) then
+ ALPHA=T
+ JJ=IP
+ endif
+ endif
+ 240 CONTINUE
+C
+C IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE THEN ALPHA WILL
+C STILL = 2. IF SO EXIT FROM SECONDARY LOOP TO MAIN LOOP.
+C
+ IF (ALPHA.EQ.TWO) GO TO 330
+C
+C OTHERWISE USE ALPHA WHICH WILL BE BETWEEN 0. AND 1. TO
+C INTERPOLATE BETWEEN THE OLD X AND THE NEW ZZ.
+C
+ DO 250 IP=1,NSETP
+ L=INDEX(IP)
+ X(L)=X(L)+ALPHA*(ZZ(IP)-X(L))
+ 250 continue
+C
+C MODIFY A AND B AND THE INDEX ARRAYS TO MOVE COEFFICIENT I
+C FROM SET P TO SET Z.
+C
+ I=INDEX(JJ)
+ 260 continue
+ X(I)=ZERO
+C
+ IF (JJ .ne. NSETP) then
+ JJ=JJ+1
+ DO 280 J=JJ,NSETP
+ II=INDEX(J)
+ INDEX(J-1)=II
+ CALL G1 (A(J-1,II),A(J,II),CC,SS,A(J-1,II))
+ A(J,II)=ZERO
+ DO 270 L=1,N
+ IF (L.NE.II) then
+c
+c Apply procedure G2 (CC,SS,A(J-1,L),A(J,L))
+c
+ TEMP = A(J-1,L)
+ A(J-1,L) = CC*TEMP + SS*A(J,L)
+ A(J,L) =-SS*TEMP + CC*A(J,L)
+ endif
+ 270 CONTINUE
+c
+c Apply procedure G2 (CC,SS,B(J-1),B(J))
+c
+ TEMP = B(J-1)
+ B(J-1) = CC*TEMP + SS*B(J)
+ B(J) =-SS*TEMP + CC*B(J)
+ 280 continue
+ endif
+c
+ NPP1=NSETP
+ NSETP=NSETP-1
+ IZ1=IZ1-1
+ INDEX(IZ1)=I
+C
+C SEE IF THE REMAINING COEFFS IN SET P ARE FEASIBLE. THEY SHOULD
+C BE BECAUSE OF THE WAY ALPHA WAS DETERMINED.
+C IF ANY ARE INFEASIBLE IT IS DUE TO ROUND-OFF ERROR. ANY
+C THAT ARE NONPOSITIVE WILL BE SET TO ZERO
+C AND MOVED FROM SET P TO SET Z.
+C
+ DO 300 JJ=1,NSETP
+ I=INDEX(JJ)
+ IF (X(I) .le. ZERO) go to 260
+ 300 CONTINUE
+C
+C COPY B( ) INTO ZZ( ). THEN SOLVE AGAIN AND LOOP BACK.
+C
+ DO 310 I=1,M
+ 310 ZZ(I)=B(I)
+ RTNKEY = 2
+ GO TO 400
+ 320 CONTINUE
+ GO TO 210
+C ****** END OF SECONDARY LOOP ******
+C
+ 330 continue
+ DO 340 IP=1,NSETP
+ I=INDEX(IP)
+ 340 X(I)=ZZ(IP)
+C ALL NEW COEFFS ARE POSITIVE. LOOP BACK TO BEGINNING.
+ GO TO 30
+C
+C ****** END OF MAIN LOOP ******
+C
+C COME TO HERE FOR TERMINATION.
+C COMPUTE THE NORM OF THE FINAL RESIDUAL VECTOR.
+C
+ 350 continue
+ SM=ZERO
+ IF (NPP1 .le. M) then
+ DO 360 I=NPP1,M
+ 360 SM=SM+B(I)**2
+ else
+ DO 380 J=1,N
+ 380 W(J)=ZERO
+ endif
+ RNORM=sqrt(SM)
+ RETURN
+C
+C THE FOLLOWING BLOCK OF CODE IS USED AS AN INTERNAL SUBROUTINE
+C TO SOLVE THE TRIANGULAR SYSTEM, PUTTING THE SOLUTION IN ZZ().
+C
+ 400 continue
+ DO 430 L=1,NSETP
+ IP=NSETP+1-L
+ IF (L .ne. 1) then
+ DO 410 II=1,IP
+ ZZ(II)=ZZ(II)-A(II,JJ)*ZZ(IP+1)
+ 410 continue
+ endif
+ JJ=INDEX(IP)
+ ZZ(IP)=ZZ(IP)/A(IP,JJ)
+ 430 continue
+ go to (200, 320), RTNKEY
+ END
diff --git a/src/nnnpls.f b/src/nnnpls.f
new file mode 100644
index 0000000..4075e1a
--- /dev/null
+++ b/src/nnnpls.f
@@ -0,0 +1,353 @@
+c Note that calls to 'write' have been commented out
+c since such call trigger warnings in R -KMM, March 2012
+c Also, made DUMMY double precision DUMMY(*)
+
+C SUBROUTINE NNNPLS (A,MDA,M,N,CON,B,X,RNORM,W,ZZ,INDEX,MODE)
+C
+C Algorithm NNNPLS: NONNEGATIVE LEAST SQUARES with trivial modifications
+C to allow for non-positive constraints as well (made by Katharine Mullen,
+C 11.2007)
+C
+c The original version of this code was developed by
+c Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
+c 1973 JUN 15, and published in the book
+c "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
+c Revised FEB 1995 to accompany reprinting of the book by SIAM.
+c
+C GIVEN AN M BY N MATRIX, A, AND AN M-VECTOR, B, COMPUTE AN
+C N-VECTOR, X, THAT SOLVES THE LEAST SQUARES PROBLEM
+C
+C A * X = B SUBJECT TO X .GE. 0
+C ------------------------------------------------------------------
+c Subroutine Arguments
+c
+C A(),MDA,M,N MDA IS THE FIRST DIMENSIONING PARAMETER FOR THE
+C ARRAY, A(). ON ENTRY A() CONTAINS THE M BY N
+C MATRIX, A. ON EXIT A() CONTAINS
+C THE PRODUCT MATRIX, Q*A , WHERE Q IS AN
+C M BY M ORTHOGONAL MATRIX GENERATED IMPLICITLY BY
+C THIS SUBROUTINE.
+C B() ON ENTRY B() CONTAINS THE M-VECTOR, B. ON EXIT B() CON-
+C TAINS Q*B.
+C CON() M-VECTOR CON, that contains -1 where B is constrained to a
+C non-positive value, 1 where B is constrained to a non-negative
+C value
+C X() ON ENTRY X() NEED NOT BE INITIALIZED. ON EXIT X() WILL
+C CONTAIN THE SOLUTION VECTOR.
+C RNORM ON EXIT RNORM CONTAINS THE EUCLIDEAN NORM OF THE
+C RESIDUAL VECTOR.
+C W() AN N-ARRAY OF WORKING SPACE. ON EXIT W() WILL CONTAIN
+C THE DUAL SOLUTION VECTOR. W WILL SATISFY W(I) = 0.
+C FOR ALL I IN SET P AND W(I) .LE. 0. FOR ALL I IN SET Z
+C ZZ() AN M-ARRAY OF WORKING SPACE.
+C INDEX() AN INTEGER WORKING ARRAY OF LENGTH AT LEAST N.
+C ON EXIT THE CONTENTS OF THIS ARRAY DEFINE THE SETS
+C P AND Z AS FOLLOWS..
+C
+C INDEX(1) THRU INDEX(NSETP) = SET P.
+C INDEX(IZ1) THRU INDEX(IZ2) = SET Z.
+C IZ1 = NSETP + 1 = NPP1
+C IZ2 = N
+C MODE THIS IS A SUCCESS-FAILURE FLAG WITH THE FOLLOWING
+C MEANINGS.
+C 1 THE SOLUTION HAS BEEN COMPUTED SUCCESSFULLY.
+C 2 THE DIMENSIONS OF THE PROBLEM ARE BAD.
+C EITHER M .LE. 0 OR N .LE. 0.
+C 3 ITERATION COUNT EXCEEDED. MORE THAN 3*N ITERATIONS.
+C
+C ------------------------------------------------------------------
+ SUBROUTINE nnnpls(A,MDA,M,N,CON,B,X,RNORM,W,ZZ,INDEX,MODE,NSETP)
+C ------------------------------------------------------------------
+ integer I, II, IP, ITER, ITMAX, IZ, IZ1, IZ2, IZMAX, J, JJ, JZ, L
+ integer M, MDA, MODE,N, NPP1, NSETP, RTNKEY
+c integer INDEX(N)
+c double precision A(MDA,N), B(M), W(N), X(N), ZZ(M)
+ integer INDEX(*)
+ integer C1, C2
+ double precision CON(N)
+ double precision A(MDA,*), B(*), W(*), X(*), ZZ(*), DUMMY(1)
+ double precision ALPHA, ASAVE, CC, DIFF, FACTOR, RNORM
+ double precision SM, SS, T, TEMP, TWO, UNORM, UP, WMAX
+ double precision ZERO, ZTEST
+ parameter(FACTOR = 0.01d0)
+ parameter(TWO = 2.0d0, ZERO = 0.0d0)
+C ------------------------------------------------------------------
+ MODE=1
+ IF (M .le. 0 .or. N .le. 0) then
+ MODE=2
+ RETURN
+ endif
+ ITER=0
+ ITMAX=3*N
+C USE CON() TO FLIP THE SIGN OF A() WHERE A CONSTRAINT TO NON-POSITIVE VALUES
+C IS DESIRED
+ DO C2 = 1,N
+ IF (CON(C2) .lt. ZERO) then
+ DO C1 = 1,M
+ A(C1,C2) = - A(C1,C2)
+ enddo
+ endif
+ enddo
+C
+C INITIALIZE THE ARRAYS INDEX() AND X().
+C
+ DO 20 I=1,N
+ X(I)=ZERO
+ 20 INDEX(I)=I
+C
+ IZ2=N
+ IZ1=1
+ NSETP=0
+ NPP1=1
+C ****** MAIN LOOP BEGINS HERE ******
+ 30 CONTINUE
+C QUIT IF ALL COEFFICIENTS ARE ALREADY IN THE SOLUTION.
+C OR IF M COLS OF A HAVE BEEN TRIANGULARIZED.
+C
+ IF (IZ1 .GT.IZ2.OR.NSETP.GE.M) GO TO 350
+C
+C COMPUTE COMPONENTS OF THE DUAL (NEGATIVE GRADIENT) VECTOR W().
+C
+ DO 50 IZ=IZ1,IZ2
+ J=INDEX(IZ)
+ SM=ZERO
+ DO 40 L=NPP1,M
+ 40 SM=SM+A(L,J)*B(L)
+ W(J)=SM
+ 50 continue
+C FIND LARGEST POSITIVE W(J).
+ 60 continue
+ WMAX=ZERO
+ DO 70 IZ=IZ1,IZ2
+ J=INDEX(IZ)
+ IF (W(J) .gt. WMAX) then
+ WMAX=W(J)
+ IZMAX=IZ
+ endif
+ 70 CONTINUE
+C
+C IF WMAX .LE. 0. GO TO TERMINATION.
+C THIS INDICATES SATISFACTION OF THE KUHN-TUCKER CONDITIONS.
+C
+ IF (WMAX .le. ZERO) go to 350
+ IZ=IZMAX
+ J=INDEX(IZ)
+C
+C THE SIGN OF W(J) IS OK FOR J TO BE MOVED TO SET P.
+C BEGIN THE TRANSFORMATION AND CHECK NEW DIAGONAL ELEMENT TO AVOID
+C NEAR LINEAR DEPENDENCE.
+C
+ ASAVE=A(NPP1,J)
+ CALL H12 (1,NPP1,NPP1+1,M,A(1,J),1,UP,DUMMY,1,1,0)
+ UNORM=ZERO
+ IF (NSETP .ne. 0) then
+ DO 90 L=1,NSETP
+ 90 UNORM=UNORM+A(L,J)**2
+ endif
+ UNORM=sqrt(UNORM)
+ IF (DIFF(UNORM+ABS(A(NPP1,J))*FACTOR,UNORM) .gt. ZERO) then
+C
+C COL J IS SUFFICIENTLY INDEPENDENT. COPY B INTO ZZ, UPDATE ZZ
+C AND SOLVE FOR ZTEST ( = PROPOSED NEW VALUE FOR X(J) ).
+C
+ DO 120 L=1,M
+ 120 ZZ(L)=B(L)
+ CALL H12 (2,NPP1,NPP1+1,M,A(1,J),1,UP,ZZ,1,1,1)
+ ZTEST=ZZ(NPP1)/A(NPP1,J)
+C
+C SEE IF ZTEST IS POSITIVE
+C
+ IF (ZTEST .gt. ZERO) go to 140
+ endif
+C
+C REJECT J AS A CANDIDATE TO BE MOVED FROM SET Z TO SET P.
+C RESTORE A(NPP1,J), SET W(J)=0., AND LOOP BACK TO TEST DUAL
+C COEFFS AGAIN.
+C
+ A(NPP1,J)=ASAVE
+ W(J)=ZERO
+ GO TO 60
+C
+C THE INDEX J=INDEX(IZ) HAS BEEN SELECTED TO BE MOVED FROM
+C SET Z TO SET P. UPDATE B, UPDATE INDICES, APPLY HOUSEHOLDER
+C TRANSFORMATIONS TO COLS IN NEW SET Z, ZERO SUBDIAGONAL ELTS IN
+C COL J, SET W(J)=0.
+C
+ 140 continue
+ DO 150 L=1,M
+ 150 B(L)=ZZ(L)
+C
+ INDEX(IZ)=INDEX(IZ1)
+ INDEX(IZ1)=J
+ IZ1=IZ1+1
+ NSETP=NPP1
+ NPP1=NPP1+1
+C
+ IF (IZ1 .le. IZ2) then
+ DO 160 JZ=IZ1,IZ2
+ JJ=INDEX(JZ)
+ CALL H12 (2,NSETP,NPP1,M,A(1,J),1,UP,A(1,JJ),1,MDA,1)
+ 160 continue
+ endif
+C
+ IF (NSETP .ne. M) then
+ DO 180 L=NPP1,M
+ 180 A(L,J)=ZERO
+ endif
+C
+ W(J)=ZERO
+C SOLVE THE TRIANGULAR SYSTEM.
+C STORE THE SOLUTION TEMPORARILY IN ZZ().
+ RTNKEY = 1
+ GO TO 400
+ 200 CONTINUE
+C
+C ****** SECONDARY LOOP BEGINS HERE ******
+C
+C ITERATION COUNTER.
+C
+ 210 continue
+ ITER=ITER+1
+ IF (ITER .gt. ITMAX) then
+ MODE=3
+c write (*,'(/a)') ' NNLS quitting on iteration count.'
+ GO TO 350
+ endif
+C
+C SEE IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE.
+C IF NOT COMPUTE ALPHA.
+C
+ ALPHA=TWO
+ DO 240 IP=1,NSETP
+ L=INDEX(IP)
+ IF (ZZ(IP) .le. ZERO) then
+ T=-X(L)/(ZZ(IP)-X(L))
+ IF (ALPHA .gt. T) then
+ ALPHA=T
+ JJ=IP
+ endif
+ endif
+ 240 CONTINUE
+C
+C IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE THEN ALPHA WILL
+C STILL = 2. IF SO EXIT FROM SECONDARY LOOP TO MAIN LOOP.
+C
+ IF (ALPHA.EQ.TWO) GO TO 330
+C
+C OTHERWISE USE ALPHA WHICH WILL BE BETWEEN 0. AND 1. TO
+C INTERPOLATE BETWEEN THE OLD X AND THE NEW ZZ.
+C
+ DO 250 IP=1,NSETP
+ L=INDEX(IP)
+ X(L)=X(L)+ALPHA*(ZZ(IP)-X(L))
+ 250 continue
+C
+C MODIFY A AND B AND THE INDEX ARRAYS TO MOVE COEFFICIENT I
+C FROM SET P TO SET Z.
+C
+ I=INDEX(JJ)
+ 260 continue
+ X(I)=ZERO
+C
+ IF (JJ .ne. NSETP) then
+ JJ=JJ+1
+ DO 280 J=JJ,NSETP
+ II=INDEX(J)
+ INDEX(J-1)=II
+ CALL G1 (A(J-1,II),A(J,II),CC,SS,A(J-1,II))
+ A(J,II)=ZERO
+ DO 270 L=1,N
+ IF (L.NE.II) then
+c
+c Apply procedure G2 (CC,SS,A(J-1,L),A(J,L))
+c
+ TEMP = A(J-1,L)
+ A(J-1,L) = CC*TEMP + SS*A(J,L)
+ A(J,L) =-SS*TEMP + CC*A(J,L)
+ endif
+ 270 CONTINUE
+c
+c Apply procedure G2 (CC,SS,B(J-1),B(J))
+c
+ TEMP = B(J-1)
+ B(J-1) = CC*TEMP + SS*B(J)
+ B(J) =-SS*TEMP + CC*B(J)
+ 280 continue
+ endif
+c
+ NPP1=NSETP
+ NSETP=NSETP-1
+ IZ1=IZ1-1
+ INDEX(IZ1)=I
+C
+C SEE IF THE REMAINING COEFFS IN SET P ARE FEASIBLE. THEY SHOULD
+C BE BECAUSE OF THE WAY ALPHA WAS DETERMINED.
+C IF ANY ARE INFEASIBLE IT IS DUE TO ROUND-OFF ERROR. ANY
+C THAT ARE NONPOSITIVE WILL BE SET TO ZERO
+C AND MOVED FROM SET P TO SET Z.
+C
+ DO 300 JJ=1,NSETP
+ I=INDEX(JJ)
+ IF (X(I) .le. ZERO) go to 260
+ 300 CONTINUE
+C
+C COPY B( ) INTO ZZ( ). THEN SOLVE AGAIN AND LOOP BACK.
+C
+ DO 310 I=1,M
+ 310 ZZ(I)=B(I)
+ RTNKEY = 2
+ GO TO 400
+ 320 CONTINUE
+ GO TO 210
+C ****** END OF SECONDARY LOOP ******
+C
+ 330 continue
+ DO 340 IP=1,NSETP
+ I=INDEX(IP)
+ 340 X(I)=ZZ(IP)
+C ALL NEW COEFFS ARE POSITIVE. LOOP BACK TO BEGINNING.
+ GO TO 30
+C
+C ****** END OF MAIN LOOP ******
+C
+C COME TO HERE FOR TERMINATION.
+C COMPUTE THE NORM OF THE FINAL RESIDUAL VECTOR.
+C
+ 350 continue
+C USE CON() TO FLIP THE SIGN OF A() WHERE A CONSTRAINT TO NON-POSITIVE VALUES
+C IS DESIRED
+ DO C2 = 1,N
+ IF (CON(C2) .lt. ZERO) then
+ DO C1 = 1,M
+ A(C1,C2) = - A(C1,C2)
+ enddo
+ x(c2)=-x(c2)
+ endif
+ enddo
+ SM=ZERO
+ IF (NPP1 .le. M) then
+ DO 360 I=NPP1,M
+ 360 SM=SM+B(I)**2
+ else
+ DO 380 J=1,N
+ 380 W(J)=ZERO
+ endif
+ RNORM=sqrt(SM)
+ RETURN
+C
+C THE FOLLOWING BLOCK OF CODE IS USED AS AN INTERNAL SUBROUTINE
+C TO SOLVE THE TRIANGULAR SYSTEM, PUTTING THE SOLUTION IN ZZ().
+C
+ 400 continue
+ DO 430 L=1,NSETP
+ IP=NSETP+1-L
+ IF (L .ne. 1) then
+ DO 410 II=1,IP
+ ZZ(II)=ZZ(II)-A(II,JJ)*ZZ(IP+1)
+ 410 continue
+ endif
+ JJ=INDEX(IP)
+ ZZ(IP)=ZZ(IP)/A(IP,JJ)
+ 430 continue
+ go to (200, 320), RTNKEY
+ END
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/r-cran-nnls.git
More information about the debian-med-commit
mailing list