https://github.com/cran/fields
Tip revision: 6c8b30169bba182a68765ee3cb9b4e2ef7d38332 authored by Doug Nychka on 16 November 2011, 00:00:00 UTC
version 6.6.3
version 6.6.3
Tip revision: 6c8b301
conjugate.gradient.r
# fields, Tools for spatial data
# Copyright 2004-2011, Institute for Mathematics Applied Geosciences
# University Corporation for Atmospheric Research
# Licensed under the GPL -- www.gpl.org/licenses/gpl.html
"conjugate.gradient" <- function(b, multAx, start,
tol = 1e-05, kmax = 25, verbose = TRUE, ...) {
call <- match.call()
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))
}