function [F,P,W,G,Q] = F_CGM(method,X,W_km1,G_km1,Q_km1)
%function to compute the "next" iternation of the CGM filter matrix
%
% method - text identifier for the selected CGM technique
% valid choices are:
% SD - steepest descent
% FR Fletcher-Reeves
% PR Polak-Ribiere
% HS Hestenes-Stiefel
% BB Barzilai-Borwein
% X - NxN sensitivity matrix for forward problem
% W_km1 - NxN matrix of "w_(k-1): vectors stored columnwise
% G_km1 - NxN matrix of "g_(k-1): vectors stored columnwise
% Q_km1 - NxN matrix of "q_(k-1): vectors stored columnwise
% for first iternation, W, G, Q, are zero matricies
%
% F - NxN filter matrix for the current ("k") iteration
% P - NxN P = F*X matrix
% W,G,Q - "k-1" matrices for the next iteration
%
N = length( X );
Ymat = eye( N,N ); % all of them have same square size
if( det( W_km1 ) == 0 ) % need to know if this is the first iteration
it = 1;
else
it = 2;
end
% F = zeros( N, N);
% P = F;
% W = F;
% G = F;
% Q = F;
% Detailed explanation goes here
for iy = 1:N
%
% load the current column for the main loop
%
Y = Ymat(:,iy);
q = Q_km1(:,iy);
w_km1 = W_km1(:,iy);
g_km1 = G_km1(:,iy);
q_km1 = Q_km1(:,iy); % Barzali and Borwein need this
% compute the gradient - same for all methods
g_k=-2*X'*(Y-X*q);
mag_g = g_k'*g_k;
if(it == 1)
w_k = -g_k;
else
switch method
case 'SD'
w_k = -g_k;
case 'BB'
w_k = -g_k;
case 'FR'
gamma = g_k'*g_k/(g_km1'*g_km1);
w_k = -g_k + gamma * w_km1;
case 'PR'
gamma = (g_k-g_km1)'*g_k/(g_km1'*g_km1);
w_k = -g_k + gamma * w_km1;
case 'HS'
gamma = (g_k-g_km1)'*g_k/((g_k-g_km1)'*w_km1);
w_k = -g_k + gamma * w_km1;
end
end
% compute the optimal step size
if( method ~= 'BB' | it <= 1 ) % use on first iteration always
alphasd=-w_k'*g_k/2/(w_k'*X'*X*w_k);
else
s_km1 = q - q_km1;
alphasd = s_km1'*s_km1/(s_km1'*(g_k-g_km1));
end
% must save the last iteration information, too!
W(:,iy) = w_k;
G(:,iy) = g_k;
Q(:,iy) = q; % Barzali and Borwein need this
% update the q vector
q=q+alphasd*w_k;
% reflect the column vector into the rows of F
% this is due to the persymmetry of the F matrix
for i = 1:N
f(i) = q(N-i+1);
end
% compute the P filter also
p = f*X;
% p = X*f';
frow = f;
prow = p;
F(N-iy+1,:) = frow; % reflection to bottom of F due to persymmetry
P(N-iy+1,:) = prow;
Q(:,iy) = q; % save this, too
end % end of the iy iterations
end