https://github.com/cran/fields
Tip revision: 6769ffc81115fbf0bf7d9c566cf7ac81be0049dc authored by Doug Nychka on 25 July 2005, 00:00:00 UTC
version 3.04
version 3.04
Tip revision: 6769ffc
conjugate.gradient.r
"conjugate.gradient" <-
function (b, multAx, start = NULL, tol = 1e-05, kmax = 25, verbose = TRUE,
...)
{
call <- match.call()
if (is.null(start)) {
x <- rep(0, length(y))
r <- b
}
else {
x <- start
r <- b - multAx(x, ...)
}
p <- r
rho <- rep(NA, kmax + 1)
rho[1 + (0)] <- sum(r^2)
test <- sqrt(sum(b^2)) * tol
niter <- 1
k <- 0
for (k in 1:kmax) {
niter <- niter + 1
if (k != 1) {
beta <- rho[1 + (k - 1)]/rho[1 + (k - 2)]
p <- r + beta * p
}
w <- multAx(p, ...)
alpha <- rho[1 + (k - 1)]/(sum(p * w))
x <- x + alpha * p
r <- r - alpha * w
rho[1 + k] <- sum(r^2)
if (verbose) {
cat("iter", k, " crit : ", signif(sqrt(rho[1 + k]),
4), signif(test, 4), fill = TRUE)
}
if (sqrt(rho[1 + k]) < test) {
niter <- k + 1
break
}
}
list(call = call, x = x, residuals = r, niter, conv = list(rho = rho,
test = test, maxiter = kmax, niter = niter))
}