echo on
%%%%% laplace_conjug_grad.m %%%%%

disp(sprintf('**** Exampel of Conjugate-gradient iteration for Laplace eq.****\n'));

%Two dimensional 3x3 elliptic boundary problem (Laplace equation)
%with boundary values 0 (left, bottom right) and 1 up.
%Natural ordering of nodes: 1 left bottom, ..., 3 right bottom,
%4 left second row (from bottom) etc. 9 right up.
%Laplace matrix L is composed of sub-matrices T and -I where
format short;
b=[0 0 0 0 0 0 1 1 1]'; %boundary conditions (upper edge = 1)
T=diag([4 4 4],0)+ diag([-1 -1],1)+ diag([-1 -1],-1)
I=eye(3)
Z=zeros(3,3)
LA=[T -I Z; -I T -I; Z -I T;]
%All diagonal elements are the same. Matrix LA can be saved by diagonals,
%or not at all because their structure is known.
%The system dimension (number of inner points, elements of vector x) is 9.
%Matrix-vector multiplication A*s_k could be improved for Laplace eq.
%Boundary conditions are incorporated in right-hand side vector b.
%Initial approximation x0
%x0=4*[ones(1,9)]';
%%%
%This method can be used for any general positive definite and symmetric system. 
%Note: Use general multiplication for m.
%LA=diag([4 4 4 4 4 4 4 4 4],0)+ diag([-1 -1 -1 -1 -1 -1 -1 -1],1)+ ...
%   diag([-1 -1 -1 -1 -1 -1 -1 -1],-1)+ diag([-1 -1 -1 -1 -1 -1 -1],2)+ ...
%   diag([-1 -1 -1 -1 -1 -1 -1],-2);
%%%
%flops(0);
%first approximation:
x0=ones(9,1);
r0=b-LA*x0;	%residual
s0=r0;		%search direction
%
   m=LA*s0;		%matrix-vector multiplication (can be improved for banded matrices)
   %m=mult(LA,s0);	%tailored multuplication of Laplace matrix for 2-D domain
	alfa=(r0'*r0)/(s0'*m);	%compute search parameter
   x1=x0+alfa*s0;		%update solution
   r1=r0-alfa*m;	%compute new residual
   beta=(r1'*r1)/(r0'*r0);
   s1=r1+beta*s0;		%compute new search direction
%other steps
echo off;
for i=1:4
	x0=x1;
	r0=r1;
	s0=s1;
   m=LA*s0;		%matrix-vector multiplication (can be improved for banded matrices)
   %m=mult(LA,s0);	%tailored multuplication of Laplace matrix for 2-D domain
   alfa=(r0'*r0)/(s0'*m);	%compute search parameter
   x1=x0+alfa*s0;		%update solution
   r1=r0-alfa*m;	%compute new residual
   beta=(r1'*r1)/(r0'*r0);
   s1=r1+beta*s0;		%compute new search direction
   i=i+1;
end
echo on;
x1	%approximate solution after 5 steps (faster convergence than 
   %Jacobi and Gauss-Seidel and SOR)
%flops
x=LA\b	%exact solution