[med-svn] [r-cran-minpack.lm] 04/09: New upstream version 1.2-1
Andreas Tille
tille at debian.org
Wed Nov 29 18:09:11 UTC 2017
This is an automated email from the git hooks/post-receive script.
tille pushed a commit to branch master
in repository r-cran-minpack.lm.
commit 79f3db1f934d55ca2efa1b7af7c5417da7987431
Author: Andreas Tille <tille at debian.org>
Date: Wed Nov 29 19:00:26 2017 +0100
New upstream version 1.2-1
---
DESCRIPTION | 13 ++
INDEX | 4 +
LICENSE.note | 52 ++++++
MD5 | 37 ++++
NAMESPACE | 15 ++
NEWS | 22 +++
R/nls.lm.R | 110 ++++++++++++
R/nlsLM.R | 275 ++++++++++++++++++++++++++++++
R/wfct.R | 57 +++++++
THANKS | 8 +
debian/changelog | 15 --
debian/compat | 1 -
debian/control | 26 ---
debian/copyright | 27 ---
debian/rules | 4 -
debian/source/format | 1 -
debian/watch | 2 -
man/nls.lm.Rd | 294 ++++++++++++++++++++++++++++++++
man/nls.lm.control.Rd | 85 ++++++++++
man/nlsLM.Rd | 165 ++++++++++++++++++
man/wfct.Rd | 69 ++++++++
src/Makevars | 2 +
src/chkder.f | 140 ++++++++++++++++
src/dogleg.f | 177 ++++++++++++++++++++
src/dpmpar.f | 177 ++++++++++++++++++++
src/enorm.f | 108 ++++++++++++
src/fcn_lmder.c | 67 ++++++++
src/fcn_lmdif.c | 57 +++++++
src/fcn_message.c | 26 +++
src/fdjac2.f | 107 ++++++++++++
src/get_element.c | 15 ++
src/lmder.f | 452 +++++++++++++++++++++++++++++++++++++++++++++++++
src/lmdif.f | 454 ++++++++++++++++++++++++++++++++++++++++++++++++++
src/lmpar.f | 264 +++++++++++++++++++++++++++++
src/minpack_lm.h | 59 +++++++
src/nls_lm.c | 274 ++++++++++++++++++++++++++++++
src/prod.c | 52 ++++++
src/qform.f | 95 +++++++++++
src/qrfac.f | 164 ++++++++++++++++++
src/qrsolv.f | 193 +++++++++++++++++++++
src/r1mpyq.f | 92 ++++++++++
src/r1updt.f | 207 +++++++++++++++++++++++
src/rwupdt.f | 113 +++++++++++++
src/transpose.c | 8 +
src/vector.c | 11 ++
45 files changed, 4520 insertions(+), 76 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
new file mode 100644
index 0000000..2268792
--- /dev/null
+++ b/DESCRIPTION
@@ -0,0 +1,13 @@
+Package: minpack.lm
+Version: 1.2-1
+Title: R Interface to the Levenberg-Marquardt Nonlinear Least-Squares
+ Algorithm Found in MINPACK, Plus Support for Bounds
+Author: Timur V. Elzhov, Katharine M. Mullen, Andrej-Nikolai Spiess, Ben Bolker
+Maintainer: Katharine M. Mullen <mullenkate at gmail.com>
+Description: The nls.lm function provides an R interface to lmder and lmdif from the MINPACK library, for solving nonlinear least-squares problems by a modification of the Levenberg-Marquardt algorithm, with support for lower and upper parameter bounds. The implementation can be used via nls-like calls using the nlsLM function.
+Suggests: MASS
+License: GPL-3
+NeedsCompilation: yes
+Packaged: 2016-11-20 02:45:43 UTC; kmm
+Repository: CRAN
+Date/Publication: 2016-11-20 16:40:02
diff --git a/INDEX b/INDEX
new file mode 100644
index 0000000..5546a5f
--- /dev/null
+++ b/INDEX
@@ -0,0 +1,4 @@
+nls.lm Addresses NLS problems with the
+ Levenberg-Marquardt algorithm
+nls.lm.control Control various aspects of the
+ Levenberg-Marquardt algorithm
diff --git a/LICENSE.note b/LICENSE.note
new file mode 100644
index 0000000..11d8a9a
--- /dev/null
+++ b/LICENSE.note
@@ -0,0 +1,52 @@
+Minpack Copyright Notice (1999) University of Chicago. All rights reserved
+
+Redistribution and use in source and binary forms, with or
+without modification, are permitted provided that the
+following conditions are met:
+
+1. Redistributions of source code must retain the above
+copyright notice, this list of conditions and the following
+disclaimer.
+
+2. Redistributions in binary form must reproduce the above
+copyright notice, this list of conditions and the following
+disclaimer in the documentation and/or other materials
+provided with the distribution.
+
+3. The end-user documentation included with the
+redistribution, if any, must include the following
+acknowledgment:
+
+ "This product includes software developed by the
+ University of Chicago, as Operator of Argonne National
+ Laboratory.
+
+Alternately, this acknowledgment may appear in the software
+itself, if and wherever such third-party acknowledgments
+normally appear.
+
+4. WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS"
+WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE
+UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND
+THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES
+OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE
+OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY
+OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR
+USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF
+THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4)
+DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION
+UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL
+BE CORRECTED.
+
+5. LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT
+HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF
+ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT,
+INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF
+ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF
+PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER
+SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT
+(INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE,
+EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE
+POSSIBILITY OF SUCH LOSS OR DAMAGES.
+
diff --git a/MD5 b/MD5
new file mode 100644
index 0000000..272859d
--- /dev/null
+++ b/MD5
@@ -0,0 +1,37 @@
+335a3bb43c6ca3b7e27a1b4878800112 *DESCRIPTION
+3677943af2b5ea8bf902576c7e2d9749 *INDEX
+5fe4603e80ef7390306f51ef74449bbd *LICENSE.note
+30c62ba66108dc517d49bd7e00a67fca *NAMESPACE
+0e3335d81b9dff79badce8cc301b5593 *NEWS
+265ff6dbf621ad27f6b6a0aa3a764aaa *R/nls.lm.R
+4cb406cb55d72e9936d3fd34e70fe392 *R/nlsLM.R
+3dbf4ac36ed7d60d91945b849218b402 *R/wfct.R
+4a6a3e370ef553a6572e95939efed59d *THANKS
+cd1248eeaf577b58ca88d49220f42256 *man/nls.lm.Rd
+aeb92193f7c054a0f1e209583f329d04 *man/nls.lm.control.Rd
+21d97cdad27d47c66c947a2019fde2a6 *man/nlsLM.Rd
+4a2664f64117659a022ac8c953b764f3 *man/wfct.Rd
+0c6eb1ce22538691b7a80032eb2a6483 *src/Makevars
+1d0363640860c4485d1bae9482975fd7 *src/chkder.f
+45628269d1f116651b9b6a9602785c85 *src/dogleg.f
+290b5ab2f116903e49c8f542d7afbac6 *src/dpmpar.f
+a63a84008c57c577e03b7d892b964bd5 *src/enorm.f
+2044df416d6ff7f47cbaeaf4eff1ab0c *src/fcn_lmder.c
+e45ad6ae5d49544a957771619c23b3ff *src/fcn_lmdif.c
+475d6c4428b6ffb8aca76a39c725583e *src/fcn_message.c
+748ae7494ea2ef9d5e30b0324386204f *src/fdjac2.f
+5a9d8a815eb087d5d78fadbb259562ff *src/get_element.c
+c85d3371f3d1da5087d23a10a8b96759 *src/lmder.f
+62371f45465dafa3a9aebb136523b1b4 *src/lmdif.f
+8ed1e10412c27be1f2490ba0f37a1f7b *src/lmpar.f
+e056dbe99f9f59baf8a7473e638d2b62 *src/minpack_lm.h
+987ccfc990378cd92b844135deddb950 *src/nls_lm.c
+e297546d0149b0fdc2f807c232c4009a *src/prod.c
+85e35e9406aaf057cec7003b9c10747b *src/qform.f
+fc1df4af782730f0fa5110dec34d4e6b *src/qrfac.f
+a74ea0548499e332f79b0226c6ae83b8 *src/qrsolv.f
+b7a12daf8df48a7ada4835cd1e3d7acc *src/r1mpyq.f
+1697d6f366ed5dfe59ca835eeb042917 *src/r1updt.f
+2c691b14b200dc2178718fa5de9d0268 *src/rwupdt.f
+0c1d0ddb4570597b918d1a85c3770263 *src/transpose.c
+a4b444caea9bdfd0cf8e11e81a135419 *src/vector.c
diff --git a/NAMESPACE b/NAMESPACE
new file mode 100644
index 0000000..7c80cf0
--- /dev/null
+++ b/NAMESPACE
@@ -0,0 +1,15 @@
+useDynLib(minpack.lm)
+export(nls.lm, nls.lm.control, nlsLM, wfct)
+S3method("print", "nls.lm")
+S3method("residuals", "nls.lm")
+S3method("coef", "nls.lm")
+S3method("deviance", "nls.lm")
+S3method("summary", "nls.lm")
+S3method("vcov", "nls.lm")
+S3method("df.residual", "nls.lm")
+S3method("print", "summary.nls.lm")
+importFrom("stats", "as.formula", "coef", "deviance", "getInitial",
+ "model.weights", "na.omit", "nls.control", "numericDeriv",
+ "printCoefmat", "pt", "qt", "resid", "residuals",
+ "setNames", "var", "vcov")
+
diff --git a/NEWS b/NEWS
new file mode 100644
index 0000000..c6c5d14
--- /dev/null
+++ b/NEWS
@@ -0,0 +1,22 @@
+Version 1.2-1
+ o In some cases 'profile.nls' would not work on models created with 'nlsLM', due to a wrong class assignment in the output ("call" instead of "numeric"). This bug was kindly made aware by Patrice Kiener and removed.
+ o License is now only GPL-3; per clause 7 we include the University of
+ Chicago license terms as a supplement.
+
+Version 1.2-0
+ o Removed a small bug in 'nlsLM' that gave an error if an external function was supplied within the formula syntax, e.g. y ~ fitModel(x, a). This was made aware by an email from C. Softley, pointing to a discussion and solution in http://stackoverflow.com/questions/30089698/error-in-evalexpr-envir-enclos-could-not-find-function-nested-functions.
+
+Version 1.1-9
+ o Modified NAMESPACE and added Suggests: MASS to DESCRIPTION file.
+Version 1.1-8
+ o The new function wfct for easy specification of weights in nls and
+ nlsLM is now included.
+Version 1.1-7
+ o changed the License to GPL-3 + the terms specified in the file
+ LICENSE, which give the original copyright distributed with
+ the minpack Fortran source. This is possible given Clause 7 in
+ GPL-3.
+Version 1.1-6
+ o Support for lower and upper bounds added
+ o Added 'nlsLM', a modified version of 'nls' that uses ‘nls.lm’ for
+ fitting.
\ No newline at end of file
diff --git a/R/nls.lm.R b/R/nls.lm.R
new file mode 100644
index 0000000..a07ae9b
--- /dev/null
+++ b/R/nls.lm.R
@@ -0,0 +1,110 @@
+nls.lm.control <- function(ftol = sqrt(.Machine$double.eps), ptol =
+ sqrt(.Machine$double.eps), gtol = 0, diag
+ = list(), epsfcn = 0, factor = 100, maxfev
+ = integer(), maxiter = 50, nprint = 0)
+ list(ftol = ftol, ptol = ptol, gtol = gtol, diag = diag, epsfcn =
+ epsfcn, factor = factor, maxfev = maxfev,
+ maxiter = maxiter, nprint = nprint)
+nls.lm <- function(par, lower=NULL, upper=NULL,fn, jac = NULL,
+ control = nls.lm.control(), ...)
+{
+ fn1 <- function(par) fn(par, ...)
+ jac1 <- if (!is.null(jac))
+ function(par) jac(par, ...)
+ if(is.null(lower))
+ lower <- rep(-Inf,length(par))
+ if(is.null(upper))
+ upper <- rep(Inf,length(par))
+ if(length(lower) != length(par))
+ stop("length(lower) must be equal to length(par)")
+ if(length(upper) != length(par))
+ stop("length(upper) must be equal to length(par)")
+ ctrl <- nls.lm.control()
+ if(!missing(control)) {
+ control <- as.list(control)
+ ctrl[names(control)] <- control
+ }
+ if(length(ctrl[["maxfev"]]) == 0)
+ ctrl[["maxfev"]] <- 100*(length(unlist(par)) + 1)
+ out <- .Call("nls_lm", par, lower, upper, fn1, jac1, ctrl, new.env(),
+ PACKAGE = "minpack.lm")
+
+ out$hessian <- matrix(out$hessian, nrow = length(unlist(par)))
+
+ names(out$par) <-
+ rownames(out$hessian) <-
+ colnames(out$hessian) <-
+ names(out$diag) <- names(par)
+ out$deviance <- sum(out$fvec^2)
+ class(out) <- "nls.lm"
+ out
+}
+print.nls.lm <- function(x, digits = max(3, getOption("digits") - 3), ...)
+{
+ cat("Nonlinear regression via the Levenberg-Marquardt algorithm\n")
+
+ cat("parameter estimates:", toString(x$par), "\n")
+ cat("residual sum-of-squares: ", format(x$deviance, digits = digits),
+ "\n", sep = '')
+ cat("reason terminated: ", x$message, "\n", sep='')
+ invisible(x)
+}
+deviance.nls.lm <- function(object, ...) object$deviance
+coef.nls.lm <- function(object, ...) unlist(object$par)
+residuals.nls.lm <- function(object, ...) object$fvec
+df.residual.nls.lm <- function(object, ...)
+ length(resid(object)) - length(coef(object))
+summary.nls.lm <- function (object, ...)
+{
+ param <- coef(object)
+ pnames <- names(param)
+ ibb <- chol(object$hessian)
+ ih <- chol2inv(ibb)
+ p <- length(param)
+ rdf <- length(object$fvec) - p
+ resvar <- deviance(object) / rdf
+ se <- sqrt(diag(ih) * resvar)
+ names(se) <- pnames
+ tval <- param/se
+
+ param <- cbind(param, se, tval, 2 * pt(abs(tval), rdf, lower.tail = FALSE))
+ dimnames(param) <- list(pnames, c("Estimate", "Std. Error", "t value", "Pr(>|t|)"))
+ ans <- list(residuals = object$fvec, sigma = sqrt(object$deviance/rdf),
+ df = c(p, rdf), cov.unscaled = ih,
+ info = object$info, niter = object$niter,
+ stopmess = object$message,
+ coefficients = param)
+ class(ans) <- "summary.nls.lm"
+ ans
+}
+print.summary.nls.lm <-
+ function (x, digits = max(3, getOption("digits") - 3), ...)
+
+{
+ df <- x$df
+ rdf <- df[2]
+ cat("\nParameters:\n")
+ printCoefmat(x$coefficients, digits = digits, ...)
+ cat("\nResidual standard error:",
+ format(signif(x$sigma, digits)), "on", rdf, "degrees of freedom\n")
+ cat("Number of iterations to termination:", x$niter, "\n")
+ cat("Reason for termination:", x$stopmess, "\n")
+ invisible(x)
+}
+
+vcov.nls.lm <- function(object,...) {
+ object$deviance/(length(object$fvec)-length(object$par))*solve(object$hessian)
+}
+
+confint.nls.lm <- function(object, parm, level = 0.95, ...) {
+ cc <- coef(object)
+ if (missing(parm)) parm <- seq_along(cc)
+ levs <- c((1-level)/2,0.5+level/2)
+ dfval <- (length(object$fvec)-length(object$par))
+ tdist <- qt(levs[2],dfval)
+ m1 <- outer(sqrt(diag(vcov(object))),c(-1,1))*tdist
+ m2 <- sweep(m1,cc,MARGIN=1,"+")
+ colnames(m2) <- paste(100*levs,"%",sep="")
+ m2[parm,]
+}
+
diff --git a/R/nlsLM.R b/R/nlsLM.R
new file mode 100644
index 0000000..cf7eaf7
--- /dev/null
+++ b/R/nlsLM.R
@@ -0,0 +1,275 @@
+nlsModel <- function (form, data, start, wts, upper = NULL)
+{
+ thisEnv <- environment()
+ env <- new.env(hash = TRUE, parent = environment(form))
+ for (i in names(data)) assign(i, data[[i]], envir = env)
+ ind <- as.list(start)
+ parLength <- 0
+ for (i in names(ind)) {
+ temp <- start[[i]]
+ storage.mode(temp) <- "double"
+ assign(i, temp, envir = env)
+ ind[[i]] <- parLength + seq_along(start[[i]])
+ parLength <- parLength + length(start[[i]])
+ }
+ getPars.noVarying <- function() unlist(setNames(lapply(names(ind),
+ get, envir = env), names(ind)))
+ getPars <- getPars.noVarying
+ internalPars <- getPars()
+ if (!is.null(upper)) upper <- rep_len(upper, parLength)
+ useParams <- rep(TRUE, parLength)
+ lhs <- eval(form[[2L]], envir = env)
+ rhs <- eval(form[[3L]], envir = env)
+ .swts <- if (!missing(wts) && length(wts)) sqrt(wts) else rep_len(1, length(rhs))
+ assign(".swts", .swts, envir = env)
+ resid <- .swts * (lhs - rhs)
+ dev <- sum(resid^2)
+ if (is.null(attr(rhs, "gradient"))) {
+ getRHS.noVarying <- function() {
+ if (is.null(upper))
+ numericDeriv(form[[3L]], names(ind), env)
+ else numericDeriv(form[[3L]], names(ind), env, ifelse(internalPars <
+ upper, 1, -1))
+ }
+ getRHS <- getRHS.noVarying
+ rhs <- getRHS()
+ }
+ else {
+ getRHS.noVarying <- function() eval(form[[3L]], envir = env)
+ getRHS <- getRHS.noVarying
+ }
+ dimGrad <- dim(attr(rhs, "gradient"))
+ marg <- length(dimGrad)
+ if (marg > 0L) {
+ gradSetArgs <- vector("list", marg + 1L)
+ for (i in 2L:marg) gradSetArgs[[i]] <- rep(TRUE, dimGrad[i -
+ 1])
+ useParams <- rep(TRUE, dimGrad[marg])
+ }
+ else {
+ gradSetArgs <- vector("list", 2L)
+ useParams <- rep(TRUE, length(attr(rhs, "gradient")))
+ }
+ npar <- length(useParams)
+ gradSetArgs[[1L]] <- (~attr(ans, "gradient"))[[2L]]
+ gradCall <- switch(length(gradSetArgs) - 1L, call("[", gradSetArgs[[1L]],
+ gradSetArgs[[2L]], drop = FALSE), call("[", gradSetArgs[[1L]],
+ gradSetArgs[[2L]], gradSetArgs[[2L]], drop = FALSE),
+ call("[", gradSetArgs[[1L]], gradSetArgs[[2L]], gradSetArgs[[2L]],
+ gradSetArgs[[3L]], drop = FALSE), call("[", gradSetArgs[[1L]],
+ gradSetArgs[[2L]], gradSetArgs[[2L]], gradSetArgs[[3L]],
+ gradSetArgs[[4L]], drop = FALSE))
+ getRHS.varying <- function() {
+ ans <- getRHS.noVarying()
+ attr(ans, "gradient") <- eval(gradCall)
+ ans
+ }
+ QR <- qr(.swts * attr(rhs, "gradient"))
+ qrDim <- min(dim(QR$qr))
+ if (QR$rank < qrDim) stop("singular gradient matrix at initial parameter estimates")
+
+ getPars.varying <- function() unlist(setNames(lapply(names(ind),
+ get, envir = env), names(ind)))[useParams]
+ setPars.noVarying <- function(newPars) {
+ assign("internalPars", newPars, envir = thisEnv)
+ for (i in names(ind)) assign(i, unname(newPars[ind[[i]]]),
+ envir = env)
+ }
+ setPars.varying <- function(newPars) {
+ internalPars[useParams] <- newPars
+ for (i in names(ind)) assign(i, unname(internalPars[ind[[i]]]),
+ envir = env)
+ }
+ setPars <- setPars.noVarying
+ on.exit(remove(i, data, parLength, start, temp, m))
+ m <- list(resid = function() resid, fitted = function() rhs,
+ formula = function() form, deviance = function() dev,
+ lhs = function() lhs, gradient = function() .swts * attr(rhs,
+ "gradient"), conv = function() {
+ if (npar == 0) return(0)
+ rr <- qr.qty(QR, resid)
+ sqrt(sum(rr[1L:npar]^2)/sum(rr[-(1L:npar)]^2))
+ }, incr = function() qr.coef(QR, resid), setVarying = function(vary = rep(TRUE,
+ length(useParams))) {
+ assign("useParams", if (is.character(vary)) {
+ temp <- logical(length(useParams))
+ temp[unlist(ind[vary])] <- TRUE
+ temp
+ } else if (is.logical(vary) && length(vary) != length(useParams)) stop("setVarying : 'vary' length must match length of parameters") else {
+ vary
+ }, envir = thisEnv)
+ gradCall[[length(gradCall) - 1L]] <<- useParams
+ if (all(useParams)) {
+ assign("setPars", setPars.noVarying, envir = thisEnv)
+ assign("getPars", getPars.noVarying, envir = thisEnv)
+ assign("getRHS", getRHS.noVarying, envir = thisEnv)
+ assign("npar", length(useParams), envir = thisEnv)
+ } else {
+ assign("setPars", setPars.varying, envir = thisEnv)
+ assign("getPars", getPars.varying, envir = thisEnv)
+ assign("getRHS", getRHS.varying, envir = thisEnv)
+ assign("npar", length(seq_along(useParams)[useParams]),
+ envir = thisEnv)
+ }
+ }, setPars = function(newPars) {
+ setPars(newPars)
+ assign("resid", .swts * (lhs - assign("rhs", getRHS(),
+ envir = thisEnv)), envir = thisEnv)
+ assign("dev", sum(resid^2), envir = thisEnv)
+ assign("QR", qr(.swts * attr(rhs, "gradient")), envir = thisEnv)
+ return(QR$rank < min(dim(QR$qr)))
+ }, getPars = function() getPars(), getAllPars = function() getPars(),
+ getEnv = function() env, trace = function() {
+ cat(format(dev), ": ", format(getPars()))
+ cat("\n")
+ }, Rmat = function() qr.R(QR), predict = function(newdata = list(),
+ qr = FALSE) eval(form[[3L]], as.list(newdata), env))
+ class(m) <- "nlsModel"
+ m
+}
+
+nlsLM <- function (formula, data = parent.frame(), start, jac = NULL,
+ algorithm = "LM", control = nls.lm.control(),
+ lower = NULL, upper = NULL, trace = FALSE, subset, weights, na.action,
+ model = FALSE, ...)
+{
+ formula <- as.formula(formula)
+ if (!is.list(data) && !is.environment(data)) stop("'data' must be a list or an environment")
+ mf <- match.call()
+ varNames <- all.vars(formula)
+
+ if (length(formula) == 2L) {
+ formula[[3L]] <- formula[[2L]]
+ formula[[2L]] <- 0
+ }
+
+ form2 <- formula
+ form2[[2L]] <- 0
+ varNamesRHS <- all.vars(form2)
+ mWeights <- missing(weights)
+
+ ## if trace = TRUE, set nls.lm.control$nprint = 1
+ if (trace) control$nprint <- 1
+
+ pnames <- if (missing(start)) {
+ if (!is.null(attr(data, "parameters"))) {
+ names(attr(data, "parameters"))
+ }
+ else {
+ cll <- formula[[length(formula)]]
+ func <- get(as.character(cll[[1L]]))
+ if (!is.null(pn <- attr(func, "pnames")))
+ as.character(as.list(match.call(func, call = cll))[-1L][pn])
+ }
+ } else names(start)
+
+ env <- environment(formula)
+ if (is.null(env)) env <- parent.frame()
+ if (length(pnames)) varNames <- varNames[is.na(match(varNames, pnames))]
+ lenVar <- function(var) tryCatch(length(eval(as.name(var), data, env)), error = function(e) -1)
+
+ if (length(varNames)) {
+ n <- sapply(varNames, lenVar)
+ if (any(not.there <- n == -1)) {
+ nnn <- names(n[not.there])
+ if (missing(start)) {
+ warning("No starting values specified for some parameters.\n",
+ "Initializing ", paste(sQuote(nnn), collapse = ", "),
+ " to '1.'.\n", "Consider specifying 'start' or using a selfStart model")
+ start <- as.list(rep(1, length(nnn)))
+ names(start) <- nnn
+ varNames <- varNames[i <- is.na(match(varNames, nnn))]
+ n <- n[i]
+ }
+ else stop("parameters without starting value in 'data': ", paste(nnn, collapse = ", "))
+ }
+ } else {
+ if (length(pnames) && any((np <- sapply(pnames, lenVar)) == -1)) {
+ message("fitting parameters ", paste(sQuote(pnames[np == -1]), collapse = ", "), " without any variables")
+ n <- integer()
+ } else stop("no parameters to fit")
+ }
+
+ respLength <- length(eval(formula[[2L]], data, env))
+
+ if (length(n) > 0L) {
+ varIndex <- n%%respLength == 0
+ if (is.list(data) && diff(range(n[names(n) %in% names(data)])) > 0) {
+ mf <- data
+ if (!missing(subset)) warning("argument 'subset' will be ignored")
+ if (!missing(na.action)) warning("argument 'na.action' will be ignored")
+ if (missing(start)) start <- getInitial(formula, mf)
+ startEnv <- new.env(hash = FALSE, parent = environment(formula))
+ for (i in names(start)) assign(i, start[[i]], envir = startEnv)
+ rhs <- eval(formula[[3L]], data, startEnv)
+ n <- NROW(rhs)
+ wts <- if (mWeights) rep(1, n) else eval(substitute(weights), data, environment(formula))
+ } else {
+ mf$formula <- as.formula(paste("~", paste(varNames[varIndex], collapse = "+")), env = environment(formula))
+ mf$start <- mf$control <- mf$algorithm <- mf$trace <- mf$model <- NULL
+ mf$lower <- mf$upper <- NULL
+ mf[[1L]] <- as.name("model.frame")
+ mf <- eval.parent(mf)
+ n <- nrow(mf)
+ mf <- as.list(mf)
+ wts <- if (!mWeights) model.weights(mf) else rep(1, n)
+ }
+ if (any(wts < 0 | is.na(wts))) stop("missing or negative weights not allowed")
+ } else {
+ varIndex <- logical()
+ mf <- list(0)
+ wts <- numeric()
+ }
+
+ if (missing(start)) start <- getInitial(formula, mf)
+ for (var in varNames[!varIndex]) mf[[var]] <- eval(as.name(var), data, env)
+ varNamesRHS <- varNamesRHS[varNamesRHS %in% varNames[varIndex]]
+
+ ## added nls.lm from 'minpack.lm' package
+ mf <- c(mf, start)
+ lhs <- eval(formula[[2L]], envir = mf)
+ m <- match(names(start), names(mf))
+ .swts <- if (!missing(wts) && length(wts)) sqrt(wts)
+
+ FCT <- function(par) {
+ mf[m] <- par
+ rhs <- eval(formula[[3L]], envir = mf, environment(formula))
+ res <- lhs - rhs
+ res <- .swts * res
+ res
+ }
+
+ NLS <- nls.lm(par = start, fn = FCT, jac = jac, control = control, lower = lower, upper = upper, ...)
+ ## previous versions before boundaries were included
+ ## NLS <- nls.lm(par = start, fn = FCT, jac = jac, control = control, ...)
+ start <- NLS$par
+
+ ## pass optimized parameters to 'nlsModel'
+ m <- nlsModel(formula, mf, start, wts)
+
+ ## => internal 'nls' iterations by 'C_nls_iter' switched off
+ ## use 'nls.lm' output for convergence info
+ if (NLS$info %in% c(1, 2, 3, 4)) isConv <- TRUE else isConv <- FALSE
+ finIter <- NLS$niter
+ finTol <- nls.lm.control()$ftol
+ convInfo <- list(isConv = isConv, finIter = finIter, finTol = finTol,
+ stopCode = NLS$info, stopMessage = NLS$message)
+ nls.out <- list(m = m, convInfo = convInfo, data = substitute(data), call = match.call())
+
+ nls.out$call$algorithm <- algorithm
+ ## need to use '$tol' parameter from nls.control to make 'predict.nls' work
+ nls.out$call$control <- nls.control()
+ nls.out$call$trace <- FALSE
+ ## lower/upper changed in version 1.2-1
+ nls.out$call$lower <- lower
+ nls.out$call$upper <- upper
+ nls.out$na.action <- attr(mf, "na.action")
+ nls.out$dataClasses <- attr(attr(mf, "terms"), "dataClasses")[varNamesRHS]
+ if (model)
+ nls.out$model <- mf
+ if (!mWeights)
+ nls.out$weights <- wts
+ nls.out$control <- control
+ class(nls.out) <- "nls"
+ nls.out
+}
diff --git a/R/wfct.R b/R/wfct.R
new file mode 100644
index 0000000..74286b9
--- /dev/null
+++ b/R/wfct.R
@@ -0,0 +1,57 @@
+wfct <- function(expr)
+{
+expr <- deparse(substitute(expr))
+
+## create new environment
+newEnv <- new.env()
+
+## get call
+mc <- sys.calls()[[1]]
+mcL <- as.list(mc)
+
+## get data and write to newEnv
+DATA <- mcL[["data"]]
+DATA <- eval(DATA)
+DATA <- as.list(DATA)
+NAMES <- names(DATA)
+for (i in 1:length(DATA)) assign(NAMES[i], DATA[[i]], envir = newEnv)
+
+## get parameter, response and predictor names
+formula <- as.formula(mcL[[2]])
+VARS <- all.vars(formula)
+RESP <- VARS[1]
+RHS <- VARS[-1]
+PRED <- match(RHS, names(DATA))
+PRED <- names(DATA)[na.omit(PRED)]
+
+## calculate variances for response values if "error" is in expression
+## and write to newEnv
+if (length(grep("error", expr)) > 0) {
+y <- DATA[[RESP]]
+x <- DATA[[PRED]]
+## test for replication
+if (!any(duplicated(x))) stop("No replicates available to calculate error from!")
+## calculate error
+error <- tapply(y, x, function(e) var(e, na.rm = TRUE))
+error <- as.numeric(sqrt(error))
+## convert to original repititions
+error <- rep(error, as.numeric(table(x)))
+assign("error", error, envir = newEnv)
+}
+
+## calculate fitted or residual values if "fitted"/"resid" is in expression
+## and write to newEnv
+if (length(grep("fitted", expr)) > 0 || length(grep("resid", expr)) > 0) {
+mc2 <- mc
+mc2$weights <- NULL
+MODEL <- eval(mc2)
+fitted <- fitted(MODEL)
+resid <- residuals(MODEL)
+assign("fitted", fitted, newEnv)
+assign("resid", resid, newEnv)
+}
+
+## return evaluation in newEnv: vector of weights
+OUT <- eval(parse(text = expr), envir = newEnv)
+return(OUT)
+}
diff --git a/THANKS b/THANKS
new file mode 100644
index 0000000..9194cac
--- /dev/null
+++ b/THANKS
@@ -0,0 +1,8 @@
+Since version 1.1-6, the following people have helped improve minpack.lm:
+
+Romain Barillot
+Daniel Kaschek
+Andrej-Nikolai Spiess
+
+
+
diff --git a/debian/changelog b/debian/changelog
deleted file mode 100644
index a71552a..0000000
--- a/debian/changelog
+++ /dev/null
@@ -1,15 +0,0 @@
-r-cran-minpack.lm (1.2-1-1) unstable; urgency=medium
-
- * New upstream version
- * Convert to dh-r
- * Canonical homepage for CRAN
- * debhelper 10
- * d/watch: version=4
-
- -- Andreas Tille <tille at debian.org> Wed, 23 Nov 2016 21:13:13 +0100
-
-r-cran-minpack.lm (1.2-0-1) unstable; urgency=low
-
- * Initial release (closes: #829574)
-
- -- Andreas Tille <tille at debian.org> Mon, 04 Jul 2016 14:06:18 +0200
diff --git a/debian/compat b/debian/compat
deleted file mode 100644
index f599e28..0000000
--- a/debian/compat
+++ /dev/null
@@ -1 +0,0 @@
-10
diff --git a/debian/control b/debian/control
deleted file mode 100644
index bb2ddda..0000000
--- a/debian/control
+++ /dev/null
@@ -1,26 +0,0 @@
-Source: r-cran-minpack.lm
-Maintainer: Debian Med Packaging Team <debian-med-packaging at lists.alioth.debian.org>
-Uploaders: Andreas Tille <tille at debian.org>
-Section: gnu-r
-Priority: optional
-Build-Depends: debhelper (>= 10),
- dh-r,
- r-base-dev
-Standards-Version: 3.9.8
-Vcs-Browser: https://anonscm.debian.org/viewvc/debian-med/trunk/packages/R/r-cran-minpack.lm/trunk/
-Vcs-Svn: svn://anonscm.debian.org/debian-med/trunk/packages/R/r-cran-minpack.lm/trunk/
-Homepage: https://cran.r-project.org/package=minpack.lm
-
-Package: r-cran-minpack.lm
-Architecture: any
-Depends: ${shlibs:Depends},
- ${misc:Depends},
- ${R:Depends}
-Recommends: ${R:Recommends}
-Suggests: ${R:Suggests}
-Description: GNU R Levenberg-Marquardt nonlinear least-squares algorithm found in MINPACK
- The nls.lm function provides an R interface to lmder and lmdif from the
- MINPACK library, for solving nonlinear least-squares problems by a
- modification of the Levenberg-Marquardt algorithm, with support for
- lower and upper parameter bounds. The implementation can be used via nls-
- like calls using the nlsLM function.
diff --git a/debian/copyright b/debian/copyright
deleted file mode 100644
index 4961b00..0000000
--- a/debian/copyright
+++ /dev/null
@@ -1,27 +0,0 @@
-Format: https://www.debian.org/doc/packaging-manuals/copyright-format/1.0/
-Upstream-Name: minpack.lm
-Upstream-Contact: Katharine M. Mullen <mullenkate at gmail.com>
-Source: https://cran.r-project.org/package=minpack.lm
-
-Files: *
-Copyright: 2010-2016 Katharine M. Mullen <mullenkate at gmail.com>,
- Timur V. Elzhov, Andrej-Nikolai Spiess, Ben Bolker
-License: GPL-3+
-
-Files: debian/*
-Copyright: 2016 Andreas Tille <tille at debian.org>
-License: GPL-3+
-
-License: GPL-3+
- This program is free software; you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation; either version 3 of the License, or
- (at your option) any later version.
- .
- This program is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
- .
- On Debian systems, the complete text of the GNU General Public
- License can be found in `/usr/share/common-licenses/GPL-3'.
diff --git a/debian/rules b/debian/rules
deleted file mode 100755
index 68d9a36..0000000
--- a/debian/rules
+++ /dev/null
@@ -1,4 +0,0 @@
-#!/usr/bin/make -f
-
-%:
- dh $@ --buildsystem R
diff --git a/debian/source/format b/debian/source/format
deleted file mode 100644
index 163aaf8..0000000
--- a/debian/source/format
+++ /dev/null
@@ -1 +0,0 @@
-3.0 (quilt)
diff --git a/debian/watch b/debian/watch
deleted file mode 100644
index 9ace895..0000000
--- a/debian/watch
+++ /dev/null
@@ -1,2 +0,0 @@
-version=4
-http://cran.r-project.org/src/contrib/minpack.lm_([-\d.]*)\.tar\.gz
diff --git a/man/nls.lm.Rd b/man/nls.lm.Rd
new file mode 100644
index 0000000..737fc64
--- /dev/null
+++ b/man/nls.lm.Rd
@@ -0,0 +1,294 @@
+\name{nls.lm}
+\alias{nls.lm}
+\encoding{UTF-8}
+\title{Addresses NLS problems with the Levenberg-Marquardt algorithm}
+\description{
+ The purpose of \code{nls.lm} is to minimize the sum square of the
+ vector returned by the function \code{fn}, by a modification of the
+ Levenberg-Marquardt algorithm. The user may also provide a
+ function \code{jac} which calculates the Jacobian.
+}
+\usage{
+nls.lm(par, lower=NULL, upper=NULL, fn, jac = NULL,
+ control = nls.lm.control(), \dots)
+}
+\arguments{
+ \item{par}{A list or numeric vector of starting estimates. If
+ \code{par} is a list, then each element must be of length 1. }
+ \item{lower}{A numeric vector of lower bounds on each parameter. If
+ not given, the default lower bound for each parameter is set to
+ \code{-Inf}. }
+ \item{upper}{A numeric vector of upper bounds on each parameter. If
+ not given, the default upper bound for each parameter is set to
+ \code{Inf}. }
+ \item{fn}{A function that returns a vector of residuals, the sum square
+ of which is to be minimized. The first argument of \code{fn} must be
+ \code{par}. }
+ \item{jac}{A function to return the Jacobian for the \code{fn} function.}
+ \item{control}{
+ An optional list of control settings. See \code{\link{nls.lm.control}} for
+ the names of the settable control values and their effect.
+}
+ \item{\dots}{Further arguments to be passed to \code{fn} and \code{jac}.}
+}
+\details{
+ Both functions \code{fn} and \code{jac} (if provided) must return
+ numeric vectors. Length of the vector returned by \code{fn} must
+ not be lower than the length of \code{par}. The vector returned by
+ \code{jac} must have length equal to
+ \eqn{length(\code{fn}(\code{par}, \dots))\cdot length(\code{par})}{%
+ length(\code{fn}(\code{par}, \dots)) * length(\code{par})}.
+
+ The \code{control} argument is a list; see \code{\link{nls.lm.control}} for
+ details.
+
+ \bold{Successful completion.}\cr
+ \cr
+ The accuracy of \code{nls.lm} is controlled by the convergence
+ parameters \code{ftol}, \code{ptol}, and \code{gtol}. These
+ parameters are used in tests which make three types of comparisons
+ between the approximation \eqn{par} and a solution
+ \eqn{par_0}{par0}. \code{nls.lm} terminates when any of the tests
+ is satisfied. If any of the convergence parameters is less than
+ the machine precision, then \code{nls.lm} only attempts to satisfy
+ the test defined by the machine precision. Further progress is not
+ usually possible.\cr
+ The tests assume that \code{fn} as well as \code{jac} are
+ reasonably well behaved. If this condition is not satisfied, then
+ \code{nls.lm} may incorrectly indicate convergence. The validity
+ of the answer can be checked, for example, by rerunning
+ \code{nls.lm} with tighter tolerances.\cr
+ \cr
+ \emph{First convergence test.}\cr
+ If \eqn{|z|} denotes the Euclidean norm of a vector \eqn{z}, then
+ this test attempts to guarantee that
+ \deqn{|fvec| < (1 + \code{ftol})\,|fvec_0|,}{%
+ |fvec| < (1 + \code{ftol})|fvec0|,}
+ where \eqn{fvec_0}{fvec0} denotes the result of \code{fn} function
+ evaluated at \eqn{par_0}{par0}. If this condition is satisfied
+ with \code{ftol} \eqn{\simeq 10^{-k}}{~ 10^(-k)}, then the final
+ residual norm \eqn{|fvec|} has \eqn{k} significant decimal digits
+ and \code{info} is set to 1 (or to 3 if the second test is also
+ satisfied). Unless high precision solutions are required, the
+ recommended value for \code{ftol} is the square root of the machine
+ precision.\cr
+ \cr
+ \emph{Second convergence test.}\cr
+ If \eqn{D} is the diagonal matrix whose entries are defined by the
+ array \code{diag}, then this test attempt to guarantee that
+ \deqn{|D\,(par - par_0)| < \code{ptol}\,|D\,par_0|,}{%
+ |D*(par - par0)| < \code{ptol}|D*par0|,}
+ If this condition is satisfied with \code{ptol} \eqn{\simeq
+ 10^{-k}}{~ 10^(-k)}, then the larger components of
+ \eqn{(D\,par)}{D*par} have \eqn{k} significant decimal digits and
+ \code{info} is set to 2 (or to 3 if the first test is also
+ satisfied). There is a danger that the smaller components of
+ \eqn{(D\,par)}{D*par} may have large relative errors, but if
+ \code{diag} is internally set, then the accuracy of the components
+ of \eqn{par} is usually related to their sensitivity. Unless high
+ precision solutions are required, the recommended value for
+ \code{ptol} is the square root of the machine precision.\cr
+ \cr
+ \emph{Third convergence test.}\cr
+ This test is satisfied when the cosine of the angle between the
+ result of \code{fn} evaluation \eqn{fvec} and any column of the
+ Jacobian at \eqn{par} is at most \code{gtol} in absolute value.
+ There is no clear relationship between this test and the accuracy
+ of \code{nls.lm}, and furthermore, the test is equally well
+ satisfied at other critical points, namely maximizers and saddle
+ points. Therefore, termination caused by this test (\code{info} =
+ 4) should be examined carefully. The recommended value for
+ \code{gtol} is zero.\cr
+ \cr
+ \bold{Unsuccessful completion.}\cr
+ \cr
+ Unsuccessful termination of \code{nls.lm} can be due to improper
+ input parameters, arithmetic interrupts, an excessive number of
+ function evaluations, or an excessive number of iterations. \cr
+ \cr
+ \emph{Improper input parameters.}\cr
+ \code{info} is set to 0 if \eqn{length(\code{par}) = 0}, or
+ \eqn{length(fvec) < length(\code{par})}, or \code{ftol} \eqn{< 0},
+ or \code{ptol} \eqn{< 0}, or \code{gtol} \eqn{< 0}, or \code{maxfev}
+ \eqn{\leq 0}{<= 0}, or \code{factor} \eqn{\leq 0}{<= 0}.\cr
+ \cr
+ \emph{Arithmetic interrupts.}\cr
+ If these interrupts occur in the \code{fn} function during an
+ early stage of the computation, they may be caused by an
+ unacceptable choice of \eqn{par} by \code{nls.lm}. In this case,
+ it may be possible to remedy the situation by rerunning
+ \code{nls.lm} with a smaller value of \code{factor}.\cr
+ \cr
+ \emph{Excessive number of function evaluations.}\cr
+ A reasonable value for \code{maxfev} is \eqn{100\cdot
+ (length(\code{par}) + 1)}{100*(length(\code{par}) + 1)}. If the
+ number of calls to \code{fn} reaches \code{maxfev}, then this
+ indicates that the routine is converging very slowly as measured
+ by the progress of \eqn{fvec} and \code{info} is set to 5. In this
+ case, it may be helpful to force \code{diag} to be internally set.
+
+ \emph{Excessive number of function iterations.}\cr
+ The allowed number of iterations defaults to 50, can be increased if
+ desired. \cr
+
+ The list returned by \code{nls.lm} has methods
+ for the generic functions \code{\link{coef}},
+ \code{\link{deviance}}, \code{\link{df.residual}},
+ \code{\link{print}}, \code{\link{residuals}}, \code{\link{summary}},
+ \code{\link{confint}},
+ and \code{\link{vcov}}.
+
+ }
+\value{
+ A list with components:
+ \item{par}{The best set of parameters found.}
+ \item{hessian}{A symmetric matrix giving an estimate of the Hessian
+ at the solution found.}
+ \item{fvec}{The result of the last \code{fn} evaluation; that is, the
+ residuals. }
+ \item{info}{\code{info} is an integer code indicating
+ the reason for termination.
+ \describe{
+ \item{0}{Improper input parameters.}
+ \item{1}{Both actual and predicted relative reductions in the
+ sum of squares are at most \code{ftol}.}
+ \item{2}{Relative error between two consecutive iterates is
+ at most \code{ptol}.}
+ \item{3}{Conditions for \code{info} = 1 and \code{info} = 2 both hold.}
+ \item{4}{The cosine of the angle between \code{fvec} and any column
+ of the Jacobian is at most \code{gtol} in absolute value.}
+ \item{5}{Number of calls to \code{fn} has reached \code{maxfev}.}
+ \item{6}{\code{ftol} is too small. No further reduction in the sum
+ of squares is possible.}
+ \item{7}{\code{ptol} is too small. No further improvement in the
+ approximate solution \code{par} is possible.}
+ \item{8}{\code{gtol} is too small. \code{fvec} is orthogonal to the
+ columns of the Jacobian to machine precision.}
+ \item{9}{The number of iterations has reached \code{maxiter}.}
+ }}
+ \item{message}{character string indicating reason for termination}.
+ \item{diag}{The result list of \code{diag}. See \bold{Details}.}
+ \item{niter}{The number of iterations completed before termination.}
+ \item{rsstrace}{The residual sum of squares at each iteration.
+ Can be used to check the progress each iteration. }
+ \item{deviance}{The sum of the squared residual vector.}
+}
+\references{
+ J.J. Moré, "The Levenberg-Marquardt algorithm: implementation and
+ theory," in \emph{Lecture Notes in Mathematics}
+ \bold{630}: Numerical Analysis, G.A. Watson (Ed.),
+ Springer-Verlag: Berlin, 1978, pp. 105-116.
+}
+\note{
+ The public domain FORTRAN sources of MINPACK package by J.J. Moré,
+ implementing the Levenberg-Marquardt algorithm were downloaded from
+ \url{http://netlib.org/minpack/}, and left unchanged.
+ The contents of this manual page are largely extracted from
+ the comments of MINPACK sources.
+}
+
+\seealso{\code{\link{optim}}, \code{\link{nls}}, \code{\link{nls.lm.control}}}
+
+\examples{
+
+###### example 1
+
+## values over which to simulate data
+x <- seq(0,5,length=100)
+
+## model based on a list of parameters
+getPred <- function(parS, xx) parS$a * exp(xx * parS$b) + parS$c
+
+## parameter values used to simulate data
+pp <- list(a=9,b=-1, c=6)
+
+## simulated data, with noise
+simDNoisy <- getPred(pp,x) + rnorm(length(x),sd=.1)
+
+## plot data
+plot(x,simDNoisy, main="data")
+
+## residual function
+residFun <- function(p, observed, xx) observed - getPred(p,xx)
+
+## starting values for parameters
+parStart <- list(a=3,b=-.001, c=1)
+
+## perform fit
+nls.out <- nls.lm(par=parStart, fn = residFun, observed = simDNoisy,
+xx = x, control = nls.lm.control(nprint=1))
+
+## plot model evaluated at final parameter estimates
+lines(x,getPred(as.list(coef(nls.out)), x), col=2, lwd=2)
+
+## summary information on parameter estimates
+summary(nls.out)
+
+###### example 2
+
+## function to simulate data
+f <- function(TT, tau, N0, a, f0) {
+ expr <- expression(N0*exp(-TT/tau)*(1 + a*cos(f0*TT)))
+ eval(expr)
+}
+
+## helper function for an analytical gradient
+j <- function(TT, tau, N0, a, f0) {
+ expr <- expression(N0*exp(-TT/tau)*(1 + a*cos(f0*TT)))
+ c(eval(D(expr, "tau")), eval(D(expr, "N0" )),
+ eval(D(expr, "a" )), eval(D(expr, "f0" )))
+}
+
+## values over which to simulate data
+TT <- seq(0, 8, length=501)
+
+## parameter values underlying simulated data
+p <- c(tau = 2.2, N0 = 1000, a = 0.25, f0 = 8)
+
+## get data
+Ndet <- do.call("f", c(list(TT = TT), as.list(p)))
+## with noise
+N <- Ndet + rnorm(length(Ndet), mean=Ndet, sd=.01*max(Ndet))
+
+## plot the data to fit
+par(mfrow=c(2,1), mar = c(3,5,2,1))
+plot(TT, N, bg = "black", cex = 0.5, main="data")
+
+## define a residual function
+fcn <- function(p, TT, N, fcall, jcall)
+ (N - do.call("fcall", c(list(TT = TT), as.list(p))))
+
+## define analytical expression for the gradient
+fcn.jac <- function(p, TT, N, fcall, jcall)
+ -do.call("jcall", c(list(TT = TT), as.list(p)))
+
+## starting values
+guess <- c(tau = 2.2, N0 = 1500, a = 0.25, f0 = 10)
+
+## to use an analytical expression for the gradient found in fcn.jac
+## uncomment jac = fcn.jac
+out <- nls.lm(par = guess, fn = fcn, jac = fcn.jac,
+ fcall = f, jcall = j,
+ TT = TT, N = N, control = nls.lm.control(nprint=1))
+
+## get the fitted values
+N1 <- do.call("f", c(list(TT = TT), out$par))
+
+## add a blue line representing the fitting values to the plot of data
+lines(TT, N1, col="blue", lwd=2)
+
+## add a plot of the log residual sum of squares as it is made to
+## decrease each iteration; note that the RSS at the starting parameter
+## values is also stored
+plot(1:(out$niter+1), log(out$rsstrace), type="b",
+main="log residual sum of squares vs. iteration number",
+xlab="iteration", ylab="log residual sum of squares", pch=21,bg=2)
+
+## get information regarding standard errors
+summary(out)
+
+}
+\keyword{nonlinear}
+\keyword{optimize}
+\keyword{regression}
+
diff --git a/man/nls.lm.control.Rd b/man/nls.lm.control.Rd
new file mode 100644
index 0000000..4241f51
--- /dev/null
+++ b/man/nls.lm.control.Rd
@@ -0,0 +1,85 @@
+\name{nls.lm.control}
+\alias{nls.lm.control}
+\encoding{UTF-8}
+\title{Control various aspects of the Levenberg-Marquardt algorithm}
+\description{
+ Allow the user to set some characteristics
+ Levenberg-Marquardt nonlinear least squares algorithm implemented
+ in \code{nls.lm}.
+}
+\usage{
+nls.lm.control(ftol = sqrt(.Machine$double.eps),
+ptol = sqrt(.Machine$double.eps), gtol = 0, diag = list(), epsfcn = 0,
+factor = 100, maxfev = integer(), maxiter = 50, nprint = 0)
+}
+\arguments{
+ \item{ftol}{non-negative numeric. Termination occurs when
+ both the actual and predicted relative reductions in the sum of
+ squares are at most \code{ftol}. Therefore, \code{ftol} measures
+ the relative error desired in the sum of squares.}
+ \item{ptol}{non-negative numeric. Termination occurs when
+ the relative error between two consecutive iterates is at most
+ \code{ptol}. Therefore, \code{ptol} measures the relative error
+ desired in the approximate solution.}
+ \item{gtol}{non-negative numeric. Termination occurs when
+ the cosine of the angle between result of \code{fn} evaluation
+ \eqn{fvec} and any column of the Jacobian is at most \code{gtol}
+ in absolute value. Therefore, \code{gtol} measures the
+ orthogonality desired between the function vector and the
+ columns of the Jacobian.}
+ \item{diag}{a list or numeric vector containing positive
+ entries that serve as multiplicative scale factors for the
+ parameters. Length of \code{diag} should be equal to that of
+ \code{par}. If not, user-provided \code{diag} is ignored and
+ \code{diag} is internally set.}
+ \item{epsfcn}{(used if \code{jac} is not provided) is a
+ numeric used in determining a suitable step for the
+ forward-difference approximation. This approximation assumes
+ that the relative errors in the functions are of the order of
+ \code{epsfcn}. If \code{epsfcn} is less than the machine
+ precision, it is assumed that the relative errors in the
+ functions are of the order of the machine precision.}
+ \item{factor}{positive numeric, used in determining the
+ initial step bound. This bound is set to the product of
+ \code{factor} and the \eqn{|\code{diag}*\code{par}|} if nonzero,
+ or else to \code{factor} itself. In most cases \code{factor}
+ should lie in the interval (0.1,100). 100 is a generally
+ recommended value.}
+ \item{maxfev}{integer; termination occurs
+ when the number of calls to \code{fn} has reached \code{maxfev}.
+ Note that \code{nls.lm} sets the value of \code{maxfev} to
+ \code{100*(length(par) + 1)} if
+ \code{maxfev = integer()}, where \code{par} is the list or
+ vector of parameters to be optimized. }
+ \item{maxiter}{positive integer. Termination occurs
+ when the number of iterations reaches \code{maxiter}.}
+ \item{nprint}{is an integer; set \code{nprint} to be positive
+ to enable printing of iterates}
+ }
+ \value{
+ A \code{list} with exactly nine components:
+ \item{ftol}{}
+ \item{ptol}{}
+ \item{gtol}{}
+ \item{diag}{}
+ \item{epsfcn}{}
+ \item{factor}{}
+ \item{maxfev}{}
+ \item{nprint}{}
+ with meanings as explained under \sQuote{Arguments}.
+}
+\references{
+ J.J. Moré, "The Levenberg-Marquardt algorithm: implementation and
+ theory," in \emph{Lecture Notes in Mathematics}
+ \bold{630}: Numerical Analysis, G.A. Watson (Ed.),
+ Springer-Verlag: Berlin, 1978, pp. 105-116.
+}
+
+\seealso{\code{\link{nls.lm}}}
+\examples{
+nls.lm.control(maxiter = 4)
+}
+\keyword{nonlinear}
+\keyword{optimize}
+\keyword{regression}
+
diff --git a/man/nlsLM.Rd b/man/nlsLM.Rd
new file mode 100644
index 0000000..ada3e43
--- /dev/null
+++ b/man/nlsLM.Rd
@@ -0,0 +1,165 @@
+\name{nlsLM}
+\alias{nlsLM}
+\encoding{UTF-8}
+\title{Standard 'nls' framework that uses 'nls.lm' for fitting}
+\description{
+ \code{nlsLM} is a modified version of \code{\link{nls}} that uses \code{nls.lm} for fitting.
+ Since an object of class 'nls' is returned, all generic functions such as \code{\link{anova}},
+ \code{\link{coef}}, \code{\link{confint}}, \code{\link{deviance}}, \code{\link{df.residual}},
+ \code{\link{fitted}}, \code{\link{formula}}, \code{\link{logLik}}, \code{\link{predict}},
+ \code{\link{print}}, \code{\link{profile}}, \code{\link{residuals}}, \code{\link{summary}},
+ \code{\link{update}}, \code{\link{vcov}} and \code{\link{weights}} are applicable.
+}
+
+\usage{
+nlsLM(formula, data = parent.frame(), start, jac = NULL,
+ algorithm = "LM", control = nls.lm.control(),
+ lower = NULL, upper = NULL, trace = FALSE, subset,
+ weights, na.action, model = FALSE, \dots)
+}
+
+\arguments{
+ \item{formula}{a nonlinear model \code{\link{formula}} including variables and parameters. Will be coerced to a formula if necessary.}
+ \item{data}{an optional data frame in which to evaluate the variables in\code{formula} and \code{weights}.
+ Can also be a list or an environment, but not a matrix.}
+ \item{start}{a named list or named numeric vector of starting estimates.}
+ \item{jac}{A function to return the Jacobian.}
+ \item{algorithm}{only method \code{"LM"} (Levenberg-Marquardt) is implemented.}
+ \item{control}{an optional list of control settings. See \code{\link{nls.lm.control}} for
+ the names of the settable control values and their effect.}
+ \item{lower}{A numeric vector of lower bounds on each parameter. If not given, the default lower bound for each parameter is set to \code{-Inf}.}
+ \item{upper}{A numeric vector of upper bounds on each parameter. If not given, the default upper bound for each parameter is set to \code{Inf}.}
+ \item{trace}{logical value indicating if a trace of the iteration progress should be printed. Default is \code{FALSE}. If \code{TRUE}, the residual (weighted) sum-of-squares and the parameter values are printed at the conclusion of each iteration.}
+ \item{subset}{an optional vector specifying a subset of observations to be used in the fitting process.}
+ \item{weights}{an optional numeric vector of (fixed) weights. When
+ present, the objective function is weighted least squares. See the
+ \code{wfct} function for options for easy specification of weighting
+ schemes. }
+ \item{na.action}{a function which indicates what should happen when the data contain \code{NA}s. The default is set by the \code{na.action} setting of \code{\link{options}}, and is \code{\link{na.fail}} if that is unset. The \sQuote{factory-fresh} default is \code{\link{na.omit}}. Value \code{\link{na.exclude}} can be useful.}
+ \item{model}{logical. If true, the model frame is returned as part of the object. Default is \code{FALSE}.}
+ \item{\dots}{Additional optional arguments. None are used at present.}
+}
+
+\details{
+The standard \code{\link{nls}} function was modified in several ways to incorporate the Levenberg-Marquardt type \code{\link{nls.lm}} fitting algorithm. The \code{formula} is transformed into a function that returns a vector of (weighted) residuals whose sum square is minimized by \code{\link{nls.lm}}. The optimized parameters are then transferred
+to \code{nlsModel} in order to obtain an object of class 'nlsModel'. The internal C function \code{C_nls_iter} and \code{nls_port_fit} were removed to avoid subsequent "Gauss-Newton", "port" or "plinear" types of optimization of \code{nlsModel}. Several other small modifications were made in order to make all generic functions work on the output.
+}
+
+\value{
+ A list of
+ \item{m}{an \code{nlsModel} object incorporating the model.}
+ \item{data}{the expression that was passed to \code{nls} as the data argument. The actual data values are present in the environment of the \code{m} component.}
+ \item{call}{the matched call.}
+ \item{convInfo}{a list with convergence information.}
+ \item{control}{the control \code{list} used, see the \code{control} argument.}
+ \item{na.action}{the \code{"na.action"} attribute (if any) of the model frame.}
+ \item{dataClasses}{the \code{"dataClasses"} attribute (if any) of the \code{"terms"} attribute of the model frame.}
+ \item{model}{if \code{model = TRUE}, the model frame.}
+ \item{weights}{if \code{weights} is supplied, the weights.}
+}
+
+\references{
+Bates, D. M. and Watts, D. G. (1988)
+\emph{Nonlinear Regression Analysis and Its Applications}, Wiley
+
+Bates, D. M. and Chambers, J. M. (1992)
+\emph{Nonlinear models.}
+Chapter 10 of \emph{Statistical Models in S} eds J. M. Chambers and T. J. Hastie, Wadsworth & Brooks/Cole.
+
+J.J. More, "The Levenberg-Marquardt algorithm: implementation and theory,"
+in \emph{Lecture Notes in Mathematics}
+\bold{630}: Numerical Analysis, G.A. Watson (Ed.), Springer-Verlag: Berlin, 1978, pp. 105-116.
+}
+
+\author{Andrej-Nikolai Spiess and Katharine M. Mullen}
+
+\seealso{\code{\link{nls.lm}}, \code{\link{nls}}, \code{\link{nls.lm.control}}, \code{\link{optim}}}
+
+\examples{
+### Examples from 'nls' doc ###
+DNase1 <- subset(DNase, Run == 1)
+## using a selfStart model
+fm1DNase1 <- nlsLM(density ~ SSlogis(log(conc), Asym, xmid, scal), DNase1)
+## using logistic formula
+fm2DNase1 <- nlsLM(density ~ Asym/(1 + exp((xmid - log(conc))/scal)),
+ data = DNase1,
+ start = list(Asym = 3, xmid = 0, scal = 1))
+
+## all generics are applicable
+coef(fm1DNase1)
+confint(fm1DNase1)
+deviance(fm1DNase1)
+df.residual(fm1DNase1)
+fitted(fm1DNase1)
+formula(fm1DNase1)
+logLik(fm1DNase1)
+predict(fm1DNase1)
+print(fm1DNase1)
+profile(fm1DNase1)
+residuals(fm1DNase1)
+summary(fm1DNase1)
+update(fm1DNase1)
+vcov(fm1DNase1)
+weights(fm1DNase1)
+
+## weighted nonlinear regression using
+## inverse squared variance of the response
+## gives same results as original 'nls' function
+Treated <- Puromycin[Puromycin$state == "treated", ]
+var.Treated <- tapply(Treated$rate, Treated$conc, var)
+var.Treated <- rep(var.Treated, each = 2)
+Pur.wt1 <- nls(rate ~ (Vm * conc)/(K + conc), data = Treated,
+ start = list(Vm = 200, K = 0.1), weights = 1/var.Treated^2)
+Pur.wt2 <- nlsLM(rate ~ (Vm * conc)/(K + conc), data = Treated,
+ start = list(Vm = 200, K = 0.1), weights = 1/var.Treated^2)
+all.equal(coef(Pur.wt1), coef(Pur.wt2))
+
+## 'nlsLM' can fit zero-noise data
+## in contrast to 'nls'
+x <- 1:10
+y <- 2*x + 3
+\dontrun{
+nls(y ~ a + b * x, start = list(a = 0.12345, b = 0.54321))
+}
+nlsLM(y ~ a + b * x, start = list(a = 0.12345, b = 0.54321))
+
+### Examples from 'nls.lm' doc
+## values over which to simulate data
+x <- seq(0,5, length = 100)
+## model based on a list of parameters
+getPred <- function(parS, xx) parS$a * exp(xx * parS$b) + parS$c
+## parameter values used to simulate data
+pp <- list(a = 9,b = -1, c = 6)
+## simulated data with noise
+simDNoisy <- getPred(pp, x) + rnorm(length(x), sd = .1)
+## make model
+mod <- nlsLM(simDNoisy ~ a * exp(b * x) + c,
+ start = c(a = 3, b = -0.001, c = 1),
+ trace = TRUE)
+## plot data
+plot(x, simDNoisy, main = "data")
+## plot fitted values
+lines(x, fitted(mod), col = 2, lwd = 2)
+
+## create declining cosine
+## with noise
+TT <- seq(0, 8, length = 501)
+tau <- 2.2
+N0 <- 1000
+a <- 0.25
+f0 <- 8
+Ndet <- N0 * exp(-TT/tau) * (1 + a * cos(f0 * TT))
+N <- Ndet + rnorm(length(Ndet), mean = Ndet, sd = .01 * max(Ndet))
+## make model
+mod <- nlsLM(N ~ N0 * exp(-TT/tau) * (1 + a * cos(f0 * TT)),
+ start = c(tau = 2.2, N0 = 1500, a = 0.25, f0 = 10),
+ trace = TRUE)
+
+## plot data
+plot(TT, N, main = "data")
+## plot fitted values
+lines(TT, fitted(mod), col = 2, lwd = 2)
+}
+\keyword{nonlinear}
+\keyword{optimize}
+\keyword{regression}
diff --git a/man/wfct.Rd b/man/wfct.Rd
new file mode 100644
index 0000000..632be63
--- /dev/null
+++ b/man/wfct.Rd
@@ -0,0 +1,69 @@
+\name{wfct}
+\alias{wfct}
+\encoding{UTF-8}
+\title{Weighting function that can be supplied to the \code{weights} argument
+ of \code{\link{nlsLM}} or \code{\link{nls}}}
+\description{
+ \code{wfct} can be supplied to the \code{weights} argument of
+ \code{\link{nlsLM}} or \code{\link{nls}}, and facilitates specification of
+ weighting schemes.
+}
+
+\usage{
+wfct(expr)
+}
+
+\arguments{
+ \item{expr}{An expression specifying the weighting scheme
+ as described in the Details section below.}
+}
+
+\details{
+ The weighting function can take 5 different variable definitions and
+ combinations thereof:
+ \itemize{
+ \item the name of the predictor (independent) variable
+ \item the name of the response (dependent) variable
+ \item error: if replicates \eqn{y_{ij}}{y_{ij}}
+ exist, the error \eqn{\sigma(y_{ij})}{sigma_{ij}}
+ \item fitted: the fitted values \eqn{\hat{y}_i}{hat{y}_i} of the model
+ \item resid: the residuals \eqn{y_i - \hat{y}_i}{y_i - hat{y}_i} of the model
+ }
+For the last two, the model is fit unweighted, fitted values and
+residuals are extracted and the model is refit by the defined weights.
+}
+
+\value{The results of evaluation of \code{expr} in a new
+ environment, yielding the vector of weights to be applied.
+
+}
+\author{Andrej-Nikolai Spiess}
+
+\seealso{\code{\link{nlsLM}}, \code{\link{nls}}}
+
+\examples{
+
+### Examples from 'nls' doc ###
+## note that 'nlsLM' below may be replaced with calls to 'nls'
+Treated <- Puromycin[Puromycin$state == "treated", ]
+
+## Weighting by inverse of response 1/y_i:
+nlsLM(rate ~ Vm * conc/(K + conc), data = Treated,
+start = c(Vm = 200, K = 0.05), weights = wfct(1/rate))
+
+## Weighting by square root of predictor \sqrt{x_i}:
+nlsLM(rate ~ Vm * conc/(K + conc), data = Treated,
+start = c(Vm = 200, K = 0.05), weights = wfct(sqrt(conc)))
+
+## Weighting by inverse square of fitted values 1/\hat{y_i}^2:
+nlsLM(rate ~ Vm * conc/(K + conc), data = Treated,
+start = c(Vm = 200, K = 0.05), weights = wfct(1/fitted^2))
+
+## Weighting by inverse variance 1/\sigma{y_i}^2:
+nlsLM(rate ~ Vm * conc/(K + conc), data = Treated,
+start = c(Vm = 200, K = 0.05), weights = wfct(1/error^2))
+
+}
+\keyword{nonlinear}
+\keyword{optimize}
+\keyword{regression}
diff --git a/src/Makevars b/src/Makevars
new file mode 100644
index 0000000..f8469b5
--- /dev/null
+++ b/src/Makevars
@@ -0,0 +1,2 @@
+PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)
+
diff --git a/src/chkder.f b/src/chkder.f
new file mode 100644
index 0000000..29578fc
--- /dev/null
+++ b/src/chkder.f
@@ -0,0 +1,140 @@
+ subroutine chkder(m,n,x,fvec,fjac,ldfjac,xp,fvecp,mode,err)
+ integer m,n,ldfjac,mode
+ double precision x(n),fvec(m),fjac(ldfjac,n),xp(n),fvecp(m),
+ * err(m)
+c **********
+c
+c subroutine chkder
+c
+c this subroutine checks the gradients of m nonlinear functions
+c in n variables, evaluated at a point x, for consistency with
+c the functions themselves. the user must call chkder twice,
+c first with mode = 1 and then with mode = 2.
+c
+c mode = 1. on input, x must contain the point of evaluation.
+c on output, xp is set to a neighboring point.
+c
+c mode = 2. on input, fvec must contain the functions and the
+c rows of fjac must contain the gradients
+c of the respective functions each evaluated
+c at x, and fvecp must contain the functions
+c evaluated at xp.
+c on output, err contains measures of correctness of
+c the respective gradients.
+c
+c the subroutine does not perform reliably if cancellation or
+c rounding errors cause a severe loss of significance in the
+c evaluation of a function. therefore, none of the components
+c of x should be unusually small (in particular, zero) or any
+c other value which may cause loss of significance.
+c
+c the subroutine statement is
+c
+c subroutine chkder(m,n,x,fvec,fjac,ldfjac,xp,fvecp,mode,err)
+c
+c where
+c
+c m is a positive integer input variable set to the number
+c of functions.
+c
+c n is a positive integer input variable set to the number
+c of variables.
+c
+c x is an input array of length n.
+c
+c fvec is an array of length m. on input when mode = 2,
+c fvec must contain the functions evaluated at x.
+c
+c fjac is an m by n array. on input when mode = 2,
+c the rows of fjac must contain the gradients of
+c the respective functions evaluated at x.
+c
+c ldfjac is a positive integer input parameter not less than m
+c which specifies the leading dimension of the array fjac.
+c
+c xp is an array of length n. on output when mode = 1,
+c xp is set to a neighboring point of x.
+c
+c fvecp is an array of length m. on input when mode = 2,
+c fvecp must contain the functions evaluated at xp.
+c
+c mode is an integer input variable set to 1 on the first call
+c and 2 on the second. other values of mode are equivalent
+c to mode = 1.
+c
+c err is an array of length m. on output when mode = 2,
+c err contains measures of correctness of the respective
+c gradients. if there is no severe loss of significance,
+c then if err(i) is 1.0 the i-th gradient is correct,
+c while if err(i) is 0.0 the i-th gradient is incorrect.
+c for values of err between 0.0 and 1.0, the categorization
+c is less certain. in general, a value of err(i) greater
+c than 0.5 indicates that the i-th gradient is probably
+c correct, while a value of err(i) less than 0.5 indicates
+c that the i-th gradient is probably incorrect.
+c
+c subprograms called
+c
+c minpack supplied ... dpmpar
+c
+c fortran supplied ... dabs,dlog10,dsqrt
+c
+c argonne national laboratory. minpack project. march 1980.
+c burton s. garbow, kenneth e. hillstrom, jorge j. more
+c
+c **********
+ integer i,j
+ double precision eps,epsf,epslog,epsmch,factor,one,temp,zero
+ double precision dpmpar
+ data factor,one,zero /1.0d2,1.0d0,0.0d0/
+c
+c epsmch is the machine precision.
+c
+ epsmch = dpmpar(1)
+c
+ eps = dsqrt(epsmch)
+c
+ if (mode .eq. 2) go to 20
+c
+c mode = 1.
+c
+ do 10 j = 1, n
+ temp = eps*dabs(x(j))
+ if (temp .eq. zero) temp = eps
+ xp(j) = x(j) + temp
+ 10 continue
+ go to 70
+ 20 continue
+c
+c mode = 2.
+c
+ epsf = factor*epsmch
+ epslog = dlog10(eps)
+ do 30 i = 1, m
+ err(i) = zero
+ 30 continue
+ do 50 j = 1, n
+ temp = dabs(x(j))
+ if (temp .eq. zero) temp = one
+ do 40 i = 1, m
+ err(i) = err(i) + temp*fjac(i,j)
+ 40 continue
+ 50 continue
+ do 60 i = 1, m
+ temp = one
+ if (fvec(i) .ne. zero .and. fvecp(i) .ne. zero
+ * .and. dabs(fvecp(i)-fvec(i)) .ge. epsf*dabs(fvec(i)))
+ * temp = eps*dabs((fvecp(i)-fvec(i))/eps-err(i))
+ * /(dabs(fvec(i)) + dabs(fvecp(i)))
+ err(i) = one
+ if (temp .gt. epsmch .and. temp .lt. eps)
+ * err(i) = (dlog10(temp) - epslog)/epslog
+ if (temp .ge. eps) err(i) = zero
+ 60 continue
+ 70 continue
+c
+ return
+c
+c last card of subroutine chkder.
+c
+ end
diff --git a/src/dogleg.f b/src/dogleg.f
new file mode 100644
index 0000000..b812f19
--- /dev/null
+++ b/src/dogleg.f
@@ -0,0 +1,177 @@
+ subroutine dogleg(n,r,lr,diag,qtb,delta,x,wa1,wa2)
+ integer n,lr
+ double precision delta
+ double precision r(lr),diag(n),qtb(n),x(n),wa1(n),wa2(n)
+c **********
+c
+c subroutine dogleg
+c
+c given an m by n matrix a, an n by n nonsingular diagonal
+c matrix d, an m-vector b, and a positive number delta, the
+c problem is to determine the convex combination x of the
+c gauss-newton and scaled gradient directions that minimizes
+c (a*x - b) in the least squares sense, subject to the
+c restriction that the euclidean norm of d*x be at most delta.
+c
+c this subroutine completes the solution of the problem
+c if it is provided with the necessary information from the
+c qr factorization of a. that is, if a = q*r, where q has
+c orthogonal columns and r is an upper triangular matrix,
+c then dogleg expects the full upper triangle of r and
+c the first n components of (q transpose)*b.
+c
+c the subroutine statement is
+c
+c subroutine dogleg(n,r,lr,diag,qtb,delta,x,wa1,wa2)
+c
+c where
+c
+c n is a positive integer input variable set to the order of r.
+c
+c r is an input array of length lr which must contain the upper
+c triangular matrix r stored by rows.
+c
+c lr is a positive integer input variable not less than
+c (n*(n+1))/2.
+c
+c diag is an input array of length n which must contain the
+c diagonal elements of the matrix d.
+c
+c qtb is an input array of length n which must contain the first
+c n elements of the vector (q transpose)*b.
+c
+c delta is a positive input variable which specifies an upper
+c bound on the euclidean norm of d*x.
+c
+c x is an output array of length n which contains the desired
+c convex combination of the gauss-newton direction and the
+c scaled gradient direction.
+c
+c wa1 and wa2 are work arrays of length n.
+c
+c subprograms called
+c
+c minpack-supplied ... dpmpar,enorm
+c
+c fortran-supplied ... dabs,dmax1,dmin1,dsqrt
+c
+c argonne national laboratory. minpack project. march 1980.
+c burton s. garbow, kenneth e. hillstrom, jorge j. more
+c
+c **********
+ integer i,j,jj,jp1,k,l
+ double precision alpha,bnorm,epsmch,gnorm,one,qnorm,sgnorm,sum,
+ * temp,zero
+ double precision dpmpar,enorm
+ data one,zero /1.0d0,0.0d0/
+c
+c epsmch is the machine precision.
+c
+ epsmch = dpmpar(1)
+c
+c first, calculate the gauss-newton direction.
+c
+ jj = (n*(n + 1))/2 + 1
+ do 50 k = 1, n
+ j = n - k + 1
+ jp1 = j + 1
+ jj = jj - k
+ l = jj + 1
+ sum = zero
+ if (n .lt. jp1) go to 20
+ do 10 i = jp1, n
+ sum = sum + r(l)*x(i)
+ l = l + 1
+ 10 continue
+ 20 continue
+ temp = r(jj)
+ if (temp .ne. zero) go to 40
+ l = j
+ do 30 i = 1, j
+ temp = dmax1(temp,dabs(r(l)))
+ l = l + n - i
+ 30 continue
+ temp = epsmch*temp
+ if (temp .eq. zero) temp = epsmch
+ 40 continue
+ x(j) = (qtb(j) - sum)/temp
+ 50 continue
+c
+c test whether the gauss-newton direction is acceptable.
+c
+ do 60 j = 1, n
+ wa1(j) = zero
+ wa2(j) = diag(j)*x(j)
+ 60 continue
+ qnorm = enorm(n,wa2)
+ if (qnorm .le. delta) go to 140
+c
+c the gauss-newton direction is not acceptable.
+c next, calculate the scaled gradient direction.
+c
+ l = 1
+ do 80 j = 1, n
+ temp = qtb(j)
+ do 70 i = j, n
+ wa1(i) = wa1(i) + r(l)*temp
+ l = l + 1
+ 70 continue
+ wa1(j) = wa1(j)/diag(j)
+ 80 continue
+c
+c calculate the norm of the scaled gradient and test for
+c the special case in which the scaled gradient is zero.
+c
+ gnorm = enorm(n,wa1)
+ sgnorm = zero
+ alpha = delta/qnorm
+ if (gnorm .eq. zero) go to 120
+c
+c calculate the point along the scaled gradient
+c at which the quadratic is minimized.
+c
+ do 90 j = 1, n
+ wa1(j) = (wa1(j)/gnorm)/diag(j)
+ 90 continue
+ l = 1
+ do 110 j = 1, n
+ sum = zero
+ do 100 i = j, n
+ sum = sum + r(l)*wa1(i)
+ l = l + 1
+ 100 continue
+ wa2(j) = sum
+ 110 continue
+ temp = enorm(n,wa2)
+ sgnorm = (gnorm/temp)/temp
+c
+c test whether the scaled gradient direction is acceptable.
+c
+ alpha = zero
+ if (sgnorm .ge. delta) go to 120
+c
+c the scaled gradient direction is not acceptable.
+c finally, calculate the point along the dogleg
+c at which the quadratic is minimized.
+c
+ bnorm = enorm(n,qtb)
+ temp = (bnorm/gnorm)*(bnorm/qnorm)*(sgnorm/delta)
+ temp = temp - (delta/qnorm)*(sgnorm/delta)**2
+ * + dsqrt((temp-(delta/qnorm))**2
+ * +(one-(delta/qnorm)**2)*(one-(sgnorm/delta)**2))
+ alpha = ((delta/qnorm)*(one - (sgnorm/delta)**2))/temp
+ 120 continue
+c
+c form appropriate convex combination of the gauss-newton
+c direction and the scaled gradient direction.
+c
+ temp = (one - alpha)*dmin1(sgnorm,delta)
+ do 130 j = 1, n
+ x(j) = temp*wa1(j) + alpha*x(j)
+ 130 continue
+ 140 continue
+ return
+c
+c last card of subroutine dogleg.
+c
+ end
diff --git a/src/dpmpar.f b/src/dpmpar.f
new file mode 100644
index 0000000..cb6545a
--- /dev/null
+++ b/src/dpmpar.f
@@ -0,0 +1,177 @@
+ double precision function dpmpar(i)
+ integer i
+c **********
+c
+c Function dpmpar
+c
+c This function provides double precision machine parameters
+c when the appropriate set of data statements is activated (by
+c removing the c from column 1) and all other data statements are
+c rendered inactive. Most of the parameter values were obtained
+c from the corresponding Bell Laboratories Port Library function.
+c
+c The function statement is
+c
+c double precision function dpmpar(i)
+c
+c where
+c
+c i is an integer input variable set to 1, 2, or 3 which
+c selects the desired machine parameter. If the machine has
+c t base b digits and its smallest and largest exponents are
+c emin and emax, respectively, then these parameters are
+c
+c dpmpar(1) = b**(1 - t), the machine precision,
+c
+c dpmpar(2) = b**(emin - 1), the smallest magnitude,
+c
+c dpmpar(3) = b**emax*(1 - b**(-t)), the largest magnitude.
+c
+c Argonne National Laboratory. MINPACK Project. November 1996.
+c Burton S. Garbow, Kenneth E. Hillstrom, Jorge J. More'
+c
+c **********
+ integer mcheps(4)
+ integer minmag(4)
+ integer maxmag(4)
+ double precision dmach(3)
+ equivalence (dmach(1),mcheps(1))
+ equivalence (dmach(2),minmag(1))
+ equivalence (dmach(3),maxmag(1))
+c
+c Machine constants for the IBM 360/370 series,
+c the Amdahl 470/V6, the ICL 2900, the Itel AS/6,
+c the Xerox Sigma 5/7/9 and the Sel systems 85/86.
+c
+c data mcheps(1),mcheps(2) / z34100000, z00000000 /
+c data minmag(1),minmag(2) / z00100000, z00000000 /
+c data maxmag(1),maxmag(2) / z7fffffff, zffffffff /
+c
+c Machine constants for the Honeywell 600/6000 series.
+c
+c data mcheps(1),mcheps(2) / o606400000000, o000000000000 /
+c data minmag(1),minmag(2) / o402400000000, o000000000000 /
+c data maxmag(1),maxmag(2) / o376777777777, o777777777777 /
+c
+c Machine constants for the CDC 6000/7000 series.
+c
+c data mcheps(1) / 15614000000000000000b /
+c data mcheps(2) / 15010000000000000000b /
+c
+c data minmag(1) / 00604000000000000000b /
+c data minmag(2) / 00000000000000000000b /
+c
+c data maxmag(1) / 37767777777777777777b /
+c data maxmag(2) / 37167777777777777777b /
+c
+c Machine constants for the PDP-10 (KA processor).
+c
+c data mcheps(1),mcheps(2) / "114400000000, "000000000000 /
+c data minmag(1),minmag(2) / "033400000000, "000000000000 /
+c data maxmag(1),maxmag(2) / "377777777777, "344777777777 /
+c
+c Machine constants for the PDP-10 (KI processor).
+c
+c data mcheps(1),mcheps(2) / "104400000000, "000000000000 /
+c data minmag(1),minmag(2) / "000400000000, "000000000000 /
+c data maxmag(1),maxmag(2) / "377777777777, "377777777777 /
+c
+c Machine constants for the PDP-11.
+c
+c data mcheps(1),mcheps(2) / 9472, 0 /
+c data mcheps(3),mcheps(4) / 0, 0 /
+c
+c data minmag(1),minmag(2) / 128, 0 /
+c data minmag(3),minmag(4) / 0, 0 /
+c
+c data maxmag(1),maxmag(2) / 32767, -1 /
+c data maxmag(3),maxmag(4) / -1, -1 /
+c
+c Machine constants for the Burroughs 6700/7700 systems.
+c
+c data mcheps(1) / o1451000000000000 /
+c data mcheps(2) / o0000000000000000 /
+c
+c data minmag(1) / o1771000000000000 /
+c data minmag(2) / o7770000000000000 /
+c
+c data maxmag(1) / o0777777777777777 /
+c data maxmag(2) / o7777777777777777 /
+c
+c Machine constants for the Burroughs 5700 system.
+c
+c data mcheps(1) / o1451000000000000 /
+c data mcheps(2) / o0000000000000000 /
+c
+c data minmag(1) / o1771000000000000 /
+c data minmag(2) / o0000000000000000 /
+c
+c data maxmag(1) / o0777777777777777 /
+c data maxmag(2) / o0007777777777777 /
+c
+c Machine constants for the Burroughs 1700 system.
+c
+c data mcheps(1) / zcc6800000 /
+c data mcheps(2) / z000000000 /
+c
+c data minmag(1) / zc00800000 /
+c data minmag(2) / z000000000 /
+c
+c data maxmag(1) / zdffffffff /
+c data maxmag(2) / zfffffffff /
+c
+c Machine constants for the Univac 1100 series.
+c
+c data mcheps(1),mcheps(2) / o170640000000, o000000000000 /
+c data minmag(1),minmag(2) / o000040000000, o000000000000 /
+c data maxmag(1),maxmag(2) / o377777777777, o777777777777 /
+c
+c Machine constants for the Data General Eclipse S/200.
+c
+c Note - it may be appropriate to include the following card -
+c static dmach(3)
+c
+c data minmag/20k,3*0/,maxmag/77777k,3*177777k/
+c data mcheps/32020k,3*0/
+c
+c Machine constants for the Harris 220.
+c
+c data mcheps(1),mcheps(2) / '20000000, '00000334 /
+c data minmag(1),minmag(2) / '20000000, '00000201 /
+c data maxmag(1),maxmag(2) / '37777777, '37777577 /
+c
+c Machine constants for the Cray-1.
+c
+c data mcheps(1) / 0376424000000000000000b /
+c data mcheps(2) / 0000000000000000000000b /
+c
+c data minmag(1) / 0200034000000000000000b /
+c data minmag(2) / 0000000000000000000000b /
+c
+c data maxmag(1) / 0577777777777777777777b /
+c data maxmag(2) / 0000007777777777777776b /
+c
+c Machine constants for the Prime 400.
+c
+c data mcheps(1),mcheps(2) / :10000000000, :00000000123 /
+c data minmag(1),minmag(2) / :10000000000, :00000100000 /
+c data maxmag(1),maxmag(2) / :17777777777, :37777677776 /
+c
+c Machine constants for the VAX-11.
+c
+c data mcheps(1),mcheps(2) / 9472, 0 /
+c data minmag(1),minmag(2) / 128, 0 /
+c data maxmag(1),maxmag(2) / -32769, -1 /
+c
+c Machine constants for IEEE machines.
+c
+ data dmach(1) /2.22044604926d-16/
+ data dmach(2) /2.22507385852d-308/
+ data dmach(3) /1.79769313485d+308/
+c
+ dpmpar = dmach(i)
+ return
+c
+c Last card of function dpmpar.
+c
+ end
diff --git a/src/enorm.f b/src/enorm.f
new file mode 100644
index 0000000..2cb5b60
--- /dev/null
+++ b/src/enorm.f
@@ -0,0 +1,108 @@
+ double precision function enorm(n,x)
+ integer n
+ double precision x(n)
+c **********
+c
+c function enorm
+c
+c given an n-vector x, this function calculates the
+c euclidean norm of x.
+c
+c the euclidean norm is computed by accumulating the sum of
+c squares in three different sums. the sums of squares for the
+c small and large components are scaled so that no overflows
+c occur. non-destructive underflows are permitted. underflows
+c and overflows do not occur in the computation of the unscaled
+c sum of squares for the intermediate components.
+c the definitions of small, intermediate and large components
+c depend on two constants, rdwarf and rgiant. the main
+c restrictions on these constants are that rdwarf**2 not
+c underflow and rgiant**2 not overflow. the constants
+c given here are suitable for every known computer.
+c
+c the function statement is
+c
+c double precision function enorm(n,x)
+c
+c where
+c
+c n is a positive integer input variable.
+c
+c x is an input array of length n.
+c
+c subprograms called
+c
+c fortran-supplied ... dabs,dsqrt
+c
+c argonne national laboratory. minpack project. march 1980.
+c burton s. garbow, kenneth e. hillstrom, jorge j. more
+c
+c **********
+ integer i
+ double precision agiant,floatn,one,rdwarf,rgiant,s1,s2,s3,xabs,
+ * x1max,x3max,zero
+ data one,zero,rdwarf,rgiant /1.0d0,0.0d0,3.834d-20,1.304d19/
+ s1 = zero
+ s2 = zero
+ s3 = zero
+ x1max = zero
+ x3max = zero
+ floatn = n
+ agiant = rgiant/floatn
+ do 90 i = 1, n
+ xabs = dabs(x(i))
+ if (xabs .gt. rdwarf .and. xabs .lt. agiant) go to 70
+ if (xabs .le. rdwarf) go to 30
+c
+c sum for large components.
+c
+ if (xabs .le. x1max) go to 10
+ s1 = one + s1*(x1max/xabs)**2
+ x1max = xabs
+ go to 20
+ 10 continue
+ s1 = s1 + (xabs/x1max)**2
+ 20 continue
+ go to 60
+ 30 continue
+c
+c sum for small components.
+c
+ if (xabs .le. x3max) go to 40
+ s3 = one + s3*(x3max/xabs)**2
+ x3max = xabs
+ go to 50
+ 40 continue
+ if (xabs .ne. zero) s3 = s3 + (xabs/x3max)**2
+ 50 continue
+ 60 continue
+ go to 80
+ 70 continue
+c
+c sum for intermediate components.
+c
+ s2 = s2 + xabs**2
+ 80 continue
+ 90 continue
+c
+c calculation of norm.
+c
+ if (s1 .eq. zero) go to 100
+ enorm = x1max*dsqrt(s1+(s2/x1max)/x1max)
+ go to 130
+ 100 continue
+ if (s2 .eq. zero) go to 110
+ if (s2 .ge. x3max)
+ * enorm = dsqrt(s2*(one+(x3max/s2)*(x3max*s3)))
+ if (s2 .lt. x3max)
+ * enorm = dsqrt(x3max*((s2/x3max)+(x3max*s3)))
+ go to 120
+ 110 continue
+ enorm = x3max*dsqrt(s3)
+ 120 continue
+ 130 continue
+ return
+c
+c last card of function enorm.
+c
+ end
diff --git a/src/fcn_lmder.c b/src/fcn_lmder.c
new file mode 100644
index 0000000..fce5f8f
--- /dev/null
+++ b/src/fcn_lmder.c
@@ -0,0 +1,67 @@
+#include <R.h>
+#include <Rdefines.h>
+#include "minpack_lm.h"
+
+
+void fcn_lmder(int *m, int *n, double *par, double *fvec, double *fjac,
+ int *ldfjac, int *iflag)
+{
+ int i, j;
+ double sumf;
+ SEXP sexp_fvec, sexp_fjac;
+
+ /* Rprintf("fcn-lmder calling...\n"); */
+ if (IS_NUMERIC(OS->par))
+ for (i = 0; i < *n; i++) {
+ if(par[i] < NUMERIC_POINTER(OS->lower)[i])
+ par[i] = NUMERIC_POINTER(OS->lower)[i];
+ if(par[i] > NUMERIC_POINTER(OS->upper)[i])
+ par[i] = NUMERIC_POINTER(OS->upper)[i];
+ NUMERIC_POINTER(OS->par)[i] = par[i];
+ }
+ else
+ for (i = 0; i < *n; i++) {
+ if(par[i] < NUMERIC_POINTER(OS->lower)[i])
+ par[i] = NUMERIC_POINTER(OS->lower)[i];
+ if(par[i] > NUMERIC_POINTER(OS->upper)[i])
+ par[i] = NUMERIC_POINTER(OS->upper)[i];
+ NUMERIC_POINTER(VECTOR_ELT(OS->par, i))[0] = par[i];
+ }
+
+ if(*iflag == 0) {
+ if(OS->nprint > 0) {
+ Rprintf("It. %4d, RSS = %10g, Par. =", OS->niter,
+ OS->rsstrace[OS->niter]);
+ for (i = 0; i < *n; i++)
+ Rprintf(" % 10g", par[i]);
+ Rprintf("\n");
+
+ }
+ OS->niter++;
+ }
+ else if (*iflag == 1) {
+ SETCADR(OS->fcall, OS->par);
+ PROTECT(sexp_fvec = eval(OS->fcall, OS->env));
+ UNPROTECT(1);
+ sumf = 0;
+ for (i = 0; i < *m; i++) {
+ fvec[i] = NUMERIC_POINTER(sexp_fvec)[i];
+ sumf += fvec[i]*fvec[i];
+ }
+
+ OS->rsstrace[OS->niter] = sumf;
+ }
+ else if (*iflag == 2) {
+ SETCADR(OS->jcall, OS->par);
+ PROTECT(sexp_fjac = eval(OS->jcall, OS->env));
+
+ for (j = 0; j < *n; j++)
+ for (i = 0; i < *m; i++)
+ fjac[(*ldfjac)*j + i] = NUMERIC_POINTER(sexp_fjac)[(*m)*j + i];
+
+ UNPROTECT(1);
+
+ }
+ if(OS->niter == OS->maxiter)
+ *iflag = -1;
+}
diff --git a/src/fcn_lmdif.c b/src/fcn_lmdif.c
new file mode 100644
index 0000000..0574d22
--- /dev/null
+++ b/src/fcn_lmdif.c
@@ -0,0 +1,57 @@
+#include <R.h>
+#include <Rdefines.h>
+#include "minpack_lm.h"
+
+
+void fcn_lmdif(int *m, int *n, double *par, double *fvec, int *iflag)
+{
+ int i;
+ double sumf;
+ SEXP sexp_fvec;
+
+ /* Rprintf("fcn-lmdif calling...\n"); */
+ if (IS_NUMERIC(OS->par))
+ for (i = 0; i < *n; i++) {
+ if(par[i] < NUMERIC_POINTER(OS->lower)[i])
+ par[i] = NUMERIC_POINTER(OS->lower)[i];
+ if(par[i] > NUMERIC_POINTER(OS->upper)[i])
+ par[i] = NUMERIC_POINTER(OS->upper)[i];
+ NUMERIC_POINTER(OS->par)[i] = par[i];
+ }
+ else
+ for (i = 0; i < *n; i++) {
+ if(par[i] < NUMERIC_POINTER(OS->lower)[i])
+ par[i] = NUMERIC_POINTER(OS->lower)[i];
+ if(par[i] > NUMERIC_POINTER(OS->upper)[i])
+ par[i] = NUMERIC_POINTER(OS->upper)[i];
+ NUMERIC_POINTER(VECTOR_ELT(OS->par, i))[0] = par[i];
+ }
+ if(*iflag == 0){
+ if(OS->nprint > 0) {
+ Rprintf("It. %4d, RSS = %10g, Par. =", OS->niter,
+ OS->rsstrace[OS->niter]);
+ for (i = 0; i < *n; i++)
+ Rprintf(" % 10g", par[i]);
+ Rprintf("\n");
+
+ }
+ OS->niter++;
+ }
+ else if (*iflag == 1 || *iflag == 2) {
+ SETCADR(OS->fcall, OS->par);
+ PROTECT(sexp_fvec = eval(OS->fcall, OS->env));
+
+ for (i = 0; i < *m; i++)
+ fvec[i] = NUMERIC_POINTER(sexp_fvec)[i];
+ UNPROTECT(1);
+ sumf = 0;
+ for (i = 0; i < *m; i++)
+ sumf += (fvec[i]*fvec[i]);
+
+ OS->rsstrace[OS->niter] = sumf;
+
+ }
+ if(OS->niter == OS->maxiter)
+ *iflag = -1;
+
+}
diff --git a/src/fcn_message.c b/src/fcn_message.c
new file mode 100644
index 0000000..8c00400
--- /dev/null
+++ b/src/fcn_message.c
@@ -0,0 +1,26 @@
+#include <R.h>
+
+char *fcn_message(char *msg, int info, int n, int nit)
+{
+ if (info == 1)
+ sprintf(msg, "Relative error in the sum of squares is at most `ftol'.");
+ else if (info == 2)
+ sprintf(msg, "Relative error between `par' and the solution is at most `ptol'.");
+ else if (info == 3)
+ sprintf(msg, "Conditions for `info = 1' and `info = 2' both hold.");
+ else if (info == 4)
+ sprintf(msg, "The cosine of the angle between `fvec' and any column of the Jacobian is at most `gtol' in absolute value.");
+ else if (info == 5)
+ sprintf(msg, "Number of calls to `fcn' has reached or exceeded `maxfev' == %d.", n);
+ else if (info == 6)
+ sprintf(msg, "`ftol' is too small. No further reduction in the sum of squares is possible.");
+ else if (info == 7)
+ sprintf(msg, "`ptol' is too small. No further improvement in the approximate solution `par' is possible.");
+ else if (info == 8)
+ sprintf(msg, "`gtol' is too small. `fvec' is orthogonal to the columns of the Jacobian to machine precision.");
+ else if (info < 0)
+ sprintf(msg, "Number of iterations has reached `maxiter' == %d.", nit);
+ else if (info == 0)
+ sprintf(msg, "Improper input parameters.");
+ return msg;
+}
diff --git a/src/fdjac2.f b/src/fdjac2.f
new file mode 100644
index 0000000..218ab94
--- /dev/null
+++ b/src/fdjac2.f
@@ -0,0 +1,107 @@
+ subroutine fdjac2(fcn,m,n,x,fvec,fjac,ldfjac,iflag,epsfcn,wa)
+ integer m,n,ldfjac,iflag
+ double precision epsfcn
+ double precision x(n),fvec(m),fjac(ldfjac,n),wa(m)
+c **********
+c
+c subroutine fdjac2
+c
+c this subroutine computes a forward-difference approximation
+c to the m by n jacobian matrix associated with a specified
+c problem of m functions in n variables.
+c
+c the subroutine statement is
+c
+c subroutine fdjac2(fcn,m,n,x,fvec,fjac,ldfjac,iflag,epsfcn,wa)
+c
+c where
+c
+c fcn is the name of the user-supplied subroutine which
+c calculates the functions. fcn must be declared
+c in an external statement in the user calling
+c program, and should be written as follows.
+c
+c subroutine fcn(m,n,x,fvec,iflag)
+c integer m,n,iflag
+c double precision x(n),fvec(m)
+c ----------
+c calculate the functions at x and
+c return this vector in fvec.
+c ----------
+c return
+c end
+c
+c the value of iflag should not be changed by fcn unless
+c the user wants to terminate execution of fdjac2.
+c in this case set iflag to a negative integer.
+c
+c m is a positive integer input variable set to the number
+c of functions.
+c
+c n is a positive integer input variable set to the number
+c of variables. n must not exceed m.
+c
+c x is an input array of length n.
+c
+c fvec is an input array of length m which must contain the
+c functions evaluated at x.
+c
+c fjac is an output m by n array which contains the
+c approximation to the jacobian matrix evaluated at x.
+c
+c ldfjac is a positive integer input variable not less than m
+c which specifies the leading dimension of the array fjac.
+c
+c iflag is an integer variable which can be used to terminate
+c the execution of fdjac2. see description of fcn.
+c
+c epsfcn is an input variable used in determining a suitable
+c step length for the forward-difference approximation. this
+c approximation assumes that the relative errors in the
+c functions are of the order of epsfcn. if epsfcn is less
+c than the machine precision, it is assumed that the relative
+c errors in the functions are of the order of the machine
+c precision.
+c
+c wa is a work array of length m.
+c
+c subprograms called
+c
+c user-supplied ...... fcn
+c
+c minpack-supplied ... dpmpar
+c
+c fortran-supplied ... dabs,dmax1,dsqrt
+c
+c argonne national laboratory. minpack project. march 1980.
+c burton s. garbow, kenneth e. hillstrom, jorge j. more
+c
+c **********
+ integer i,j
+ double precision eps,epsmch,h,temp,zero
+ double precision dpmpar
+ data zero /0.0d0/
+c
+c epsmch is the machine precision.
+c
+ epsmch = dpmpar(1)
+c
+ eps = dsqrt(dmax1(epsfcn,epsmch))
+ do 20 j = 1, n
+ temp = x(j)
+ h = eps*dabs(temp)
+ if (h .eq. zero) h = eps
+ x(j) = temp + h
+ call fcn(m,n,x,wa,iflag)
+ if (iflag .lt. 0) go to 30
+ x(j) = temp
+ do 10 i = 1, m
+ fjac(i,j) = (wa(i) - fvec(i))/h
+ 10 continue
+ 20 continue
+ 30 continue
+ return
+c
+c last card of subroutine fdjac2.
+c
+ end
diff --git a/src/get_element.c b/src/get_element.c
new file mode 100644
index 0000000..8e3b1ea
--- /dev/null
+++ b/src/get_element.c
@@ -0,0 +1,15 @@
+#include <R.h>
+#include <Rinternals.h>
+
+SEXP getListElement(SEXP list, char *str)
+{
+ SEXP elmt = R_NilValue, names = getAttrib(list, R_NamesSymbol);
+ int i;
+
+ for (i = 0; i < length(list); i++)
+ if (strcmp(CHAR(STRING_ELT(names, i)), str) == 0) {
+ elmt = VECTOR_ELT(list, i);
+ break;
+ }
+ return elmt;
+}
diff --git a/src/lmder.f b/src/lmder.f
new file mode 100644
index 0000000..8797d8b
--- /dev/null
+++ b/src/lmder.f
@@ -0,0 +1,452 @@
+ subroutine lmder(fcn,m,n,x,fvec,fjac,ldfjac,ftol,xtol,gtol,
+ * maxfev,diag,mode,factor,nprint,info,nfev,njev,
+ * ipvt,qtf,wa1,wa2,wa3,wa4)
+ integer m,n,ldfjac,maxfev,mode,nprint,info,nfev,njev
+ integer ipvt(n)
+ double precision ftol,xtol,gtol,factor
+ double precision x(n),fvec(m),fjac(ldfjac,n),diag(n),qtf(n),
+ * wa1(n),wa2(n),wa3(n),wa4(m)
+c **********
+c
+c subroutine lmder
+c
+c the purpose of lmder is to minimize the sum of the squares of
+c m nonlinear functions in n variables by a modification of
+c the levenberg-marquardt algorithm. the user must provide a
+c subroutine which calculates the functions and the jacobian.
+c
+c the subroutine statement is
+c
+c subroutine lmder(fcn,m,n,x,fvec,fjac,ldfjac,ftol,xtol,gtol,
+c maxfev,diag,mode,factor,nprint,info,nfev,
+c njev,ipvt,qtf,wa1,wa2,wa3,wa4)
+c
+c where
+c
+c fcn is the name of the user-supplied subroutine which
+c calculates the functions and the jacobian. fcn must
+c be declared in an external statement in the user
+c calling program, and should be written as follows.
+c
+c subroutine fcn(m,n,x,fvec,fjac,ldfjac,iflag)
+c integer m,n,ldfjac,iflag
+c double precision x(n),fvec(m),fjac(ldfjac,n)
+c ----------
+c if iflag = 1 calculate the functions at x and
+c return this vector in fvec. do not alter fjac.
+c if iflag = 2 calculate the jacobian at x and
+c return this matrix in fjac. do not alter fvec.
+c ----------
+c return
+c end
+c
+c the value of iflag should not be changed by fcn unless
+c the user wants to terminate execution of lmder.
+c in this case set iflag to a negative integer.
+c
+c m is a positive integer input variable set to the number
+c of functions.
+c
+c n is a positive integer input variable set to the number
+c of variables. n must not exceed m.
+c
+c x is an array of length n. on input x must contain
+c an initial estimate of the solution vector. on output x
+c contains the final estimate of the solution vector.
+c
+c fvec is an output array of length m which contains
+c the functions evaluated at the output x.
+c
+c fjac is an output m by n array. the upper n by n submatrix
+c of fjac contains an upper triangular matrix r with
+c diagonal elements of nonincreasing magnitude such that
+c
+c t t t
+c p *(jac *jac)*p = r *r,
+c
+c where p is a permutation matrix and jac is the final
+c calculated jacobian. column j of p is column ipvt(j)
+c (see below) of the identity matrix. the lower trapezoidal
+c part of fjac contains information generated during
+c the computation of r.
+c
+c ldfjac is a positive integer input variable not less than m
+c which specifies the leading dimension of the array fjac.
+c
+c ftol is a nonnegative input variable. termination
+c occurs when both the actual and predicted relative
+c reductions in the sum of squares are at most ftol.
+c therefore, ftol measures the relative error desired
+c in the sum of squares.
+c
+c xtol is a nonnegative input variable. termination
+c occurs when the relative error between two consecutive
+c iterates is at most xtol. therefore, xtol measures the
+c relative error desired in the approximate solution.
+c
+c gtol is a nonnegative input variable. termination
+c occurs when the cosine of the angle between fvec and
+c any column of the jacobian is at most gtol in absolute
+c value. therefore, gtol measures the orthogonality
+c desired between the function vector and the columns
+c of the jacobian.
+c
+c maxfev is a positive integer input variable. termination
+c occurs when the number of calls to fcn with iflag = 1
+c has reached maxfev.
+c
+c diag is an array of length n. if mode = 1 (see
+c below), diag is internally set. if mode = 2, diag
+c must contain positive entries that serve as
+c multiplicative scale factors for the variables.
+c
+c mode is an integer input variable. if mode = 1, the
+c variables will be scaled internally. if mode = 2,
+c the scaling is specified by the input diag. other
+c values of mode are equivalent to mode = 1.
+c
+c factor is a positive input variable used in determining the
+c initial step bound. this bound is set to the product of
+c factor and the euclidean norm of diag*x if nonzero, or else
+c to factor itself. in most cases factor should lie in the
+c interval (.1,100.).100. is a generally recommended value.
+c
+c nprint is an integer input variable that enables controlled
+c printing of iterates if it is positive. in this case,
+c fcn is called with iflag = 0 at the beginning of the first
+c iteration and every nprint iterations thereafter and
+c immediately prior to return, with x, fvec, and fjac
+c available for printing. fvec and fjac should not be
+c altered. if nprint is not positive, no special calls
+c of fcn with iflag = 0 are made.
+c
+c info is an integer output variable. if the user has
+c terminated execution, info is set to the (negative)
+c value of iflag. see description of fcn. otherwise,
+c info is set as follows.
+c
+c info = 0 improper input parameters.
+c
+c info = 1 both actual and predicted relative reductions
+c in the sum of squares are at most ftol.
+c
+c info = 2 relative error between two consecutive iterates
+c is at most xtol.
+c
+c info = 3 conditions for info = 1 and info = 2 both hold.
+c
+c info = 4 the cosine of the angle between fvec and any
+c column of the jacobian is at most gtol in
+c absolute value.
+c
+c info = 5 number of calls to fcn with iflag = 1 has
+c reached maxfev.
+c
+c info = 6 ftol is too small. no further reduction in
+c the sum of squares is possible.
+c
+c info = 7 xtol is too small. no further improvement in
+c the approximate solution x is possible.
+c
+c info = 8 gtol is too small. fvec is orthogonal to the
+c columns of the jacobian to machine precision.
+c
+c nfev is an integer output variable set to the number of
+c calls to fcn with iflag = 1.
+c
+c njev is an integer output variable set to the number of
+c calls to fcn with iflag = 2.
+c
+c ipvt is an integer output array of length n. ipvt
+c defines a permutation matrix p such that jac*p = q*r,
+c where jac is the final calculated jacobian, q is
+c orthogonal (not stored), and r is upper triangular
+c with diagonal elements of nonincreasing magnitude.
+c column j of p is column ipvt(j) of the identity matrix.
+c
+c qtf is an output array of length n which contains
+c the first n elements of the vector (q transpose)*fvec.
+c
+c wa1, wa2, and wa3 are work arrays of length n.
+c
+c wa4 is a work array of length m.
+c
+c subprograms called
+c
+c user-supplied ...... fcn
+c
+c minpack-supplied ... dpmpar,enorm,lmpar,qrfac
+c
+c fortran-supplied ... dabs,dmax1,dmin1,dsqrt,mod
+c
+c argonne national laboratory. minpack project. march 1980.
+c burton s. garbow, kenneth e. hillstrom, jorge j. more
+c
+c **********
+ integer i,iflag,iter,j,l
+ double precision actred,delta,dirder,epsmch,fnorm,fnorm1,gnorm,
+ * one,par,pnorm,prered,p1,p5,p25,p75,p0001,ratio,
+ * sum,temp,temp1,temp2,xnorm,zero
+ double precision dpmpar,enorm
+ data one,p1,p5,p25,p75,p0001,zero
+ * /1.0d0,1.0d-1,5.0d-1,2.5d-1,7.5d-1,1.0d-4,0.0d0/
+c
+c epsmch is the machine precision.
+c
+ epsmch = dpmpar(1)
+c
+ info = 0
+ iflag = 0
+ nfev = 0
+ njev = 0
+c
+c check the input parameters for errors.
+c
+ if (n .le. 0 .or. m .lt. n .or. ldfjac .lt. m
+ * .or. ftol .lt. zero .or. xtol .lt. zero .or. gtol .lt. zero
+ * .or. maxfev .le. 0 .or. factor .le. zero) go to 300
+ if (mode .ne. 2) go to 20
+ do 10 j = 1, n
+ if (diag(j) .le. zero) go to 300
+ 10 continue
+ 20 continue
+c
+c evaluate the function at the starting point
+c and calculate its norm.
+c
+ iflag = 1
+ call fcn(m,n,x,fvec,fjac,ldfjac,iflag)
+ nfev = 1
+ if (iflag .lt. 0) go to 300
+ fnorm = enorm(m,fvec)
+c
+c initialize levenberg-marquardt parameter and iteration counter.
+c
+ par = zero
+ iter = 1
+c
+c beginning of the outer loop.
+c
+ 30 continue
+c
+c calculate the jacobian matrix.
+c
+ iflag = 2
+ call fcn(m,n,x,fvec,fjac,ldfjac,iflag)
+ njev = njev + 1
+ if (iflag .lt. 0) go to 300
+c
+c if requested, call fcn to enable printing of iterates.
+c
+ if (nprint .le. 0) go to 40
+ iflag = 0
+ if (mod(iter-1,nprint) .eq. 0)
+ * call fcn(m,n,x,fvec,fjac,ldfjac,iflag)
+ if (iflag .lt. 0) go to 300
+ 40 continue
+c
+c compute the qr factorization of the jacobian.
+c
+ call qrfac(m,n,fjac,ldfjac,.true.,ipvt,n,wa1,wa2,wa3)
+c
+c on the first iteration and if mode is 1, scale according
+c to the norms of the columns of the initial jacobian.
+c
+ if (iter .ne. 1) go to 80
+ if (mode .eq. 2) go to 60
+ do 50 j = 1, n
+ diag(j) = wa2(j)
+ if (wa2(j) .eq. zero) diag(j) = one
+ 50 continue
+ 60 continue
+c
+c on the first iteration, calculate the norm of the scaled x
+c and initialize the step bound delta.
+c
+ do 70 j = 1, n
+ wa3(j) = diag(j)*x(j)
+ 70 continue
+ xnorm = enorm(n,wa3)
+ delta = factor*xnorm
+ if (delta .eq. zero) delta = factor
+ 80 continue
+c
+c form (q transpose)*fvec and store the first n components in
+c qtf.
+c
+ do 90 i = 1, m
+ wa4(i) = fvec(i)
+ 90 continue
+ do 130 j = 1, n
+ if (fjac(j,j) .eq. zero) go to 120
+ sum = zero
+ do 100 i = j, m
+ sum = sum + fjac(i,j)*wa4(i)
+ 100 continue
+ temp = -sum/fjac(j,j)
+ do 110 i = j, m
+ wa4(i) = wa4(i) + fjac(i,j)*temp
+ 110 continue
+ 120 continue
+ fjac(j,j) = wa1(j)
+ qtf(j) = wa4(j)
+ 130 continue
+c
+c compute the norm of the scaled gradient.
+c
+ gnorm = zero
+ if (fnorm .eq. zero) go to 170
+ do 160 j = 1, n
+ l = ipvt(j)
+ if (wa2(l) .eq. zero) go to 150
+ sum = zero
+ do 140 i = 1, j
+ sum = sum + fjac(i,j)*(qtf(i)/fnorm)
+ 140 continue
+ gnorm = dmax1(gnorm,dabs(sum/wa2(l)))
+ 150 continue
+ 160 continue
+ 170 continue
+c
+c test for convergence of the gradient norm.
+c
+ if (gnorm .le. gtol) info = 4
+ if (info .ne. 0) go to 300
+c
+c rescale if necessary.
+c
+ if (mode .eq. 2) go to 190
+ do 180 j = 1, n
+ diag(j) = dmax1(diag(j),wa2(j))
+ 180 continue
+ 190 continue
+c
+c beginning of the inner loop.
+c
+ 200 continue
+c
+c determine the levenberg-marquardt parameter.
+c
+ call lmpar(n,fjac,ldfjac,ipvt,diag,qtf,delta,par,wa1,wa2,
+ * wa3,wa4)
+c
+c store the direction p and x + p. calculate the norm of p.
+c
+ do 210 j = 1, n
+ wa1(j) = -wa1(j)
+ wa2(j) = x(j) + wa1(j)
+ wa3(j) = diag(j)*wa1(j)
+ 210 continue
+ pnorm = enorm(n,wa3)
+c
+c on the first iteration, adjust the initial step bound.
+c
+ if (iter .eq. 1) delta = dmin1(delta,pnorm)
+c
+c evaluate the function at x + p and calculate its norm.
+c
+ iflag = 1
+ call fcn(m,n,wa2,wa4,fjac,ldfjac,iflag)
+ nfev = nfev + 1
+ if (iflag .lt. 0) go to 300
+ fnorm1 = enorm(m,wa4)
+c
+c compute the scaled actual reduction.
+c
+ actred = -one
+ if (p1*fnorm1 .lt. fnorm) actred = one - (fnorm1/fnorm)**2
+c
+c compute the scaled predicted reduction and
+c the scaled directional derivative.
+c
+ do 230 j = 1, n
+ wa3(j) = zero
+ l = ipvt(j)
+ temp = wa1(l)
+ do 220 i = 1, j
+ wa3(i) = wa3(i) + fjac(i,j)*temp
+ 220 continue
+ 230 continue
+ temp1 = enorm(n,wa3)/fnorm
+ temp2 = (dsqrt(par)*pnorm)/fnorm
+ prered = temp1**2 + temp2**2/p5
+ dirder = -(temp1**2 + temp2**2)
+c
+c compute the ratio of the actual to the predicted
+c reduction.
+c
+ ratio = zero
+ if (prered .ne. zero) ratio = actred/prered
+c
+c update the step bound.
+c
+ if (ratio .gt. p25) go to 240
+ if (actred .ge. zero) temp = p5
+ if (actred .lt. zero)
+ * temp = p5*dirder/(dirder + p5*actred)
+ if (p1*fnorm1 .ge. fnorm .or. temp .lt. p1) temp = p1
+ delta = temp*dmin1(delta,pnorm/p1)
+ par = par/temp
+ go to 260
+ 240 continue
+ if (par .ne. zero .and. ratio .lt. p75) go to 250
+ delta = pnorm/p5
+ par = p5*par
+ 250 continue
+ 260 continue
+c
+c test for successful iteration.
+c
+ if (ratio .lt. p0001) go to 290
+c
+c successful iteration. update x, fvec, and their norms.
+c
+ do 270 j = 1, n
+ x(j) = wa2(j)
+ wa2(j) = diag(j)*x(j)
+ 270 continue
+ do 280 i = 1, m
+ fvec(i) = wa4(i)
+ 280 continue
+ xnorm = enorm(n,wa2)
+ fnorm = fnorm1
+ iter = iter + 1
+ 290 continue
+c
+c tests for convergence.
+c
+ if (dabs(actred) .le. ftol .and. prered .le. ftol
+ * .and. p5*ratio .le. one) info = 1
+ if (delta .le. xtol*xnorm) info = 2
+ if (dabs(actred) .le. ftol .and. prered .le. ftol
+ * .and. p5*ratio .le. one .and. info .eq. 2) info = 3
+ if (info .ne. 0) go to 300
+c
+c tests for termination and stringent tolerances.
+c
+ if (nfev .ge. maxfev) info = 5
+ if (dabs(actred) .le. epsmch .and. prered .le. epsmch
+ * .and. p5*ratio .le. one) info = 6
+ if (delta .le. epsmch*xnorm) info = 7
+ if (gnorm .le. epsmch) info = 8
+ if (info .ne. 0) go to 300
+c
+c end of the inner loop. repeat if iteration unsuccessful.
+c
+ if (ratio .lt. p0001) go to 200
+c
+c end of the outer loop.
+c
+ go to 30
+ 300 continue
+c
+c termination, either normal or user imposed.
+c
+ if (iflag .lt. 0) info = iflag
+ iflag = 0
+ if (nprint .gt. 0) call fcn(m,n,x,fvec,fjac,ldfjac,iflag)
+ return
+c
+c last card of subroutine lmder.
+c
+ end
diff --git a/src/lmdif.f b/src/lmdif.f
new file mode 100644
index 0000000..dd3d4ee
--- /dev/null
+++ b/src/lmdif.f
@@ -0,0 +1,454 @@
+ subroutine lmdif(fcn,m,n,x,fvec,ftol,xtol,gtol,maxfev,epsfcn,
+ * diag,mode,factor,nprint,info,nfev,fjac,ldfjac,
+ * ipvt,qtf,wa1,wa2,wa3,wa4)
+ integer m,n,maxfev,mode,nprint,info,nfev,ldfjac
+ integer ipvt(n)
+ double precision ftol,xtol,gtol,epsfcn,factor
+ double precision x(n),fvec(m),diag(n),fjac(ldfjac,n),qtf(n),
+ * wa1(n),wa2(n),wa3(n),wa4(m)
+ external fcn
+c **********
+c
+c subroutine lmdif
+c
+c the purpose of lmdif is to minimize the sum of the squares of
+c m nonlinear functions in n variables by a modification of
+c the levenberg-marquardt algorithm. the user must provide a
+c subroutine which calculates the functions. the jacobian is
+c then calculated by a forward-difference approximation.
+c
+c the subroutine statement is
+c
+c subroutine lmdif(fcn,m,n,x,fvec,ftol,xtol,gtol,maxfev,epsfcn,
+c diag,mode,factor,nprint,info,nfev,fjac,
+c ldfjac,ipvt,qtf,wa1,wa2,wa3,wa4)
+c
+c where
+c
+c fcn is the name of the user-supplied subroutine which
+c calculates the functions. fcn must be declared
+c in an external statement in the user calling
+c program, and should be written as follows.
+c
+c subroutine fcn(m,n,x,fvec,iflag)
+c integer m,n,iflag
+c double precision x(n),fvec(m)
+c ----------
+c calculate the functions at x and
+c return this vector in fvec.
+c ----------
+c return
+c end
+c
+c the value of iflag should not be changed by fcn unless
+c the user wants to terminate execution of lmdif.
+c in this case set iflag to a negative integer.
+c
+c m is a positive integer input variable set to the number
+c of functions.
+c
+c n is a positive integer input variable set to the number
+c of variables. n must not exceed m.
+c
+c x is an array of length n. on input x must contain
+c an initial estimate of the solution vector. on output x
+c contains the final estimate of the solution vector.
+c
+c fvec is an output array of length m which contains
+c the functions evaluated at the output x.
+c
+c ftol is a nonnegative input variable. termination
+c occurs when both the actual and predicted relative
+c reductions in the sum of squares are at most ftol.
+c therefore, ftol measures the relative error desired
+c in the sum of squares.
+c
+c xtol is a nonnegative input variable. termination
+c occurs when the relative error between two consecutive
+c iterates is at most xtol. therefore, xtol measures the
+c relative error desired in the approximate solution.
+c
+c gtol is a nonnegative input variable. termination
+c occurs when the cosine of the angle between fvec and
+c any column of the jacobian is at most gtol in absolute
+c value. therefore, gtol measures the orthogonality
+c desired between the function vector and the columns
+c of the jacobian.
+c
+c maxfev is a positive integer input variable. termination
+c occurs when the number of calls to fcn is at least
+c maxfev by the end of an iteration.
+c
+c epsfcn is an input variable used in determining a suitable
+c step length for the forward-difference approximation. this
+c approximation assumes that the relative errors in the
+c functions are of the order of epsfcn. if epsfcn is less
+c than the machine precision, it is assumed that the relative
+c errors in the functions are of the order of the machine
+c precision.
+c
+c diag is an array of length n. if mode = 1 (see
+c below), diag is internally set. if mode = 2, diag
+c must contain positive entries that serve as
+c multiplicative scale factors for the variables.
+c
+c mode is an integer input variable. if mode = 1, the
+c variables will be scaled internally. if mode = 2,
+c the scaling is specified by the input diag. other
+c values of mode are equivalent to mode = 1.
+c
+c factor is a positive input variable used in determining the
+c initial step bound. this bound is set to the product of
+c factor and the euclidean norm of diag*x if nonzero, or else
+c to factor itself. in most cases factor should lie in the
+c interval (.1,100.). 100. is a generally recommended value.
+c
+c nprint is an integer input variable that enables controlled
+c printing of iterates if it is positive. in this case,
+c fcn is called with iflag = 0 at the beginning of the first
+c iteration and every nprint iterations thereafter and
+c immediately prior to return, with x and fvec available
+c for printing. if nprint is not positive, no special calls
+c of fcn with iflag = 0 are made.
+c
+c info is an integer output variable. if the user has
+c terminated execution, info is set to the (negative)
+c value of iflag. see description of fcn. otherwise,
+c info is set as follows.
+c
+c info = 0 improper input parameters.
+c
+c info = 1 both actual and predicted relative reductions
+c in the sum of squares are at most ftol.
+c
+c info = 2 relative error between two consecutive iterates
+c is at most xtol.
+c
+c info = 3 conditions for info = 1 and info = 2 both hold.
+c
+c info = 4 the cosine of the angle between fvec and any
+c column of the jacobian is at most gtol in
+c absolute value.
+c
+c info = 5 number of calls to fcn has reached or
+c exceeded maxfev.
+c
+c info = 6 ftol is too small. no further reduction in
+c the sum of squares is possible.
+c
+c info = 7 xtol is too small. no further improvement in
+c the approximate solution x is possible.
+c
+c info = 8 gtol is too small. fvec is orthogonal to the
+c columns of the jacobian to machine precision.
+c
+c nfev is an integer output variable set to the number of
+c calls to fcn.
+c
+c fjac is an output m by n array. the upper n by n submatrix
+c of fjac contains an upper triangular matrix r with
+c diagonal elements of nonincreasing magnitude such that
+c
+c t t t
+c p *(jac *jac)*p = r *r,
+c
+c where p is a permutation matrix and jac is the final
+c calculated jacobian. column j of p is column ipvt(j)
+c (see below) of the identity matrix. the lower trapezoidal
+c part of fjac contains information generated during
+c the computation of r.
+c
+c ldfjac is a positive integer input variable not less than m
+c which specifies the leading dimension of the array fjac.
+c
+c ipvt is an integer output array of length n. ipvt
+c defines a permutation matrix p such that jac*p = q*r,
+c where jac is the final calculated jacobian, q is
+c orthogonal (not stored), and r is upper triangular
+c with diagonal elements of nonincreasing magnitude.
+c column j of p is column ipvt(j) of the identity matrix.
+c
+c qtf is an output array of length n which contains
+c the first n elements of the vector (q transpose)*fvec.
+c
+c wa1, wa2, and wa3 are work arrays of length n.
+c
+c wa4 is a work array of length m.
+c
+c subprograms called
+c
+c user-supplied ...... fcn
+c
+c minpack-supplied ... dpmpar,enorm,fdjac2,lmpar,qrfac
+c
+c fortran-supplied ... dabs,dmax1,dmin1,dsqrt,mod
+c
+c argonne national laboratory. minpack project. march 1980.
+c burton s. garbow, kenneth e. hillstrom, jorge j. more
+c
+c **********
+ integer i,iflag,iter,j,l
+ double precision actred,delta,dirder,epsmch,fnorm,fnorm1,gnorm,
+ * one,par,pnorm,prered,p1,p5,p25,p75,p0001,ratio,
+ * sum,temp,temp1,temp2,xnorm,zero
+ double precision dpmpar,enorm
+ data one,p1,p5,p25,p75,p0001,zero
+ * /1.0d0,1.0d-1,5.0d-1,2.5d-1,7.5d-1,1.0d-4,0.0d0/
+c
+c epsmch is the machine precision.
+c
+ epsmch = dpmpar(1)
+c
+ info = 0
+ iflag = 0
+ nfev = 0
+c
+c check the input parameters for errors.
+c
+ if (n .le. 0 .or. m .lt. n .or. ldfjac .lt. m
+ * .or. ftol .lt. zero .or. xtol .lt. zero .or. gtol .lt. zero
+ * .or. maxfev .le. 0 .or. factor .le. zero) go to 300
+ if (mode .ne. 2) go to 20
+ do 10 j = 1, n
+ if (diag(j) .le. zero) go to 300
+ 10 continue
+ 20 continue
+c
+c evaluate the function at the starting point
+c and calculate its norm.
+c
+ iflag = 1
+ call fcn(m,n,x,fvec,iflag)
+ nfev = 1
+ if (iflag .lt. 0) go to 300
+ fnorm = enorm(m,fvec)
+c
+c initialize levenberg-marquardt parameter and iteration counter.
+c
+ par = zero
+ iter = 1
+c
+c beginning of the outer loop.
+c
+ 30 continue
+c
+c calculate the jacobian matrix.
+c
+ iflag = 2
+ call fdjac2(fcn,m,n,x,fvec,fjac,ldfjac,iflag,epsfcn,wa4)
+ nfev = nfev + n
+ if (iflag .lt. 0) go to 300
+c
+c if requested, call fcn to enable printing of iterates.
+c
+ if (nprint .le. 0) go to 40
+ iflag = 0
+ if (mod(iter-1,nprint) .eq. 0) call fcn(m,n,x,fvec,iflag)
+ if (iflag .lt. 0) go to 300
+ 40 continue
+c
+c compute the qr factorization of the jacobian.
+c
+ call qrfac(m,n,fjac,ldfjac,.true.,ipvt,n,wa1,wa2,wa3)
+c
+c on the first iteration and if mode is 1, scale according
+c to the norms of the columns of the initial jacobian.
+c
+ if (iter .ne. 1) go to 80
+ if (mode .eq. 2) go to 60
+ do 50 j = 1, n
+ diag(j) = wa2(j)
+ if (wa2(j) .eq. zero) diag(j) = one
+ 50 continue
+ 60 continue
+c
+c on the first iteration, calculate the norm of the scaled x
+c and initialize the step bound delta.
+c
+ do 70 j = 1, n
+ wa3(j) = diag(j)*x(j)
+ 70 continue
+ xnorm = enorm(n,wa3)
+ delta = factor*xnorm
+ if (delta .eq. zero) delta = factor
+ 80 continue
+c
+c form (q transpose)*fvec and store the first n components in
+c qtf.
+c
+ do 90 i = 1, m
+ wa4(i) = fvec(i)
+ 90 continue
+ do 130 j = 1, n
+ if (fjac(j,j) .eq. zero) go to 120
+ sum = zero
+ do 100 i = j, m
+ sum = sum + fjac(i,j)*wa4(i)
+ 100 continue
+ temp = -sum/fjac(j,j)
+ do 110 i = j, m
+ wa4(i) = wa4(i) + fjac(i,j)*temp
+ 110 continue
+ 120 continue
+ fjac(j,j) = wa1(j)
+ qtf(j) = wa4(j)
+ 130 continue
+c
+c compute the norm of the scaled gradient.
+c
+ gnorm = zero
+ if (fnorm .eq. zero) go to 170
+ do 160 j = 1, n
+ l = ipvt(j)
+ if (wa2(l) .eq. zero) go to 150
+ sum = zero
+ do 140 i = 1, j
+ sum = sum + fjac(i,j)*(qtf(i)/fnorm)
+ 140 continue
+ gnorm = dmax1(gnorm,dabs(sum/wa2(l)))
+ 150 continue
+ 160 continue
+ 170 continue
+c
+c test for convergence of the gradient norm.
+c
+ if (gnorm .le. gtol) info = 4
+ if (info .ne. 0) go to 300
+c
+c rescale if necessary.
+c
+ if (mode .eq. 2) go to 190
+ do 180 j = 1, n
+ diag(j) = dmax1(diag(j),wa2(j))
+ 180 continue
+ 190 continue
+c
+c beginning of the inner loop.
+c
+ 200 continue
+c
+c determine the levenberg-marquardt parameter.
+c
+ call lmpar(n,fjac,ldfjac,ipvt,diag,qtf,delta,par,wa1,wa2,
+ * wa3,wa4)
+c
+c store the direction p and x + p. calculate the norm of p.
+c
+ do 210 j = 1, n
+ wa1(j) = -wa1(j)
+ wa2(j) = x(j) + wa1(j)
+ wa3(j) = diag(j)*wa1(j)
+ 210 continue
+ pnorm = enorm(n,wa3)
+c
+c on the first iteration, adjust the initial step bound.
+c
+ if (iter .eq. 1) delta = dmin1(delta,pnorm)
+c
+c evaluate the function at x + p and calculate its norm.
+c
+ iflag = 1
+ call fcn(m,n,wa2,wa4,iflag)
+ nfev = nfev + 1
+ if (iflag .lt. 0) go to 300
+ fnorm1 = enorm(m,wa4)
+c
+c compute the scaled actual reduction.
+c
+ actred = -one
+ if (p1*fnorm1 .lt. fnorm) actred = one - (fnorm1/fnorm)**2
+c
+c compute the scaled predicted reduction and
+c the scaled directional derivative.
+c
+ do 230 j = 1, n
+ wa3(j) = zero
+ l = ipvt(j)
+ temp = wa1(l)
+ do 220 i = 1, j
+ wa3(i) = wa3(i) + fjac(i,j)*temp
+ 220 continue
+ 230 continue
+ temp1 = enorm(n,wa3)/fnorm
+ temp2 = (dsqrt(par)*pnorm)/fnorm
+ prered = temp1**2 + temp2**2/p5
+ dirder = -(temp1**2 + temp2**2)
+c
+c compute the ratio of the actual to the predicted
+c reduction.
+c
+ ratio = zero
+ if (prered .ne. zero) ratio = actred/prered
+c
+c update the step bound.
+c
+ if (ratio .gt. p25) go to 240
+ if (actred .ge. zero) temp = p5
+ if (actred .lt. zero)
+ * temp = p5*dirder/(dirder + p5*actred)
+ if (p1*fnorm1 .ge. fnorm .or. temp .lt. p1) temp = p1
+ delta = temp*dmin1(delta,pnorm/p1)
+ par = par/temp
+ go to 260
+ 240 continue
+ if (par .ne. zero .and. ratio .lt. p75) go to 250
+ delta = pnorm/p5
+ par = p5*par
+ 250 continue
+ 260 continue
+c
+c test for successful iteration.
+c
+ if (ratio .lt. p0001) go to 290
+c
+c successful iteration. update x, fvec, and their norms.
+c
+ do 270 j = 1, n
+ x(j) = wa2(j)
+ wa2(j) = diag(j)*x(j)
+ 270 continue
+ do 280 i = 1, m
+ fvec(i) = wa4(i)
+ 280 continue
+ xnorm = enorm(n,wa2)
+ fnorm = fnorm1
+ iter = iter + 1
+ 290 continue
+c
+c tests for convergence.
+c
+ if (dabs(actred) .le. ftol .and. prered .le. ftol
+ * .and. p5*ratio .le. one) info = 1
+ if (delta .le. xtol*xnorm) info = 2
+ if (dabs(actred) .le. ftol .and. prered .le. ftol
+ * .and. p5*ratio .le. one .and. info .eq. 2) info = 3
+ if (info .ne. 0) go to 300
+c
+c tests for termination and stringent tolerances.
+c
+ if (nfev .ge. maxfev) info = 5
+ if (dabs(actred) .le. epsmch .and. prered .le. epsmch
+ * .and. p5*ratio .le. one) info = 6
+ if (delta .le. epsmch*xnorm) info = 7
+ if (gnorm .le. epsmch) info = 8
+ if (info .ne. 0) go to 300
+c
+c end of the inner loop. repeat if iteration unsuccessful.
+c
+ if (ratio .lt. p0001) go to 200
+c
+c end of the outer loop.
+c
+ go to 30
+ 300 continue
+c
+c termination, either normal or user imposed.
+c
+ if (iflag .lt. 0) info = iflag
+ iflag = 0
+ if (nprint .gt. 0) call fcn(m,n,x,fvec,iflag)
+ return
+c
+c last card of subroutine lmdif.
+c
+ end
diff --git a/src/lmpar.f b/src/lmpar.f
new file mode 100644
index 0000000..26c422a
--- /dev/null
+++ b/src/lmpar.f
@@ -0,0 +1,264 @@
+ subroutine lmpar(n,r,ldr,ipvt,diag,qtb,delta,par,x,sdiag,wa1,
+ * wa2)
+ integer n,ldr
+ integer ipvt(n)
+ double precision delta,par
+ double precision r(ldr,n),diag(n),qtb(n),x(n),sdiag(n),wa1(n),
+ * wa2(n)
+c **********
+c
+c subroutine lmpar
+c
+c given an m by n matrix a, an n by n nonsingular diagonal
+c matrix d, an m-vector b, and a positive number delta,
+c the problem is to determine a value for the parameter
+c par such that if x solves the system
+c
+c a*x = b , sqrt(par)*d*x = 0 ,
+c
+c in the least squares sense, and dxnorm is the euclidean
+c norm of d*x, then either par is zero and
+c
+c (dxnorm-delta) .le. 0.1*delta ,
+c
+c or par is positive and
+c
+c abs(dxnorm-delta) .le. 0.1*delta .
+c
+c this subroutine completes the solution of the problem
+c if it is provided with the necessary information from the
+c qr factorization, with column pivoting, of a. that is, if
+c a*p = q*r, where p is a permutation matrix, q has orthogonal
+c columns, and r is an upper triangular matrix with diagonal
+c elements of nonincreasing magnitude, then lmpar expects
+c the full upper triangle of r, the permutation matrix p,
+c and the first n components of (q transpose)*b. on output
+c lmpar also provides an upper triangular matrix s such that
+c
+c t t t
+c p *(a *a + par*d*d)*p = s *s .
+c
+c s is employed within lmpar and may be of separate interest.
+c
+c only a few iterations are generally needed for convergence
+c of the algorithm. if, however, the limit of 10 iterations
+c is reached, then the output par will contain the best
+c value obtained so far.
+c
+c the subroutine statement is
+c
+c subroutine lmpar(n,r,ldr,ipvt,diag,qtb,delta,par,x,sdiag,
+c wa1,wa2)
+c
+c where
+c
+c n is a positive integer input variable set to the order of r.
+c
+c r is an n by n array. on input the full upper triangle
+c must contain the full upper triangle of the matrix r.
+c on output the full upper triangle is unaltered, and the
+c strict lower triangle contains the strict upper triangle
+c (transposed) of the upper triangular matrix s.
+c
+c ldr is a positive integer input variable not less than n
+c which specifies the leading dimension of the array r.
+c
+c ipvt is an integer input array of length n which defines the
+c permutation matrix p such that a*p = q*r. column j of p
+c is column ipvt(j) of the identity matrix.
+c
+c diag is an input array of length n which must contain the
+c diagonal elements of the matrix d.
+c
+c qtb is an input array of length n which must contain the first
+c n elements of the vector (q transpose)*b.
+c
+c delta is a positive input variable which specifies an upper
+c bound on the euclidean norm of d*x.
+c
+c par is a nonnegative variable. on input par contains an
+c initial estimate of the levenberg-marquardt parameter.
+c on output par contains the final estimate.
+c
+c x is an output array of length n which contains the least
+c squares solution of the system a*x = b, sqrt(par)*d*x = 0,
+c for the output par.
+c
+c sdiag is an output array of length n which contains the
+c diagonal elements of the upper triangular matrix s.
+c
+c wa1 and wa2 are work arrays of length n.
+c
+c subprograms called
+c
+c minpack-supplied ... dpmpar,enorm,qrsolv
+c
+c fortran-supplied ... dabs,dmax1,dmin1,dsqrt
+c
+c argonne national laboratory. minpack project. march 1980.
+c burton s. garbow, kenneth e. hillstrom, jorge j. more
+c
+c **********
+ integer i,iter,j,jm1,jp1,k,l,nsing
+ double precision dxnorm,dwarf,fp,gnorm,parc,parl,paru,p1,p001,
+ * sum,temp,zero
+ double precision dpmpar,enorm
+ data p1,p001,zero /1.0d-1,1.0d-3,0.0d0/
+c
+c dwarf is the smallest positive magnitude.
+c
+ dwarf = dpmpar(2)
+c
+c compute and store in x the gauss-newton direction. if the
+c jacobian is rank-deficient, obtain a least squares solution.
+c
+ nsing = n
+ do 10 j = 1, n
+ wa1(j) = qtb(j)
+ if (r(j,j) .eq. zero .and. nsing .eq. n) nsing = j - 1
+ if (nsing .lt. n) wa1(j) = zero
+ 10 continue
+ if (nsing .lt. 1) go to 50
+ do 40 k = 1, nsing
+ j = nsing - k + 1
+ wa1(j) = wa1(j)/r(j,j)
+ temp = wa1(j)
+ jm1 = j - 1
+ if (jm1 .lt. 1) go to 30
+ do 20 i = 1, jm1
+ wa1(i) = wa1(i) - r(i,j)*temp
+ 20 continue
+ 30 continue
+ 40 continue
+ 50 continue
+ do 60 j = 1, n
+ l = ipvt(j)
+ x(l) = wa1(j)
+ 60 continue
+c
+c initialize the iteration counter.
+c evaluate the function at the origin, and test
+c for acceptance of the gauss-newton direction.
+c
+ iter = 0
+ do 70 j = 1, n
+ wa2(j) = diag(j)*x(j)
+ 70 continue
+ dxnorm = enorm(n,wa2)
+ fp = dxnorm - delta
+ if (fp .le. p1*delta) go to 220
+c
+c if the jacobian is not rank deficient, the newton
+c step provides a lower bound, parl, for the zero of
+c the function. otherwise set this bound to zero.
+c
+ parl = zero
+ if (nsing .lt. n) go to 120
+ do 80 j = 1, n
+ l = ipvt(j)
+ wa1(j) = diag(l)*(wa2(l)/dxnorm)
+ 80 continue
+ do 110 j = 1, n
+ sum = zero
+ jm1 = j - 1
+ if (jm1 .lt. 1) go to 100
+ do 90 i = 1, jm1
+ sum = sum + r(i,j)*wa1(i)
+ 90 continue
+ 100 continue
+ wa1(j) = (wa1(j) - sum)/r(j,j)
+ 110 continue
+ temp = enorm(n,wa1)
+ parl = ((fp/delta)/temp)/temp
+ 120 continue
+c
+c calculate an upper bound, paru, for the zero of the function.
+c
+ do 140 j = 1, n
+ sum = zero
+ do 130 i = 1, j
+ sum = sum + r(i,j)*qtb(i)
+ 130 continue
+ l = ipvt(j)
+ wa1(j) = sum/diag(l)
+ 140 continue
+ gnorm = enorm(n,wa1)
+ paru = gnorm/delta
+ if (paru .eq. zero) paru = dwarf/dmin1(delta,p1)
+c
+c if the input par lies outside of the interval (parl,paru),
+c set par to the closer endpoint.
+c
+ par = dmax1(par,parl)
+ par = dmin1(par,paru)
+ if (par .eq. zero) par = gnorm/dxnorm
+c
+c beginning of an iteration.
+c
+ 150 continue
+ iter = iter + 1
+c
+c evaluate the function at the current value of par.
+c
+ if (par .eq. zero) par = dmax1(dwarf,p001*paru)
+ temp = dsqrt(par)
+ do 160 j = 1, n
+ wa1(j) = temp*diag(j)
+ 160 continue
+ call qrsolv(n,r,ldr,ipvt,wa1,qtb,x,sdiag,wa2)
+ do 170 j = 1, n
+ wa2(j) = diag(j)*x(j)
+ 170 continue
+ dxnorm = enorm(n,wa2)
+ temp = fp
+ fp = dxnorm - delta
+c
+c if the function is small enough, accept the current value
+c of par. also test for the exceptional cases where parl
+c is zero or the number of iterations has reached 10.
+c
+ if (dabs(fp) .le. p1*delta
+ * .or. parl .eq. zero .and. fp .le. temp
+ * .and. temp .lt. zero .or. iter .eq. 10) go to 220
+c
+c compute the newton correction.
+c
+ do 180 j = 1, n
+ l = ipvt(j)
+ wa1(j) = diag(l)*(wa2(l)/dxnorm)
+ 180 continue
+ do 210 j = 1, n
+ wa1(j) = wa1(j)/sdiag(j)
+ temp = wa1(j)
+ jp1 = j + 1
+ if (n .lt. jp1) go to 200
+ do 190 i = jp1, n
+ wa1(i) = wa1(i) - r(i,j)*temp
+ 190 continue
+ 200 continue
+ 210 continue
+ temp = enorm(n,wa1)
+ parc = ((fp/delta)/temp)/temp
+c
+c depending on the sign of the function, update parl or paru.
+c
+ if (fp .gt. zero) parl = dmax1(parl,par)
+ if (fp .lt. zero) paru = dmin1(paru,par)
+c
+c compute an improved estimate for par.
+c
+ par = dmax1(parl,par+parc)
+c
+c end of an iteration.
+c
+ go to 150
+ 220 continue
+c
+c termination.
+c
+ if (iter .eq. 0) par = zero
+ return
+c
+c last card of subroutine lmpar.
+c
+ end
diff --git a/src/minpack_lm.h b/src/minpack_lm.h
new file mode 100644
index 0000000..5505942
--- /dev/null
+++ b/src/minpack_lm.h
@@ -0,0 +1,59 @@
+#ifndef MINPACK_LM_HEADER_INCLUDED
+#define MINPACK_LM_HEADER_INCLUDED 1
+
+typedef struct opt_struct {
+ SEXP par;
+ SEXP lower;
+ SEXP upper;
+ SEXP fcall;
+ SEXP jcall;
+ SEXP env;
+ double ftol;
+ double ptol;
+ double gtol;
+ double epsfcn;
+ double *diag;
+ double factor;
+ int nprint;
+ int maxiter;
+ int niter;
+ double rsstrace[1024];
+} opt_struct, *OptStruct;
+
+void fcn_lmdif(int *m, int *n, double *par, double *fvec, int *iflag);
+void fcn_lmder(int *m, int *n, double *par, double *fvec, double
+ *fjac, int *ldfjac, int *iflag);
+void F77_NAME(lmdif)(void (*fcn_lmdif)(int *m, int *n, double *par,
+ double *fvec, int *iflag),
+ int *m, int *n, double *par, double *fvec,
+ double *ftol, double *ptol, double *gtol, int *maxfev,
+ double *epsfcn, double *diag, int *mode, double *factor,
+ int *nprint, int *info, int *nfev, double *fjac,
+ int *ldfjac, int *ipvt, double *qtf,
+ double *wa1, double *wa2, double *wa3, double *wa4);
+void F77_NAME(lmder)(void (*fcn_lmder)(int *m, int *n, double *par,
+ double *fvec, double *fjac, int *ldfjac,
+ int *iflag),
+ int *m, int *n, double *par, double *fvec,
+ double *fjac, int *ldfjac,
+ double *ftol, double *ptol, double *gtol,
+ int *maxfev, double *diag, int *mode,
+ double *factor, int *nprint, int *info,
+ int *nfev, int *njev, int *ipvt,
+ double *qtf, double *wa1, double *wa2,
+ double *wa3, double *wa4);
+ char *fcn_message(char*, int, int, int);
+
+void transpose(double *x, int nrx, int ncx, double *y);
+void matprod (double *x, int nrx, int ncx, double *y, int nry, int ncy, double *z);
+void crossprod(double *x, int nrx, int ncx, double *y, int nry, int ncy, double *z);
+SEXP getListElement(SEXP list, char *str);
+double *real_vector(int n);
+int *int_vector(int n);
+
+SEXP nls_lm(SEXP par_arg, SEXP lower_arg, SEXP upper_arg,
+ SEXP fn, SEXP jac, SEXP control, SEXP rho);
+
+extern OptStruct OS;
+
+#endif /* MINPACK_LM_HEADER_INCLUDED */
diff --git a/src/nls_lm.c b/src/nls_lm.c
new file mode 100644
index 0000000..49b1dd7
--- /dev/null
+++ b/src/nls_lm.c
@@ -0,0 +1,274 @@
+#include <R.h>
+#include <Rdefines.h>
+#include "minpack_lm.h"
+
+OptStruct OS;
+
+SEXP nls_lm(SEXP par_arg, SEXP lower_arg, SEXP upper_arg, SEXP fn, SEXP jac,
+ SEXP control, SEXP rho)
+{
+ int i, j;
+ int n, m, ldfjac;
+ int info, nfev, njev;
+
+ double *par, *fvec, *fjac, *qtf, *wa1, *wa2, *wa3, *wa4,
+ *perm, *perm_t, *r, *r2, *r2_x_perm_t, *hess;
+ int *ipvt;
+
+ SEXP eval_test;
+ SEXP sexp_diag, sexp_hess, sexp_fvec, sexp_info, sexp_niter,
+ sexp_message, sexp_rsstrace;
+ SEXP out, out_names;
+
+ char lmfun_name[8], message[256];
+
+ int maxfev;
+ int mode;
+ int npr = 1;
+
+ PROTECT_INDEX ipx;
+
+ OS = (OptStruct) R_alloc(1, sizeof(opt_struct));
+
+ PROTECT(OS->par = duplicate(par_arg));
+ PROTECT(OS->lower = duplicate(lower_arg));
+ PROTECT(OS->upper = duplicate(upper_arg));
+ n = length(OS->par);
+
+ switch (TYPEOF(OS->par)) {
+ case REALSXP:
+ break;
+ case VECSXP:
+ for (i = 0; i < n; i++)
+ SET_VECTOR_ELT(OS->par, i, AS_NUMERIC(VECTOR_ELT(OS->par, i)));
+ break;
+ default:
+ error("`par' that you provided is non-list and non-numeric!");
+ }
+ switch (TYPEOF(OS->lower)) {
+ case REALSXP:
+ break;
+ case VECSXP:
+ for (i = 0; i < n; i++)
+ SET_VECTOR_ELT(OS->lower, i, AS_NUMERIC(VECTOR_ELT(OS->lower, i)));
+ break;
+ default:
+ error("`lower' that you provided is non-list and non-numeric!");
+ }
+ switch (TYPEOF(OS->upper)) {
+ case REALSXP:
+ break;
+ case VECSXP:
+ for (i = 0; i < n; i++)
+ SET_VECTOR_ELT(OS->upper, i, AS_NUMERIC(VECTOR_ELT(OS->upper, i)));
+ break;
+ default:
+ error("`upper' that you provided is non-list and non-numeric!");
+ }
+
+ if (!isFunction(fn)) error("fn is not a function!");
+ PROTECT(OS->fcall = lang2(fn, OS->par));
+
+ if (!isEnvironment(rho)) error("rho is not an environment!");
+ OS->env = rho;
+
+ PROTECT(eval_test = eval(OS->fcall, OS->env));
+ if (!IS_NUMERIC(eval_test) || length(eval_test) == 0)
+ error("evaluation of fn function returns non-sensible value!");
+ m = length(eval_test);
+ UNPROTECT(1);
+
+ ldfjac = m;
+
+ par = real_vector(n);
+ fvec = real_vector(m);
+ fjac = real_vector(ldfjac * n);
+ qtf = real_vector(n);
+ wa1 = real_vector(n);
+ wa2 = real_vector(n);
+ wa3 = real_vector(n);
+ wa4 = real_vector(m);
+ ipvt = int_vector(n);
+ perm = real_vector(n * n);
+ perm_t = real_vector(n * n);
+ r = real_vector(n * n);
+ r2 = real_vector(n * n);
+ r2_x_perm_t = real_vector(n * n);
+ hess = real_vector(n * n);
+
+ OS->ftol = NUMERIC_VALUE(getListElement(control, "ftol"));
+ OS->ptol = NUMERIC_VALUE(getListElement(control, "ptol"));
+ OS->gtol = NUMERIC_VALUE(getListElement(control, "gtol"));
+ OS->epsfcn = NUMERIC_VALUE(getListElement(control, "epsfcn"));
+ OS->factor = NUMERIC_VALUE(getListElement(control, "factor"));
+ OS->diag = real_vector(n);
+
+ PROTECT_WITH_INDEX(sexp_diag = getListElement(control, "diag"), &ipx);
+ switch (TYPEOF(sexp_diag)) {
+ case REALSXP:
+ if (length(sexp_diag) == n) {
+ REPROTECT(sexp_diag = duplicate(sexp_diag), ipx);
+ for (i = 0; i < n; i++)
+ OS->diag[i] = NUMERIC_POINTER(sexp_diag)[i];
+ mode = 2;
+ }
+ else {
+ REPROTECT(sexp_diag = NEW_NUMERIC(n), ipx);
+ mode = 1;
+ }
+ break;
+ case VECSXP:
+ if (length(sexp_diag) == n) {
+ REPROTECT(sexp_diag = duplicate(sexp_diag), ipx);
+ for (i = 0; i < n; i++) {
+ SET_VECTOR_ELT(sexp_diag, i, AS_NUMERIC(VECTOR_ELT(sexp_diag, i)));
+ OS->diag[i] = NUMERIC_VALUE(VECTOR_ELT(sexp_diag, i));
+ }
+ mode = 2;
+ }
+ else {
+ REPROTECT(sexp_diag = NEW_LIST(n), ipx);
+ for (i = 0; i < n; i++)
+ SET_VECTOR_ELT(sexp_diag, i, NEW_NUMERIC(1));
+ mode = 1;
+ }
+ break;
+ default:
+ error("`diag' that you provided is non-list and non-numeric!");
+ }
+ maxfev = INTEGER_VALUE(getListElement(control, "maxfev"));
+ OS->maxiter = INTEGER_VALUE(getListElement(control, "maxiter"));
+ if(OS->maxiter > 1024) {
+ OS->maxiter = 1024;
+ warning("resetting `maxiter' to 1024!");
+ }
+ OS->nprint = INTEGER_VALUE(getListElement(control, "nprint"));
+ if(OS->nprint > 0)
+ OS->nprint = 1;
+ if (IS_NUMERIC(OS->par)) {
+ for (i = 0; i < n; i++)
+ if(R_FINITE(NUMERIC_POINTER(OS->par)[i]))
+ par[i] = NUMERIC_POINTER(OS->par)[i];
+ else
+ error("Non-finite (or null) value for a parameter specified!");
+ }
+ else {
+ for (i = 0; i < n; i++)
+ if(R_FINITE(NUMERIC_VALUE(VECTOR_ELT(OS->par, i))))
+ par[i] = NUMERIC_VALUE(VECTOR_ELT(OS->par, i));
+ else
+ error("Non-finite (or null) value for a parameter specified!");
+ }
+
+ OS->niter = 0;
+
+/*========================================================================*/
+
+ if (isNull(jac)) {
+ F77_CALL(lmdif)(&fcn_lmdif, &m, &n, par, fvec,
+ &OS->ftol, &OS->ptol, &OS->gtol,
+ &maxfev, &OS->epsfcn, OS->diag, &mode,
+ &OS->factor, &npr, &info, &nfev, fjac, &ldfjac,
+ ipvt, qtf, wa1, wa2, wa3, wa4);
+ strcpy(lmfun_name, "lmdif");
+ }
+ else {
+ if (!isFunction(jac))
+ error("jac is not a function!");
+ PROTECT(OS->jcall = lang2(jac, OS->par));
+ PROTECT(eval_test = eval(OS->jcall, OS->env));
+ if (!IS_NUMERIC(eval_test))
+ error("evaluation of jac function returns non-numeric vector!");
+ if (length(eval_test) != n*m)
+ error("jac function must return numeric vector with length"
+ " == length(par) * length(fn(par, ...)). Your function"
+ " returns one with length %d while %d expected.",
+ length(eval_test), n*m);
+ UNPROTECT(1);
+ F77_CALL(lmder)(&fcn_lmder, &m, &n, par, fvec,
+ fjac, &ldfjac,
+ &OS->ftol, &OS->ptol, &OS->gtol,
+ &maxfev, OS->diag, &mode,
+ &OS->factor, &npr, &info, &nfev, &njev,
+ ipvt, qtf, wa1, wa2, wa3, wa4);
+ UNPROTECT(1);
+ strcpy(lmfun_name, "lmder");
+ }
+
+/*========================================================================*/
+
+ fcn_message(message, info, maxfev, OS->maxiter);
+ if (info < 1 || 9 < info)
+ warning("%s: info = %d. %s\n\n", lmfun_name, info, message);
+
+ PROTECT(sexp_hess = NEW_NUMERIC(n*n));
+ for (j = 0; j < n; j++)
+ for (i = 0; i < n; i++) {
+ perm[j*n + i] = (i + 1 == ipvt[j]) ? 1 : 0;
+ r [j*n + i] = (i <= j) ? fjac[j*ldfjac + i] : 0;
+ }
+
+ /* perm %*% t(r) %*% r %*% t(perm) *
+ * | |___r2__| | *
+ * | |_r2_x_perm_t_| *
+ * |_______hess_______| */
+
+ transpose(perm, n, n, perm_t);
+ crossprod(r, n, n, r, n, n, r2);
+ matprod (r2, n, n, perm_t, n, n, r2_x_perm_t);
+ matprod (perm, n, n, r2_x_perm_t, n, n, hess);
+
+ for (i = 0; i < n*n; i++)
+ NUMERIC_POINTER(sexp_hess)[i] = hess[i];
+
+ PROTECT(sexp_fvec = NEW_NUMERIC(m));
+ for (i = 0; i < m; i++)
+ NUMERIC_POINTER(sexp_fvec)[i] = fvec[i];
+
+ PROTECT(sexp_rsstrace = NEW_NUMERIC(OS->niter));
+ for (i = 0; i < OS->niter; i++)
+ NUMERIC_POINTER(sexp_rsstrace)[i] = OS->rsstrace[i];
+
+ PROTECT(sexp_info = NEW_INTEGER(1));
+ INTEGER_POINTER(sexp_info)[0] = info;
+
+ PROTECT(sexp_niter = NEW_INTEGER(1));
+ INTEGER_POINTER(sexp_niter)[0] = OS->niter-1;
+
+ PROTECT(sexp_message = NEW_STRING(1));
+ SET_STRING_ELT(sexp_message, 0, mkChar(message));
+ if (IS_NUMERIC(sexp_diag)) {
+ for (i = 0; i < n; i++)
+ NUMERIC_POINTER(sexp_diag)[i] = OS->diag[i];
+ }
+ else {
+ for (i = 0; i < n; i++)
+ NUMERIC_POINTER(VECTOR_ELT(sexp_diag, i))[0] = OS->diag[i];
+ }
+
+ PROTECT(out = NEW_LIST(8));
+ SET_VECTOR_ELT(out, 0, OS->par);
+ SET_VECTOR_ELT(out, 1, sexp_hess);
+ SET_VECTOR_ELT(out, 2, sexp_fvec);
+ SET_VECTOR_ELT(out, 3, sexp_info);
+ SET_VECTOR_ELT(out, 4, sexp_message);
+ SET_VECTOR_ELT(out, 5, sexp_diag);
+ SET_VECTOR_ELT(out, 6, sexp_niter);
+ SET_VECTOR_ELT(out, 7, sexp_rsstrace);
+
+ PROTECT(out_names = NEW_STRING(8));
+ SET_STRING_ELT(out_names, 0, mkChar("par"));
+ SET_STRING_ELT(out_names, 1, mkChar("hessian"));
+ SET_STRING_ELT(out_names, 2, mkChar("fvec"));
+ SET_STRING_ELT(out_names, 3, mkChar("info"));
+ SET_STRING_ELT(out_names, 4, mkChar("message"));
+ SET_STRING_ELT(out_names, 5, mkChar("diag"));
+ SET_STRING_ELT(out_names, 6, mkChar("niter"));
+ SET_STRING_ELT(out_names, 7, mkChar("rsstrace"));
+
+ SET_NAMES(out, out_names);
+
+ UNPROTECT(13);
+
+ return out;
+}
diff --git a/src/prod.c b/src/prod.c
new file mode 100644
index 0000000..3df0314
--- /dev/null
+++ b/src/prod.c
@@ -0,0 +1,52 @@
+/* The following are two functions from `src/main/array.c' file
+ * from R-2.0.1 sources */
+
+#include <R.h>
+#include <R_ext/Applic.h> /* for dgemm */
+
+void matprod(double *x, int nrx, int ncx,
+ double *y, int nry, int ncy, double *z)
+{
+ char *transa = "N", *transb = "N";
+ int i, j, k;
+ double one = 1.0, zero = 0.0, sum;
+ Rboolean have_na = FALSE;
+
+ if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
+ /* Don't trust the BLAS to handle NA/NaNs correctly: PR#4582
+ * The test is only O(n) here
+ */
+ for (i = 0; i < nrx*ncx; i++)
+ if (ISNAN(x[i])) {have_na = TRUE; break;}
+ if (!have_na)
+ for (i = 0; i < nry*ncy; i++)
+ if (ISNAN(y[i])) {have_na = TRUE; break;}
+ if (have_na) {
+ for (i = 0; i < nrx; i++)
+ for (k = 0; k < ncy; k++) {
+ sum = 0.0;
+ for (j = 0; j < ncx; j++)
+ sum += x[i + j * nrx] * y[j + k * nry];
+ z[i + k * nrx] = sum;
+ }
+ } else
+ F77_CALL(dgemm)(transa, transb, &nrx, &ncy, &ncx, &one,
+ x, &nrx, y, &nry, &zero, z, &nrx);
+ } else /* zero-extent operations should return zeroes */
+ for(i = 0; i < nrx*ncy; i++) z[i] = 0;
+}
+
+void crossprod(double *x, int nrx, int ncx,
+ double *y, int nry, int ncy, double *z)
+{
+ char *transa = "T", *transb = "N";
+ double one = 1.0, zero = 0.0;
+ if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
+ F77_CALL(dgemm)(transa, transb, &ncx, &ncy, &nrx, &one,
+ x, &nrx, y, &nry, &zero, z, &ncx);
+ }
+ else { /* zero-extent operations should return zeroes */
+ int i;
+ for(i = 0; i < ncx*ncy; i++) z[i] = 0;
+ }
+}
diff --git a/src/qform.f b/src/qform.f
new file mode 100644
index 0000000..087b247
--- /dev/null
+++ b/src/qform.f
@@ -0,0 +1,95 @@
+ subroutine qform(m,n,q,ldq,wa)
+ integer m,n,ldq
+ double precision q(ldq,m),wa(m)
+c **********
+c
+c subroutine qform
+c
+c this subroutine proceeds from the computed qr factorization of
+c an m by n matrix a to accumulate the m by m orthogonal matrix
+c q from its factored form.
+c
+c the subroutine statement is
+c
+c subroutine qform(m,n,q,ldq,wa)
+c
+c where
+c
+c m is a positive integer input variable set to the number
+c of rows of a and the order of q.
+c
+c n is a positive integer input variable set to the number
+c of columns of a.
+c
+c q is an m by m array. on input the full lower trapezoid in
+c the first min(m,n) columns of q contains the factored form.
+c on output q has been accumulated into a square matrix.
+c
+c ldq is a positive integer input variable not less than m
+c which specifies the leading dimension of the array q.
+c
+c wa is a work array of length m.
+c
+c subprograms called
+c
+c fortran-supplied ... min0
+c
+c argonne national laboratory. minpack project. march 1980.
+c burton s. garbow, kenneth e. hillstrom, jorge j. more
+c
+c **********
+ integer i,j,jm1,k,l,minmn,np1
+ double precision one,sum,temp,zero
+ data one,zero /1.0d0,0.0d0/
+c
+c zero out upper triangle of q in the first min(m,n) columns.
+c
+ minmn = min0(m,n)
+ if (minmn .lt. 2) go to 30
+ do 20 j = 2, minmn
+ jm1 = j - 1
+ do 10 i = 1, jm1
+ q(i,j) = zero
+ 10 continue
+ 20 continue
+ 30 continue
+c
+c initialize remaining columns to those of the identity matrix.
+c
+ np1 = n + 1
+ if (m .lt. np1) go to 60
+ do 50 j = np1, m
+ do 40 i = 1, m
+ q(i,j) = zero
+ 40 continue
+ q(j,j) = one
+ 50 continue
+ 60 continue
+c
+c accumulate q from its factored form.
+c
+ do 120 l = 1, minmn
+ k = minmn - l + 1
+ do 70 i = k, m
+ wa(i) = q(i,k)
+ q(i,k) = zero
+ 70 continue
+ q(k,k) = one
+ if (wa(k) .eq. zero) go to 110
+ do 100 j = k, m
+ sum = zero
+ do 80 i = k, m
+ sum = sum + q(i,j)*wa(i)
+ 80 continue
+ temp = sum/wa(k)
+ do 90 i = k, m
+ q(i,j) = q(i,j) - temp*wa(i)
+ 90 continue
+ 100 continue
+ 110 continue
+ 120 continue
+ return
+c
+c last card of subroutine qform.
+c
+ end
diff --git a/src/qrfac.f b/src/qrfac.f
new file mode 100644
index 0000000..cb68608
--- /dev/null
+++ b/src/qrfac.f
@@ -0,0 +1,164 @@
+ subroutine qrfac(m,n,a,lda,pivot,ipvt,lipvt,rdiag,acnorm,wa)
+ integer m,n,lda,lipvt
+ integer ipvt(lipvt)
+ logical pivot
+ double precision a(lda,n),rdiag(n),acnorm(n),wa(n)
+c **********
+c
+c subroutine qrfac
+c
+c this subroutine uses householder transformations with column
+c pivoting (optional) to compute a qr factorization of the
+c m by n matrix a. that is, qrfac determines an orthogonal
+c matrix q, a permutation matrix p, and an upper trapezoidal
+c matrix r with diagonal elements of nonincreasing magnitude,
+c such that a*p = q*r. the householder transformation for
+c column k, k = 1,2,...,min(m,n), is of the form
+c
+c t
+c i - (1/u(k))*u*u
+c
+c where u has zeros in the first k-1 positions. the form of
+c this transformation and the method of pivoting first
+c appeared in the corresponding linpack subroutine.
+c
+c the subroutine statement is
+c
+c subroutine qrfac(m,n,a,lda,pivot,ipvt,lipvt,rdiag,acnorm,wa)
+c
+c where
+c
+c m is a positive integer input variable set to the number
+c of rows of a.
+c
+c n is a positive integer input variable set to the number
+c of columns of a.
+c
+c a is an m by n array. on input a contains the matrix for
+c which the qr factorization is to be computed. on output
+c the strict upper trapezoidal part of a contains the strict
+c upper trapezoidal part of r, and the lower trapezoidal
+c part of a contains a factored form of q (the non-trivial
+c elements of the u vectors described above).
+c
+c lda is a positive integer input variable not less than m
+c which specifies the leading dimension of the array a.
+c
+c pivot is a logical input variable. if pivot is set true,
+c then column pivoting is enforced. if pivot is set false,
+c then no column pivoting is done.
+c
+c ipvt is an integer output array of length lipvt. ipvt
+c defines the permutation matrix p such that a*p = q*r.
+c column j of p is column ipvt(j) of the identity matrix.
+c if pivot is false, ipvt is not referenced.
+c
+c lipvt is a positive integer input variable. if pivot is false,
+c then lipvt may be as small as 1. if pivot is true, then
+c lipvt must be at least n.
+c
+c rdiag is an output array of length n which contains the
+c diagonal elements of r.
+c
+c acnorm is an output array of length n which contains the
+c norms of the corresponding columns of the input matrix a.
+c if this information is not needed, then acnorm can coincide
+c with rdiag.
+c
+c wa is a work array of length n. if pivot is false, then wa
+c can coincide with rdiag.
+c
+c subprograms called
+c
+c minpack-supplied ... dpmpar,enorm
+c
+c fortran-supplied ... dmax1,dsqrt,min0
+c
+c argonne national laboratory. minpack project. march 1980.
+c burton s. garbow, kenneth e. hillstrom, jorge j. more
+c
+c **********
+ integer i,j,jp1,k,kmax,minmn
+ double precision ajnorm,epsmch,one,p05,sum,temp,zero
+ double precision dpmpar,enorm
+ data one,p05,zero /1.0d0,5.0d-2,0.0d0/
+c
+c epsmch is the machine precision.
+c
+ epsmch = dpmpar(1)
+c
+c compute the initial column norms and initialize several arrays.
+c
+ do 10 j = 1, n
+ acnorm(j) = enorm(m,a(1,j))
+ rdiag(j) = acnorm(j)
+ wa(j) = rdiag(j)
+ if (pivot) ipvt(j) = j
+ 10 continue
+c
+c reduce a to r with householder transformations.
+c
+ minmn = min0(m,n)
+ do 110 j = 1, minmn
+ if (.not.pivot) go to 40
+c
+c bring the column of largest norm into the pivot position.
+c
+ kmax = j
+ do 20 k = j, n
+ if (rdiag(k) .gt. rdiag(kmax)) kmax = k
+ 20 continue
+ if (kmax .eq. j) go to 40
+ do 30 i = 1, m
+ temp = a(i,j)
+ a(i,j) = a(i,kmax)
+ a(i,kmax) = temp
+ 30 continue
+ rdiag(kmax) = rdiag(j)
+ wa(kmax) = wa(j)
+ k = ipvt(j)
+ ipvt(j) = ipvt(kmax)
+ ipvt(kmax) = k
+ 40 continue
+c
+c compute the householder transformation to reduce the
+c j-th column of a to a multiple of the j-th unit vector.
+c
+ ajnorm = enorm(m-j+1,a(j,j))
+ if (ajnorm .eq. zero) go to 100
+ if (a(j,j) .lt. zero) ajnorm = -ajnorm
+ do 50 i = j, m
+ a(i,j) = a(i,j)/ajnorm
+ 50 continue
+ a(j,j) = a(j,j) + one
+c
+c apply the transformation to the remaining columns
+c and update the norms.
+c
+ jp1 = j + 1
+ if (n .lt. jp1) go to 100
+ do 90 k = jp1, n
+ sum = zero
+ do 60 i = j, m
+ sum = sum + a(i,j)*a(i,k)
+ 60 continue
+ temp = sum/a(j,j)
+ do 70 i = j, m
+ a(i,k) = a(i,k) - temp*a(i,j)
+ 70 continue
+ if (.not.pivot .or. rdiag(k) .eq. zero) go to 80
+ temp = a(j,k)/rdiag(k)
+ rdiag(k) = rdiag(k)*dsqrt(dmax1(zero,one-temp**2))
+ if (p05*(rdiag(k)/wa(k))**2 .gt. epsmch) go to 80
+ rdiag(k) = enorm(m-j,a(jp1,k))
+ wa(k) = rdiag(k)
+ 80 continue
+ 90 continue
+ 100 continue
+ rdiag(j) = -ajnorm
+ 110 continue
+ return
+c
+c last card of subroutine qrfac.
+c
+ end
diff --git a/src/qrsolv.f b/src/qrsolv.f
new file mode 100644
index 0000000..f48954b
--- /dev/null
+++ b/src/qrsolv.f
@@ -0,0 +1,193 @@
+ subroutine qrsolv(n,r,ldr,ipvt,diag,qtb,x,sdiag,wa)
+ integer n,ldr
+ integer ipvt(n)
+ double precision r(ldr,n),diag(n),qtb(n),x(n),sdiag(n),wa(n)
+c **********
+c
+c subroutine qrsolv
+c
+c given an m by n matrix a, an n by n diagonal matrix d,
+c and an m-vector b, the problem is to determine an x which
+c solves the system
+c
+c a*x = b , d*x = 0 ,
+c
+c in the least squares sense.
+c
+c this subroutine completes the solution of the problem
+c if it is provided with the necessary information from the
+c qr factorization, with column pivoting, of a. that is, if
+c a*p = q*r, where p is a permutation matrix, q has orthogonal
+c columns, and r is an upper triangular matrix with diagonal
+c elements of nonincreasing magnitude, then qrsolv expects
+c the full upper triangle of r, the permutation matrix p,
+c and the first n components of (q transpose)*b. the system
+c a*x = b, d*x = 0, is then equivalent to
+c
+c t t
+c r*z = q *b , p *d*p*z = 0 ,
+c
+c where x = p*z. if this system does not have full rank,
+c then a least squares solution is obtained. on output qrsolv
+c also provides an upper triangular matrix s such that
+c
+c t t t
+c p *(a *a + d*d)*p = s *s .
+c
+c s is computed within qrsolv and may be of separate interest.
+c
+c the subroutine statement is
+c
+c subroutine qrsolv(n,r,ldr,ipvt,diag,qtb,x,sdiag,wa)
+c
+c where
+c
+c n is a positive integer input variable set to the order of r.
+c
+c r is an n by n array. on input the full upper triangle
+c must contain the full upper triangle of the matrix r.
+c on output the full upper triangle is unaltered, and the
+c strict lower triangle contains the strict upper triangle
+c (transposed) of the upper triangular matrix s.
+c
+c ldr is a positive integer input variable not less than n
+c which specifies the leading dimension of the array r.
+c
+c ipvt is an integer input array of length n which defines the
+c permutation matrix p such that a*p = q*r. column j of p
+c is column ipvt(j) of the identity matrix.
+c
+c diag is an input array of length n which must contain the
+c diagonal elements of the matrix d.
+c
+c qtb is an input array of length n which must contain the first
+c n elements of the vector (q transpose)*b.
+c
+c x is an output array of length n which contains the least
+c squares solution of the system a*x = b, d*x = 0.
+c
+c sdiag is an output array of length n which contains the
+c diagonal elements of the upper triangular matrix s.
+c
+c wa is a work array of length n.
+c
+c subprograms called
+c
+c fortran-supplied ... dabs,dsqrt
+c
+c argonne national laboratory. minpack project. march 1980.
+c burton s. garbow, kenneth e. hillstrom, jorge j. more
+c
+c **********
+ integer i,j,jp1,k,kp1,l,nsing
+ double precision cos,cotan,p5,p25,qtbpj,sin,sum,tan,temp,zero
+ data p5,p25,zero /5.0d-1,2.5d-1,0.0d0/
+c
+c copy r and (q transpose)*b to preserve input and initialize s.
+c in particular, save the diagonal elements of r in x.
+c
+ do 20 j = 1, n
+ do 10 i = j, n
+ r(i,j) = r(j,i)
+ 10 continue
+ x(j) = r(j,j)
+ wa(j) = qtb(j)
+ 20 continue
+c
+c eliminate the diagonal matrix d using a givens rotation.
+c
+ do 100 j = 1, n
+c
+c prepare the row of d to be eliminated, locating the
+c diagonal element using p from the qr factorization.
+c
+ l = ipvt(j)
+ if (diag(l) .eq. zero) go to 90
+ do 30 k = j, n
+ sdiag(k) = zero
+ 30 continue
+ sdiag(j) = diag(l)
+c
+c the transformations to eliminate the row of d
+c modify only a single element of (q transpose)*b
+c beyond the first n, which is initially zero.
+c
+ qtbpj = zero
+ do 80 k = j, n
+c
+c determine a givens rotation which eliminates the
+c appropriate element in the current row of d.
+c
+ if (sdiag(k) .eq. zero) go to 70
+ if (dabs(r(k,k)) .ge. dabs(sdiag(k))) go to 40
+ cotan = r(k,k)/sdiag(k)
+ sin = p5/dsqrt(p25+p25*cotan**2)
+ cos = sin*cotan
+ go to 50
+ 40 continue
+ tan = sdiag(k)/r(k,k)
+ cos = p5/dsqrt(p25+p25*tan**2)
+ sin = cos*tan
+ 50 continue
+c
+c compute the modified diagonal element of r and
+c the modified element of ((q transpose)*b,0).
+c
+ r(k,k) = cos*r(k,k) + sin*sdiag(k)
+ temp = cos*wa(k) + sin*qtbpj
+ qtbpj = -sin*wa(k) + cos*qtbpj
+ wa(k) = temp
+c
+c accumulate the tranformation in the row of s.
+c
+ kp1 = k + 1
+ if (n .lt. kp1) go to 70
+ do 60 i = kp1, n
+ temp = cos*r(i,k) + sin*sdiag(i)
+ sdiag(i) = -sin*r(i,k) + cos*sdiag(i)
+ r(i,k) = temp
+ 60 continue
+ 70 continue
+ 80 continue
+ 90 continue
+c
+c store the diagonal element of s and restore
+c the corresponding diagonal element of r.
+c
+ sdiag(j) = r(j,j)
+ r(j,j) = x(j)
+ 100 continue
+c
+c solve the triangular system for z. if the system is
+c singular, then obtain a least squares solution.
+c
+ nsing = n
+ do 110 j = 1, n
+ if (sdiag(j) .eq. zero .and. nsing .eq. n) nsing = j - 1
+ if (nsing .lt. n) wa(j) = zero
+ 110 continue
+ if (nsing .lt. 1) go to 150
+ do 140 k = 1, nsing
+ j = nsing - k + 1
+ sum = zero
+ jp1 = j + 1
+ if (nsing .lt. jp1) go to 130
+ do 120 i = jp1, nsing
+ sum = sum + r(i,j)*wa(i)
+ 120 continue
+ 130 continue
+ wa(j) = (wa(j) - sum)/sdiag(j)
+ 140 continue
+ 150 continue
+c
+c permute the components of z back to components of x.
+c
+ do 160 j = 1, n
+ l = ipvt(j)
+ x(l) = wa(j)
+ 160 continue
+ return
+c
+c last card of subroutine qrsolv.
+c
+ end
diff --git a/src/r1mpyq.f b/src/r1mpyq.f
new file mode 100644
index 0000000..ec99b96
--- /dev/null
+++ b/src/r1mpyq.f
@@ -0,0 +1,92 @@
+ subroutine r1mpyq(m,n,a,lda,v,w)
+ integer m,n,lda
+ double precision a(lda,n),v(n),w(n)
+c **********
+c
+c subroutine r1mpyq
+c
+c given an m by n matrix a, this subroutine computes a*q where
+c q is the product of 2*(n - 1) transformations
+c
+c gv(n-1)*...*gv(1)*gw(1)*...*gw(n-1)
+c
+c and gv(i), gw(i) are givens rotations in the (i,n) plane which
+c eliminate elements in the i-th and n-th planes, respectively.
+c q itself is not given, rather the information to recover the
+c gv, gw rotations is supplied.
+c
+c the subroutine statement is
+c
+c subroutine r1mpyq(m,n,a,lda,v,w)
+c
+c where
+c
+c m is a positive integer input variable set to the number
+c of rows of a.
+c
+c n is a positive integer input variable set to the number
+c of columns of a.
+c
+c a is an m by n array. on input a must contain the matrix
+c to be postmultiplied by the orthogonal matrix q
+c described above. on output a*q has replaced a.
+c
+c lda is a positive integer input variable not less than m
+c which specifies the leading dimension of the array a.
+c
+c v is an input array of length n. v(i) must contain the
+c information necessary to recover the givens rotation gv(i)
+c described above.
+c
+c w is an input array of length n. w(i) must contain the
+c information necessary to recover the givens rotation gw(i)
+c described above.
+c
+c subroutines called
+c
+c fortran-supplied ... dabs,dsqrt
+c
+c argonne national laboratory. minpack project. march 1980.
+c burton s. garbow, kenneth e. hillstrom, jorge j. more
+c
+c **********
+ integer i,j,nmj,nm1
+ double precision cos,one,sin,temp
+ data one /1.0d0/
+c
+c apply the first set of givens rotations to a.
+c
+ nm1 = n - 1
+ if (nm1 .lt. 1) go to 50
+ do 20 nmj = 1, nm1
+ j = n - nmj
+ if (dabs(v(j)) .gt. one) cos = one/v(j)
+ if (dabs(v(j)) .gt. one) sin = dsqrt(one-cos**2)
+ if (dabs(v(j)) .le. one) sin = v(j)
+ if (dabs(v(j)) .le. one) cos = dsqrt(one-sin**2)
+ do 10 i = 1, m
+ temp = cos*a(i,j) - sin*a(i,n)
+ a(i,n) = sin*a(i,j) + cos*a(i,n)
+ a(i,j) = temp
+ 10 continue
+ 20 continue
+c
+c apply the second set of givens rotations to a.
+c
+ do 40 j = 1, nm1
+ if (dabs(w(j)) .gt. one) cos = one/w(j)
+ if (dabs(w(j)) .gt. one) sin = dsqrt(one-cos**2)
+ if (dabs(w(j)) .le. one) sin = w(j)
+ if (dabs(w(j)) .le. one) cos = dsqrt(one-sin**2)
+ do 30 i = 1, m
+ temp = cos*a(i,j) + sin*a(i,n)
+ a(i,n) = -sin*a(i,j) + cos*a(i,n)
+ a(i,j) = temp
+ 30 continue
+ 40 continue
+ 50 continue
+ return
+c
+c last card of subroutine r1mpyq.
+c
+ end
diff --git a/src/r1updt.f b/src/r1updt.f
new file mode 100644
index 0000000..e034973
--- /dev/null
+++ b/src/r1updt.f
@@ -0,0 +1,207 @@
+ subroutine r1updt(m,n,s,ls,u,v,w,sing)
+ integer m,n,ls
+ logical sing
+ double precision s(ls),u(m),v(n),w(m)
+c **********
+c
+c subroutine r1updt
+c
+c given an m by n lower trapezoidal matrix s, an m-vector u,
+c and an n-vector v, the problem is to determine an
+c orthogonal matrix q such that
+c
+c t
+c (s + u*v )*q
+c
+c is again lower trapezoidal.
+c
+c this subroutine determines q as the product of 2*(n - 1)
+c transformations
+c
+c gv(n-1)*...*gv(1)*gw(1)*...*gw(n-1)
+c
+c where gv(i), gw(i) are givens rotations in the (i,n) plane
+c which eliminate elements in the i-th and n-th planes,
+c respectively. q itself is not accumulated, rather the
+c information to recover the gv, gw rotations is returned.
+c
+c the subroutine statement is
+c
+c subroutine r1updt(m,n,s,ls,u,v,w,sing)
+c
+c where
+c
+c m is a positive integer input variable set to the number
+c of rows of s.
+c
+c n is a positive integer input variable set to the number
+c of columns of s. n must not exceed m.
+c
+c s is an array of length ls. on input s must contain the lower
+c trapezoidal matrix s stored by columns. on output s contains
+c the lower trapezoidal matrix produced as described above.
+c
+c ls is a positive integer input variable not less than
+c (n*(2*m-n+1))/2.
+c
+c u is an input array of length m which must contain the
+c vector u.
+c
+c v is an array of length n. on input v must contain the vector
+c v. on output v(i) contains the information necessary to
+c recover the givens rotation gv(i) described above.
+c
+c w is an output array of length m. w(i) contains information
+c necessary to recover the givens rotation gw(i) described
+c above.
+c
+c sing is a logical output variable. sing is set true if any
+c of the diagonal elements of the output s are zero. otherwise
+c sing is set false.
+c
+c subprograms called
+c
+c minpack-supplied ... dpmpar
+c
+c fortran-supplied ... dabs,dsqrt
+c
+c argonne national laboratory. minpack project. march 1980.
+c burton s. garbow, kenneth e. hillstrom, jorge j. more,
+c john l. nazareth
+c
+c **********
+ integer i,j,jj,l,nmj,nm1
+ double precision cos,cotan,giant,one,p5,p25,sin,tan,tau,temp,
+ * zero
+ double precision dpmpar
+ data one,p5,p25,zero /1.0d0,5.0d-1,2.5d-1,0.0d0/
+c
+c giant is the largest magnitude.
+c
+ giant = dpmpar(3)
+c
+c initialize the diagonal element pointer.
+c
+ jj = (n*(2*m - n + 1))/2 - (m - n)
+c
+c move the nontrivial part of the last column of s into w.
+c
+ l = jj
+ do 10 i = n, m
+ w(i) = s(l)
+ l = l + 1
+ 10 continue
+c
+c rotate the vector v into a multiple of the n-th unit vector
+c in such a way that a spike is introduced into w.
+c
+ nm1 = n - 1
+ if (nm1 .lt. 1) go to 70
+ do 60 nmj = 1, nm1
+ j = n - nmj
+ jj = jj - (m - j + 1)
+ w(j) = zero
+ if (v(j) .eq. zero) go to 50
+c
+c determine a givens rotation which eliminates the
+c j-th element of v.
+c
+ if (dabs(v(n)) .ge. dabs(v(j))) go to 20
+ cotan = v(n)/v(j)
+ sin = p5/dsqrt(p25+p25*cotan**2)
+ cos = sin*cotan
+ tau = one
+ if (dabs(cos)*giant .gt. one) tau = one/cos
+ go to 30
+ 20 continue
+ tan = v(j)/v(n)
+ cos = p5/dsqrt(p25+p25*tan**2)
+ sin = cos*tan
+ tau = sin
+ 30 continue
+c
+c apply the transformation to v and store the information
+c necessary to recover the givens rotation.
+c
+ v(n) = sin*v(j) + cos*v(n)
+ v(j) = tau
+c
+c apply the transformation to s and extend the spike in w.
+c
+ l = jj
+ do 40 i = j, m
+ temp = cos*s(l) - sin*w(i)
+ w(i) = sin*s(l) + cos*w(i)
+ s(l) = temp
+ l = l + 1
+ 40 continue
+ 50 continue
+ 60 continue
+ 70 continue
+c
+c add the spike from the rank 1 update to w.
+c
+ do 80 i = 1, m
+ w(i) = w(i) + v(n)*u(i)
+ 80 continue
+c
+c eliminate the spike.
+c
+ sing = .false.
+ if (nm1 .lt. 1) go to 140
+ do 130 j = 1, nm1
+ if (w(j) .eq. zero) go to 120
+c
+c determine a givens rotation which eliminates the
+c j-th element of the spike.
+c
+ if (dabs(s(jj)) .ge. dabs(w(j))) go to 90
+ cotan = s(jj)/w(j)
+ sin = p5/dsqrt(p25+p25*cotan**2)
+ cos = sin*cotan
+ tau = one
+ if (dabs(cos)*giant .gt. one) tau = one/cos
+ go to 100
+ 90 continue
+ tan = w(j)/s(jj)
+ cos = p5/dsqrt(p25+p25*tan**2)
+ sin = cos*tan
+ tau = sin
+ 100 continue
+c
+c apply the transformation to s and reduce the spike in w.
+c
+ l = jj
+ do 110 i = j, m
+ temp = cos*s(l) + sin*w(i)
+ w(i) = -sin*s(l) + cos*w(i)
+ s(l) = temp
+ l = l + 1
+ 110 continue
+c
+c store the information necessary to recover the
+c givens rotation.
+c
+ w(j) = tau
+ 120 continue
+c
+c test for zero diagonal elements in the output s.
+c
+ if (s(jj) .eq. zero) sing = .true.
+ jj = jj + (m - j + 1)
+ 130 continue
+ 140 continue
+c
+c move w back into the last column of the output s.
+c
+ l = jj
+ do 150 i = n, m
+ s(l) = w(i)
+ l = l + 1
+ 150 continue
+ if (s(jj) .eq. zero) sing = .true.
+ return
+c
+c last card of subroutine r1updt.
+c
+ end
diff --git a/src/rwupdt.f b/src/rwupdt.f
new file mode 100644
index 0000000..05282b5
--- /dev/null
+++ b/src/rwupdt.f
@@ -0,0 +1,113 @@
+ subroutine rwupdt(n,r,ldr,w,b,alpha,cos,sin)
+ integer n,ldr
+ double precision alpha
+ double precision r(ldr,n),w(n),b(n),cos(n),sin(n)
+c **********
+c
+c subroutine rwupdt
+c
+c given an n by n upper triangular matrix r, this subroutine
+c computes the qr decomposition of the matrix formed when a row
+c is added to r. if the row is specified by the vector w, then
+c rwupdt determines an orthogonal matrix q such that when the
+c n+1 by n matrix composed of r augmented by w is premultiplied
+c by (q transpose), the resulting matrix is upper trapezoidal.
+c the matrix (q transpose) is the product of n transformations
+c
+c g(n)*g(n-1)* ... *g(1)
+c
+c where g(i) is a givens rotation in the (i,n+1) plane which
+c eliminates elements in the (n+1)-st plane. rwupdt also
+c computes the product (q transpose)*c where c is the
+c (n+1)-vector (b,alpha). q itself is not accumulated, rather
+c the information to recover the g rotations is supplied.
+c
+c the subroutine statement is
+c
+c subroutine rwupdt(n,r,ldr,w,b,alpha,cos,sin)
+c
+c where
+c
+c n is a positive integer input variable set to the order of r.
+c
+c r is an n by n array. on input the upper triangular part of
+c r must contain the matrix to be updated. on output r
+c contains the updated triangular matrix.
+c
+c ldr is a positive integer input variable not less than n
+c which specifies the leading dimension of the array r.
+c
+c w is an input array of length n which must contain the row
+c vector to be added to r.
+c
+c b is an array of length n. on input b must contain the
+c first n elements of the vector c. on output b contains
+c the first n elements of the vector (q transpose)*c.
+c
+c alpha is a variable. on input alpha must contain the
+c (n+1)-st element of the vector c. on output alpha contains
+c the (n+1)-st element of the vector (q transpose)*c.
+c
+c cos is an output array of length n which contains the
+c cosines of the transforming givens rotations.
+c
+c sin is an output array of length n which contains the
+c sines of the transforming givens rotations.
+c
+c subprograms called
+c
+c fortran-supplied ... dabs,dsqrt
+c
+c argonne national laboratory. minpack project. march 1980.
+c burton s. garbow, dudley v. goetschel, kenneth e. hillstrom,
+c jorge j. more
+c
+c **********
+ integer i,j,jm1
+ double precision cotan,one,p5,p25,rowj,tan,temp,zero
+ data one,p5,p25,zero /1.0d0,5.0d-1,2.5d-1,0.0d0/
+c
+ do 60 j = 1, n
+ rowj = w(j)
+ jm1 = j - 1
+c
+c apply the previous transformations to
+c r(i,j), i=1,2,...,j-1, and to w(j).
+c
+ if (jm1 .lt. 1) go to 20
+ do 10 i = 1, jm1
+ temp = cos(i)*r(i,j) + sin(i)*rowj
+ rowj = -sin(i)*r(i,j) + cos(i)*rowj
+ r(i,j) = temp
+ 10 continue
+ 20 continue
+c
+c determine a givens rotation which eliminates w(j).
+c
+ cos(j) = one
+ sin(j) = zero
+ if (rowj .eq. zero) go to 50
+ if (dabs(r(j,j)) .ge. dabs(rowj)) go to 30
+ cotan = r(j,j)/rowj
+ sin(j) = p5/dsqrt(p25+p25*cotan**2)
+ cos(j) = sin(j)*cotan
+ go to 40
+ 30 continue
+ tan = rowj/r(j,j)
+ cos(j) = p5/dsqrt(p25+p25*tan**2)
+ sin(j) = cos(j)*tan
+ 40 continue
+c
+c apply the current transformation to r(j,j), b(j), and alpha.
+c
+ r(j,j) = cos(j)*r(j,j) + sin(j)*rowj
+ temp = cos(j)*b(j) + sin(j)*alpha
+ alpha = -sin(j)*b(j) + cos(j)*alpha
+ b(j) = temp
+ 50 continue
+ 60 continue
+ return
+c
+c last card of subroutine rwupdt.
+c
+ end
diff --git a/src/transpose.c b/src/transpose.c
new file mode 100644
index 0000000..ca1dfec
--- /dev/null
+++ b/src/transpose.c
@@ -0,0 +1,8 @@
+void transpose(double *x, int nrx, int ncx, double *y)
+{
+ int i, j;
+
+ for (i = 0; i < ncx; i++)
+ for (j = 0; j < nrx; j++)
+ y[j*ncx + i] = x[i*nrx + j];
+}
diff --git a/src/vector.c b/src/vector.c
new file mode 100644
index 0000000..f70f94a
--- /dev/null
+++ b/src/vector.c
@@ -0,0 +1,11 @@
+#include <R.h>
+
+double *real_vector(int n)
+{
+ return (double*) R_alloc(n, sizeof(double));
+}
+
+int *int_vector(int n)
+{
+ return (int*) R_alloc(n, sizeof(int));
+}
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/r-cran-minpack.lm.git
More information about the debian-med-commit
mailing list