%% Implement Hardware-Efficient Complex Burst QR Decomposition
% This example demonstrates how to compute the QR decomposition of complex
% matrices using hardware-efficient MATLAB(R) code in Simulink(R). This
% model shares computational resources across steps of the QR Decomposition
% algorithm. It thus uses fewer on chip resources than a fully pipelined
% approach, while sacrificing some total throughput.

% Copyright 2018-2019 The MathWorks, Inc.
%% Solving Matrix Equations with the QR Decomposition
%
% When solving matrix equations, it is seldom, if ever, necessary to
% compute the inverse of a matrix [1][3]. The Complex Burst QR
% Decomposition block provides a hardware efficient way of solving the
% equation
%
% $$ Ax = B, $$
%
% where $A$ is a complex m x n matrix, $x$ is a complex n x p matrix, and
% $B$ is a complex m x p matrix. This equation can be converted to the form
%
% $$ Rx = Q^H B $$
%
% through a series of orthogonal transformations. Here, $R$ is a complex m
% x n upper triangular matrix. $Q$ is a complex m x m orthogonal matrix
% such that $QR = A$. The Complex Burst QR Decomposition retains only the
% non-zero rows of $R$, along with the corresponding rows of  $Q^H B$.
%
% The QR Decomposition is well suited to fixed-point architectures because
% it can be entirely performed with Givens rotations. These have an
% efficient fixed-point implementation in terms of the CORDIC Algorithm .
% For more details on the algorithm for fixed-point QR decomposition, see
% <matlab:showdemo('cordicqr_demo'); Perform QR Factorization Using
% CORDIC>. The Simulink model used in this example is:

fxpdemo_complex_burst_QR_model
%%
% This model works with fixed-point, double, single, and scaled double data
% types.
%
%% Input Parameters for Complex Burst QR Decomposition
%
% The Complex Burst QR Decomposition block takes four input parameters. The
% mask for this block is shown below:
%
% <<complex_burst_qr_mask.png>>
%
% You must assign input data to variables in the MATLAB workspace. To
% customize this algorithm, change the data according to the procedures
% used in the following sections.
%
%% Define Matrix Dimensions
% The model's architecture allocates the minimal memory necessary to store
% the data needed to perform the QR decomposition. Therefore, the size of
% the input matrices must be known at compile time. Here, |m| is the number
% of rows in the complex matrices |A| and |B|, |n| is the number of columns
% in |A|, and |p| is the number of columns in |B|.
%
% Additionally, since the block can handle decomposition of an indefinite
% number of matrices in series, a number of samples is specified. The model
% terminates when all of the samples have been used.
m = 4;
n = 4;
p = 4;
num_samples = 100;
%% Define Input Datatypes
% Before generating input data, it is important to specify the data type of
% the matrix data as shown below. This example uses signed fixed-point
% types with 16 bit word lengths and fractional types. For a complex 4 x 4
% matrix, the output type needs an additional 3 bits to accommodate the
% data growth in the integer part of the input data; see
% <matlab:showdemo('fi_complex_mldivide_systolic_array_example'); Use
% Hardware-Efficient Algorithm to Solve Systems of Complex-Valued Linear
% Equations> for further details.
word_length = 16;
qr_growth_bits = 3;
fraction_length = word_length - 1;
nt = numerictype(1,word_length + qr_growth_bits,fraction_length)
%% Define the Number of CORDIC Iterations
% The CORDIC approximation for Givens rotations is an iterative algorithm
% that gives an additional bit of accuracy per iteration. For signed
% fixed-point datatypes, the greatest possible accuracy is achieved when
% the number of iterations performed is one fewer than the wordlength.
NumberOfCORDICIterations = nt.WordLength - 1
%% Initialize Random Data
% The data handler in this example takes complex matrices |A| and |B| as
% inputs. In this example, |A| and |B| are defined as random matrices
% elements drawn from a uniform distribution from -1 to 1. Note that |A|
% and |B| are each defined as three-dimensional arrays, where each sample
% matrix is stored in the first two dimensions.
A = fi(complex(2*rand(m,n,num_samples) - 1, 2*rand(m,n,num_samples) - 1),nt);
B = fi(complex(2*rand(m,p,num_samples) - 1, 2*rand(m,p,num_samples) - 1),nt);

%% Transferring Data to Complex Burst QR Decomposition
% The |ready| port triggers the Data Handler subsystem. When |ready| is
% high, the block asserts |validIn| and sends the next row of |A| and |B|
% to |aIn| and |bIn|. The protocol for transferring data allows data to be
% sent whenever |ready| is high, ensuring all data is processed. If data is
% sent when |ready| is not high, it will not be processed.
%% Handling Outputs
% Complex Burst QR Decomposition outputs data one row at a time. Whenever
% the block outputs a result row, it asserts |validOut|. Note that |rOut|
% outputs the rows of $R$, while |cOut| outputs the rows of $Q^HB$.  The
% rows are output in reverse order, as this is the natural order to consume
% them for back substitution. In this example, they are saved as |(m *
% num_samples)| x |n| matrices, with the rows ordered as they were received.
%% Simulate
sim fxpdemo_complex_burst_QR_model
%% Reconstruct the Solutions from the Output Data
% Because the data are output in reverse order, you must reconstruct the
% data to interpret the results. The following code puts the data in the
% correct order.

C = cell(1,num_samples);
R = cell(1,num_samples);
for ii = 1:num_samples
    C{ii} = flipud(C_Out((ii-1)*m + 1:ii*m,:));
    R{ii} = flipud(R_Out((ii-1)*m + 1:ii*m,:));
end
%% Evaluate the Accuracy of the Result
% To evaluate the accuracy of the Burst QR Decomposition, examine the
% magnitude of the difference between the solution to the matrix equation
% using the 'Complex QR Burst Decomposition' block, and that obtained using
% MATLAB(R)'s built-in backsolve for doubles. Plot the absolute error and
% condition number for each sample. The data show that the accuracy of the
% solution tracks the condition number, as expected [2].
xAct = cell(1,num_samples);
xExp = cell(1,num_samples);
xErr = zeros(1,num_samples);
condNumber = zeros(1,num_samples);
for ii = 1:num_samples
    xAct{ii} = double(R{ii})\double(C{ii});
    xExp{ii} = double(A(:,:,ii))\double(B(:,:,ii));
    xErr(ii) = norm(xAct{ii} - xExp{ii});
    condNumber(ii) = cond(double(A(:,:,ii)));
end
figure(1)
clf
h1 = subplot(2,1,1);
hold on;
h1.YScale = 'log';
plot(xErr)
grid on
title('Error of ''Complex Burst QR Decomposition''');
ylabel('norm(R\\C - A\\B)');
h2 = subplot(2,1,2);
hold on;
h2.YScale = 'log';
plot(condNumber);
grid on
title('Condition Number of Samples');
ylabel('cond(A)');
xlabel('Test Point');
linkaxes([h1,h2],'x');

%% Close the model

close_system fxpdemo_complex_burst_QR_model

%% References
%
% [1] George E. Forsythe, M.A. Malcom and Cleve B. Moler. Computer Methods
% for Mathematical Computations. Englewood Cliffs, N.J.: Prentice-Hall,
% 1977.
%
% [2] Cleve B. Moler. Cleve's Corner:
% <matlab:web('https://blogs.mathworks.com/cleve/2017/07/17/what-is-the-condition-number-of-a-matrix');
% What is the Condition Number of a Matrix?>, The MathWorks, Inc. 2017.
%
% [3] Cleve B. Moler. Numerical Computing with MATLAB. SIAM, 2004. isbn:
% 978-0-898716-60-3.
%

%#ok<*NASGU,*NOPTS>
