%% QR Decomposition on NVIDIA GPU Using cuSOLVER Libraries
% This example shows how to create a standalone CUDA(R) executable 
% that leverages the CUDA Solver library (cuSOLVER). The example uses a 
% curve fitting application that mimics automatic lane tracking on a road 
% to illustrate:
%
% * Fitting an arbitrary-order polynomial to noisy data by using matrix QR
% factorization.
% * Using the <docid:coder_ref#bu7pod8-1 coder.LAPACKCallback> class to 
% provide the LAPACK library information for the code generator when 
% generating standalone executables.

% Copyright 2019 The MathWorks, Inc.

%% Prerequisites
% * CUDA enabled NVIDIA(R) GPU with compute capability 3.2 or higher.
% * NVIDIA CUDA toolkit and driver.
% * LAPACK library that is optimized for your execution environment. For more information, see 
% <https://www.netlib.org/lapack/faq.html#_what_and_where_are_the_lapack_vendors_implementations
% LAPACK vendors implementations>. This example uses the |mwlapack|
% libraries that MATLAB(R) provides in _matlabroot/extern/lib_.
% * Environment variables for the compilers and libraries. For information 
% on the supported versions of the compilers and libraries, see 
% <docid:gpucoder_gs#mw_aa8b0a39-45ea-4295-b244-52d6e6907bff
% Third-Party Products>. For setting up the environment variables, see 
% <docid:gpucoder_gs#mw_453fbbd7-fd08-44a8-9113-a132ed383275
% Environment Variables>.

%% Create a Folder and Copy Relevant Files
% The following line of code creates a folder in your current working 
% folder (pwd) and copies all the relevant files into this folder. If you 
% do not want to perform this operation or if you cannot generate files in 
% this folder, change your current working folder.
gpucoderdemo_setup('gpucoderdemo_qr_decomposition');

%% Verify GPU Environment
% To verify that the compilers and libraries necessary for running this example
% are set up correctly, use the <docid:gpucoder_ref#mw_0957d820-192f-400a-8045-0bb746a75278 coder.checkGpuInstall>
% function.
envCfg = coder.gpuEnvConfig('host');
envCfg.BasicCodegen = 1;
envCfg.Quiet = 1;
coder.checkGpuInstall(envCfg);

%% Solve a Linear System by Using Matrix Factorization
% In curve fitting applications, the objective is to estimate the
% coefficients of a low-order polynomial. The polynomial is then used as a 
% model for observed noisy data, which in this example represents the lane 
% boundary of the road ahead of a vehicle. For example, when using a quadratic 
% polynomial, there are three coefficients ($a$, $b$, and $c$) to estimate:
% 
% $$ax^2 + bx + c$$
% 
% The polynomial that fits best is defined as the one that minimizes
% the sum of the squared errors between itself and the noisy
% data. To solve this least-squares problem, you get and solve an 
% overdetermined linear system. An explicit matrix inverse is not required
% to solve the system.
%
% In this example, the unknowns are the coefficients of each term in the 
% polynomial. Because the polynomial you use as a model always starts from 
% the current position on the road, the constant term in the polynomial is 
% assumed to be zero. Estimate the coefficients for the linear and 
% higher-order terms. Set up a matrix equation _Ax=y_ such that:
%
% * $y$ contains the sensor outputs.
% * $x$ contains the polynomial coefficients that we need to obtain.
% * $A$ is a constant matrix related to the order of the polynomial and the
% locations of the sensors.
%
% Solve the equation using the _QR_ factorization of $A$:
%
% $$Ax = QRx=y$$
%
% and
%
% $$x = pinv(A) * y = R^{-1}Q^T*y$$
% 
% where _pinv()_ represents pseudo-inverse. Given the matrix $A$, you can 
% use the following code to implement a solution of this matrix equation.
% Factoring $A$ allows for an easier solution of the system.
%
%  [Q,R,P] = qr(A);
%  z = Q' * y;
%  x = R \ z;
%  yhat = A * x;
%
% Use the <matlab:edit(fullfile(matlabroot,'toolbox','gpucoder','gpucoderdemos','gpucoderdemo_qr_decomposition','linsolveQR.m')) linsolveQR> 
% function to solve the equation using the _QR_ factorization of $A$.

type linsolveQR.m

%% Signal Model for the Road
% To test the algorithm, a continuously curving road model is
% used, that is, a sinusoid that is contaminated with additive noise. By 
% varying the frequency of the sinusoid in the model, you can stress the
% algorithm by different amounts. This code simulates noisy sensor outputs 
% using our road model:

Duration = 2;                                                       % Distance that we look ahead
N = 25;                                                             % Total number of sensors providing estimates of road boundary
Ts = Duration / N;                                                  % Sample interval
FracPeriod = 0.5;                                                   % Fraction of period of sinusoid to match
y = sin(2*pi* (0:N-1)' * (FracPeriod/N)) + sqrt(0.3) * randn(N,1);  % y will contain the simulated sensor outputs

%%
% Use this code to form the Vandermonde matrix $A$:

Npoly = 3;                  % Order of polynomial to use in fitting
v = (0:Ts:((N-1)*Ts))';
A = zeros(length(v), Npoly);
for i = Npoly : -1 : 1
    A(:,i) = v.^i;
end
 
%% 
% The Vandermonde matrix $A$ and sensor outputs matrix $y$ are passed as input
% parameters to the |linsolveQR| entry-point function. These inputs are
% written to comma-separated files and are read from the handwritten main 
% <matlab:edit(fullfile(matlabroot,'toolbox','gpucoder','gpucoderdemos','gpucoderdemo_qr_decomposition','qrmain.cu')) qrmain.cu.> 

 writematrix(reshape(A, 1, 75), 'inputA.csv');
 writematrix(reshape(y, 1, 25), 'inputY.csv');

%% Custom Callback Class for Standalone Code Generation
% The |qr| function is only partially supported in the |cuSOLVER| library.
% In such cases, GPU Coder(TM) uses the |LAPACK| library for 
% certain linear algebra function calls. |LAPACK| is an external software library 
% for numeric linear algebra. For MEX targets, the code generator uses 
% the |LAPACK| library included in MATLAB. 
%
% For standalone targets, you must define a custom 
% <docid:coder_ref#bu7pod8-1 coder.LAPACKCallback> class that 
% specifies the LAPACK libraries along with the header files 
% to use for linear algebra calls in the generated code. In this example, the 
% <matlab:edit(fullfile(matlabroot,'toolbox','gpucoder','gpucoderdemos','gpucoderdemo_qr_decomposition','lapackCallback.m')) lapackCallback>
% callback class specifies the paths to these libraries in 
% |updateBuildInfo| method. You must modify this file with the
% library names and paths for the custom LAPACK installation on your computer.

type lapackCallback.m


%% Standalone Code Generation
% Generate a standalone executable by the specifying |CustomLAPACKCallback| 
% property in the code configuration object and using a handwritten main
% |qrmain.cu|.

cfg = coder.gpuConfig('exe');
cfg.GpuConfig.EnableCUSOLVER = 1;
cfg.CustomLAPACKCallback = 'lapackCallback';
cfg.CustomSource = 'qrmain.cu';
cfg.CustomInclude = '.';
codegen -config cfg -args {A,y} linsolveQR -report

%% Standalone Code Execution
% When you execute the generated standalone executable, the outputs $yhat$ and
% $x$ are computed and written to comma-separated files. Read these outputs 
% back in MATLAB and use the <docid:matlab_ref#buggqqy-1 plot> function to 
% visualize the sensor data and fitted curve.
%
if ispc
    system('linsolveQR.exe');
else
    system('./linsolveQR');
end
yhat = reshape(readmatrix('outputYhat.csv'), 25, 1);
x = reshape(readmatrix('outputX.csv'), 3, 1);
figure
plot(v, y, 'k.', v, yhat, 'r')
axis([0 N*Ts -3 3]);
grid;
xlabel('Distance Ahead of the Vehicle');
legend('Sensor data','Curve fit');
title('Estimate the Lane Boundary Ahead of the Vehicle');

%% Cleanup
% Remove the generated files and return to the original folder.
cleanup

displayEndOfDemoMessage(mfilename)
