#############################################################
########### definition of CGEMevalOnGrid ####################
#############################################################
CGEMevalOnGrid <- function(y,w,candidateRanges.Grid,nu)
{
CGEMvaluesOnthegrid<- matrix(NA,length(candidateRanges.Grid))
trace.term<- matrix(NA,length(candidateRanges.Grid))
niter <- matrix(NA,length(candidateRanges.Grid),2)
n=length(y)
n1=round(sqrt(n))
##########################
# the matrice-vector product required by conjugate.gradient:
bEV<-sum(y**2)/n -1.
prod.bCorrelPlusI.Timesx <- function( x) {
   junk<- matrix(x, n1,n1)
   result<- stationary.image.cov(Y=junk,cov.obj=cov.obj)
   result<- bEV*result + junk
   result
}
##########################
coefProvFory<-y
coefProvForw<-w
# loop over ranges
for(k in 1:length(candidateRanges.Grid)){
   rangeCand<- candidateRanges.Grid[k]
   cov.obj<- stationary.image.cov( setup=TRUE, 
            grid=list(x=1:n1,y=1:n1), 
            range= rangeCand*n1,smoothness=nu)
   #
   out2<- conjugate.gradient(y, 
            prod.bCorrelPlusI.Timesx,start=coefProvFory,
            tol=1e-03,    kmax=200,verbose=FALSE)
   niter[k,1]<-out2$conv$niter
   # coefProv<- solve( bCorrelPlusI.mat,y)
   coefProvFory<- out2$x
   #
   # bCorrelTimesCoef <- bEV*( Correl.mat %*% coefProv)
   bCorrelTimesCoef <- (prod.bCorrelPlusI.Timesx(coefProvFory)
              - coefProvFory)
   cgem.quadratic.term <- sum(bCorrelTimesCoef *  coefProvFory) /n
   #
   out3<- conjugate.gradient(w, 
            prod.bCorrelPlusI.Timesx,start= coefProvForw,
            tol=1e-03,      kmax=200,verbose=FALSE)
   niter[k,2]<-out3$conv$niter
   #
   coefProvForw <- out3$x  # that is (I-A).w
   cgem.trace.term <- 1-sum( w * (w-coefProvForw)) /n
   CGEMvaluesOnthegrid[k] <- bEV*(cgem.quadratic.term 
            + cgem.trace.term) 
   trace.term[k] <- cgem.trace.term         
   }
list(values = CGEMvaluesOnthegrid,  
            niterForY=niter[,1],   niterForW=niter[,2],
            trace.term=trace.term)
}
########### end of CGEMevalOnGrid ###########################
#############################################################