I don't know the specification of the return value of the dgeqr function for lapack.

#include<array>
# include <cassert>
# include <iostream>
# include <vector>

using namespace std;

US>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&lwork, int&info);
}

void print_matrix(intm, intn, double*a){
  // column-major
  for(inti=0;i<m;i++){
    for(int j=0;j<n;j++){
      cout<<a[i+j*m]<";
    }
    cout<<endl;
  }
}

int main() {
  // example
  constintm=3;//#ofrows
  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_query[5];
  double work_query[1];
  int info;

  dgeqr_(m, n, a.data(), m, t_query, -1, work_query, -1, info);
  assert(info==0);

  intsize=(int)t_query[0];
  intlwork=(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 intrm = min(m,n);
  const intrn = n;
  array<double,(rm*rn)>r={};
  for(inti=0;i<rm;i++){
    for(int j=0;j<rn;j++){
      if(i<=j){
        r[i+j*rm] = a[i+j*m]; // diagnostic and above
      }
    }
  }

  // workspace query
  const int qm = m;
  const int qn = min(m,n);
  array<double,(qm*qn)>q={};
  for(inti=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());
}

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