SVD_divide_conquer() [1/4]

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

static

Definition at line 1561 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 = min * min + 2 * min + max;
 
  double *workspace = new double[2 * wsize];
  int *iworkspace = new int[8 * min];
  double *rwork = new double[5 * min * min + 7 * min]; // documentation not clear
 
  std::vector<double> copyA = make_lapack_array(*A);
  std::vector<double> u(2 * rows * rows);
  std::vector<double> vt(2 * cols * cols);
  std::vector<double> sigma(2 * min);
 
  zgesdd_(&doit,
          &rows,
          &cols,
          copyA.data(),
          &rows,
          sigma.data(),
          u.data(),
          &rows,
          vt.data(),
          &cols,
          workspace,
          &wsize,
          rwork,
          iworkspace,
          &info);
 
  if (info < 0)
    {
      ERROR("argument passed to zgesdd had an illegal value");
      ret = false;
    }
  else if (info > 0)
    {
      ERROR("zgesdd did not converge");
      ret = false;
    }
  else
    {
      U->resize(rows, rows);
      fill_from_lapack_array(u, *U);
      VT->resize(cols, cols);
      fill_from_lapack_array(vt, *VT);
      Sigma->resize(min, 1);
      fill_from_lapack_array(sigma, *Sigma);
    }
 
  delete[] workspace;
  delete[] iworkspace;
  delete[] rwork;  
 
  return ret;
}

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