QR() [2/2]

bool Lapack::QR ( const DMatRR * A, DMatRR * Q, DMatRR * R, bool return_QR )

static

Definition at line 898 of file lapack.cpp.

{
  // sizes:
  //  input: A[m,n]
  //  output for returnQR==true:
  //   case m >= n:  Q[m,n], R[n,n]
  //   case m < n: Q[m,m], R[m,n]
  //  output for returnQR==false:
  //   Q[m,n], R[1,min(m,n)]
  bool ret = true;
  int rows = static_cast<int>(A->numRows());
  int cols = static_cast<int>(A->numColumns());
  int info1 = 0, info2 = 0, info3 = 0, info4 = 0;
  int min = (rows <= cols) ? rows : cols;
 
  if (min == 0)
    {
      ERROR("expected a matrix with positive dimensions");
      return false;
    }
 
  std::vector<double> copyA = make_lapack_array(*A); 
  std::vector<double> tau (min); // TODO: set to 0??
  double workspace_size[1];
  int work_size = -1;
  // find optimal workspace
  dgeqrf_(&rows, &cols, copyA.data(), &rows, tau.data(), workspace_size, &work_size, &info1);
  work_size = static_cast<int>(workspace_size[0]);
  // std::cout << "work size for QR:  " << work_size << std::endl;
  double *workspace = new double[work_size];
 
  dgeqrf_(&rows, &cols, copyA.data(), &rows, tau.data(), workspace, &work_size, &info2);
 
  if (info1 < 0 or info2 < 0)
    {
      ERROR("argument passed to dgeqrf had an illegal value");
      ret = false;
    }
  if (info1 > 0 or info2 > 0)
    {
      ERROR("can this happen?");
      ret = false;
    }
 
  if (ret)
    {
      if (return_QR)
        {
          Q->resize(rows, cols);
          R->resize(cols, cols);
          fill_from_lapack_upper(copyA, rows, cols, *R);
 
          // Reset Q, R, with their values.
          int orgqr_work_size = -1;
          dorgqr_(&rows,
                  &cols,
                  &min,
                  copyA.data(),
                  &rows,  // lda?
                  tau.data(),
                  workspace_size,
                  &orgqr_work_size,
                  &info3);
          orgqr_work_size = static_cast<int>(workspace_size[0]);
          if (orgqr_work_size > work_size)
            {
              delete[] workspace;
              work_size = orgqr_work_size;
              workspace = new double[work_size];
              std::cout << "work size increased to: " << work_size << std::endl;
            }
          dorgqr_(&rows,
                  &cols,
                  &min,
                  copyA.data(),
                  &rows,  // lda?
                  tau.data(),
                  workspace,
                  &work_size,
                  &info4);
          if (info3 < 0 or info4 < 0)
            {
              ERROR("argument passed to dorgqr or dorgqr had an illegal value");
              ret = false;
            }
          else if (info3 > 0 or info4 > 0)
            {
              ERROR("can this happen?");
              ret = false;
            }
          else
            {
              fill_from_lapack_array(copyA, *Q);
            }
        }
      else
        {
          // Return the raw values for QR: the "A" matrix encodes R and the
          // Householders, and tau has the multipliers.
          Q->resize(rows, cols);
          R->resize(1, min);
          fill_from_lapack_array(copyA, *Q);
          fill_from_lapack_array(tau, *R);
        }
    }
 
  delete[] workspace;
 
  return ret;
}

References dgeqrf_(), dorgqr_(), ERROR, fill_from_lapack_array(), fill_from_lapack_upper(), make_lapack_array(), DMat< ACoeffRing >::numColumns(), DMat< ACoeffRing >::numRows(), and DMat< ACoeffRing >::resize().

Referenced by MatrixOps::QR(), and MatrixOps::QR().