%%
%% Listing 7.2 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
%% (C) Piotr Krzyżanowski, 2011
%%
%% Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
%% zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
%%

function x = solvegr(e, W, D, b, tol, maxit)
%
% Rozwiazuje (I - e*W*D^{-1})x = (1-e)b
% 	0 <= e <= 1, 
% 	W - symetryczna
% 	D - diagonalna, dod. okresl.

if nargin < 5
	tol = 1e-8;
end
if nargin < 6
	maxit = 100;
end

switch e
case 0
	x = b;
case 1
	opts.maxit = maxit; opts.tol = tol; k = 16;
	[V, EV, info] = eigs(W, D, k, 'lm', opts);
	EV = diag(EV);
	nof1 = sum(abs(EV-1) < 1e-6); % liczba wartości własnych bliskich "1"
	if nof1 > 1
		warning(['Wykryto ', num2str(nof1), ' wartosci wlasnych bliskich 1']);
		disp(sort(EV, 'descend'));
	end
	[m, i] = max(EV);
	y = V(:,i);
	x = D*y;
	x = x/sum(x);
otherwise % 0 < e < 1
	DW = @(x) (D*x - e*W*x);
	b = (1-e)*b;
	if exist('OCTAVE_VERSION')
	% Octave
		[y, flag, relres, iter, resid] = pcr(DW, b,  tol, maxit, D);
		fprintf(stderr, 'PCR: \n');
	else
	% MATLAB
		[y, flag, relres, iter, resid] = minres(DW, b, tol, maxit, D);
		stderr = 2; % deskryptor stderr w MATLABie
		fprintf(stderr, 'MINRES: \n');
	end
	semilogy(abs(resid));
	fprintf(stderr, '\tKod zakonczenia: %d\n\tRelRes: %e\n\tIteracji: %d\n', flag, relres, prod(iter));
	x = D*y;
end
end
