[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