#################################################################
## definition of CGEMbisectionLogScaleSearch ####################
#################################################################
CGEMbisectionLogScaleSearch <- function(y,w,nu,
	candidateRanges.LBound,candidateRanges.UBound,tol)
{
n=length(y)
n1=round(sqrt(n))
bEV<-sum(y**2)/n -1.
#
#############################################
CGEM.eval <- function(y,w, coefProvFory, coefProvForw, rangeCand)
   {
      ##########################
      # la fonction necessaire au conjugate.gradient
      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
      }
      ##########################
   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)
   # coefProvFory <- 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)
   #
   coefProvForw <- out3$x  # that is (I-A).w
   cgem.trace.term <- sum( w * coefProvForw) /n
   CGEMvalue <- bEV*(cgem.quadratic.term + cgem.trace.term) 
   list(value = CGEMvalue,  niterForY=out2$conv$niter,   
                 niterForW=out3$conv$niter, 
                 coefForY=coefProvFory, coefForW= coefProvForw)
   }
#
## end of CGEM.eval ####################
#
niterBisectionMax <- 20
niterCGiterationsHistory<-matrix(0., niterBisectionMax, 2)
#
x1 <- candidateRanges.LBound
x2 <- candidateRanges.UBound
out1 <- CGEM.eval(y,w,  y,  w, candidateRanges.LBound)
out2 <- CGEM.eval(y,w,  y,  w, candidateRanges.UBound)
startForY  <- out2$coefForY
startForW  <- out2$coefForW
#      
f1 <- out1$value  - bEV
f2 <- out2$value  - bEV
#listSuccessiveValues<-c(x1,f1,x2,f2)
if (f1 > f2) 
        stop(" f1 must be < f2 ")
        #
        for (k in 1:niterBisectionMax) {
        xm <- sqrt(x1 * x2)     ################ instead of the mean
        outm <- CGEM.eval(y,w,  startForY,  startForW, xm)
        startForY  <- outm$coefForY
        startForW  <- outm$coefForW
        niterCGiterationsHistory[k,1] <- outm$niterForY 
        niterCGiterationsHistory[k,2] <- outm$niterForW 
        fm <- outm$value  - bEV
        if (fm < 0) {
            x1 <- xm
            f1 <- fm
        }
        else {
            x2 <- xm
            f2 <- fm
        }
        if (abs(1- x2/x1) < tol) {
            break
        }
       #listSuccessiveValues<-append(listSuccessiveValues,xm,fm) 
    }
    xm <- sqrt(x1 * x2)         ################ instead of the mean
#
list(root=xm, niterCGiterationsHistory = niterCGiterationsHistory)
}
## end of CGEMbisectionLogScaleSearch ####################
##########################################################