Commit 7aa70fc3 authored by Dirk Eddelbuettel's avatar Dirk Eddelbuettel

Import Upstream version 1.7-13

parent d3eb4372
Package: mgcv
Version: 1.7-12
Version: 1.7-13
Author: Simon Wood <simon.wood@r-project.org>
Maintainer: Simon Wood <simon.wood@r-project.org>
Title: GAMs with GCV/AIC/REML smoothness estimation and GAMMs by PQL
Description: Routines for GAMs and other generalized ridge regression
with multiple smoothing parameter selection by GCV, REML
or UBRE/AIC. Also GAMMs by REML or PQL. Includes a gam()
function.
Description: Routines for GAMs and other generalized ridge regression
with multiple smoothing parameter selection by GCV, REML or
UBRE/AIC. Also GAMMs by REML or PQL. Includes a gam() function.
Priority: recommended
Depends: R (>= 2.14.0)
Imports: graphics, stats, nlme, Matrix
Depends: R (>= 2.14.0), stats, graphics
Imports: nlme, methods, Matrix
Suggests: nlme (>= 3.1-64), splines, Matrix, parallel
LazyLoad: yes
ByteCompile: yes
License: GPL (>= 2)
Packaged: 2011-12-11 21:23:23 UTC; sw283
Packaged: 2012-01-21 16:33:06 UTC; sw283
Repository: CRAN
Date/Publication: 2011-12-12 12:21:51
Date/Publication: 2012-01-22 09:59:42
af5be24b4a8a4b948f42d61631c198fc *DESCRIPTION
abc5e760b0f5dd22c81d7614c0b304f0 *NAMESPACE
5d5d72ee6d284c96b4525e9eb748bc0f *DESCRIPTION
29e13076f8e7c500f10e2b64b0821984 *NAMESPACE
ecfb144fb5214dde68dffac22f219a1f *R/bam.r
b160632e8f38fa99470e2f8cba82a495 *R/gam.fit3.r
902657a0ee2dedc3fdfa501bf3b37c5b *R/gam.sim.r
e137c06cabb48551c18cf0cc3512d297 *R/gamm.r
da8dd27eb3b4b59a46ef58c789a38e56 *R/mgcv.r
c61836edb704dbd7b718c754d714a291 *R/mgcv.r
bf70158e37e33ea1136efdeab97569f8 *R/plots.r
8e2c7eddd2a206d05d708d4496f7eb31 *R/smooth.r
d20083082ba7e1bf361ac7d404efd8a3 *R/smooth.r
fe9745f610246ee1f31eb915ca0d76a9 *R/sparse.r
c5051a5b85fa7e9fa8ffdf9d52339408 *changeLog
76637934ae66a4b74a0637e698f71469 *changeLog
e468195a83fab90da8e760c2c3884bd3 *data/columb.polys.rda
40874e3ced720a596750f499ded8a60a *data/columb.rda
88d77139cc983317b6acd8c5f1252ab9 *gnugpl2.txt
......@@ -18,7 +18,7 @@ f693920e12f8a6f2b6cab93648628150 *index
c51c9b8c9c73f81895176ded39b91394 *man/Predict.matrix.cr.smooth.Rd
612ab6354541ebe38a242634d73b66ba *man/Tweedie.Rd
6d711de718a09e1e1ae2a6967abade33 *man/anova.gam.Rd
03fadfd98a6f40d6d517ec895859cff3 *man/bam.Rd
c4d1ad309698994e7c1c75f7db294a58 *man/bam.Rd
b385d6d5419d0d6aefe03af1a79d5c4e *man/bam.update.Rd
4e925cb579f4693d1b8ec2d5092c0b37 *man/cSplineDes.Rd
9753a8051d9b495d9855537f3f26f491 *man/choose.k.Rd
......@@ -40,7 +40,7 @@ fd98327327ba74bb1a61a6519f12e936 *man/gam.convergence.Rd
dd35a8a851460c2d2106c03d544c8241 *man/gam.models.Rd
468d116a2ef9e60f683af48f4f100ef5 *man/gam.outer.Rd
7e5ba69a44bc937ddca04e4f153c7975 *man/gam.selection.Rd
71e12bd2f7a166ca573aaf27fd984b25 *man/gam.side.Rd
76651917bd61fc6bc447bbb40b887236 *man/gam.side.Rd
78588cf8ed0af8eca70bba3bbed64dbe *man/gam.vcomp.Rd
278e0b3aa7baa44dfb96e235ceb07f4c *man/gam2objective.Rd
4d5b3b1266edc31ce3b0e6be11ee9166 *man/gamObject.Rd
......@@ -85,14 +85,14 @@ f77ca1471881d2f93c74864d076c0a0e *man/rTweedie.Rd
472cc8de5077f5511fe7eb2ad454767e *man/rig.Rd
7258dfc0149fff020054117fd2ee6bd8 *man/s.Rd
690dddb35fd5986a6eeb39dd79fa33f9 *man/slanczos.Rd
47f5ed6af7c784f26e0fe8e027ee1d90 *man/smooth.construct.Rd
3026df475af053acd2b742af875f626c *man/smooth.construct.Rd
10a456bc3d7151f60be24a7e3bba3d30 *man/smooth.construct.ad.smooth.spec.Rd
9f05449c0114748608f99693ec25cf3f *man/smooth.construct.cr.smooth.spec.Rd
131c3f2131138ba5d6f6bcde4be5ac31 *man/smooth.construct.ds.smooth.spec.Rd
1bb6748d4d2934e48f0572bc5114ffcb *man/smooth.construct.fs.smooth.spec.Rd
e0c041d7ec47290741e82fd71c408995 *man/smooth.construct.mrf.smooth.spec.Rd
78688702b6396f24692b74c554483428 *man/smooth.construct.ps.smooth.spec.Rd
0b134a171b9c5f8194787bf243f9c962 *man/smooth.construct.re.smooth.spec.Rd
d202c6718fb1138fdd99e6102250aedf *man/smooth.construct.re.smooth.spec.Rd
53d7986dd7a54c1edb13ce84dbfe34a2 *man/smooth.construct.sos.smooth.spec.Rd
8672633a1fad8df3cb1f53d7fa883620 *man/smooth.construct.tensor.smooth.spec.Rd
4b9bd43c3acbab6ab0159d59967e19db *man/smooth.construct.tp.smooth.spec.Rd
......@@ -123,7 +123,7 @@ d40012dcda1a10ee535a9b3de9b46c19 *src/gdi.c
49af97195accb65adc75620183d39a4c *src/general.h
da280ee5538a828afde0a4f6c7b8328a *src/init.c
8b37eb0db498a3867dc83364dc65f146 *src/magic.c
ebd3647ff460932bada8f49159f40ac8 *src/mat.c
066af9db587e5fe6e5cc4ff8c09ae9c2 *src/mat.c
d21847ac9a1f91ee9446c70bd93a490a *src/matrix.c
54ce9309b17024ca524e279612a869d6 *src/matrix.h
08c94a2af4cd047ecd79871ecbafe33a *src/mgcv.c
......
......@@ -53,15 +53,22 @@ export(anova.gam, bam, bam.update, concurvity, cSplineDes,
Tweedie,uniquecombs, vcov.gam, vis.gam)
importFrom(grDevices,cm.colors,gray,heat.colors,terrain.colors,topo.colors)
importFrom(graphics,axis,box,contour,hist,lines,mtext, par, persp,plot,points,polygon,strheight,strwidth,text)
importFrom(stats,.checkMFClasses,.getXlevels, as.formula, anova,anova.glmlist,coef,cooks.distance,cor,
delete.response,family,fitted,formula,
gaussian,glm,influence,lm,logLik,median,model.frame,model.matrix,model.offset,na.pass,napredict,naresid,nlm,optim,
pchisq,pf,pnorm,pt,
predict,printCoefmat ,quantile,qqnorm, reformulate,residuals,rnorm,runif,termplot,terms,terms.formula, uniroot,
var,vcov)
importFrom(graphics,axis,box,contour,hist,lines,mtext, par, persp,plot,points,
polygon,strheight,strwidth,text)
#importFrom(stats,.checkMFClasses,.getXlevels, as.formula, anova,anova.glmlist,
# coef,cooks.distance,cor,delete.response,family,fitted,formula,
# gaussian,glm,influence,lm,logLik,median,model.frame,model.matrix,
# model.offset,na.pass,napredict,naresid,nlm,optim,pchisq,pf,pnorm,pt,
# predict,printCoefmat ,quantile,qqnorm, reformulate,residuals,rnorm,
# runif,termplot,terms,terms.formula, uniroot,
# var,vcov)
importFrom(stats,anova,influence,cooks.distance,logLik,vcov,residuals,predict,
model.matrix)
importFrom(nlme,Dim,corMatrix,logDet,pdConstruct,pdFactor,pdMatrix)
importFrom(Matrix,t,mean,colMeans,colSums,chol,solve,diag)
#importFrom(Matrix,t,mean,colMeans,colSums,chol,solve,diag)
import(Matrix)
importFrom(methods,cbind2)
S3method(anova, gam)
S3method(influence, gam)
......
......@@ -374,7 +374,24 @@ fixDependence <- function(X1,X2,tol=.Machine$double.eps^.5,rank.def=0)
}
gam.side <- function(sm,Xp,tol=.Machine$double.eps^.5)
augment.smX <- function(sm,nobs,np) {
## augments a smooth model matrix with a square root penalty matrix for
## identifiability constraint purposes.
ns <- length(sm$S) ## number of penalty matrices
if (ns==0) { ## nothing to do
return(rbind(sm$X),matrix(0,np,np))
}
St <- sm$S[[1]] ## total scaled penalty
St <- St/sqrt(sum(St^2))
if (ns>1) for (i in 2:ns) St <- St + sm$S[[i]]/sqrt(sum(sm$S[[i]]^2))
rS <- mroot(St,rank=ncol(St)) ## get sqrt of penalty
rS <- rS/sqrt(mean(rS^2))
X <- rbind(sm$X/sqrt(mean(sm$X^2)),matrix(0,np,ncol(sm$X))) ## create augmented model matrix
X[nobs+sm$p.ind,] <- t(rS) ## add in
X ## scaled augmented model matrix
}
gam.side <- function(sm,Xp,tol=.Machine$double.eps^.5,with.pen=FALSE)
# works through a list of smooths, sm, aiming to identify nested or partially
# nested terms, and impose identifiability constraints on them.
# Xp is the parametric model matrix. It is needed in order to check whether
......@@ -402,7 +419,9 @@ gam.side <- function(sm,Xp,tol=.Machine$double.eps^.5)
## Need to test for intercept or equivalent in Xp
intercept <- FALSE
if (ncol(Xp)) {
## first check columns directly...
if (sum(apply(Xp,2,sd)<.Machine$double.eps^.75)>0) intercept <- TRUE else {
## no constant column, so need to check span of Xp...
f <- rep(1,nrow(Xp))
ff <- qr.fitted(qr(Xp),f)
if (max(abs(ff-f))<.Machine$double.eps^.75) intercept <- TRUE
......@@ -428,19 +447,41 @@ gam.side <- function(sm,Xp,tol=.Machine$double.eps^.5)
## the smooths of which it is an argument, arranged in ascending
## order of dimension.
if (maxDim==1) stop("model has repeated 1-d smooths of same variable.")
## Now set things up to enable term specific model matrices to be
## augmented with square root penalties, on the fly...
if (with.pen) {
k <- 1
for (i in 1:m) { ## create parameter indices for each term
k1 <- k + ncol(sm[[i]]$X) - 1
sm[[i]]$p.ind <- k:k1
k <- k1 + 1
}
np <- k-1 ## number of penalized parameters
}
nobs <- nrow(sm[[1]]$X) ## number of observations
for (d in 2:maxDim) { ## work up through dimensions
for (i in 1:m) { ## work through smooths
if (sm[[i]]$dim == d) { ## check for nesting
X1 <- matrix(1,nrow(sm[[i]]$X),as.integer(intercept))
if (sm[[i]]$dim == d&&sm[[i]]$side.constrain) { ## check for nesting
if (with.pen) X1 <- matrix(c(rep(1,nobs),rep(0,np)),nobs+np,as.integer(intercept)) else
X1 <- matrix(1,nobs,as.integer(intercept))
for (j in 1:d) { ## work through variables
b <- sm.id[[sm[[i]]$vn[j]]] # list of smooths dependent on this variable
k <- (1:length(b))[b==i] ## locate current smooth in list
if (k>1) for (l in 1:(k-1)) { ## collect X columns
X1 <- cbind(X1,sm[[b[l]]]$X)
if (with.pen) { ## need to use augmented model matrix in testing
if (is.null(sm[[b[l]]]$Xa)) sm[[b[l]]]$Xa <- augment.smX(sm[[b[l]]],nobs,np)
X1 <- cbind(X1,sm[[b[l]]]$Xa)
} else X1 <- cbind(X1,sm[[b[l]]]$X) ## penalties not considered
}
} ## Now X1 contains columns for all lower dimensional terms
if (ncol(X1)==as.integer(intercept)) ind <- NULL else
ind <- fixDependence(X1,sm[[i]]$X,tol=tol)
if (ncol(X1)==as.integer(intercept)) ind <- NULL else {
if (with.pen) {
if (is.null(sm[[i]]$Xa)) sm[[i]]$Xa <- augment.smX(sm[[i]],nobs,np)
ind <- fixDependence(X1,sm[[i]]$Xa,tol=tol)
} else ind <- fixDependence(X1,sm[[i]]$X,tol=tol)
}
## ... the columns to zero to ensure independence
if (!is.null(ind)) {
sm[[i]]$X <- sm[[i]]$X[,-ind]
......@@ -464,18 +505,28 @@ gam.side <- function(sm,Xp,tol=.Machine$double.eps^.5)
sm[[i]]$df <- ncol(sm[[i]]$X)
attr(sm[[i]],"del.index") <- ind
## Now deal with case in which prediction constraints differ from fit constraints
if (!is.null(sm[[i]]$Xp)) { ## need to get deletion indeces under prediction parameterization
if (!is.null(sm[[i]]$Xp)) { ## need to get deletion indices under prediction parameterization
## Note that: i) it doesn't matter what the identifiability con on X1 is
## ii) the degree of rank deficiency can't be changed by an identifiability con
ind <- fixDependence(X1,sm[[i]]$Xp,rank.def=length(ind))
if (with.pen) {
smi <- sm[[i]] ## clone smooth
smi$X <- smi$Xp ## copy prediction Xp to X slot.
smi$S <- smi$Sp ## and make sure penalty parameterization matches.
Xpa <- augment.smX(smi,nobs,np)
ind <- fixDependence(X1,Xpa,rank.def=length(ind))
} else ind <- fixDependence(X1,sm[[i]]$Xp,rank.def=length(ind))
sm[[i]]$Xp <- sm[[i]]$Xp[,-ind]
attr(sm[[i]],"del.index") <- ind ## over-writes original
}
}
sm[[i]]$vn <- NULL
} ## end if
} ## end if (sm[[i]]$dim == d)
} ## end i in 1:m loop
}
if (with.pen) for (i in 1:m) { ## remove working variables that were added
sm[[i]]$p.ind <- NULL
if (!is.null(sm[[i]]$Xa)) sm[[i]]$Xa <- NULL
}
sm
}
......@@ -3466,6 +3517,10 @@ set.mgcv.options <- function()
options(mgcv.vc.logrange=25)
}
.onLoad <- function(...) {
set.mgcv.options()
}
.onAttach <- function(...) {
print.mgcv.version()
set.mgcv.options()
......@@ -3473,11 +3528,12 @@ set.mgcv.options <- function()
.onUnload <- function(libpath) library.dynam.unload("mgcv", libpath)
.First.lib <- function(lib, pkg) {
library.dynam("mgcv", pkg, lib)
print.mgcv.version()
set.mgcv.options()
}
#.First.lib <- function(lib, pkg) {
## defunct
# library.dynam("mgcv", pkg, lib)
# print.mgcv.version()
# set.mgcv.options()
#}
###############################################################################
......
......@@ -1806,7 +1806,7 @@ smooth.construct.ad.smooth.spec<-function(object,data,knots)
########################################################
smooth.construct.re.smooth.spec<-function(object,data,knots)
smooth.construct.re.smooth.spec <- function(object,data,knots)
## a simple random effects constructor method function
## basic idea is that s(x,f,z,...,bs="re") generates model matrix
## corresponding to ~ x:f:z: ... - 1. Corresponding coefficients
......@@ -1832,6 +1832,7 @@ smooth.construct.re.smooth.spec<-function(object,data,knots)
## need to store formula (levels taken care of by calling function)
object$form <- form
object$side.constrain <- FALSE ## don't apply side constraints
object$plot.me <- TRUE ## "re" terms can be plotted by plot.gam
object$te.ok <- 2 ## these terms are suitable as te marginals, but
## can not be plotted
......@@ -2671,6 +2672,11 @@ smoothCon <- function(object,data,knots,absorb.cons=FALSE,scale.penalty=TRUE,n=n
## plot.me tells `plot.gam' whether or not to plot the term
if (is.null(sm$plot.me)) sm$plot.me <- TRUE
## add side.constrain indicator if missing
## `side.constrain' tells gam.side, whether term should be constrained
## as a result of any nesting detected...
if (is.null(sm$side.constrain)) sm$side.constrain <- TRUE
## automatically produce centering constraint...
## must be done here on original model matrix to ensure same
## basis for all `id' linked terms
......@@ -2865,6 +2871,13 @@ smoothCon <- function(object,data,knots,absorb.cons=FALSE,scale.penalty=TRUE,n=n
for (i in 1:length(sml)) { ## loop through smooth list
sml[[i]]$Xp <- t(qr.qty(qrcp,t(sml[[i]]$X))[(pj+1):k,]) ## form XZ
sml[[i]]$Cp <- NULL
if (length(sml[[i]]$S)) { ## gam.side requires penalties in prediction para
sml[[i]]$Sp <- sml[[i]]$S ## penalties in prediction parameterization
for (l in 1:length(sml[[i]]$S)) { # some smooths have > 1 penalty
ZSZ <- qr.qty(qrcp,sml[[i]]$S[[l]])[(pj+1):k,]
sml[[i]]$Sp[[l]]<-t(qr.qty(qrcp,t(ZSZ))[(pj+1):k,]) ## Z'SZ
}
}
}
} else qrcp <- NULL ## rest of Cp processing is after C processing
......@@ -2955,7 +2968,7 @@ smoothCon <- function(object,data,knots,absorb.cons=FALSE,scale.penalty=TRUE,n=n
} ## end smooth list loop
}
## finish of treatment of case where prediction constraints are different
## finish off treatment of case where prediction constraints are different
if (!is.null(qrcp)) {
for (i in 1:length(sml)) { ## loop through smooth list
attr(sml[[i]],"qrc") <- qrcp
......
** denotes quite substantial/important changes
*** denotes really big changes
1.7-13
** The Lanczos routine in mat.c was using a stupidly inefficient check for
convergence of the largest magnitude eigenvectors. This resulted in
far too many Lanczos steps being used in setting up thin plate regression
splines, and a noticeable speed penalty. This is now fixed, with many thanks
David Shavlik for reporting the slow down.
* Namespace modified to import from methods. Dependency on stats and graphics
made explicit.
* "re" smooths are no longer subject to side constraint under nesting (since
this is almost alway un-necessary and undesirable, and often unexpected).
* side.con modified to allow smooths to be excluded and to allow side
constraint computation to take account of penalties (unused at present).
1.7-12
* bam can now compute the leading order QR decomposition on a cluster
set up using the parallel package.
* Default k for "tp" and "ds" modified so that it doesn't exceed the 100 +
* Default k for "tp" and "ds" modified so that it doesn't exceed 100 +
the null space dimension (to avoid complaints from users smoothing in
quite alot of dimensions). Also defualt sub-sample size reduced to 2000.
quite alot of dimensions). Also default sub-sample size reduced to 2000.
* Greater use of BLAS routines in the underlying method code. In particular
all leading order operations count steps for gam fitting now use BLAS.
......@@ -143,7 +160,7 @@
coefficients if the original start values lead to an infinite deviance.
* Duchon spline bug fix (could fail to create model matrix if
number of data was one greateer than number of unique data).
number of data was one greater than number of unique data).
* fix so that 'main' is not ignored by plot.gam (got broken in 1.7-0
object orientation of smooth plotting)
......
......@@ -174,7 +174,7 @@ You must have more unique combinations of covariates than the model has total
parameters. (Total parameters is sum of basis dimensions plus sum of non-spline
terms less the number of spline terms).
This routine is less stable than `gam' for a the same dataset.
This routine is less stable than `gam' for the same dataset.
The negbin family is only supported for the *known theta* case.
......
......@@ -14,7 +14,7 @@
applied to models constructed using any smooths, including user defined smooths.
}
\usage{
gam.side(sm,Xp,tol=.Machine$double.eps^.5)
gam.side(sm,Xp,tol=.Machine$double.eps^.5,with.pen=FALSE)
}
%- maybe also `usage' for other objects documented here.
\arguments{
......@@ -22,6 +22,8 @@ gam.side(sm,Xp,tol=.Machine$double.eps^.5)
\code{\link{smooth.construct}}.}
\item{Xp}{The model matrix for the strictly parametric model components.}
\item{tol}{The tolerance to use when assessing linear dependence of smooths.}
\item{with.pen}{Should the computation of dependence consider the penalties or not.
Doing so will lead to fewer constraints.}
}
\details{ Models such as \code{y~s(x)+s(z)+s(x,z)} can be estimated by
\code{\link{gam}}, but require identifiability constraints to be applied, to
......
......@@ -125,6 +125,8 @@ it can be used and plotted, and \code{2} is it can be used but not plotted. Set
\item{plot.me}{Set to \code{FALSE} if this smooth should not be plotted by \code{\link{plot.gam}}. Set to \code{TRUE} if \code{NULL}.}
\item{side.constrain}{Set to \code{FALSE} to ensure that the smooth is never subject to side constraints as a result of nesting. }
\item{L}{smooths may depend on fewer `underlying' smoothing parameters than there are elements of
\code{S}. In this case \code{L} is the matrix mapping the vector of underlying log smoothing
parameters to the vector of logs of the smoothing parameters actually multiplying the \code{S[[i]]}.
......@@ -161,7 +163,10 @@ covariate combinations; i.e. that the basis is the same whether generated from
the full set of covariates, or just the unique combinations of covariates.
Plotting of smooths is handled by plot methods for smooth objects. A default \code{mgcv.smooth} method
is used if there is no more specific method available. Plot methods can
is used if there is no more specific method available. Plot methods can be added for specific smooth classes, see
source code for \code{mgcv:::plot.sos.smooth}, \code{mgcv:::plot.random.effect}, \code{mgcv:::plot.mgcv.smooth}
for example code.
}
......@@ -178,7 +183,7 @@ Ruppert, D., M.P. Wand and R.J. Carroll (2003) Semiparametric Regression. Cambri
University Press.
However if you want p-splines, rather than splines with derivative based penalties,
then the built in "ps" class is probably a better bet. It's based on
then the built in "ps" class is probably a marginally better bet. It's based on
Eilers, P.H.C. and B.D. Marx (1996) Flexible Smoothing with B-splines and Penalties.
Statistical Science, 11(2):89-121
......
......@@ -47,8 +47,9 @@ while \code{s(x,v,w,bs="re")} would result in \code{model.matrix(~x:v:w-1)} be
matrix. In both cases the corresponding model coefficients are assumed i.i.d. normal, and are hence
subject to ridge penalties.
Note that smooth \code{id}s are not supported for random effect terms, and that \code{"re"} terms are
not suitable for use as marginal smooths in a tensor product smooth.
Note that smooth \code{id}s are not supported for random effect terms. Unlike most smooth terms, side
conditions are never applied to random effect terms in the event of nesting (since they are identifiable
without side conditions).
Random effects implemented in this way do not exploit the sparse structure of many random effects, and
may therefore be relatively inefficient for models with large numbers of random effects, when \code{gamm4}
......
......@@ -729,7 +729,7 @@ void mgcv_trisymeig(double *d,double *g,double *v,int *n,int getvec,int descendi
void Rlanczos(double *A,double *U,double *D,int *n, int *m, int *lm,double *tol) {
/* Prototype faster lanczos_spd for calling from R.
/* Faster version of lanczos_spd for calling from R.
A is n by n symmetric matrix. Let k = m + max(0,lm).
U is n by k and D is a k-vector.
m is the number of upper eigenvalues required and lm the number of lower.
......@@ -738,18 +738,16 @@ void Rlanczos(double *A,double *U,double *D,int *n, int *m, int *lm,double *tol)
Matrices are stored in R (and LAPACK) format (1 column after another).
ISSUE: 1. eps_stop tolerance is set *very* tight.
2. Currently all eigenvectors of Tj are found, although only the next unconverged one
ISSUE: 1. Currently all eigenvectors of Tj are found, although only the next unconverged one
is really needed. Might be better to be more selective using dstein from LAPACK.
3. Basing whole thing on dstevx might be faster
4. Is random start vector really best? convergence seems very slow. Might be better to
use e.g. a sine wave, and simply change its frequency if it seems to be in null space.
Demmel (1997) suggests using a random vector, to avoid any chance of orthogonality with
an eigenvector!
2. Basing whole thing on dstevx might be faster
3. Is random start vector really best? Actually Demmel (1997) suggests using a random vector,
to avoid any chance of orthogonality with an eigenvector!
4. Could use selective orthogonalization, but cost of full orth is only 2nj, while n^2 of method is
unavoidable, so probably not worth it.
*/
int biggest=0,f_check,i,k,kk,ok,l,j,vlength=0,neg_conv,pos_conv,ni,pi,neg_closed,pos_closed,converged,incx=1;
double **q,*v=NULL,bt,xx,yy,*a,*b,*d,*g,*z,*err,*p0,*p1,*zp,*qp,normTj,eps_stop=DOUBLE_EPS,max_err,alpha=1.0,beta=0.0;
int biggest=0,f_check,i,k,kk,ok,l,j,vlength=0,ni,pi,converged,incx=1;
double **q,*v=NULL,bt,xx,yy,*a,*b,*d,*g,*z,*err,*p0,*p1,*zp,*qp,normTj,eps_stop,max_err,alpha=1.0,beta=0.0;
unsigned long jran=1,ia=106,ic=1283,im=6075; /* simple RNG constants */
const char uplo='U';
eps_stop = *tol;
......@@ -806,8 +804,6 @@ void Rlanczos(double *A,double *U,double *D,int *n, int *m, int *lm,double *tol)
/* Now stabilize by full re-orthogonalization.... */
for (i=0;i<=j;i++)
{ /* form xx= z'q[i] */
/*for (xx=0.0,qp=q[i],p0=qp + *n,zp=z;qp<p0;zp++,qp++) xx += *zp * *qp;*/
......@@ -867,34 +863,23 @@ void Rlanczos(double *A,double *U,double *D,int *n, int *m, int *lm,double *tol)
if (j >= *m + *lm)
{ max_err=normTj*eps_stop;
if (biggest) { /* getting m largest magnitude eigen values */
/* Finished only when the smallest element of the positive converged set and the
smallest element of the negative converged set are both not in the largest magnitude
set, or these sets are finished, and the largest magnitude set is of size *m, or when the total number
converged equals the matrix dimension */
pos_closed=neg_closed=0;
neg_conv=0;i=j; while (i>=0&&err[i]<max_err&&d[i]<0.0) { neg_conv++;i--;}
if (i>=0&&err[i]<max_err) neg_closed=1; /* all negatives found */
pos_conv=0;i=0; while (i<=j&&err[i]<max_err&&d[i]>=0.0) { i++;pos_conv++;}
if (i<=j&&err[i]<max_err) pos_closed=1; /* all positives found */
if (neg_conv+pos_conv >= *m) { /* some chance of having finished */
pi=0;ni=0; /* counters for how many of neg and pos converged to include in largest set */
while (pi+ni < *m) {
if ((d[pi] > -d[j-ni]&&pi<pos_conv)||ni>=neg_conv) pi++; else ni++;
/* only one convergence test is sane here:
1. Find the *m largest magnitude elements of d. (*lm is 0)
2. When all these have converged, we are done.
*/
pi=ni=0;converged=1;
while (pi+ni < *m) if (fabs(d[pi])>= fabs(d[j-ni])) { /* include d[pi] in largest set */
if (err[pi]>max_err) {converged=0;break;} else pi++;
} else { /* include d[j-ni] in largest set */
if (err[ni]>max_err) {converged=0;break;} else ni++;
}
/* now pi and ni are the number of terms to include in the largest m
from the +ve and -ve converged sets */
converged=1;
if (!pos_closed&&pi==pos_conv) converged=0; /* don't know that there is not a larger value to come */
if (!neg_closed&&ni==neg_conv) converged=0; /* ditto */
} else converged = 0;
if (converged) {
*m = pi;
*lm = ni;
j++;break;
}
} else
} else /* number of largest and smallest supplied */
{ ok=1;
for (i=0;i < *m;i++) if (err[i]>max_err) ok=0;
for (i=j;i > j - *lm;i--) if (err[i]>max_err) ok=0;
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment