%% Using Function Handle Instead of Numeric Matrix
% Solve a linear system by providing bicg with a function handle that computes A*x and A'*x in place of the coefficient matrix A.
%
% Create an nonsymmetric tridiagonal matrix. Preview the matrix.

A = gallery('wilk',21) + diag(ones(20,1),1)

%
% Since this tridiagonal matrix has a special structure, you can represent the operation A*x with a function handle. When A multiplies a vector, most of the elements in the resulting vector are zeros. The nonzero elements in the result correspond with the nonzero tridiagonal elements of A.
% The expression A\;x becomes:
% A\;x=\left\lbrack \begin{array}{ccccccc}
% 10 & 2 & 0 & \cdots  &  & \cdots  & 0\\
% 1 & 9 & 2 & 0 &  &  & \vdots \\
% 0 & 1 & \ddots  & 2 & 0 &  & \\
% \vdots  & 0 & 1 & 0 & \ddots  & \ddots  & \vdots \\
%  &  & 0 & \ddots  & 1 & \ddots  & 0\\
% \vdots  &  &  & \ddots  & \ddots  & \ddots  & 2\\
% 0 & \cdots  &  & \cdots  & 0 & 1 & 10
% \end{array}\right\rbrack \left\lbrack \begin{array}{c}
% x_1 \\
% x_2 \\
% x_3 \\
% \vdots \\
% \vdots \\
% x_{21}
% \end{array}\right\rbrack =\left\lbrack \begin{array}{c}
% {10x}_1 +2x_2 \\
% x_1 +9x_2 +2x_3 \\
% \vdots \\
% \vdots \\
% x_{19} +9x_{20} +2x_{21} \\
% x_{20} +10x_{21}
% \end{array}\right\rbrack.
% The resulting vector can be written as the sum of three vectors:
% A\;x=\left\lbrack \begin{array}{c}
% {10x}_1 +2x_2 \\
% x_1 +9x_2 +2x_3 \\
% \vdots \\
% \vdots \\
% x_{19} +9x_{20} +2x_{21} \\
% x_{20} +10x_{21}
% \end{array}\right\rbrack=\left\lbrack \begin{array}{c}
% 0\\
% x_1 \\
% x_2 \\
% \vdots \\
% x_{20}
% \end{array}\right\rbrack +\left\lbrack \begin{array}{c}
% {10x}_1 \\
% {9x}_2 \\
% \vdots \\
% 9x_{20} \\
% 10x_{21}
% \end{array}\right\rbrack +2\cdot \left\lbrack \begin{array}{c}
% x_2 \\
% x_3 \\
% \vdots \\
% x_{21} \\
% 0
% \end{array}\right\rbrack.
% Likewise, the expression for A^T \;x becomes:
% A^T \;x=\left\lbrack \begin{array}{ccccccc}
% 10 & 1 & 0 & \cdots  &  & \cdots  & 0\\
% 2 & 9 & 1 & 0 &  &  & \vdots \\
% 0 & 2 & \ddots  & 1 & 0 &  & \\
% \vdots  & 0 & 2 & 0 & \ddots  & \ddots  & \vdots \\
%  &  & 0 & \ddots  & 1 & \ddots  & 0\\
% \vdots  &  &  & \ddots  & \ddots  & \ddots  & 1\\
% 0 & \cdots  &  & \cdots  & 0 & 2 & 10
% \end{array}\right\rbrack \left\lbrack \begin{array}{c}
% x_1 \\
% x_2 \\
% x_3 \\
% \vdots \\
% \vdots \\
% x_{21}
% \end{array}\right\rbrack =\left\lbrack \begin{array}{c}
% {10x}_1 +x_2 \\
% {2x}_1 +9x_2 +x_3 \\
% \vdots \\
% \vdots \\
% {2x}_{19} +9x_{20} +x_{21} \\
% {2x}_{20} +10x_{21}
% \end{array}\right\rbrack.
% A^T \;x=\left\lbrack \begin{array}{c}
% {10x}_1 +x_2 \\
% {2x}_1 +9x_2 +x_3 \\
% \vdots \\
% \vdots \\
% {2x}_{19} +9x_{20} +x_{21} \\
% {2x}_{20} +10x_{21}
% \end{array}\right\rbrack =2\cdot \left\lbrack \begin{array}{c}
% 0\\
% x_1 \\
% x_2 \\
% \vdots \\
% x_{20}
% \end{array}\right\rbrack +\left\lbrack \begin{array}{c}
% {10x}_1 \\
% {9x}_2 \\
% \vdots \\
% 9x_{20} \\
% 10x_{21}
% \end{array}\right\rbrack +\left\lbrack \begin{array}{c}
% x_2 \\
% x_3 \\
% \vdots \\
% x_{21} \\
% 0
% \end{array}\right\rbrack.
% In MATLAB®, write a function that creates these vectors and adds them together, giving the value of A*x or A'*x, depending on the flag input:
% function y = afun(x,flag)
% if strcmp(flag,'notransp') % Compute A*x
%     y = [0; x(1:20)] ...
%         + [(10:-1:0)'; (1:10)'].*x ...
%         + 2*[x(2:end); 0];
% elseif strcmp(flag,'transp') % Compute A'*x
%     y = 2*[0; x(1:20)] ...
%         + [(10:-1:0)'; (1:10)'].*x ...
%         + [x(2:end); 0];
% end
% end
% (This function is saved as a local function at the end of the example.)
%
% Now, solve the linear system \textrm{Ax}=b by providing bicg with the function handle that calculates A*x and A'*x. Use a tolerance of 1e-6 and 25 iterations. Specify b as the row sums of A so that the true solution for x is a vector of ones.

b = full(sum(A,2));
tol = 1e-6;  
maxit = 25;
x1 = bicg(@afun,b,tol,maxit)

%
%% Local Functions

function y = afun(x,flag)
if strcmp(flag,'notransp') % Compute A*x
    y = [0; x(1:20)] ...
        + [(10:-1:0)'; (1:10)'].*x ...
        + 2*[x(2:end); 0];
elseif strcmp(flag,'transp') % Compute A'*x
    y = 2*[0; x(1:20)] ...
        + [(10:-1:0)'; (1:10)'].*x ...
        + [x(2:end); 0];
end
end

% Copyright 2019 The MathWorks, Inc.