solver for system of linear equation

[Cplusc]

Hello everyone.
I write a conjugate gradient and preconditioned conjugate gradient solver and they are working fine. I am using these solvers in finite difference method. In order to test my code I just run it for a simple mesh grid which is a 2d square with four cells. In total, There will be 9 points in mesh in which 8 of them are boundary points and are being set previously and I'm going to find the solution for center point.
The equation for this point will be
unew[i][j] = 0.25 * ( u[i - 1][j] + u[i][j + 1] + u[i][j - 1] + u[i + 1][j] + f[i][j] * dx * dy);

consider that I'm solving the Poisson equation and f will be the right hand side function in Poisson equation.

Now my problem is I don't understand why the boundary point are not set in the first iteration of my solver but for example after 1000 iteration everything is fine. I know I don't need any solver for this simple equation but I wonder if it should not be working for all cases?

[lastchance]

I don't understand why the boundary point are not set in the first iteration of my solver

Your boundary points should correspond to boundary conditions of your differential equation. They are either fixed (Dirichlet boundary conditions), in which case they will be defined before you start iterating and they will never change, or they are specified (usually-zero) gradient (Neumann boundary conditions), in which case they are updated after each iteration. Regardless, your code as written will certainly use those boundary points, so you had better make sure that they are initialised before use.

Without knowing what physical system you are trying to solve for it's impossible to comment further.

[Cplusc]

Thanks for your answer.
Let's say my A matrix is like

1 0 0 0 0 0 0 0 0
0 1 0 0 0 0 0 0 0
0 0 1 0 0 0 0 0 0 
0 1 0 1 -4 1 0 1 0
0 0 0 0 0 1 0 0 0
0 0 0 0 0 0 1 0 0
0 0 0 0 0 0 0 1 0
0 0 0 0 0 0 0 0 1

and this is when I applied the Dirichlet boundary condition.

and the x vector containing the boundary point (it is the initial guess for r=Ax-b) is like

x={0}

and the right hand side is

rhs={0.5,0.25,0.5,0.25,1,0.25,0.5,0.25,0.5}

and this is my solver

void Blas::MPI_CG(SparseMatrix& A, VecDbl_t& x, VecDbl_t& b, size_t iter_max) {

	double TOLERANCE = 1.0e-10;
	const double NEARZERO = 1.0e-10;
	int n = A.Size_;
	x.resize(n);
	VecDbl_t Ap(n);
	VecDbl_t r = b;
	VecDbl_t p = r;
	int k = 0;

	while (k < 1000)
	{
		VecDbl_t rold = r;                                        
		MxV(A, p, Ap);

		double alpha = dotproduct(r,r) / dotproduct(p, Ap);
		x = vectorCombination(1.0, x, alpha, p);            
		r = vectorCombination(1.0, r, -alpha, Ap);          

		if (vectorNorm(r) < TOLERANCE) break;       

		double beta = dotproduct(r, r) / dotproduct(rold, rold);
		p = vectorCombination(1.0, r, beta, p);             
		k++;
	}

}

The dot product and MxV are working properly. and the Dirichlet boundary condition already has been applied on A matrix.
I don't undrestand why in the first iteration the boundary are not set?

[Cplusc]

Also I used this solver http://www.cplusplus.com/forum/general/222617/ but again the same problem.