cuSOLVER ライブラリを使用した NVIDIA GPU での QR 分解

ライブスクリプトを開く

この例では、CUDA Solver (cuSOLVER) ライブラリを使用するスタンドアロンの CUDA® 実行可能ファイルを作成する方法を示します。この例では、路上での自動車線トラッキングを再現する曲線近似アプリケーションを使用して、以下の点について説明します。

行列の QR 分解を使用して、任意次数の多項式でノイズデータへの近似を行う。
coder.LAPACKCallbackクラスを使用して、スタンドアロンの実行可能ファイルを生成する際に、コードジェネレーターに LAPACK ライブラリ情報を提供する。

必要条件

CUDA 対応 NVIDIA® GPU。
NVIDIA CUDA Toolkit およびドライバー。
実行環境に最適化された LAPACK ライブラリ。詳細については、LAPACK ベンダーの実装を参照してください。この例では、matlabroot/extern/lib に用意されている MATLAB® の mwlapack ライブラリを使用します。
コンパイラおよびライブラリの環境変数。サポートされているコンパイラおよびライブラリのバージョンの詳細については、サードパーティハードウェアを参照してください。環境変数の設定は、前提条件となる製品の設定を参照してください。

GPU 環境の検証

この例を実行するのに必要なコンパイラおよびライブラリが正しく設定されていることを検証するために、関数coder.checkGpuInstallを使用します。

envCfg = coder.gpuEnvConfig('host');
envCfg.BasicCodegen = 1;
envCfg.Quiet = 1;
coder.checkGpuInstall(envCfg);

行列因数分解を使用して線形システムを解く

曲線近似を利用する場合、その目的は低次多項式の係数を推定することです。この多項式は、この例では車両の前方の道路の車線境界線を表す観察されたノイズデータのモデルとして使用されます。たとえば、2 次多項式を使用する場合は、3 つの係数 ( $a$ 、 $b$ および $c$ ) を推定します。

$a x^{2} + b x + c$

最も近似する多項式が、その多項式とノイズデータの間の 2 乗誤差和を最小にする多項式として定義されます。この最小二乗問題を解くには、過決定線形システムを取得して解きます。明示的な逆行列は、システムを解くのに必要ありません。

この例では、多項式内の各項の係数が不明です。モデルとして使用する多項式は、常に路上の現在の位置を起点とするため、多項式内の定数項はゼロと仮定されます。線形項および高次項の係数を推定します。行列方程式 Ax=y を以下のように設定します。

$y$ はセンサー出力を含む。
$x$ は取得する多項式の係数を含む。
$A$ は多項式の次数とセンサーの位置に関連した定数行列。

以下のように $A$ の "QR" 分解を使用して方程式を解きます。

$A x = Q R x = y$

と

$x = p i n v (A) * y = R^{- 1} Q^{T} * y$

ここで、pinv() は疑似逆行列を表します。行列 $A$ については、次のコードを使用してこの行列方程式の解を実装できます。 $A$ を分解することによって、系の求解がより簡単になります。

[Q,R,P] = qr(A);
z = Q' * y;
x = R \ z;
yhat = A * x;

$A$ の "QR" 分解を使用して方程式を解くには、関数 linsolveQR を使用します。

type linsolveQR.m

function [yhat,x] = linsolveQR(A,y)
%#codegen

%   Copyright 2019 The MathWorks, Inc.

[Q,R,P] = qr(A);
z = Q' * y;
x = R \ z;
yhat = A * x;

end

道路用の信号モデル

アルゴリズムをテストするには、連続した曲線道路のモデル、つまり、加法性ノイズによって歪んだ正弦波を使用します。モデルの正弦波の周波数を変化させることで、アルゴリズムにさまざまな量の応力を加えます。このコードは、道路モデルを使用してノイズを含むセンサー出力をシミュレートします。

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

このコードを使用して、ヴァンデルモンド行列 $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

ヴァンデルモンド行列 $A$ とセンサー出力行列 $y$ はエントリポイント関数 linsolveQR への入力パラメーターとして渡されます。これらの入力はコンマ区切りファイルに書き込まれ、カスタムのメインファイル qrmain.cu から読み取られます。

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

スタンドアロンのコード生成用のカスタムコールバッククラス

関数 qr は cuSOLVER ライブラリで部分的にのみサポートされています。このような場合、GPU Coder™ は特定の線形代数の関数呼び出しに LAPACK ライブラリを使用します。LAPACK は、数値の線形代数のための外部ソフトウェアライブラリです。MEX ターゲットの場合、コードジェネレーターは MATLAB に付属の LAPACK ライブラリを使用します。

スタンドアロンのターゲットの場合、生成コードでの線形代数呼び出しのため、使用するヘッダーファイルと共に LAPACK ライブラリを指定するカスタムのcoder.LAPACKCallbackクラスを定義しなければなりません。この例では、lapackCallback コールバッククラスが updateBuildInfo メソッドでこれらのライブラリへのパスを指定します。ライブラリ名とコンピューター上のカスタム LAPACK インストールのパスに合わせて、このファイルを編集しなければなりません。

type lapackCallback.m

classdef lapackCallback < coder.LAPACKCallback
%

%   Copyright 2019 The MathWorks, Inc.

    methods (Static)
        function hn = getHeaderFilename()
            hn = 'lapacke.h';
        end

        function updateBuildInfo(buildInfo, buildctx)
            [~, libExt] = buildctx.getStdLibInfo();          
          
            % Specify path to LAPACK library
            if ispc
                lapackLocation = [matlabroot,'\extern'];
                libName = ['libmwlapack' libExt];
                buildInfo.addIncludePaths([lapackLocation,'\include']);            
                libPath = [lapackLocation,'\lib\win64\microsoft\'];
            else
                lapackLocation = [matlabroot];
                libName = ['libmwlapack' libExt];
                buildInfo.addIncludePaths([lapackLocation,'/extern/include']);            
                libPath = [lapackLocation,'/bin/glnxa64'];
            end
            
            % Add include path and LAPACK library for linking
            buildInfo.addLinkObjects(libName, libPath, 1000, true, true);

            buildInfo.addDefines('HAVE_LAPACK_CONFIG_H');
            buildInfo.addDefines('LAPACK_COMPLEX_STRUCTURE');
        end
    end
end

スタンドアロンコードの生成

コード構成オブジェクトで CustomLAPACKCallback プロパティを指定して、また、カスタムのメインファイル 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

Code generation successful: View report

スタンドアロンコードの実行

生成されたスタンドアロン実行可能ファイルを実行すると、出力 $y h a t$ および $x$ が計算されてコンマ区切りのファイルに書き込まれます。これらの出力を MATLAB に読み戻し、関数plotを使用してセンサーデータと近似曲線を可視化します。

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');