least_squares_deficient() [1/4]

bool Lapack::least_squares_deficient ( const DMatCC * A, const DMatCC * b, DMatCC * x )

static

Definition at line 1728 of file lapack.cpp.

{
  bool ret = true;
  int rows = static_cast<int>(A->numRows());
  int cols = static_cast<int>(A->numColumns());
  int brows = static_cast<int>(b->numRows());
  int bcols = static_cast<int>(b->numColumns());
  double rcond = -1.0;
  int rank, info;
  int min = (rows < cols) ? rows : cols;
  int max = (rows > cols) ? rows : cols;
  int wsize = 2 * min + ((bcols > max) ? bcols : max);
 
  if (brows != rows)
    {
      ERROR("expected compatible right hand side");
      return false;
    }
 
  if (min == 0 || bcols == 0)
    {
      x->resize(cols, bcols);
      return true;
    }
 
  std::vector<double> copyA = make_lapack_array(*A);
  std::vector<double> copyb = make_lapack_array(*b);
  double *workspace = new double[2 * wsize];
  double *sing = new double[min];
  double *rwork = new double[5 * min];
 
  if (rows < cols)
    {
      // Make 'b' (copyb) into a cols x bcols matrix, by adding a zero block at
      // the bottom
      std::vector<double> copyb2(2 * cols * bcols);
      for (auto& a : copyb2) a = 0;
      int bloc = 0;
      int copyloc = 0;
      for (int j = 0; j < bcols; j++)
        {
          copyloc = 2 * j * cols;
          for (int i = 0; i < 2 * brows; i++)
            {
              copyb2[copyloc++] = copyb[bloc++];
            }
        }
      std::swap(copyb, copyb2);
    }
  
  zgelss_(&rows,
          &cols,
          &bcols,
          copyA.data(),
          &rows,
          copyb.data(),
          &max,
          sing,
          &rcond,
          &rank,
          workspace,
          &wsize,
          rwork,
          &info);
 
  if (info != 0)
    {
      ERROR("argument passed to zgelss had an illegal value");
      ret = false;
    }
  else
    {
      x->resize(cols, bcols);
      if (rows > cols)
        {
          // We only need the first 'cols' rows of copyb
          int copyloc = 0;
          for (int j = 0; j < bcols; j++)
            {
              copyloc = 2 * j * rows;
              for (int i = 0; i < cols; i++)
                {
                  double re = copyb[copyloc++];
                  double im = copyb[copyloc++];
                  x->ring().set_from_doubles(x->entry(i,j), re, im);
                }
            }
        }
      else
        {
          fill_from_lapack_array(copyb, *x);
        }
    }
 
  delete[] workspace;
  delete[] sing;
  delete[] rwork;
 
  return ret;
}

References ERROR, fill_from_lapack_array(), make_lapack_array(), max, DMat< ACoeffRing >::numColumns(), DMat< ACoeffRing >::numRows(), std::swap(), x, and zgelss_().