# The standard Conjugate Gradient implementation
# with a matrix A. No preconditioning.
# Parameters: zero initial guess; b-right_hand_side; y-the solution
# ---------------------------------------------------------------
 y <- cg(A,b,eps_cg){
 
     it <- 1  
     n <- nrow(A)  
     y = array(0, dim=c(n,1)) # creates a zero initia guess
     resvec = array(0, dim=c(n+1,1)) # convergence history
     g = -b  # g = A%*%y - b
qq0 <- t(g)%*%g;
cat("___________cg_delta0:",qq0,"\n")
ceps  <-  eps_cg*qq0
resvec[it,1] = qq0
qq1 <- 1e+12
d = -g
    while (qq1>ceps){
        h <- A%*%d
        tau <- qq0/(t(d)%*%h)
        y <- y+d%*%tau 
        g <- g+h%*%tau
        qq1 <- t(g)%*%g
        cat("______________________________cg_delta1:",qq1,"\n")
        beta <- qq1/qq0
        d <- -g+d%*%beta
        it <- it+1
        qq0 <- qq1
	resvec[it,1] = qq0
	y
    }
it <- it-1
cat("___________ Total number of iterations:",it,"\n")
cat("___________ Plot of the resudual curve:","\n")
plot(1:it+1,resvec[1:it+1,1],type="pl",ylog=TRUE)
return(y)
# ---------------------------------------------------------------
}
#	p1=proc.time()
#	Sys.sleep(2)
#	proc.time()-p1