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
fields.diagonalize2.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
"fields.diagonalize2" <- function(A,B, verbose=FALSE){
M<- nrow(A)
hold.AB <- eigen( A+B, symmetric=TRUE)
if( verbose){
cat("log 10 condition number of A +B in fields.diagonlize2", fill=TRUE)
print( log10( max(hold.AB$values)/ min(hold.AB$values)))}
# inverse square root of A+B
hold.AB<- (t(hold.AB$vectors)*(1/sqrt( hold.AB$values)))
hold.B <- eigen( hold.AB %*% A %*% t(hold.AB), symmetric = TRUE)
G <- t(hold.B$vectors)%*%hold.AB
D.A<- hold.B$values
# remove some large temporary matrices.
remove( hold.AB)
remove( hold.B)
# crank on finding G and D.
G<- (1/sqrt(D.A))*G
D<- colSums( t(G) * (B) %*% t(G))
# sort from largest to smallest and take transpose ---
# this will now matches old version in fields.diagonalize
D<- D[M:1]
G<- t(G[M:1,])
# to test:
# test.for.zero( t(G) %*% (A) %*% (G), diag(1,M), tag="A test" )
# test.for.zero( t(G) %*% (B) %*% (G), diag(D,M), tag="B test" )
list(G = G, D = D)
}