solve_SD.m
View Options

			function [X, varargout] = solve_SD(A, b, varargin)
			[Ni, Nj] = size(A);
			assert(Ni == Nj \|\| Ni == 0, 'Matrix must be square or non-zero')
			X = rand(Ni,1);
			X = X / norm(X);
			r = b - A * X;

			if length(varargin) > 0
			limit = varargin{1};
			else
			limit = 1000000;
			end

			X0 = X(:); % force copy
			solution = A \ b;
			errors = zeros(limit,1);
			k = 1;
			xnorm = norm(solution);

			while norm(r) > eps && k < limit
			errors(k) = err_A(A, b, X, X0, solution);
			alpha = (r' * r) / (r' * A * r);
			X = X + alpha * r;
			r = b - A * X;
			k = k + 1;
			end

			errors(k) = err_A(A, b, X, X0, solution);
			errors = errors(1:k);

			if length(varargin) > 0
			varargout{1} = errors;
			end
			end

			function err = err_A(A, b, Xi, X0, X)
			err = sqrt(abs(diff_A(A, b, Xi, X) / diff_A(A, b, X0, X)));
			end

			function err = diff_A(A, b, Xi, X)
			err = 0.5 * ( Xi' * A * Xi - X' * A * X ) - b' * (Xi - X);
			end

compphys/exo4/solve_SD.m33bba04159d6master