RTK測位アルゴリズムの基礎とプログラミング

RTK-GPS測位の基礎と
プログラミング (3)
Basics of RTK-GPS Positioning and Its Programing
東京海洋大産学官連携研究員/技術コンサルタント
高須知二
Tomoji TAKASU
内容
• 行列演算ライブラリ
• BLASS
• LAPACK
• C→Fortran呼び出し法
• 可変長行列・ベクトル
• 最小二乗法解法
• 行列内部表現
http://gpspp.sakura.ne.jp
行列演算ライブラリ (1)
• メリット
実行速度が速い (数倍～数10倍)
信頼性が高い
コードを書かなくて良い
• デメリット
使うのが面倒 (使いこなしが大変)
make時にライブラリが必要(可搬性)
http://gpspp.sakura.ne.jp
行列演算ライブラリ (2)
• BLAS
基礎行列演算: 代入、内積、ノルム、
加減算、乗算、除算
• LAPACK
線型方程式(連立一次方程式)
最小二乗、固有値分解、特異値分解
http://gpspp.sakura.ne.jp
BLAS (1)
• Level1 : y←ax+y...(ベクトル間)
• Level2 : y←aAx+by...(行列-ベクトル)
• Level3 : C←aAB+bC...(行列-行列)
• 関数名 : xAAAA
(x : S=単精度実数, D=倍精度実数
C=単精度複素数, Z=倍精度複素数)
http://gpspp.sakura.ne.jp
BLAS (2)
• 主要ルーチン
DGEMM : C←aA*B*+bC
• プロトタイプ
DGEMM(char *transa, char *transb,
int* m, int *n, int *k, double *alpha,
double *a, int *lda, double *b, int *ldb,
double * beta, double *c, int *ldc);
http://gpspp.sakura.ne.jp
LAPACK (1)
• ドライバルーチン
• 個別ルーチン
• 補助ルーチン
http://gpspp.sakura.ne.jp
LAPACK (2)
• DGESV:線形方程式
y  Ax  x
• DGELS:最小二乗
min y  Ax
x
2
x
http://gpspp.sakura.ne.jp
C→Fortran呼び出し方法
• 関
数
名
置
換
GEMM→gemm_() (GCC)
• 引数: すべてポインタ(参照)渡し。
(値渡しはできない)
• リンク時オプション: ライブラリ指定
-llapack -lblas (GCC)
http://gpspp.sakura.ne.jp
可変長行列/ベクトル
• 標
{double
...;
free(a);}
• ALLOCA
C(89)
準
*a=malloc(sizeof(double)*n*m);
b=a[i+n*j];
...
(
非
標
準
)
{double
*a=alloca(sizeof(double)*n*m);
...; b=a[i+n*j]; ...}
• C99
{double
...; a[i][j]; ...} (Row-Major)
a[n][m];
http://gpspp.sakura.ne.jp
最小二乗法解法 (1)
• 正
規x
 ( A A)
T方
1
T程
A y
式
• 計画行列QR分解
A  QR
Q Q  I ,R 
T
AT Ax  AT y
R Q QRx  R Q y
T
T
T
T
Rx  Q T y
http://gpspp.sakura.ne.jp
最小二乗法解法 (2)
• QR分解
修正Gram-Shmidt
Householder変換
Givens-Rotation
• 逐次改良
http://gpspp.sakura.ne.jp
行列内部表現 (1)
FORTRAN/MATLAB
(BLAS/LAPACK)
C
1-ORIGIN
0-ORIGIN
Index
Memory
Order
n
n
m
Column-Major
Access to
Member
A(i,j)
m
Row-Major
A[i][j]
(a[i+j*n]:Column-Major)
(i=1,2,...n,j=1,2,...m) (i=0,1,...n-1,j=0,1,...m-1)
http://gpspp.sakura.ne.jp
行列内部表現 (2)
• BLAS/LAPACK互換性
• 標準C(89)サポート (可搬性)
• メモリ確保用ライブラリ : A=mat(n,m)
• Column-Major Order (Fortran形式)
• 行列要素参照 : A[i+j*n]
http://gpspp.sakura.ne.jp

Download Report