

c++でopenblasからlapackを使っています。  
QR分解のためのルーチンとしてdgeqr関数があるのは調べたのですが、返り値の仕様が理解できていません。  
目的としては、行列AをQR分解して、直行行列Qと上三角行列Rを得たいです。  
どのように呼び出して、上記の2つの値を得ればいいのでしょうか？


CC BY-SA 4.0
2022-11-30T05:49:18.390
92511
2022-12-03T06:01:52.570
40768
question
c++
lapackのdgeqr関数の返り値の仕様がわからない


`dgeqr`を呼び出した後、`R`は（書き換えられた）`A`の中に含まれています。`Q`を得るには`dgemqr`を呼び出す必要があります。LAPACKはFortranのサブルーチンなので列優先であることに注意してください。以下はC++14での例です。

[code]

    #include <array>
    #include <cassert>
    #include <iostream>
    #include <vector>
    
    using namespace std;
    
    extern "C" {
    void dgeqr_(const int &m, const int &n, double *a, const int &lda, double *t,
                const int &tsize, double *work, const int &lwork, int &info);
    void dgemqr_(const char &side, const char &trans, const int &m, const int &n,
                 const int &k, const double *a, const int &lda, const double *t,
                 const int &tsize, double *c, const int &ldc, double *work,
                 const int &lwork, int &info);
    }
    
    void print_matrix(int m, int n, double *a) {
      // column-major
      for (int i = 0; i < m; i++) {
        for (int j = 0; j < n; j++) {
          cout << a[i + j * m] << " ";
        }
        cout << endl;
      }
    }
    
    int main() {
      // example
      const int m = 3; // # of rows
      const int n = 3; // # of columns
      array<double, (m * n)> a = {
          12,  6,   -4,  // 1st column
          -51, 167, 24,  // 2nd column
          4,   -68, -41, // 3rd column
      };
    
      // print A
      cout << "A =" << endl;
      print_matrix(m, n, a.data());
    
      // workspace query
      double t_query[5];
      double work_query[1];
      int info;
    
      dgeqr_(m, n, a.data(), m, t_query, -1, work_query, -1, info);
      assert(info == 0);
    
      int tsize = (int)t_query[0];
      int lwork = (int)work_query[0];
    
      // perform QR factorization
      vector<double> t(max(5, tsize));
      vector<double> work(max(1, lwork));
    
      dgeqr_(m, n, a.data(), m, t.data(), tsize, work.data(), lwork, info);
      assert(info == 0);
    
      // extract R
      const int rm = min(m, n);
      const int rn = n;
      array<double, (rm * rn)> r = {};
      for (int i = 0; i < rm; i++) {
        for (int j = 0; j < rn; j++) {
          if (i <= j) {
            r[i + j * rm] = a[i + j * m];  // diagonal and above
          }
        }
      }
    
      // workspace query
      const int qm = m;
      const int qn = min(m, n);
      array<double, (qm * qn)> q = {};
      for (int i = 0; i < min(qm, qn); i++) {
        q[i + i * qm] = 1;  // (possibly truncated) identity matrix
      }
      dgemqr_('L', 'N', qm, qn, qn, a.data(), m, t.data(), tsize, q.data(), qm,
              work_query, -1, info);
      assert(info == 0);
    
      lwork = (int)work_query[0];
    
      // reconstruct Q
      work.reserve(max(1, lwork));
    
      dgemqr_('L', 'N', qm, qn, qn, a.data(), m, t.data(), tsize, q.data(), qm,
              work.data(), lwork, info);
      assert(info == 0);
    
      // print Q
      cout << "Q =" << endl;
      print_matrix(qm, qn, q.data());
    
      // print R
      cout << "R =" << endl;
      print_matrix(rm, rn, r.data());
    }
    
[/code]

結果：

[code]

    A =
    12 -51 4
    6 167 -68
    -4 24 -41
    Q =
    -0.857143 0.394286 0.331429
    -0.428571 -0.902857 -0.0342857
    0.285714 -0.171429 0.942857
    R =
    -14 -21 14
    0 -175 70
    0 0 -35
    
[/code]


CC BY-SA 4.0
2022-12-02T13:53:47.333
92541
2022-12-03T06:01:52.570
2022-12-03T06:01:52.570
36010
36010
92511
answer