eigenvectors() [3/4]

bool Lapack::eigenvectors ( const DMatRR * A, DMatCC * eigenvals, DMatCC * eigenvecs )

static

Definition at line 365 of file lapack.cpp.

{
  int size = static_cast<int>(A->numRows());
  if (size != static_cast<int>(A->numColumns()))
    {
      ERROR("expected a square matrix");
      return false;
    }
 
  if (size == 0)
    {
      eigvals->resize(0, 1); // TODO: is this 1 correct?
      eigvecs->resize(0, 0);
      return true;
    }
 
  bool ret = true;
  char dont = 'N';
  char doit = 'V';
  int wsize = 4 * size;
 
  double *workspace = new double[2*wsize]; 
  int info;
 
  std::vector<double> copyA = make_lapack_array(*A);
  double *real = new double[size]; // real components of eigvals
  double *imag = new double[size];  // imaginary components
  double *eigen = new double[size * size]; // eigvecs
 
  dgeev_(&dont, /* left e-vectors */
         &doit, /* right e-vectors */
         &size,
         copyA.data(),
         &size,
         real,
         imag,
         static_cast<double *>(nullptr),
         &size, /* left eigvecs */
         eigen,
         &size, /* right eigvecs */
         workspace,
         &wsize,
         &info);
 
  if (info < 0)
    {
      ERROR("argument passed to dgeev had an illegal value");
      ret = false;
    }
  else if (info > 0)
    {
      ERROR("the QR algorithm in dgeev failed to compute all eigvals");
      ret = false;
    }
  else
    {
      // Make the complex arrays of eigvals and eigvecs
      eigvals->resize(size, 1);
      eigvecs->resize(size, size);
      //      DMatCC::ElementType *elems = eigvals->rowMajorArray();
      double* eigenLoc = eigen; // current row (eigenvector) in the eigen array
      for (int j = 0; j < size; j++, eigenLoc += size)
        {
          eigvals->ring().set_from_doubles(eigvals->entry(j,0), real[j], imag[j]);
 
          // now set j-th column of eigvecs
          if (imag[j] == 0)
            {
              for (int i = 0; i < size; ++i)
                {
                  eigvecs->entry(i,j).re = eigenLoc[i];
                  eigvecs->entry(i,j).im = 0;
                }
            }
          else if (imag[j] > 0)
            {
              for (int i = 0; i < size; ++i)
                {
                  eigvecs->ring().set_from_doubles(eigvecs->entry(i,j),
                                                   eigenLoc[i], eigenLoc[size + i]);
                  eigvecs->ring().set_from_doubles(eigvecs->entry(i,j+1),
                                                   eigenLoc[i], - eigenLoc[size + i]);
                }
            }
        }
    }
 
  delete [] workspace;
  delete [] real;
  delete [] imag;
  delete [] eigen;
 
  return ret;
}

References dgeev_(), DMat< ACoeffRing >::entry(), ERROR, cc_doubles_struct::im, make_lapack_array(), DMat< ACoeffRing >::numColumns(), DMat< ACoeffRing >::numRows(), cc_doubles_struct::re, DMat< ACoeffRing >::resize(), DMat< ACoeffRing >::ring(), and M2::ARingCC::set_from_doubles().