[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

65. lapack


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

65.1 Introduction to lapack

lapackは SLATECプロジェクトから得られるような Fortranライブラリ LAPACKの(プログラム f2clを介した) Common Lisp翻訳です。

(訳者注意書き: lapackを使用するには、 load("lapack"); load("eigensys");を実行してください。 load("lapack")には、初回だけコンパイルで時間がかかるかもしれません。)

Categories:  Numerical methods · Share packages · Package lapack


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

65.2 Functions and Variables for lapack

関数: dgeev  
    dgeev (A)  
    dgeev (A, right_p, left_p)

行列 Aの固有値と、オプションで固有ベクトルを計算します。 Aの要素はすべて整数か浮動小数点数でなければいけません。 Aは平方(行と列が同じ数)でなければいけません。 Aは対称であってもなくてもいいです。

dgeev(A)Aの固有値だけを計算します。 dgeev(A, right_p, left_p)Aの固有値と、 right_p = trueの時、右固有ベクトル、 left_p = trueの時、左固有ベクトルを計算します。

3項目のリストが返されます。 最初の項目は固有値のリストです。 二番目の項目は falseか右固有ベクトルの行列です。 三番目の項目は falseか左固有ベクトルの行列です。

右固有ベクトル v(j) (右固有ベクトル行列の j番目の列)は

A . v(j) = lambda(j) . v(j)

を満たします。ここで lambda(j)は対応する固有値です。

左固有ベクトル v(j) (左固有ベクトル行列の j番目の列)は

u(j)**H . A = lambda(j) . u(j)**H

を満たします。ここで u(j)**Hu(j)の共役転置を意味します。 Maxima関数 ctransposeが共役転置を計算します。

計算された固有ベクトルは、 Euclideanノルムが 1に等しく、最大成分の虚部が 0になるように規格化されます。

例:

 
(%i1) load ("lapack")$
(%i2) fpprintprec : 6;
(%o2)                           6
(%i3) M : matrix ([9.5, 1.75], [3.25, 10.45]);
                         [ 9.5   1.75  ]
(%o3)                    [             ]
                         [ 3.25  10.45 ]
(%i4) dgeev (M);
(%o4)          [[7.54331, 12.4067], false, false]
(%i5) [L, v, u] : dgeev (M, true, true);
                           [ - .666642  - .515792 ]
(%o5) [[7.54331, 12.4067], [                      ],
                           [  .745378   - .856714 ]
                                        [ - .856714  - .745378 ]
                                        [                      ]]
                                        [  .515792   - .666642 ]
(%i6) D : apply (diag_matrix, L);
                      [ 7.54331     0    ]
(%o6)                 [                  ]
                      [    0     12.4067 ]
(%i7) M . v - v . D;
                [      0.0       - 8.88178E-16 ]
(%o7)           [                              ]
                [ - 8.88178E-16       0.0      ]
(%i8) transpose (u) . M - D . transpose (u);
                     [ 0.0  - 4.44089E-16 ]
(%o8)                [                    ]
                     [ 0.0       0.0      ]

Categories:  Package lapack

関数: dgeqrf (A)

行列 Aの QR分解します。 Aのすべての要素は整数か浮動小数点数でなければいけません。 Aは行と列の数は同じかもしれませんし違うかもしれません。

2つの項目のリストを返します。 一番目の項目は行列 Qで、それは Aと同じ行数を持つ平方正規直交行列です。 二番目の項目は行列 Rで、それは Aと同じサイズで、 対角以下のすべての要素が零に等しいものです。 積 Q . Rは(浮動小数点の丸め誤差を除いて) Aに等しい。 ここで "."は非可換乗算演算子です。

 
(%i1) load ("lapack") $
(%i2) fpprintprec : 6 $
(%i3) M : matrix ([1, -3.2, 8], [-11, 2.7, 5.9]) $
(%i4) [q, r] : dgeqrf (M);
       [ - .0905357  .995893  ]
(%o4) [[                      ],
       [  .995893    .0905357 ]
                               [ - 11.0454   2.97863   5.15148 ]
                               [                               ]]
                               [     0      - 2.94241  8.50131 ]
(%i5) q . r - M;
         [ - 7.77156E-16   1.77636E-15   - 8.88178E-16 ]
(%o5)    [                                             ]
         [      0.0       - 1.33227E-15   8.88178E-16  ]
(%i6) mat_norm (%, 1);
(%o6)                      3.10862E-15

Categories:  Package lapack

関数: dgesv (A, b)

線形方程式 A x = bの解 xを計算します。 ここで、 Aは平方行列、 bAと同じ数の行と任意の長さの列を持つ行列です。 戻り値 xbと同じサイズです。

Abの要素は floatを介して実の浮動小数点数に評価されなければいけません; 従って、要素は任意の数値型か、数値定数のシンボルか、浮動小数点に評価される式であり得ます。 xの要素はいつも浮動小数点数です。 すべての算術は浮動小数演算として実行されます。

dgesvAの LU分解を介して解を計算します。

例:

dgesvは線形方程式 A x = bの解を計算します。

 
(%i1) A : matrix ([1, -2.5], [0.375, 5]);
                               [   1    - 2.5 ]
(%o1)                          [              ]
                               [ 0.375    5   ]
(%i2) b : matrix ([1.75], [-0.625]);
                                  [  1.75   ]
(%o2)                             [         ]
                                  [ - 0.625 ]
(%i3) x : dgesv (A, b);
                            [  1.210526315789474  ]
(%o3)                       [                     ]
                            [ - 0.215789473684211 ]
(%i4) dlange (inf_norm, b - A.x);
(%o4)                                 0.0

bAと同じ数の行と任意の長さの列を持つ行列です。 xbと同じサイズです。

 
(%i1) A : matrix ([1, -0.15], [1.82, 2]);
                               [  1    - 0.15 ]
(%o1)                          [              ]
                               [ 1.82    2    ]
(%i2) b : matrix ([3.7, 1, 8], [-2.3, 5, -3.9]);
                              [  3.7   1    8   ]
(%o2)                         [                 ]
                              [ - 2.3  5  - 3.9 ]
(%i3) x : dgesv (A, b);
      [  3.103827540695117   1.20985481742191    6.781786185657722  ]
(%o3) [                                                             ]
      [ - 3.974483062032557  1.399032116146062  - 8.121425428948527 ]
(%i4) dlange (inf_norm, b - A . x);
(%o4)                       1.1102230246251565E-15

Abの要素は実の浮動小数点数に評価されなければいけません;

 
(%i1) A : matrix ([5, -%pi], [1b0, 11/17]);
                               [   5    - %pi ]
                               [              ]
(%o1)                          [         11   ]
                               [ 1.0b0   --   ]
                               [         17   ]
(%i2) b : matrix ([%e], [sin(1)]);
                                  [   %e   ]
(%o2)                             [        ]
                                  [ sin(1) ]
(%i3) x : dgesv (A, b);
                             [ 0.690375643155986 ]
(%o3)                        [                   ]
                             [ 0.233510982552952 ]
(%i4) dlange (inf_norm, b - A . x);
(%o4)                        2.220446049250313E-16

Categories:  Package lapack · Linear equations

関数: dgesvd  
    dgesvd (A)  
    dgesvd (A, left_p, right_p)

特異値から成る行列 Aの特異値分解 (SVD)を計算します。 オプションで左および右特異ベクトルを取ります。

Aの要素はすべて整数か浮動小数点数でなければいけません。 Aは(行と列が同じ数の)平方かもしれませんし、そうでないかもしれません。

mAの行数、nを列数とします。 Aの特異値分解は A = U . Sigma . V^Tのような 3つの行列 U, Sigma, V^Tから構成されます。 ここで、 Um掛けmのユニタリ行列、 Sigmam掛けnの対角行列、 V^Tn掛けnのユニタリ行列です。

sigma[i]Sigmaの対角要素、すなわち、 Sigma[i, i] = sigma[i]とします。 要素 sigma[i]Aのいわゆる特異値です; これらは実数で、非負で、降順で返されます。 UVの最初の min(m, n)列は Aの左と右特異ベクトルです。 dgesvdV自身ではなく Vの転置を返すことに注意してください。

dgesvd(A)Aの特異値だけを計算します。 dgesvd(A, left_p, right_p)Aの特異値と、 left_p = trueの時、左特異ベクトル、 right_p = trueの時、右特異ベクトルを計算します。

3つの項目のリストを返します。 一つ目の項目は特異値のリストです。 二つ目の項目は falseか、左特異ベクトルの行列です。 三つ目の項目は falseか、右特異ベクトルの行列です。

例:

 
(%i1) load ("lapack")$
(%i2) fpprintprec : 6;
(%o2)                           6
(%i3) M: matrix([1, 2, 3], [3.5, 0.5, 8], [-1, 2, -3], [4, 9, 7]);
                        [  1    2    3  ]
                        [               ]
                        [ 3.5  0.5   8  ]
(%o3)                   [               ]
                        [ - 1   2   - 3 ]
                        [               ]
                        [  4    9    7  ]
(%i4) dgesvd (M);
(%o4)      [[14.4744, 6.38637, .452547], false, false]
(%i5) [sigma, U, VT] : dgesvd (M, true, true);
(%o5) [[14.4744, 6.38637, .452547],
[ - .256731  .00816168   .959029    - .119523 ]
[                                             ]
[ - .526456   .672116   - .206236   - .478091 ]
[                                             ],
[  .107997   - .532278  - .0708315  - 0.83666 ]
[                                             ]
[ - .803287  - .514659  - .180867    .239046  ]
[ - .374486  - .538209  - .755044 ]
[                                 ]
[  .130623   - .836799   0.5317   ]]
[                                 ]
[ - .917986   .100488    .383672  ]
(%i6) m : length (U);
(%o6)                           4
(%i7) n : length (VT);
(%o7)                           3
(%i8) Sigma:
        genmatrix(lambda ([i, j], if i=j then sigma[i] else 0),
                  m, n);
                  [ 14.4744     0        0    ]
                  [                           ]
                  [    0     6.38637     0    ]
(%o8)             [                           ]
                  [    0        0     .452547 ]
                  [                           ]
                  [    0        0        0    ]
(%i9) U . Sigma . VT - M;
          [  1.11022E-15        0.0       1.77636E-15 ]
          [                                           ]
          [  1.33227E-15    1.66533E-15       0.0     ]
(%o9)     [                                           ]
          [ - 4.44089E-16  - 8.88178E-16  4.44089E-16 ]
          [                                           ]
          [  8.88178E-16    1.77636E-15   8.88178E-16 ]
(%i10) transpose (U) . U;
       [     1.0      5.55112E-17    2.498E-16     2.77556E-17  ]
       [                                                        ]
       [ 5.55112E-17      1.0       5.55112E-17    4.16334E-17  ]
(%o10) [                                                        ]
       [  2.498E-16   5.55112E-17       1.0       - 2.08167E-16 ]
       [                                                        ]
       [ 2.77556E-17  4.16334E-17  - 2.08167E-16       1.0      ]
(%i11) VT . transpose (VT);
          [      1.0           0.0      - 5.55112E-17 ]
          [                                           ]
(%o11)    [      0.0           1.0       5.55112E-17  ]
          [                                           ]
          [ - 5.55112E-17  5.55112E-17       1.0      ]

Categories:  Package lapack

関数: dlange (norm, A)
関数: zlange (norm, A)

行列 Aのノルムもしくはノルムのような関数を計算します。

max

max(abs(A(i, j)))を計算します。 ここで ijはそれぞれ行と列を行き渡ります。 この関数は適切な行列ノルムではないことに注意してください。

one_norm

AL[1]ノルム、すなわち、それぞれの列の要素の絶対値の和の最大値を計算します。

inf_norm

AL[inf]ノルム、すなわち、それぞれの行の要素の絶対値の和の最大値を計算します。

frobenius

Aの Frobeniusノルム、すなわち、行列要素の平方の和の平方根を計算します。

Categories:  Package lapack

関数: dgemm  
    dgemm (A, B)  
    dgemm (A, B, options)

2つの行列の積を計算します。オプションで積を三つ目の行列に足し算します。

最も簡単な形式では、 dgemm(A, B)は 2つの実行列 ABの積を計算します。

二番目の形式では、 dgemmalpha * A * B + beta * Cを計算します。 ここで A, B, Cは適当なサイズの実行列であり、 alphabetaは実数です。 オプションで、 Aと/もしくは Bは積を計算する前に転置を取ることができます。 追加のパラメータはオプションのキーワード引数で指定できます: キーワード引数はオプションで、どんな順番でも指定できます。 それらはすべて、形式 key=valを取ります。 キーワード引数は以下の通りです:

C

足すべき行列 C。 デフォルトは falseであり、行列を足さないことを意味します。

alpha

ABの積をこの値に掛けます。 デフォルトは 1です。

beta

もし行列 Cが与えられたら、 この値を、足される前にCに掛けます。 デフォルト値は 0で、これはたとえCが与えられても Cが足されないことを意味します。 故に、必ず betaに零でない値を指定してください。

transpose_a

もし trueなら、 Aの代わりに Aの転置を積に使います。 デフォルトは falseです。

transpose_b

もし trueなら Bの代わりに Bの転置を積に使います。 デフォルトは falseです。

 
(%i1) load ("lapack")$
(%i2) A : matrix([1,2,3],[4,5,6],[7,8,9]);
                                  [ 1  2  3 ]
                                  [         ]
(%o2)                             [ 4  5  6 ]
                                  [         ]
                                  [ 7  8  9 ]
(%i3) B : matrix([-1,-2,-3],[-4,-5,-6],[-7,-8,-9]);
                               [ - 1  - 2  - 3 ]
                               [               ]
(%o3)                          [ - 4  - 5  - 6 ]
                               [               ]
                               [ - 7  - 8  - 9 ]
(%i4) C : matrix([3,2,1],[6,5,4],[9,8,7]);
                                  [ 3  2  1 ]
                                  [         ]
(%o4)                             [ 6  5  4 ]
                                  [         ]
                                  [ 9  8  7 ]
(%i5) dgemm(A,B);
                         [ - 30.0   - 36.0   - 42.0  ]
                         [                           ]
(%o5)                    [ - 66.0   - 81.0   - 96.0  ]
                         [                           ]
                         [ - 102.0  - 126.0  - 150.0 ]
(%i6) A . B;
                            [ - 30   - 36   - 42  ]
                            [                     ]
(%o6)                       [ - 66   - 81   - 96  ]
                            [                     ]
                            [ - 102  - 126  - 150 ]
(%i7) dgemm(A,B,transpose_a=true);
                         [ - 66.0  - 78.0   - 90.0  ]
                         [                          ]
(%o7)                    [ - 78.0  - 93.0   - 108.0 ]
                         [                          ]
                         [ - 90.0  - 108.0  - 126.0 ]
(%i8) transpose(A) . B;
                           [ - 66  - 78   - 90  ]
                           [                    ]
(%o8)                      [ - 78  - 93   - 108 ]
                           [                    ]
                           [ - 90  - 108  - 126 ]
(%i9) dgemm(A,B,c=C,beta=1);
                         [ - 27.0  - 34.0   - 41.0  ]
                         [                          ]
(%o9)                    [ - 60.0  - 76.0   - 92.0  ]
                         [                          ]
                         [ - 93.0  - 118.0  - 143.0 ]
(%i10) A . B + C;
                            [ - 27  - 34   - 41  ]
                            [                    ]
(%o10)                      [ - 60  - 76   - 92  ]
                            [                    ]
                            [ - 93  - 118  - 143 ]
(%i11) dgemm(A,B,c=C,beta=1, alpha=-1);
                            [ 33.0   38.0   43.0  ]
                            [                     ]
(%o11)                      [ 72.0   86.0   100.0 ]
                            [                     ]
                            [ 111.0  134.0  157.0 ]
(%i12) -A . B + C;
                               [ 33   38   43  ]
                               [               ]
(%o12)                         [ 72   86   100 ]
                               [               ]
                               [ 111  134  157 ]

Categories:  Package lapack

関数: zgeev  
    zgeev (A)  
    zgeev (A, right_p, left_p)

dgeev同様、しかし行列 Aは複素行列です。

Categories:  Package lapack

関数: zheev  
    zheev (A)  
    zheev (A, eigvec_p)

zheev同様、しかし行列 Aは平方複素エルミート行列と仮定されます。 もし eigvec_ptrueなら、 行列の固有ベクトルも計算します。

実際には行列 Aがエルミートであることはチェックしません。

dgeevと同じように2つの項目のリストを返します: 固有値のリストと、 falseまたは固有ベクトルの行列。

Categories:  Package lapack


[ << ] [ >> ]           [Top] [Contents] [Index] [ ? ]

This document was generated by 市川雄二 on October, 5 2017 using texi2html 1.76.