SVD_divide_conquer() [3/4]

bool Lapack::SVD_divide_conquer ( const DMatRR * A, DMatRR * Sigma, DMatRR * U, DMatRR * VT )

static

Definition at line 636 of file lapack.cpp.

{
  bool ret = true;
  char doit = 'A';  // other options are 'S' and 'O' for singular vectors only
  int rows = static_cast<int>(A->numRows());
  int cols = static_cast<int>(A->numColumns());
  int info;
  int min = (rows <= cols) ? rows : cols;
 
  if (min == 0)
    {
      ERROR("expected a matrix with positive dimensions");
      return false;
    }
 
  int max = (rows >= cols) ? rows : cols;
  int wsize = 4 * min * min + max + 9 * min;
 
  double *workspace = new double[wsize];
  int* iworkspace = new int[8 * min];
  
  std::vector<double> copyA = make_lapack_array(*A);
  std::vector<double> u(rows * rows);
  std::vector<double> vt(cols * cols);
  std::vector<double> sigma(min);
  
  dgesdd_(&doit,
          &rows,
          &cols,
          copyA.data(),
          &rows,
          sigma.data(),
          u.data(),
          &rows,
          vt.data(),
          &cols,
          workspace,
          &wsize,
          iworkspace,
          &info);
 
  if (info < 0)
    {
      ERROR("argument passed to dgesdd had an illegal value");
      ret = false;
    }
  else if (info > 0)
    {
      ERROR("dgesdd did not converge");
      ret = false;
    }
  else
    {
      U->resize(rows, rows);
      VT->resize(cols, cols);
      Sigma->resize(min, 1);
      fill_from_lapack_array(u, *U);
      fill_from_lapack_array(vt, *VT);
      fill_from_lapack_array(sigma, *Sigma);
    }
 
  delete [] workspace;
  delete [] iworkspace;
 
  return ret;
}

References dgesdd_(), ERROR, fill_from_lapack_array(), make_lapack_array(), max, DMat< ACoeffRing >::numColumns(), DMat< ACoeffRing >::numRows(), and DMat< ACoeffRing >::resize().