// ------------------------------------------------------------------
// qpatchmath.cpp
//
// This file contains routines for math with patches including:
// converting parametric to real coordinates, computing directional
// derivatives, computing normals, and intersecting patches with
// lines.
// ------------------------------------------------------------------
// Copyright (c) 1999, 2000 by Cornell University.  All rights reserved.
// 
// See the accompanying file 'Copyright' for authorship information,
// the terms of the license governing this software, and disclaimers
// concerning this software.
// ------------------------------------------------------------------
// This file is part of the QMG software.  
// Version 2.0 of QMG, release date September 3, 1999
// Version 2.0.1 -- rewritten to estimate generalized eigenvector
//  sensitivity, October 29, 1999.  Steve Vavasis
// Version 2.0.2 -- rewritten with a new way to handle
//  degenerate problems using doubledouble and corrections.
//  Released March 31, 2000.
// ------------------------------------------------------------------
#ifdef _MSC_VER
#if _MSC_VER == 1200
#pragma warning (disable: 4786) 
#pragma warning (disable: 4788)
#endif
#endif




#include "qpatchmath.h"
#include "qmatvec.h"



#ifdef DEBUGGING
namespace QMG {
  namespace Local_to_PatchMath_CPP {
    using namespace QMG;
    ostream* debug_str;
    bool dump_everything = false;
    int invocation_count = 0;
  }
}
#endif


// Types used throughout this file.

// Types passed to Lapack routines

#ifdef LAPACK_DOUBLECOMPLEX_T
typedef LAPACK_DOUBLECOMPLEX_T Lapack_doublecomplex;
#else
typedef double Lapack_doublecomplex;
#endif

#ifdef LAPACK_INT_T
typedef LAPACK_INT_T Lapack_int;
#else
typedef long int Lapack_int;
#endif

typedef std::complex<double> Complex;


// LAPACK routine for complex standard nonsymmetric eigenvalues.


extern "C"
int zgeev_(char *jobvl, 
           char *jobvr, 
           Lapack_int *n,
           Lapack_doublecomplex *a, 
           Lapack_int *lda, 
           Lapack_doublecomplex *w, 
           Lapack_doublecomplex *vl,
           Lapack_int *ldvl, 
           Lapack_doublecomplex *vr, 
           Lapack_int *ldvr, 
           Lapack_doublecomplex* work,
           Lapack_int *lwork, 
           double *rwork, 
           Lapack_int *info);

// LAPACK routine for qr factorization

extern "C"
int zgeqrf_(Lapack_int *m, 
            Lapack_int *n, 
            Lapack_doublecomplex *a,
            Lapack_int *lda, 
            Lapack_doublecomplex *tau, 
            Lapack_doublecomplex *work, 
            Lapack_int *lwork,
            Lapack_int *info);

// LAPACK routine for applying Housholder reflections to a matrix.

extern "C"
int zunmqr_(char *side, 
            char *trans, 
            Lapack_int *m, 
            Lapack_int *n,
            Lapack_int *k, 
            Lapack_doublecomplex *a, 
            Lapack_int *lda, 
            Lapack_doublecomplex *tau,
            Lapack_doublecomplex *c, 
            Lapack_int *ldc, 
            Lapack_doublecomplex *work, 
            Lapack_int *lwork,
            Lapack_int *info);

// LAPACK routine for multiplying out a product of Householder reflections.

extern "C"
int zungqr_(Lapack_int *m, 
            Lapack_int *n, 
            Lapack_int *k,
            Lapack_doublecomplex *a, 
            Lapack_int *lda, 
            Lapack_doublecomplex *tau, 
            Lapack_doublecomplex *work, 
            Lapack_int *lwork, 
            Lapack_int *info);

// Lapack routine for pencil Hessenberg reduction.

extern "C" 
int zgghrd_(char *compq, 
            char *compz, 
            Lapack_int *n, 
            Lapack_int *ilo, 
            Lapack_int *ihi, 
            Lapack_doublecomplex *a, 
            Lapack_int *lda, 
            Lapack_doublecomplex *b,
            Lapack_int *ldb, 
            Lapack_doublecomplex *q, 
            Lapack_int *ldq, 
            Lapack_doublecomplex *z,
            Lapack_int *ldz, 
            Lapack_int *info);

// Lapack routine for QZ factorization.


extern "C" 
int zhgeqz_(char *job, 
            char *compq, 
            char *compz, 
            Lapack_int *n,
            Lapack_int *ilo, 
            Lapack_int *ihi, 
            Lapack_doublecomplex *a, 
            Lapack_int *lda, 
            Lapack_doublecomplex *b, 
            Lapack_int *ldb, 
            Lapack_doublecomplex *alpha, 
            Lapack_doublecomplex *beta, 
            Lapack_doublecomplex *q, 
            Lapack_int *ldq, 
            Lapack_doublecomplex *z, 
            Lapack_int *ldz,
            Lapack_doublecomplex *work, 
            Lapack_int *lwork, 
            double *rwork, 
            Lapack_int *info);

namespace QMG {
  extern void QR_Factor_with_col_pivoting_quadprec(double* A, 
    int m,
    int n,
    double* select_col,
    double* hhmat,
    double* betavec,
    doubledouble* wksp);
}

// Helper class that makes private routines public for this file only.

class QMG::PatchMath::Bezier_Helper {
public:
  static Real control_point_coord(const PatchMath& pm, int cpnum, int d) {
    return pm.control_point_coord_(cpnum,d);
  }
};

// Class Workspace is a workspace for math operations.
// Used so that we don't have to keep calling 'new' to
// get workspace for Lapack.
// Used as follows: Set a marker.  Then allocate vectors
// and matrices.  When the marker is destroyed, the space
// is deallocated.

class QMG::PatchMath::Workspace::DoubleDoubleVector {
private:
  Workspace& wkspa_;
  int base_;
  int sz_;
  // no copying, no assignment.
  DoubleDoubleVector(const DoubleDoubleVector& o) : wkspa_(o.wkspa_)  { }
  void operator=(const DoubleDoubleVector&) { }
public:
  DoubleDoubleVector(Workspace& wkspa, int sz) : wkspa_(wkspa),
    base_(wkspa.ddwkspa_.size()),
    sz_(sz) {
#ifdef RANGECHECK
    if (!wkspa.marker_active_)
      throw_error("Attempt to allocate in inactive workspace");
#endif
    wkspa.ddwkspa_.resize(base_ + sz_);
  }
  ~DoubleDoubleVector() { }
  doubledouble& operator[](int i) {
#ifdef RANGECHECK
    if (i < 0 || i >= sz_)
      throw_error("Range err");
#endif
    return wkspa_.ddwkspa_[base_ + i];
  }
};


class QMG::PatchMath::Workspace::RealVector {
private:
  Workspace& wkspa_;
  int base_;
  int sz_;
  // no copying, no assignment.
  RealVector(const RealVector& o) : wkspa_(o.wkspa_)  { }
  void operator=(const RealVector&) { }
public:
  RealVector(Workspace& wkspa, int sz) : wkspa_(wkspa),
    base_(wkspa.rwkspa_.size()),
    sz_(sz) {
#ifdef RANGECHECK
    if (!wkspa.marker_active_)
      throw_error("Attempt to allocate in inactive workspace");
#endif
    wkspa.rwkspa_.resize(base_ + sz_);
  }
  ~RealVector() { }
  Real& operator[](int i) {
#ifdef RANGECHECK
    if (i < 0 || i >= sz_)
      throw_error("Range err");
#endif
    return wkspa_.rwkspa_[base_ + i];
  }
  const Real& operator[](int i) const {
#ifdef RANGECHECK
    if (i < 0 || i >= sz_)
      throw_error("Range err");
#endif
    return wkspa_.rwkspa_[base_ + i];
  }
  double* make_fortran_real() {
    return &(wkspa_.rwkspa_[base_]);
  }
};


// Matrices stored by column for fortran compatibility.

class QMG::PatchMath::Workspace::ComplexMatrix {
private:
  Workspace& wkspa_;
  int base_;
  int nr_;
  int nc_;
  // no copying, no assignment
  ComplexMatrix(const ComplexMatrix& o) : wkspa_(o.wkspa_) { }
  void operator=(const ComplexMatrix&) { }
public:
  ComplexMatrix(Workspace& wkspa, int nr, int nc) : wkspa_(wkspa),
    base_(wkspa.cwkspa_.size()),
    nr_(nr),
    nc_(nc) {
    wkspa.cwkspa_.resize(base_ + nr * nc);
#ifdef RANGECHECK
    if (!wkspa.marker_active_)
      throw_error("Attempt to allocate in inactive workspace");
#endif
  }
  ~ComplexMatrix() { }
  Complex& operator()(int i, int j) {
#ifdef RANGECHECK
    if (i < 0 || i >= nr_ || j < 0 || j >= nc_)
      throw_error("Range err");
#endif
    return wkspa_.cwkspa_[base_ + i + j * nr_];
  }
  const Complex& operator()(int i, int j) const {
#ifdef RANGECHECK
    if (i < 0 || i >= nr_ || j < 0 || j >= nc_)
      throw_error("Range err");
#endif
    return wkspa_.cwkspa_[base_ + i + j * nr_];
  }
  Lapack_doublecomplex* make_fortran_complex() {
    return reinterpret_cast<Lapack_doublecomplex*>(
      static_cast<Complex*>(&(wkspa_.cwkspa_[base_])));
  }
  int numrow() const { return nr_; }
  int numcol() const { return nc_; }
};



class QMG::PatchMath::Workspace::ComplexVector {
private:
  Workspace& wkspa_;
  int base_;
  int sz_;
  // no copying, no assignment
  ComplexVector(const ComplexVector& o) : wkspa_(o.wkspa_) { }
  void operator=(const ComplexVector&) { }
public:
  ComplexVector(Workspace& wkspa, int sz) : wkspa_(wkspa),
    base_(wkspa.cwkspa_.size()),
    sz_(sz) {
    wkspa.cwkspa_.resize(base_ + sz_, Complex(0.0,0.0));
#ifdef RANGECHECK
    if (!wkspa.marker_active_)
      throw_error("Attempt to allocate in inactive workspace");
#endif
  }
  ~ComplexVector() { }
  Complex& operator[](int i) {
#ifdef RANGECHECK
    if (i < 0 || i >= sz_)
      throw_error("Range err");
#endif
    return wkspa_.cwkspa_[base_ + i];
  }
  const Complex& operator[](int i) const {
#ifdef RANGECHECK
    if (i < 0 || i >= sz_)
      throw_error("Range err");
#endif
    return wkspa_.cwkspa_[base_ + i];
  }
  Lapack_doublecomplex* make_fortran_complex() {
    return reinterpret_cast<Lapack_doublecomplex*>(
      static_cast<Complex*>(&(wkspa_.cwkspa_[base_])));
  }

};

// Real matrix stored by column for fortran compatibility.

class QMG::PatchMath::Workspace::RealMatrix {
private:
  Workspace& wkspa_;
  int base_;
  int nr_;
  int nc_;
  int nc_max_;
  // no copying, no assignment
  RealMatrix(const RealMatrix& o) : wkspa_(o.wkspa_)  { }
  void operator=(const RealMatrix&) { }
public:
  RealMatrix(Workspace& wkspa, int nr, int nc) : wkspa_(wkspa),
    base_(wkspa.rwkspa_.size()),
    nr_(nr),
    nc_(nc),
    nc_max_(nc) {
    wkspa.rwkspa_.resize(base_ + nr * nc);
#ifdef RANGECHECK
    if (!wkspa.marker_active_)
      throw_error("Attempt to allocate in inactive workspace");
#endif
  }
  ~RealMatrix() { }
  Real& operator()(int i, int j) {
#ifdef RANGECHECK
    if (i < 0 || i >= nr_ || j < 0 || j >= nc_)
      throw_error("Range err");
#endif
    return wkspa_.rwkspa_[base_ + i + j * nr_];
  }
  const Real& operator()(int i, int j) const {
#ifdef RANGECHECK
    if (i < 0 || i >= nr_ || j < 0 || j >= nc_)
      throw_error("Range err");
#endif
    return wkspa_.rwkspa_[base_ + i + j * nr_];
  }
  void change_numcol(int newcol) {
#ifdef RANGECHECK
    if (newcol < 0 || newcol > nc_max_)
      throw_error("newcol out of range");
#endif
    nc_ = newcol;
  }

  Real* make_fortran_real() {
    return &(wkspa_.rwkspa_[base_]);
  }
  int numrow() const { return nr_; }
  int numcol() const { return nc_; }
};


class QMG::PatchMath::Workspace::StartPosMarker {
private:
  Workspace& wkspa_;
  int rpos_;
  int cpos_;
  int ddpos_;
  StartPosMarker(const StartPosMarker& o) : wkspa_(o.wkspa_) { }
  void operator=(const StartPosMarker&) { }
public:
  explicit StartPosMarker(Workspace& wkspa) : wkspa_(wkspa),
    rpos_(wkspa.rwkspa_.size()),
    cpos_(wkspa.cwkspa_.size()),
  ddpos_(wkspa.ddwkspa_.size()) {
    ++wkspa.marker_active_;
  }
  ~StartPosMarker() {
    wkspa_.rwkspa_.resize(rpos_);
    wkspa_.cwkspa_.resize(cpos_);
    wkspa_.ddwkspa_.resize(ddpos_);
    --wkspa_.marker_active_;
  }
};



namespace {
  using namespace QMG;


  // Class for printing out a complex number.

  class Format_Complex;
  ostream& operator<<(ostream&, const Format_Complex&);

  class Format_Complex {
  private:
    Complex c;
  public:
    explicit Format_Complex(const Complex& c1) : c(c1) { }
    friend ostream& operator<<(ostream&, const Format_Complex&);
  };

  
  ostream& operator<<(ostream& os, const Format_Complex& fc) {
    if (fc.c.imag() >= 0)
      os << std::setprecision(17) << fc.c.real() << "+" << fc.c.imag() << "i";
    else
      os << std::setprecision(17) << fc.c.real() << "-" << -fc.c.imag() << "i";
    return os;
  }




  // QR factorization with column pivoting.
  // Return the permutation in select_col.
  // Return the Q factor in hhmat and betavec.
  // Return the R factor in A.

  void QR_Factor_with_col_pivoting(PatchMath::Workspace::RealMatrix& A, 
    PatchMath::Workspace::RealVector& select_col, 
    PatchMath::Workspace::RealMatrix& hhmat,
    PatchMath::Workspace::RealVector& betavec) {

    int m = A.numrow();
    int n = A.numcol();
    int nstep = (m < n)? m : n;
    {
      for (int i = 0; i < n; ++i)
        select_col[i] = static_cast<Real>(i);
    }
    // speed optimization

    Real* Abase = &A(0,0);

    for (int k = 0; k < nstep; ++k) {
      int select_col1;
      {
        Real mx1 = 0.0;
        for (int j = k; j < n; ++j) {
          
          // Find the norm of column j residual entries.
          Real nr1 = 0.0;
          {
            for (int i = k; i < m; ++i) {
              nr1 += Abase[i + j * m] * Abase[i + j * m];
            }
          }
          if (nr1 >= mx1) {
            mx1 = nr1;
            select_col1 = j;
          }
        }
      }
      Real tmp = select_col[k];
      select_col[k] = select_col[select_col1];
      select_col[select_col1] = tmp;
      {
        for (int i = 0; i < m; ++i) {
          Real tmp2 = Abase[i + k * m];
          Abase[i + k * m] = Abase[i + select_col1 * m];
          Abase[i + select_col1 * m] = tmp2;
        }
      }
      if (k == m - 1)
        break;
      Real* hhmat_col = &hhmat(0,k);
      Real beta = 
        Matrix::compute_hh_transform(&hhmat_col[k], &A(k, k), 1, m - k);
      betavec[k] = beta;
      {
        for (int j = 0; j < n; ++j) {
          Real ip1 = 0.0;
          {
            for (int i = k; i < m; ++i)
              ip1 += Abase[i + j * m] * hhmat_col[i];
          }
          ip1 *= beta;
          {
            for (int i = k; i < m; ++i)
              Abase[i + j * m] -= ip1 * hhmat_col[i];
          }
        }
      }
    }
  }


  // Multiply HH transforms together to get an explicit
  // orthogonal matrix Q.  Return the transpose of Q.
  // Must have Q.numcol == hhmat.numrow.
  // Must have Q.numrow <= Q.numcol.

  void formQ(const PatchMath::Workspace::RealMatrix& hhmat,
    const PatchMath::Workspace::RealVector& betavec,
    PatchMath::Workspace::RealMatrix& Q) {

    int m = hhmat.numrow();
    int n = hhmat.numcol();

    int m1 = Q.numrow();

    // Initialize Q to be the identity.

    Real* Qbase = &Q(0,0);
    {
      for (int j = 0; j < m; ++j) {
        for (int i = 0; i < m1; ++i) {
          Qbase[i + j * m1] = 0.0;
        }
        Qbase[j + j * m1] = 1.0;
      }
    }

    // Loop over HH transforms in reverse order 

    for (int k = n - 1; k >= 0; --k) {
      const Real* hhtrans = &hhmat(0,k);
      
      // Apply the kth HH transform to Q.
      // i loops over rows of Q.
      for (int i = k; i < m1; ++i) {
        Real ip1 = 0.0;
        {
          for (int j = k; j < m; ++j)
            ip1 += hhtrans[j] * Qbase[i + j * m1];
        }
        ip1 *= betavec[k];
        {
          for (int j = k; j < m; ++j)
            Qbase[i + j * m1] -= ip1 * hhtrans[j];
        }
      }
    }
  }


  inline void compute_givens(const Complex& a, 
    const Complex& b, 
    Complex& c, 
    Complex& s,
    Complex& cbar,
    Complex& sbar) {

    using std::conj;
    using std::norm;
    Real denom = sqrt(norm(a) + norm(b));
    if (denom == 0.0) {
      cbar = Complex(1.0,0.0);
      sbar = Complex(0.0,0.0);
    }
    else {
      cbar = a / denom;
      sbar = b / denom;
    }
    c = conj(cbar);
    s = conj(sbar);
  }



  // Generalized eigenvalue solver.  Solve A*x-lambda*B*x=0. 
  // Use LAPACK QZ algorithm.  Input variables are A and B.
  // (They are overwritten).  Output variable V is the eigenvectors. 
  // and eveccond, the vector of eigenvector condition numbers
  // and lambda, the eigenvectors.
  // The return value is nonzero if the QZ iteration didn't converge.

  int generalized_eigensolver(PatchMath::Workspace::ComplexMatrix& A,
    PatchMath::Workspace::ComplexMatrix& B,
    PatchMath::Workspace::ComplexMatrix& V,
    PatchMath::Workspace::RealVector& eveccond,
    PatchMath::Workspace::ComplexVector& lambda,

    // Workspace:

    PatchMath::Workspace& wksp) {

#ifdef DEBUGGING
    
    using QMG::Local_to_PatchMath_CPP::debug_str;
    using QMG::Local_to_PatchMath_CPP::dump_everything;
#endif

    PatchMath::Workspace::StartPosMarker startmark_(wksp);


    const int blksize = 64;

    Lapack_int n = static_cast<Lapack_int>(B.numrow());
    Lapack_int lwork = static_cast<Lapack_int>(n * (n + blksize));

    PatchMath::Workspace::ComplexVector cwksp(wksp, lwork);
    PatchMath::Workspace::ComplexVector tau(wksp, n);
    PatchMath::Workspace::ComplexVector alpha(wksp, n);
    PatchMath::Workspace::ComplexVector beta(wksp, n);
    PatchMath::Workspace::RealVector rwork(wksp, n);
    PatchMath::Workspace::ComplexMatrix Acopy(wksp, n, n);
    PatchMath::Workspace::ComplexMatrix Bcopy(wksp, n, n);


    Lapack_int info = 89304;

    
    // Factor B = QR.

    zgeqrf_(&n, &n, B.make_fortran_complex(), &n, tau.make_fortran_complex(),
      cwksp.make_fortran_complex(), &lwork, &info);
    if (info != 0)
      throw_error("zgeqrf failure");

    // Apply the HH transforms to A also.

    zunmqr_("L", "C", &n, &n, &n, B.make_fortran_complex(),
      &n, tau.make_fortran_complex(), A.make_fortran_complex(),
      &n, cwksp.make_fortran_complex(), &lwork, &info);

    if (info != 0)
      throw_error("zgunmqr failure");


    // Zero out B below the diagonal.

    {
      for (int j = 0; j < n; ++j)
        for (int i = j + 1; i < n; ++i)
          B(i,j) = Complex(0.0,0.0);
    }


    // Reduce to Hessenberg form.  Save Z.


    Lapack_int one = 1;


    PatchMath::Workspace::ComplexMatrix Z(wksp, n, n);
    zgghrd_("N", "I", &n, &one, &n, A.make_fortran_complex(),
      &n, B.make_fortran_complex(), &n, 0, &one, Z.make_fortran_complex(),
      &n, &info);

    if (info != 0)
      throw_error("zgghrd failure");

    // Make copies of A and B
    
    {
      Complex* A1 = &A(0,0);
      Complex* Acopy1 = &Acopy(0,0);
      for (int j = 0; j < n*n; ++j)
        Acopy1[j] = A1[j];
    }

    
    {
      Complex* B1 = &B(0,0);
      Complex* Bcopy1 = &Bcopy(0,0);
      for (int j = 0; j < n*n; ++j)
        Bcopy1[j] = B1[j];
    }


    // QZ iteration on A and B.


    zhgeqz_("S", "N", "V", &n, &one, &n, A.make_fortran_complex(),
      &n, B.make_fortran_complex(), &n, alpha.make_fortran_complex(),
      beta.make_fortran_complex(), 0, &one, Z.make_fortran_complex(),
      &n, cwksp.make_fortran_complex(), &lwork, rwork.make_fortran_real(), &info);

    if (info != 0) {
#ifdef DEBUGGING
      *debug_str << "%%%% zhgeqz_ FAILURE \n input_A = [\n";
      {
        for (int i = 0; i < n; ++i) {
          {
            for (int j = 0; j < n; ++j) {
              *debug_str << Format_Complex(Acopy(i,j));
              if (j < n - 1)
                *debug_str << ",";
            }
          }
          *debug_str << '\n';
        }
      }
      *debug_str << "];\n input_B = [\n";
      {
        for (int i = 0; i < n; ++i) {
          {
            for (int j = 0; j < n; ++j) {
              *debug_str << Format_Complex(Bcopy(i,j));
              if (j < n - 1)
                *debug_str << ",";
            }
          }
          *debug_str << '\n';
        }
      }

      *debug_str << "];\n output_A = [\n";
      {
        for (int i = 0; i < n; ++i) {
          {
            for (int j = 0; j < n; ++j) {
              *debug_str << Format_Complex(A(i,j));
              if (j < n - 1)
                *debug_str << ",";
            }
          }
          *debug_str << '\n';
        }
      }
      *debug_str << "];\n output_B = [\n";
      {
        for (int i = 0; i < n; ++i) {
          {
            for (int j = 0; j < n; ++j) {
              *debug_str << Format_Complex(A(i,j));
              if (j < n - 1)
                *debug_str << ",";
            }
          }
          *debug_str << '\n';
        }
      }
      *debug_str << "];\n";
#endif
      return -1;
    }

#ifdef DEBUGGING
    if (dump_everything) {
      *debug_str << "% After qz factorization\n AA = [\n";
      {
        for (int i = 0; i < n; ++i) {
          {
            for (int j = 0; j < i; ++j)
              *debug_str << "0,";
          }
          {
            for (int j = i; j < n; ++j) {
              *debug_str << Format_Complex(A(i,j));
              if (j < n - 1)
                *debug_str << ",";
            }
          }
          *debug_str << '\n';
        }
      }

      *debug_str << "];\n% After qz factorization\n BB = [\n";
      {
        for (int i = 0; i < n; ++i) {
          {
            for (int j = 0; j < i; ++j)
              *debug_str << "0,";
          }
          {
            for (int j = i; j < n; ++j) {
              *debug_str << Format_Complex(B(i,j));
              if (j < n - 1)
                *debug_str << ",";
            }
          }
          *debug_str << '\n';
        }
      }

      *debug_str << "];\n% After qz factorization\n Z = [\n";
      {
        for (int i = 0; i < n; ++i) {
          for (int j = 0; j < n; ++j) {
            *debug_str << Format_Complex(Z(i,j));
            if (j < n - 1)
              *debug_str << ",";
          }
          *debug_str << '\n';
        }
      }
      *debug_str << "];\n";
    }
#endif



    // Solve for right eigenvectors.  Estimate condition numbers
    // of eigenvectors (nonstandard definition of eigenvector
    // condition!) for eigenvectors whose eigenvalue is between
    // lb and ub.


    PatchMath::Workspace::ComplexMatrix T(wksp, n, n);
    PatchMath::Workspace::ComplexVector x(wksp, n - 1);
    PatchMath::Workspace::ComplexVector y(wksp, n - 1);
    PatchMath::Workspace::ComplexVector xi(wksp, n - 1);
    PatchMath::Workspace::ComplexVector z(wksp, n - 1);
    PatchMath::Workspace::ComplexMatrix savegiv(wksp, n, 2);
    PatchMath::Workspace::ComplexVector eigenvec(wksp, n);

    // Use C-style pointers to vectors for better
    // performance in forward/backsolve.

    Complex* xi1 = &xi[0];
    Complex* x1 = &x[0];
    Complex* y1 = &y[0];
    Complex* z1 = &z[0];


    using std::conj;
    using std::polar;
    using std::arg;
    using std::abs;

    for (int eigidx = 0; eigidx < n; ++eigidx) {
      Complex beta1 = B(eigidx, eigidx);
      Complex alpha1 = A(eigidx, eigidx);

      lambda[eigidx] = (abs(beta1) == 0.0)? Complex(1.0/MACHINE_EPS,0) : alpha1/beta1;

      // Set T = B(i,i)*A - A(i,i) * B.
      {
        for (int j = 0; j < n; ++j)
          for (int i = 0; i <= j; ++i)
            T(i,j) = beta1 * A(i,j) - alpha1 * B(i,j);
      }
      // Zero out diagonal entries below the pivot using Givens rotations.

      bool T_singular = false;

      {
        for (int p = eigidx; p < n - 1; ++p) {
          
          // Compute complex Givens rotation.
          Complex c, s, cbar, sbar;
          compute_givens(T(p,p+1),T(p+1,p+1), c, s, cbar, sbar);
          // Apply Givens rotations.
          {
            for (int j = p+1; j < n; ++j) {
              Complex tmp = T(p,j);
              T(p,j) = tmp * c + T(p+1,j) * s;
              T(p+1,j) = T(p+1,j) * cbar - tmp * sbar;
            }
          }
          if (T(p,p+1) == Complex(0.0,0.0))
            T_singular = true;
        }
      }
      
      // Zero out diagonal entries above the pivot using Givens rotations.
      // Save the rotations.
      
      {
        for (int p = eigidx - 1; p >= 0; --p) {

          // Compute Givens rotation.

          Complex c, s, cbar, sbar;
          compute_givens(T(p,p+1), T(p,p), c, s, cbar, sbar);
          savegiv(p,0) = c;
          savegiv(p,1) = s;
          // Apply Givens rotations.
          {
            for (int i = 0; i <= p; ++i) {
              Complex tmp = T(i,p+1);
              T(i,p+1) = tmp * c + T(i,p) * s;
              T(i,p) = T(i,p) * cbar - tmp * sbar;
            }
          }
          if (T(p,p+1) == Complex(0.0,0.0))
            T_singular = true;
        }
      }

      // Apply the transposed Givens rotations to [1;0;...0] to get the eigenvector.
      { 
        for (int i = 1; i < n; ++i)
          eigenvec[i] = Complex(0.0,0.0);
        eigenvec[0] = Complex(1.0,0.0);
      }


      {
        for (int p = 0; p < eigidx; ++p) {
          Complex c = savegiv(p,0);
          Complex s = savegiv(p,1);
          Complex tmp = eigenvec[p+1];
          eigenvec[p+1] = tmp * c - eigenvec[p] * conj(s);
          eigenvec[p] = eigenvec[p] * conj(c) + tmp * s;
        }
      }

      // Multiply Z times the eigenvector; store the answer as a column of V.

      {
        for (int i = 0; i < n; ++i) {
          Complex t(0.0, 0.0);
          for (int j = 0; j < n; ++j) {
            t += Z(i,j) * eigenvec[j];
          }
          V(i,eigidx) = t;
        }
      }
     



#ifdef DEBUGGING
      if (dump_everything) {
        *debug_str << "eigenvec(:," << eigidx+1 << ")= [\n";
        for (int j = 0;  j < n; ++j) {
          *debug_str << Format_Complex(V(j,eigidx)) << '\n';
        }
        *debug_str << "];\n";
        if (eigidx == 1) {
          *debug_str << "reducedT = [\n ";
          for (int i = 0; i < n; ++i) {
            for (int j = 0; j < n; ++j) {
              if (j < i)
                *debug_str << "0 ";
              else
                *debug_str << Format_Complex(T(i,j)) << " ";
            }
            *debug_str << '\n';
          }
          *debug_str << "];\n";
        }
      }
      
#endif

      // Compute the condition number of the reduced matrix T to get a
      // condition estimation.  Use Hager's/Higham's method.

      Real conde = 0.0;

      if (T_singular) {
        conde = 1.0 /MACHINE_EPS;
      }
      else {
        const int HAGER_TRIAL = 2;
        const int MAX_HAGER_IT = 6;
        
        for (int tr = 0; tr < HAGER_TRIAL; ++tr) {
          // Select a 'random' x with 1-norm equal to 1. 
          Real nrm = 0.0;
          {
            for (int j = 0; j < n - 1; ++j) {
              x[j] = Complex(cos(static_cast<Real>(j * 17 * tr)), 0.0);
              nrm += abs(x[j]);
            }
          }
          {
            for (int j = 0; j < n - 1; ++j) {
              x[j] /= nrm;
            }
          }
          
          for (int it = 0; true; ++it) {
            
            // Backsolve: y = T \ x.
            // Use saxpy and C-style pointers for better performance.
            {
              for (int i = 0; i < n - 1; ++i)
                y1[i] = x1[i];
            }
            {
              for (int j = n - 2; j >= 0; --j) {
                const Complex* thiscol = &T(0,j + 1);
                y1[j] /= thiscol[j];
                Complex t = y1[j];
                for (int i = 0; i < j; ++i)
                  y1[i] -= thiscol[i] * t;
              }
            }
            
            // xi = sgn(y)
            
            {
              for (int i = 0; i < n - 1; ++i)
                if (y[i]== Complex(0.0, 0.0))
                  xi[i] = Complex(1.0,0.0);
                else
                  xi[i] = polar(1.0, arg(y[i]));
            }
            // Forward sub: z = T' \ xi
            // Use innerprod algorithm and C-style pointers.
            {
              for (int i = 0; i < n - 1; ++i) {
                const Complex* thiscol = &T(0,i + 1);
                Complex t = xi1[i];
                for (int j = 0; j < i; ++j) 
                  t -= conj(thiscol[j]) * z1[j];
                z1[i] = t / conj(thiscol[i]);
              }
            }
            Real normzinf = -1.0;
            Complex ipxz(0.0, 0.0);
            int maxp;
            {
              for (int i = 0; i < n - 1; ++i) {
                Real e1 = abs(z[i]);
                if (e1 > normzinf) {
                  normzinf = e1;
                  maxp = i;
                }
                ipxz += conj(z[i]) * x[i];
              }
            }
            if (abs(ipxz) >= normzinf || it == MAX_HAGER_IT) {
              Real normy = 0.0;
              for (int i = 0; i < n - 1; ++i)
                normy += abs(y[i]);
#ifdef DEBUGGING
              if (dump_everything) {
                *debug_str << "normest = " << normy << ";\n";
              }
#endif
              if (normy >= conde)
                conde = normy;
              break;
            }
            {
              for (int i = 0; i < n - 1; ++i)
                x[i] = Complex(0.0,0.0);
            }
            x[maxp] = Complex(1.0,0.0);
          }
        }
      }
        
      eveccond[eigidx] = conde;
    }
    return 0;
  }


        
  // This is a helper class for math on 
  // a one-variable Bezier.

  class Solve_OneVar_Bezier {
  public:
    virtual Real get_bez_coeff(int j) const = 0;
    Solve_OneVar_Bezier() { }
    virtual ~Solve_OneVar_Bezier() { }
    void solve(int degree,
      vector<PatchMath::IntersectRecord>& output_stack,
      PatchMath::Workspace& workspace, 
      double scfac, 
      double tol) const;
    Real eval(Real param,
      int degree) const;
    Real eval_deriv(Real param,
      int degree) const;
  };

  // Derived class for math on a one-variable Bezier that arises
  // from intersecting a line with a curve in 2D.

  class Solve_OneVar_Bezier1 : public Solve_OneVar_Bezier {
  private: 
    const Matrix& lineeq_;
    const Point& linerhs_;
    const PatchMath& pm_;
  public:
    Solve_OneVar_Bezier1(const Matrix& lineeq,
      const Point& linerhs,
      const PatchMath& pm) : lineeq_(lineeq), linerhs_(linerhs), pm_(pm) { }
    Real get_bez_coeff(int cpnum) const {
      return lineeq_(0,0) * PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum, 0) +
        lineeq_(0,1) * PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum, 1) -
        linerhs_[0];
    }
  };

  // Derived class for math on a one-variable Bezier that
  // arises from intersecting a curve with a plane in 3D.

  class Solve_OneVar_Bezier2 : public Solve_OneVar_Bezier {
  private:
    const Point& plane_normal_;
    Real plane_rhs_;
    const PatchMath& pm_;
  public:
    Solve_OneVar_Bezier2(const Point& plane_normal,
      Real plane_rhs,
      const PatchMath& pm) : plane_normal_(plane_normal), plane_rhs_(plane_rhs), 
      pm_(pm) { }
    Real get_bez_coeff(int cpnum) const {
      return plane_normal_[0] * PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum, 0) +
        plane_normal_[1] * PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum,1) +
        plane_normal_[2] * PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum,2) -
        plane_rhs_;
    }
  };

  // Derived class for evaluating one dimension of a Bezier curve.

  class Solve_OneVar_Bezier3 : public Solve_OneVar_Bezier {
  private:
    int select_dim_;
    const PatchMath& pm_;
  public:
    Solve_OneVar_Bezier3(int select_dim, 
      const PatchMath& pm) : select_dim_(select_dim), pm_(pm) 
    { }
    Real get_bez_coeff(int cpnum) const {
      return PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum, select_dim_);
    }
  };


  // Derived class for evaluating the derivative of a Bezier curve.


  class Solve_OneVar_Bezier4 : public Solve_OneVar_Bezier {
  private:
    const Point& direc_;
    int di_;
    const PatchMath& pm_;
  public:
    Solve_OneVar_Bezier4(const Point& direc, int di,
      const PatchMath& pm) : direc_(direc), di_(di), pm_(pm) 
    { }
    Real get_bez_coeff(int cpnum) const {
      Real t = 0.0;
      for (int j = 0; j < di_; ++j) {
        t += (PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum + 1, j) -
          PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum, j)) *
          direc_[j];
      }
      return t;
    }
  };


  // Evaluate a one-variable Bezier curve using Horner algorithm.
  

  Real Solve_OneVar_Bezier::eval(Real param,
    int degree) const {
    Real u = 1.0 - param;
    Real pwrv = 1.0; 
    Real result = 0.0;
    int binom = 1;
    for (int j = 0; j < degree + 1; ++j) {
      result *= u;
      result += binom * pwrv * get_bez_coeff(j);      
      pwrv *= param;
      binom *= (degree - j);
      binom /= (j + 1);
    }
    return result;
  }
  
  // Evaluate derivative a one-variable Bezier curve using Horner algorithm.
  

  Real Solve_OneVar_Bezier::eval_deriv(Real param,
    int degree) const {
    Real u = 1.0 - param;
    Real pwrv = 1.0; 
    Real result = 0.0;
    int binom = 1;
    Real prevbez = get_bez_coeff(0);
    for (int j = 0; j < degree; ++j) {
      result *= u;
      Real nextbez = get_bez_coeff(j + 1);
      result += binom * pwrv * (nextbez - prevbez);
      prevbez = nextbez;
      pwrv *= param;
      binom *= (degree - j - 1);
      binom /= (j + 1);
    }
    return result * degree;
  }


  // solve: find zeros of a one-variable Bezier.
  // Transforms to an eigenvalue problem and invokes LAPACK.
  // Doesn't set output_stack real coord fields -- caller must do this.
 
  void Solve_OneVar_Bezier::solve(int degree,
    vector<PatchMath::IntersectRecord>& output_stack,
    PatchMath::Workspace& workspace, 
    double scfac, 
    double tol) const {


#ifdef DEBUGGING
    using QMG::Local_to_PatchMath_CPP::debug_str;
    using QMG::Local_to_PatchMath_CPP::dump_everything;

    if (dump_everything) {
      *debug_str << " In solve bezier degree = " << degree 
        << " scfac = " << scfac << " tol = " << tol << "\n bezcoeff = [";
      for (int j = 0; j < degree + 1; ++j) 
        *debug_str << std::setprecision(17) << get_bez_coeff(j) << " ";
      *debug_str << "]\n";
    }
#endif
    if (degree == 1) {

      // In the linear case just solve the one-variable equation.

      Real denom = get_bez_coeff(0) - get_bez_coeff(1);
      if (denom == 0) {
        return;
      }
      Real sol = get_bez_coeff(0) / denom;
      Real uncer = tol * scfac / fabs(denom);
      if (sol < -uncer || sol > 1.0 + uncer)
        return;
      PatchMath::IntersectRecord result;
      result.paramcoord[0] = sol;
      result.param_uncertainty = uncer;
      result.degen = fabs(sol) <= uncer || fabs(sol - 1.0) <= uncer;
      output_stack.push_back(result);
      return;
    }
    else {
      using std::abs;
      using std::setprecision;
      PatchMath::Workspace::StartPosMarker startposmarker_(workspace);
      
      int degree1 = degree + 1;

      // Multiply by binomial coefs.

      PatchMath::Workspace::RealVector bvec(workspace, degree1);
     
      Real bezmax = 0;
      {
        int binom = 1;
        for (int i = 0; i < degree + 1; ++i) {
          bvec[i] = get_bez_coeff(i) * binom;
          binom *= (degree - i);
          binom /= (i + 1);
          if (fabs(bvec[i]) > bezmax)
            bezmax = fabs(bvec[i]);
        }
      }

      if (bezmax == 0) return;

      Real sctol = scfac / bezmax;
      if (sctol < 1.0)
        sctol = 1.0;

      
      // grab space in the workspace.
      
      Lapack_int lapack_cworkspace_size = degree * 8 + 
        ((degree < 64)? 64 : degree) * degree + 256;
      
      PatchMath::Workspace::ComplexVector powerv(workspace, degree1);
      PatchMath::Workspace::ComplexVector shiftpoly(workspace, degree1);
      PatchMath::Workspace::ComplexMatrix A(workspace, degree, degree);
      PatchMath::Workspace::ComplexMatrix A0(workspace, degree, degree);
      PatchMath::Workspace::ComplexVector eigs(workspace, degree);
      PatchMath::Workspace::RealVector 
        lapack_rworkspace(workspace, degree * 2);
      PatchMath::Workspace::ComplexVector lapack_cworkspace(workspace, 
        lapack_cworkspace_size);

      // Make a substitution [u;v] = [1, zs ; -conj(zs), 1] * [u'; v']
      // where a is a random complex number.  This is to make sure
      // the leading coef does not vanish.

      Complex zs;
      Complex conjzs;

      const int NUM_SHIFT_TRY = 30; 

      // Random numbers:

      Real shift_try_r[NUM_SHIFT_TRY] = {      
       -0.38447962916636347,
       -2.11814173148805510,
       -0.79520541715653292,
        2.24417573027436390,
       -0.27771464583718358,
        0.78432463643823636,
        0.12363602396816228,
       -0.49993560972966694,
        0.42383056156273607,
       -0.43734187908789357,
        0.14121474285739941,
        0.19504210910863520,
       -1.11493224993303720,
       -0.48146184807137726,
        1.27000255609278660,
        0.33380977209691776,
        0.70210522613535886,
        2.24583720987292910,
        0.64112905641480467,
       -0.76345380712690925,
        0.67204244835124261,
        0.31748025174110595,
       -0.14116397416090734,
        0.52010663229446119,
        0.48698772195500972,
       -0.00463937484632777,
        0.93132783995359203,
        0.16364139439790015,
        2.95373272264897270,
        1.34077075576002240};

 
      Real shift_try_i[NUM_SHIFT_TRY] = { 
       -0.21055404442240019,
       -1.21741933580593490,
        0.69369406932457345,
       -0.56288076417876365,
       -1.19213020429639020,
       -0.09169712921712286,
        0.54552052527194383,
       -0.25357824432146714,
        1.05656634425974280,
        0.03185761355971936,
        0.03147825098886285,
       -0.18952988332153939,
       -0.69275829345894058,
       -0.06717677344855039,
        1.37267097105617640,
       -0.51058644303146949,
        0.98979840686694620,
       -0.81148252620916717,
        0.15095047247765064,
        0.43669407222819562,
        0.64430185026083198,
       -0.31044769697761937,
       -0.32251973074141427,
        0.32345159861044098,
       -1.92435501248309440,
       -1.46073346575500910,
       -0.40130019911383691,
       -0.12211218475088281,
        0.76970408299206805,
        1.75167038966833900};


      for (int trycount = 0; true; ++trycount) {
        if (trycount == NUM_SHIFT_TRY)
          throw_error("Repeated failure of solve onevar bez");


        // Choose the shift.

        zs = Complex(shift_try_r[trycount], shift_try_i[trycount]);
        conjzs = Complex(zs.real(), -zs.imag());


#ifdef DEBUGGING
        if (dump_everything) {
          *debug_str << "trycount = " << trycount 
                     << "\n zs = " << Format_Complex(zs)
                     << "\n conjzs = " << Format_Complex(conjzs)
                     << '\n';
        }
#endif

        // To make the substitution, use a variant of Horner's rule
        //   powerv is a polynomial that holds ( -conj(zs)*u'+ v') ^k
        // for increasing powers of k.
        //   shifpoly is the result of the substitution.

        {
          for (int j = 0; j < degree1; ++j)
            powerv[j] = Complex(0.0, 0.0);
          powerv[0] = Complex(1.0, 0.0);
        }

        {
          for (int j = 0; j < degree1; ++j)
            shiftpoly[j] = Complex(0.0, 0.0);
          shiftpoly[0] = Complex(bvec[degree], 0.0);
        }
        // Now repeatedly multiply shiftpoly by u, and then add v^k, for
        // increasing values of k.

        {
          for (int k = 0; k < degree; ++k) {
            
            // Multiply shiftpoly by u, which is u' + zs*v'
            {
              for (int l = k + 1; l > 0; --l) {
                shiftpoly[l] = shiftpoly[l - 1] + shiftpoly[l] * zs;
              }
              shiftpoly[0] = shiftpoly[0] * zs;
            }
          

            // Multiply powerv by v, which is  -conj(zs)  * u' + v'

            {
              for (int l = k + 1; l > 0; --l) {
                powerv[l] -= powerv[l - 1] * conjzs;
              }
            }
            // Add bez[k+1] * powerv to shiftpoly.

            {
              for (int l = 0; l <= k + 1; ++l) {
                shiftpoly[l] += bvec[degree - k - 1] * powerv[l];
              }
            }
          }
        }

        // If the leading coef of shiftpoly is too small, must try again.

        Complex ldgcoef = shiftpoly[degree];

#ifdef DEBUGGING
        if (dump_everything) {
          *debug_str << "after loop\npowerv = [\n";
          for (int l = 0; l <= degree; ++l)
            *debug_str << " " << Format_Complex(powerv[l]) << '\n';
          *debug_str << "]\nshiftpoly = [\n";
          for (int l2 = 0; l2 <= degree; ++l2)
            *debug_str << " " << Format_Complex(shiftpoly[l2]) << '\n';
          *debug_str << "]\n bezmax = " << bezmax 
                     << "\n ldgcoef = " << Format_Complex(ldgcoef)
                     << '\n';
        }
#endif

        if (abs(ldgcoef) < bezmax / (50 * degree)) continue;

        // Set up the companion matrix for shiftpoly.
        {
          for (int i = 0; i < degree - 1; ++i) {
            for (int j = 0; j < degree; ++j) {
              A(i,j) = Complex(0.0,0.0);
              A0(i,j) = Complex(0.0,0.0);
            }
            A(i,i+1) = Complex(1.0,0.0);
            A0(i,i+1) = Complex(1.0,0.0);
          }
        }
        {
          for (int j = 0; j < degree; ++j) {
            A(degree - 1, j) = -shiftpoly[j] / ldgcoef;
            A0(degree - 1, j) = -shiftpoly[j] / ldgcoef;
          }
        }



#ifdef DEBUGGING
        if (dump_everything) {
          *debug_str << "A = [";
          for (int i = 0; i < degree; ++i) {
            for (int j = 0; j < degree; ++j) {
              *debug_str << "  " << Format_Complex(A(i,j));
            }
            *debug_str << '\n';
          }
          *debug_str << "]\n";
        }
#endif


        // Call lapack.

        Lapack_int lapack_info = 100000;
        Lapack_int degree_l = degree;
        Lapack_int lda = degree;
        Lapack_int ldvl = degree;
        zgeev_("N", "N", &degree_l, A.make_fortran_complex(), &lda,
          eigs.make_fortran_complex(), 0, &degree_l, 0, &ldvl, 
          lapack_cworkspace.make_fortran_complex(), 
          &lapack_cworkspace_size, 
          lapack_rworkspace.make_fortran_real(), 
          &lapack_info);
        // done with loop if lapack successful.

        if (lapack_info == 0)  break;

        // If lapack failure, print some info.
        
#ifdef DEBUGGING
        *debug_str << "Lapack zgeev failure code = " << lapack_info << '\n';
        *debug_str << "initial A =\n";
        {
          for (int i = 0; i < degree; ++i) {
            for (int j = 0; j < degree; ++j) {
              *debug_str << Format_Complex(A0(i,j)) << '\n';
            }
          }
        }
        *debug_str << "final A =\n";
        {
          for (int i = 0; i < degree; ++i) {
            for (int j = 0; j < degree; ++j) {
              *debug_str << Format_Complex(A(i,j)) << '\n';
            }
          }
        }
#endif
      }

      // Eigenvals successfully computed and stored in eig.  Now
      // output the roots in [0,1].


      for (int rootidx = 0; rootidx < degree; ++rootidx) {
        Complex sol0 = eigs[rootidx];
        // Apply i transformation to [sol0;1] to get u,v.
        Complex u = sol0 + zs;
        Complex v = -conjzs * sol0 + Complex(1.0,0.0);

        // If way outside the unit interval, skip it.

        if (abs(u) >= 3 * abs(u + v)) continue;

        Complex sol = u / (u + v);
        Complex sol2 = Complex(1.0,0.0) - sol;

#ifdef DEBUGGING
        if (dump_everything)
          *debug_str << " sol0 = " << Format_Complex(sol0)
                     << "\n u = " << Format_Complex(u) 
                     <<  "\n v = " << Format_Complex(v) 
                     << "\n sol = " << Format_Complex(sol) 
                     << "\n sol2 = " << Format_Complex(sol2) 
                     << '\n';
#endif



        // Compute the residual, derivative, and
        // absolute evaluation of the polynomial at
        // sol using Horners rule.  This is for
        // estimating the uncertainty in the ansewr.

        Complex pwrv(1.0,0.0);
        Complex resid(0.0,0.0);
        Complex deriv(0.0,0.0);
        {
          Real prevbez = get_bez_coeff(0);
          Real abssol2 = abs(sol2);
          int binom = 1;
          int binom2 = 1;
          for (int j = 0; j < degree + 1; ++j) {
            resid *= sol2;
            resid += pwrv * get_bez_coeff(j) * static_cast<Real>(binom);
            if (j < degree) {
              deriv *= sol2;
              Real nextbez = get_bez_coeff(j + 1);
              deriv += pwrv * (nextbez - prevbez) * static_cast<Real>(binom2);
              prevbez = nextbez;
              binom2 *= (degree - j - 1);
              binom2 /= (j + 1);
            }
            pwrv *= sol;
            binom *= (degree - j);
            binom /= (j + 1);
          }
          deriv *= static_cast<Real>(degree);
        }
        Real trueres_bound = abs(resid) + sctol * tol;
        Real absderiv = abs(deriv);
        Real param_uncertainty = (absderiv == 0.0)? 1.0 :
          trueres_bound / absderiv;

        // Check if sol is within param_uncertainty_width of the unit interval.

        Real solreal = sol.real();
        Real solimag = sol.imag();
#ifdef DEBUGGING
        if (dump_everything) {
          *debug_str << "deriv = " << Format_Complex(deriv) << "absderiv = " << absderiv
            << "\n resid = " << Format_Complex(resid) << "trueres_bound = " << trueres_bound
            << "\n param_uncertainty = " << param_uncertainty << '\n';
        }
#endif

        Real delta1 = param_uncertainty;
        if (delta1 > 0.01)
          delta1 = 0.01;


        if (solreal < -delta1 || solreal > 1.0 + delta1
          || fabs(solimag) > delta1) continue;
        PatchMath::IntersectRecord result;
        result.paramcoord[0] = solreal;
        result.degen = (fabs(solreal) <= param_uncertainty || 
          fabs(solreal - 1.0) <= param_uncertainty);
#ifdef DEBUGGING
        if (dump_everything) 
          *debug_str << "keeping, degen = " << result.degen << "\n";
#endif
        result.param_uncertainty = param_uncertainty;
        output_stack.push_back(result);
      }
    }
  }

  // A point in R^2, used for representing Bezier parameters. 
  
  struct R2Point {
    Real coord[2];
    R2Point() { }
    R2Point(Real x, Real y) {
      coord[0] = x;
      coord[1] = y;
    }
    Real l2norm() const {
      return sqrt(coord[0] * coord[0] + coord[1] * coord[1]);
    }
  };

  // A point in C^2, used for representing complex Bezier parameters.

  struct C2Point {
    Complex coord[2];
    C2Point() { }
    C2Point(Complex x, Complex y) {
      coord[0] = x;
      coord[1] = y;
    }
    Real l2norm() const {
      return sqrt(std::norm(coord[0]) + std::norm(coord[1]));
    }
  };

  /*


  // Compute real coord from param coord using Horner's algorithm
  // of a triangular Bezier polynomial, except without binomial coefs.

  Real pl_evalr(const PatchMath::Workspace::RealMatrix& coefs,
    const R2Point& param,
    int eqnum,
    int degree) {
    Real p0 = param.coord[0];
    Real p1 = param.coord[1];
    Real p2 = 1.0 - p0 - p1;
    Real pwrv = 1.0;
    Real result = 0.0;
    for (int row = degree; row >= 0; --row) {
      int base_idx = row * (row + 1) / 2;
      Real pwru = 1.0;
      Real row_result = 0.0;
      for (int col = 0; col <= row; ++col) {
        row_result *= p2;
        row_result += pwru * coefs(eqnum, base_idx + col);
        pwru *= p0;
      }
      result += row_result * pwrv;
      pwrv *= p1;
    }
    return result;
  }
  
  // Compute real coord from param coord using Horner's algorithm
  // of a triangular Bezier polynomial, except without binomial weights.

  Complex pl_eval(const PatchMath::Workspace::RealMatrix& coefs,
    const C2Point& param,
    int eqnum,
    int degree) {
    Complex p0 = param.coord[0];
    Complex p1 = param.coord[1];
    Complex p2 = Complex(1.0,0.0) - p0 - p1;
    Complex pwrv(1.0,0.0);
    Complex result(0.0,0.0);
    for (int row = degree; row >= 0; --row) {
      int base_idx = row * (row + 1) / 2;
      Complex pwru(1.0, 0.0);
      Complex row_result(0.0, 0.0);
      for (int col = 0; col <= row; ++col) {
        row_result *= p2;
        row_result += pwru * coefs(eqnum, base_idx + col);
        pwru *= p0;
      }
      result += row_result * pwrv;
      pwrv *= p1;
    }
    return result;
  }

  */


  
  // Compute real coord from param coord using Horner's algorithm
  // of a triangular polynomial (not bezier).

  Real pl_evalr(const PatchMath::Workspace::RealMatrix& coefs,
    const R2Point& param,
    int eqnum,
    int degree) {
    Real p0 = param.coord[0];
    Real p1 = param.coord[1];
    Real pwrv = 1.0;
    Real result = 0.0;
    for (int row = degree; row >= 0; --row) {
      int base_idx = row * (row + 1) / 2;
      Real pwru = 1.0;
      Real row_result = 0.0;
      for (int col = 0; col <= row; ++col) {
        row_result += pwru * coefs(eqnum, base_idx + col);
        pwru *= p0;
      }
      result += row_result * pwrv;
      pwrv *= p1;
    }
    return result;
  }
  
  // Compute real coord from param coord using Horner's algorithm
  // of a triangular polynomial (not Bezier)

  Complex pl_eval(const PatchMath::Workspace::RealMatrix& coefs,
    const C2Point& param,
    int eqnum,
    int degree) {
    Complex p0 = param.coord[0];
    Complex p1 = param.coord[1];
    Complex pwrv(1.0,0.0);
    Complex result(0.0,0.0);
    for (int row = degree; row >= 0; --row) {
      int base_idx = row * (row + 1) / 2;
      Complex pwru(1.0, 0.0);
      Complex row_result(0.0, 0.0);
      for (int col = 0; col <= row; ++col) {
        row_result += pwru * coefs(eqnum, base_idx + col);
        pwru *= p0;
      }
      result += row_result * pwrv;
      pwrv *= p1;
    }
    return result;
  }



  // Class for math with bezier triangles.

  class Solve_BezTri {
  public:
    virtual Real get_bez_coeff(int eqnum, int j) const = 0;
    Solve_BezTri() { }
    virtual ~Solve_BezTri() { }



    bool solve(int degree,
      vector<PatchMath::IntersectRecord>& output_stack,
      PatchMath::Workspace& workspace,
      double scfac, 
      double tol,
      bool newton_improve) const;
    Real eval(const R2Point& param,
      int eqnum,
      int degree) const;
    Real eval_direc_deriv(const R2Point& param,
      const R2Point& direc,
      int eqnum,
      int degree) const;
    Complex eval_complex(const C2Point& param,
      int eqnum,
      int degree) const;
    Complex eval_direc_deriv_complex(const C2Point& param,
      const C2Point& direc,
      int eqnum,
      int degree) const;
  };

  // Derived class for intersecting a line with a triangle.

  class Solve_BezTri1 : public Solve_BezTri {
  private: 
    const Matrix& lineeq_;
    const Point& linerhs_;
    const PatchMath& pm_;
  public:
    Solve_BezTri1(const Matrix& lineeq,
      const Point& linerhs,
      const PatchMath& pm) : lineeq_(lineeq), linerhs_(linerhs), pm_(pm) { }
    Real get_bez_coeff(int eqnum, int cpnum) const {
      return lineeq_(eqnum,0) * PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum, 0) +
        lineeq_(eqnum,1) * PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum, 1) +
        lineeq_(eqnum,2) * PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum, 2) -
        linerhs_[eqnum];
    }
  };

  // Derived class for evaluating directional deriv.

  class Solve_BezTri2 : public Solve_BezTri {
  private:
    const Point& direc_;
    int di_;
    const PatchMath& pm_;
  public:
    Solve_BezTri2(const Point& direc, 
      int di,
      const PatchMath& pm) : direc_(direc), di_(di), pm_(pm) 
    { }
    Real get_bez_coeff(int eqnum, int cpnum) const;
  };
  
  

  // Derived class for evaluating one coordinate of a bezier triangle.

  class Solve_BezTri3 : public Solve_BezTri {
  private:
    int select_dim_;
    const PatchMath& pm_;
  public:
    Solve_BezTri3(int select_dim, 
      const PatchMath& pm) : select_dim_(select_dim), pm_(pm) 
    { }
    Real get_bez_coeff(int, int cpnum) const {
      return PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum, select_dim_);
    }
  };
  
  
  Real Solve_BezTri2::get_bez_coeff(int eqnum, int cpnum) const {
    Real t = 0.0;
    int offset = static_cast<int>(sqrt(2 * static_cast<double>(cpnum)));
    if (offset * (offset + 1) / 2 > cpnum)
      --offset;
    int bcp = cpnum + offset + 1;
    int ecp = (eqnum == 0)? bcp + 1 : cpnum;
    for (int j = 0; j < di_; ++j) {
      t += (PatchMath::Bezier_Helper::control_point_coord(pm_, ecp, j) -
        PatchMath::Bezier_Helper::control_point_coord(pm_, bcp, j)) *
        direc_[j];
    }
    return t;
  }
    
  // Compute real coord from param coord using Horner's algorithm.

  Real Solve_BezTri::eval(const R2Point& param,
    int eqnum,
    int degree) const {
    Real p0 = param.coord[0];
    Real p1 = param.coord[1];
    Real p2 = 1.0 - p0 - p1;
    int outer_binom = 1;
    Real pwrv = 1.0;
    Real result = 0.0;
    for (int row = degree; row >= 0; --row) {
      int base_idx = row * (row + 1) / 2;
      int inner_binom = 1;
      Real pwru = 1.0;
      Real row_result = 0.0;
      for (int col = 0; col <= row; ++col) {
        row_result *= p2;
        row_result += inner_binom * 
          get_bez_coeff(eqnum, base_idx + col) * pwru;
        pwru *= p0;
        inner_binom *= (row - col);
        inner_binom /= (col + 1);
      }
      result += row_result * pwrv * outer_binom;
      outer_binom *= row;
      outer_binom /= (degree - row + 1);
      pwrv *= p1;
    }
    return result;
  }

  // Compute complex coord from param coord using Horner's algorithm.

  Complex Solve_BezTri::eval_complex(const C2Point& param,
    int eqnum,
    int degree) const {
    Complex p0 = param.coord[0];
    Complex p1 = param.coord[1];
    Complex p2 = Complex(1.0,0.0) - p0 - p1;
    int outer_binom = 1;
    Complex pwrv(1.0,0.0);
    Complex result(0.0,0.0);
    for (int row = degree; row >= 0; --row) {
      int base_idx = row * (row + 1) / 2;
      int inner_binom = 1;
      Complex pwru(1.0,0.0);
      Complex row_result(0.0,0.0);
      for (int col = 0; col <= row; ++col) {
        row_result *= p2;
        row_result += Complex(static_cast<Real>(inner_binom), 0.0) * 
          Complex(get_bez_coeff(eqnum, base_idx + col),0.0) * pwru;
        pwru *= p0;
        inner_binom *= (row - col);
        inner_binom /= (col + 1);
      }
      result += row_result * pwrv * Complex(static_cast<Real>(outer_binom), 0);
      outer_binom *= row;
      outer_binom /= (degree - row + 1);
      pwrv *= p1;
    }
    return result;
  }
    
  // Compute directional derivative using Horner's rule.
  
  Real Solve_BezTri::eval_direc_deriv(const R2Point& param,
    const R2Point& direc,
    int eqnum,
    int degree) const {   
    Real p0 = param.coord[0];
    Real p1 = param.coord[1];
    Real p2 = 1.0 - p0 - p1;
    Real d0 = direc.coord[0];
    Real d1 = direc.coord[1];
    Real d2 = -d0 - d1;    
    int outer_binom = 1;
    Real pwrv = 1.0;
    Real result = 0.0;
    for (int row = degree - 1; row >= 0; --row) {
      int base_idx = row * (row + 1) / 2;
      int inner_binom = 1;
      Real pwru = 1.0;
      Real row_result = 0.0;
      for (int col = 0; col <= row; ++col) {
        row_result *= p2;
        row_result += inner_binom * 
          (d1 * get_bez_coeff(eqnum, base_idx + col) +
          d2 * get_bez_coeff(eqnum, base_idx + col + row + 1) +
          d0 * get_bez_coeff(eqnum, base_idx + col + row + 2))
          * pwru;
        pwru *= p0;
        inner_binom *= (row - col);
        inner_binom /= (col + 1);
      }
      result += row_result * pwrv * outer_binom;
      outer_binom *= row;
      outer_binom /= (degree - row);
      pwrv *= p1;
    }
    return degree * result;
  }

  // Compute directional derivative using Horner's rule.
  
  Complex Solve_BezTri::eval_direc_deriv_complex(const C2Point& param,
    const C2Point& direc,
    int eqnum,
    int degree) const {   
    Complex p0 = param.coord[0];
    Complex p1 = param.coord[1];
    Complex p2 = Complex(1.0,0.0) - p0 - p1;
    Complex d0 = direc.coord[0];
    Complex d1 = direc.coord[1];
    Complex d2 = -d0 - d1;    
    int outer_binom = 1;
    Complex pwrv(1.0,0.0);
    Complex result(0.0,0.0);
    for (int row = degree - 1; row >= 0; --row) {
      int base_idx = row * (row + 1) / 2;
      int inner_binom = 1;
      Complex pwru(1.0,0.0);
      Complex row_result(0.0,0.0);
      for (int col = 0; col <= row; ++col) {
        row_result *= p2;
        row_result += Complex(static_cast<Real>(inner_binom), 0.0) * 
          (d1 * get_bez_coeff(eqnum, base_idx + col) +
          d2 * get_bez_coeff(eqnum, base_idx + col + row + 1) +
          d0 * get_bez_coeff(eqnum, base_idx + col + row + 2))
          * pwru;
        pwru *= p0;
        inner_binom *= (row - col);
        inner_binom /= (col + 1);
      }
      result += row_result * pwrv * Complex(static_cast<Real>(outer_binom), 0.0);
      outer_binom *= row;
      outer_binom /= (degree - row);
      pwrv *= p1;
    }
    return result * static_cast<Real>(degree);
  }

  // Function to convert a monomial degree pair into a column number.

  inline int deg_pair_to_colnum(int u_deg, int v_deg, int total_deg) {
#ifdef RANGECHECK
    if (u_deg < 0 || v_deg < 0 || u_deg + v_deg > total_deg)
      throw_error("degrees out of range");
#endif
    return (total_deg - v_deg) * (total_deg - v_deg + 1) / 2 + u_deg;
  }

  // Inverse function -- convert a column number to a monomial degree pair.

  inline pair<int,int> colnum_to_deg_pair(int colnum, int total_deg) {
    int k = static_cast<int>(sqrt(2 * static_cast<double>(colnum)));
    if (k * (k + 1) / 2 > colnum)
      --k;
    return pair<int,int>(colnum -  k * (k + 1) / 2, total_deg - k);
  }

  
  
  Real linear_solve_C2(PatchMath::Workspace::ComplexMatrix& A,
    C2Point& rhs,
    C2Point& sol) {
    using std::abs;
    using std::norm;
    Complex det = A(0,0) * A(1,1) - A(1,0) * A(0,1);
    if (abs(det) == 0.0) {
      return BIG_REAL;
    }
    sol.coord[1] = (A(0,0) * rhs.coord[1] - A(1,0) * rhs.coord[0]) / det;
    if (abs(A(0,0)) >= abs(A(1,0))) {
      sol.coord[0] = (rhs.coord[0] - A(0,1) * sol.coord[1]) / A(0,0);
    }
    else {
      sol.coord[0] = (rhs.coord[1] - A(1,1) * sol.coord[1]) / A(1,0);
    }
    return sqrt(norm(A(0,0)) + norm(A(0,1)) + norm(A(1,0))
      + norm(A(1,1))) / abs(det);
  }

  

  // This routine applies a fractional linear transformation
  // to two triangular bezier polynomials.
  // That is, if one of the polynomials is sum_{i=0}^d sum_{j=0}^{d-i} a_{i,j} x^i*y^j,
  // then this routine computes the coefficients of the new polynomial resulting
  // from the substitution x=(a*u+b*v+c)/(g*u+h*v+i), y=(d*u+e*v+f)/(g*u+h*v+k) (after
  // clearing denominators)
  // where a,d,g,b,e,h,c,f,k (i.e., [a,b,c;d,e,f;g,h,k] transposed) 
  // are the nine given entries of the array shiftQ.
  // This routine can be used, among other things, to change
  // a bezier polynomial to a standard form polynomial.
  // In that case, the coefficients need to be premultiplied by
  // binomial factors.  The input polynomials are in coef, and
  // the output polynomials are in shiftpoly.

  void tri_poly_frac_lin_trans(int degree,
    const Real shiftQ[9],
    const PatchMath::Workspace::RealMatrix& coef,
    PatchMath::Workspace::RealMatrix& shiftpoly,
    PatchMath::Workspace& workspace) {
    
    int num_control_pt = (degree + 1) * (degree + 2) / 2;

    PatchMath::Workspace::StartPosMarker unnamed_(workspace);
    PatchMath::Workspace::RealVector pwrv(workspace, num_control_pt);
    PatchMath::Workspace::RealVector pwruv(workspace, num_control_pt);      
    PatchMath::Workspace::RealVector temp_poly(workspace, num_control_pt);    
    PatchMath::Workspace::RealMatrix shiftpoly_thisrow(workspace, 2, num_control_pt);
    
    // Use horner's rule 'symbolically';
    // Substitute Q(0,:)(u',v',w') for u, Q(1,:)(u',v',w') for v,
    // and Q(2,:)(u',v',w') for w.
    
    pwrv[0] = 1.0;
    {
      for (int i = 0; i < num_control_pt; ++i)
        for (int eqnum = 0; eqnum < 2; ++eqnum)
          shiftpoly(eqnum, i) = 0.0;
    }
    {
      for (int vdeg = 0; vdeg <= degree; ++vdeg) {
        int base_idx = (degree - vdeg) * (degree - vdeg + 1) / 2;
        {
          int vdeg_numcoef = (vdeg + 1) * (vdeg + 2) / 2;
          for (int i = 0; i < vdeg_numcoef; ++i)
            pwruv[i] = pwrv[i];
        }
        {
          for (int i = 0; i < num_control_pt; ++i)
            for (int eqnum = 0; eqnum < 2; ++eqnum)
              shiftpoly_thisrow(eqnum, i) = 0.0;
        }
        for (int udeg = 0; udeg <= degree - vdeg; ++udeg) {
          int shiftpolydeg = udeg + vdeg;

          
          // temp_poly = shiftpoly_thisrow * w symbolically
          
          
          int shiftpoly_numcoef = 
            (shiftpolydeg + 1) * (shiftpolydeg + 2) / 2;
          for (int eqnum = 0; eqnum < 2; ++eqnum) {
            {
              for (int j = 0; j < shiftpoly_numcoef; ++j) {
                pair<int,int> rval = colnum_to_deg_pair(j, shiftpolydeg);
                int ud1 = rval.first;
                int vd1 = rval.second;
                int wd1 = shiftpolydeg - ud1 - vd1;
                Real t = 0.0;
                if (ud1 > 0) {
                  int subp = deg_pair_to_colnum(ud1 - 1, vd1, shiftpolydeg - 1);
                  t += shiftpoly_thisrow(eqnum,subp) * shiftQ[2];
                }
                if (vd1 > 0) {
                  int subp = deg_pair_to_colnum(ud1, vd1 - 1, shiftpolydeg - 1);
                  t += shiftpoly_thisrow(eqnum,subp) * shiftQ[5];
                }
                if (wd1 > 0) {
                  int subp = deg_pair_to_colnum(ud1, vd1, shiftpolydeg - 1);
                  t += shiftpoly_thisrow(eqnum,subp) * shiftQ[8];
                }
                temp_poly[j] = t;
              }
            }           
            
            // Symbolic computation of 
            //  shiftpoly_thisrow = temp_poly + get_bez_coef(eqnum, base_idx + udeg) * pwruv;
            
            Real g = coef(eqnum, base_idx + udeg);
            {
              for (int j = 0; j < shiftpoly_numcoef; ++j) {
                shiftpoly_thisrow(eqnum, j) = temp_poly[j] + g * pwruv[j];
              }
            }
          }
          
          // Symbolic computation of pwruv *= u;
          
          if (shiftpolydeg < degree) {
            int shiftpolydeg1 = shiftpolydeg + 1;
            int new_numcoef = (shiftpolydeg1 + 1) * (shiftpolydeg1 + 2) / 2;
            for (int j = 0; j < new_numcoef; ++j) {
              pair<int,int> rval = colnum_to_deg_pair(j, shiftpolydeg1);
              int ud1 = rval.first;
              int vd1 = rval.second;
              int wd1 = shiftpolydeg1 - ud1 - vd1;
              Real t = 0.0;
              if (ud1 > 0) {
                int subp = deg_pair_to_colnum(ud1 - 1, vd1, shiftpolydeg);
                t += pwruv[subp] * shiftQ[0];
              }
              if (vd1 > 0) {
                int subp = deg_pair_to_colnum(ud1, vd1 - 1, shiftpolydeg);
                t += pwruv[subp] * shiftQ[3];
              }
              if (wd1 > 0) {
                int subp = deg_pair_to_colnum(ud1, vd1, shiftpolydeg);
                t += pwruv[subp] * shiftQ[6];
              }
              temp_poly[j] = t;
            }
            // Transfer result back to pwruv.
            
            for (int j2 = 0; j2 < new_numcoef; ++j2) 
              pwruv[j2] = temp_poly[j2];
          }
        }
        // Symbolic implementation of
        // shiftpoly += shiftpoly_thisrow
        
        {
          for (int j = 0; j < num_control_pt; ++j)
            for (int eqnum = 0; eqnum < 2; ++eqnum)
              shiftpoly(eqnum,j) += shiftpoly_thisrow(eqnum,j);
        }
        // Symbolic implementation of pwrv *= v;              
        
        if (vdeg < degree) {
          int vdeg1 = vdeg + 1;
          int new_numcoef = (vdeg1 + 1) * (vdeg1 + 2) / 2;
          for (int j = 0; j < new_numcoef; ++j) {
            pair<int,int> rval = colnum_to_deg_pair(j, vdeg1);
            int ud1 = rval.first;
            int vd1 = rval.second;
            int wd1 = vdeg1 - ud1 - vd1;
            Real t = 0.0;
            if (ud1 > 0) {
              int subp = deg_pair_to_colnum(ud1 - 1, vd1, vdeg);
              t += pwrv[subp] * shiftQ[1];
            }
            if (vd1 > 0) {
              int subp = deg_pair_to_colnum(ud1, vd1 - 1, vdeg);
              t += pwrv[subp] * shiftQ[4];
            }
            if (wd1 > 0) {
              int subp = deg_pair_to_colnum(ud1, vd1, vdeg);
              t += pwrv[subp] * shiftQ[7];
            }
            temp_poly[j] = t;
          }
          // Transfer result back to pwrv.
          
          for (int j2 = 0; j2 < new_numcoef; ++j2) 
            pwrv[j2] = temp_poly[j2];
        }
      }
    }
  }


  // Routine 0 for two polynomials in Bezier triangular form.
  // This routine takes a polynomial in coef and returns
  // shiftpoly, the polynomial shifted by a good shift
  // to make the leading coefficients nonzero.
  // It also returns shiftQ, the value of the shift.


  void solve_bez_tri_main0(int degree,
    const PatchMath::Workspace::RealMatrix& coef,
    PatchMath::Workspace::RealMatrix& shiftpoly,
    Real shiftQ[9],
    PatchMath::Workspace& workspace) {
    

    
    // Find a good shift.
    
    const int NUM_SHIFT_TRY = 20;
    const int LENGTH_TRY_SHIFT = NUM_SHIFT_TRY * 9;
    
    // 20 "random" 3-by-3 orthogonal matrices, listed in column major form
    // Each one is 9 consecutive entries of this array.  They have been chosen
    // so that the third column of each is close to sqrt(1/3),sqrt(1/3),sqrt(1/3)

    
    Real try_shift[LENGTH_TRY_SHIFT] = {     
      -3.9986339446162122e-001,
        7.6382965190072771e-001,
        -5.0662957735099035e-001,
        7.2155670164805497e-001,
        -7.8525628097725897e-002,
        -6.8788781937074406e-001,
        5.6521251939107797e-001,
        6.4062312525311704e-001,
        5.1974688004308345e-001,
        -4.1825926107564787e-001,
        8.0614344638659730e-001,
        -4.1855935585337523e-001,
        7.0333170126697042e-001,
        -4.1652109298009332e-003,
        -7.1084961068486163e-001,
        5.7479014302383624e-001,
        5.9170549673452388e-001,
        5.6524410356658139e-001,
        -4.1754470253415343e-001,
        8.4518400224866741e-001,
        -3.3364715453393295e-001,
        6.8435254603522711e-001,
        5.0951279808783695e-002,
        -7.2736892965052435e-001,
        5.9776083354529641e-001,
        5.3204132304621954e-001,
        5.9967826077861153e-001,
        -4.1322122997254057e-001,
        8.1075959608827297e-001,
        -4.1462886109237596e-001,
        7.0534880295277169e-001,
        -3.0171925798126442e-003,
        -7.0885397842011755e-001,
        5.7596118035252397e-001,
        5.8537148367486824e-001,
        5.7062154255444186e-001,
        -5.1303816334327279e-001,
        7.9463802706285336e-001,
        -3.2456470679822441e-001,
        6.6649222578619360e-001,
        1.3050648390060410e-001,
        -7.3400011623055428e-001,
        5.4090660554288617e-001,
        5.9288992537026419e-001,
        5.9657487415622490e-001,
        -4.4139938204522616e-001,
        7.9216459475709800e-001,
        -4.2147578856135365e-001,
        6.8324097060709998e-001,
        -7.7749871530766734e-003,
        -7.3015157718013413e-001,
        5.8167719708955445e-001,
        6.1025798183071067e-001,
        5.3780743207764214e-001,
        -3.4852485358457685e-001,
        8.0077054720415264e-001,
        -4.8713135514377592e-001,
        7.4771353416147268e-001,
        -7.5865445345176696e-002,
        -6.5967333206242174e-001,
        5.6520341229119397e-001,
        5.9414725862612561e-001,
        5.7230598267137678e-001,
        -3.9980815347937887e-001,
        8.2557180940612807e-001,
        -3.9822685485198106e-001,
        7.1220024753346789e-001,
        6.3094700785800706e-003,
        -7.0194800234817589e-001,
        5.7699588198242546e-001,
        5.6426179925741837e-001,
        5.9049502459724790e-001,
        -3.1603452538387367e-001,
        8.3314656341951498e-001,
        -4.5386009146833090e-001,
        7.5495770207820501e-001,
        -6.8883645244124769e-002,
        -6.5214562138427168e-001,
        5.7459642083656060e-001,
        5.4874570365528474e-001,
        6.0722080488211250e-001,
        -4.9177227463208201e-001,
        7.8858817058334107e-001,
        -3.6917303140831903e-001,
        7.1389245679307800e-001,
        1.2241900690711530e-001,
        -6.8947164327608057e-001,
        4.9851538595827849e-001,
        6.0261288068203045e-001,
        6.2316941998061293e-001,
        -4.6470623770685437e-001,
        8.4187825276568429e-001,
        -2.7438862978727663e-001,
        6.2464801025972061e-001,
        9.2056074942759733e-002,
        -7.7546150281281212e-001,
        6.2758503480796346e-001,
        5.3175810909315868e-001,
        5.6865660771572490e-001,
        -4.1262012930320602e-001,
        8.3240171385626938e-001,
        -3.6993515061825516e-001,
        6.8315908850609952e-001,
        1.4156105554046561e-002,
        -7.3013236092304390e-001,
        6.0252658753396204e-001,
        5.5399186947525525e-001,
        5.7450389021310455e-001,
        -4.2648059414665868e-001,
        7.6843095789576410e-001,
        -4.7710393601783618e-001,
        7.0298563506306244e-001,
        -5.0308395543764040e-002,
        -7.0942248500649807e-001,
        5.6914453323500869e-001,
        6.3795213635915560e-001,
        5.1874037051640154e-001,
        -4.6769016857136836e-001,
        7.5585587771445006e-001,
        -4.5820060928178968e-001,
        6.9393546209957990e-001,
        -7.0895782591096634e-003,
        -7.2000230021906908e-001,
        5.4746641966638254e-001,
        6.5469964869749719e-001,
        5.2119947173135761e-001,
        -3.9947369500154239e-001,
        8.1931782638932393e-001,
        -4.1126520198345096e-001,
        6.9345901185666037e-001,
        -2.3361541934381141e-002,
        -7.2011723853358234e-001,
        5.9961267988310496e-001,
        5.7286345468979083e-001,
        5.5883100880702663e-001,
        -4.4020687415549065e-001,
        7.7952095520351683e-001,
        -4.4560631542298534e-001,
        7.4275989726320368e-001,
        3.7294693410166113e-002,
        -6.6851839231317844e-001,
        5.0450534483175424e-001,
        6.2526489285904951e-001,
        5.9541428501015514e-001,
        -3.7795514074289221e-001,
        8.0257031865086748e-001,
        -4.6155259202680810e-001,
        7.5392627356590181e-001,
        -2.2542906909108040e-002,
        -6.5657215245174716e-001,
        5.3735005872620256e-001,
        5.9613144604921831e-001,
        5.9655864206146114e-001,
        -5.2456334170418273e-001,
        8.0783048438715832e-001,
        -2.6878096851330491e-001,
        6.3366002548591127e-001,
        1.5960809671463022e-001,
        -7.5696778502412343e-001,
        5.6860203362398831e-001,
        5.6739330623293627e-001,
        5.9561444190073354e-001,
        -4.5463314186292464e-001,
        8.1288382651297808e-001,
        -3.6404476498565952e-001,
        6.7501888220407247e-001,
        4.7802534365183252e-002,
        -7.3625024711454601e-001,
        5.8108365575690568e-001,
        5.8046085337590492e-001,
        5.7044454832201341e-001,
        -3.8969863538245664e-001,
        8.1931787828332880e-001,
        -4.2053916334433777e-001,
        7.3567796492208859e-001,
        2.2578150302527411e-003,
        -6.7732771551102444e-001,
        5.5399720713117784e-001,
        5.7333508230056329e-001,
        6.0363397675599839e-001
      };
    
    int bestshift;
    Real bestlead = 0.0;
    
    // Find the best shift to make the leading coefs nondegenerate.
    
    for (int shifttry = 0; shifttry < NUM_SHIFT_TRY; ++shifttry) {
      const Real* Qb = &try_shift[shifttry * 9];
      Real worst_mxlead = BIG_REAL;
      for (int whichld = 0; whichld < 3; ++whichld) {

        Real w[3];
        {
          for (int i = 0; i < 3; ++i) 
            w[i] = Qb[i + 3 * whichld];
     
        }
        if (w[2] == 0.0) continue;
        Real mxlead = 0.0;
        for (int eqnum = 0; eqnum < 2; ++eqnum) {
          Real val1 = fabs(pl_evalr(coef, R2Point(w[0]/w[2], w[1]/w[2]), eqnum,
            degree) * pow(w[2],degree));
          if (val1 > mxlead)
            mxlead = val1;
        }
        if (mxlead < worst_mxlead)
          worst_mxlead = mxlead;
      }
      if (worst_mxlead > bestlead) {
        bestlead = worst_mxlead;
        bestshift = shifttry;
      }
    }

    // shiftQ is the shift we will use -- store it in the return variable.
    
    {
      const Real* shiftQ1 = &try_shift[bestshift * 9];
      for (int i = 0; i < 9; ++i)
        shiftQ[i] = shiftQ1[i];
    }
  
    tri_poly_frac_lin_trans(degree, shiftQ, coef, shiftpoly, workspace);

    
#if 0
    // Check that the shift was done correctly.
    Real minlead = BIG_REAL;
    for (int wl2 = 0; wl2 < 3; ++wl2) {
      int col;
      if (wl2 == 0)
        col = deg_pair_to_colnum(0,0,degree);
      else if (wl2 == 1)
        col = deg_pair_to_colnum(degree,0,degree);
      else
        col = deg_pair_to_colnum(0,degree,degree);
      Real mxlead = 0.0;
      for (int eqnum = 0; eqnum < 2; ++eqnum) {
        Real f = fabs(shiftpoly(eqnum,col));
        if (f > mxlead)
          mxlead = f;
      }
      if (mxlead < minlead) {
        minlead = mxlead;
      }
    }
    if (fabs(bestlead - minlead) > 1e-10)
      throw_error("Problem computing shift");
#endif
  }

  // Part 1 of solving a system of two triangular bezier
  // polynomials.  This routine takes the coefficients
  // and computes a QR factorization of the top part of
  // of the Macaulay matrix.  These are returned in hhmat_tmp
  // and betavec_tmp.  Also compute use_monom, indicating
  // which monomials to use for the bottom part of the matrix.
  // This routine may alter coef by a perturbation if
  // necessary.  Return variable is the minimum singular value
  // estimate.

  Real solve_bez_tri_main1(int degree,
    PatchMath::Workspace::RealMatrix& coef,
    PatchMath::Workspace::RealMatrix& hhmat_tmp,
    PatchMath::Workspace::RealVector& betavec_tmp,
    PatchMath::Workspace::RealVector& use_monom,
    PatchMath::Workspace& workspace) {
    
#ifdef DEBUGGING
    using QMG::Local_to_PatchMath_CPP::dump_everything;
    using QMG::Local_to_PatchMath_CPP::debug_str;
#endif

    int tmrow = degree * (degree + 1);
    int tmcol = degree * (2 * degree + 1);

    int monomcount = (2 * degree - 1) * degree;
    int num_row_monom = degree * (degree + 1) / 2;
    int num_control_pt = (degree + 1) * (degree + 2) / 2;
   
    
    // Now we must factor the top part of the matrix.  If
    // factorization detects a singularity, try again in
    // double-double.  If there is still a singularity,
    // perturb by a small amount and try again.
    
    // Compute the top part of the matrix,
    // which contains shifted copies of the bezier coefficients.
    
    // Note, in this code, the top part is transposed, so it's
    // actually the 'left part'.

    PatchMath::Workspace::StartPosMarker unnamed_(workspace);
    PatchMath::Workspace::RealMatrix Tm(workspace, tmcol, tmrow);
    PatchMath::Workspace::RealMatrix Tmp(workspace, tmcol, tmrow);
    PatchMath::Workspace::RealVector select_col_tmp(workspace, tmrow);
    PatchMath::Workspace::DoubleDoubleVector ddworkspace(workspace, tmcol * (tmrow + 1));
    
    Real ptb_size = MACHINE_EPS * tmcol;
    int try_shiftp1 = 0;
    Real minsv_est;
    Real rminsv_est;
    const Real MIN_SV_ALLOWABLE = 1.0e-3;
    int which_factorization;
    
    for (int top_part_trycount = 0; true; ++top_part_trycount) {
      if (top_part_trycount == 30)
        throw_error("Unable to find good factorization");
      
      for (int row_monom_idx = 0; row_monom_idx < num_row_monom; ++row_monom_idx) {
        pair<int,int> row_idx = colnum_to_deg_pair(row_monom_idx, degree - 1);
        int u_row_deg = row_idx.first;
        int v_row_deg = row_idx.second;
        for (int eqnum = 0; eqnum < 2; ++eqnum) {
          int mtxrow = row_monom_idx * 2 + eqnum;
          for (int j = 0; j < tmcol; ++j)
            Tm(j, mtxrow) = 0.0;
          for (int poly_monom_idx = 0; poly_monom_idx < num_control_pt; ++poly_monom_idx) {
            pair<int,int> poly_idx = colnum_to_deg_pair(poly_monom_idx, degree);
            int u_poly_deg = poly_idx.first;
            int v_poly_deg = poly_idx.second;
            int u_col_deg = u_poly_deg + u_row_deg;
            int v_col_deg = v_poly_deg + v_row_deg;
            int mtxcol = deg_pair_to_colnum(u_col_deg,v_col_deg, 2*degree - 1);
            Tm(mtxcol, mtxrow) = coef(eqnum, poly_monom_idx);
          }
        }
      }
      
      {
        for (int i = 0; i < tmcol; ++i) {
          for (int j = 0; j < tmrow; ++j) {
            Tmp(i,j) = Tm(i,j);
          }
        }
      }
      
      
      if (top_part_trycount == 0) {
        which_factorization = 1;
        QR_Factor_with_col_pivoting(Tmp, select_col_tmp, hhmat_tmp, betavec_tmp);
      }
      else {
        which_factorization = 2;
        QR_Factor_with_col_pivoting_quadprec(&Tmp(0,0), tmcol, tmrow, &select_col_tmp[0], 
          &hhmat_tmp(0,0), &betavec_tmp[0], &ddworkspace[0]);
      }
      
      // Estimate the minimum singular value.
      
      minsv_est = fabs(Tmp(tmrow - 1, tmrow - 1));
      
#ifdef DEBUGGING
      if (dump_everything) {
        *debug_str << "% factorize top part trycount = " << top_part_trycount
          << " minsv = " << minsv_est << " ptb_size = " << ptb_size << '\n';
      }
#endif
      
      if (minsv_est >= MIN_SV_ALLOWABLE)
        break;
      if (minsv_est > MACHINE_EPS && which_factorization == 2)
        break;
      if (top_part_trycount < 2)
        rminsv_est = minsv_est;
      
      // Perturb the matrix.
      
      if (top_part_trycount >= 1) {
        ptb_size *= 10.0;
        for (int i = 0; i < 2; ++i) {
          for (int j = 0; j < num_control_pt; ++j) {
            coef(i,j) += cos(static_cast<Real>(try_shiftp1++)) * ptb_size;
          }
        }
        ++try_shiftp1;
      }
    }

    
    // Figure out which rows from the bottom part should be dropped.
    // Make Tmsmall, which is like Tm except one lower degree monomials multiply the rows.
    
    int tmsmall_row = degree * (degree - 1);
    int tmsmall_col = monomcount;
    
    PatchMath::Workspace::RealMatrix Tmsmall(workspace, tmsmall_col, tmsmall_row);
    
    {
      int num_row_monom = degree * (degree - 1) / 2;
      for (int row_monom_idx = 0; row_monom_idx < num_row_monom; ++row_monom_idx) {
        pair<int,int> row_idx = colnum_to_deg_pair(row_monom_idx, degree - 2);
        int u_row_deg = row_idx.first;
        int v_row_deg = row_idx.second;
        for (int eqnum = 0; eqnum < 2; ++eqnum) {
          int mtxrow = row_monom_idx * 2 + eqnum;
          for (int j = 0; j < tmsmall_col; ++j)
            Tmsmall(j,mtxrow) = 0.0;
          for (int poly_monom_idx = 0; poly_monom_idx < num_control_pt; ++poly_monom_idx) {
            pair<int,int> poly_idx = colnum_to_deg_pair(poly_monom_idx, degree);
            int u_poly_deg = poly_idx.first;
            int v_poly_deg = poly_idx.second;
            int u_col_deg = u_poly_deg + u_row_deg;
            int v_col_deg = v_poly_deg + v_row_deg;
            int mtxcol = deg_pair_to_colnum(u_col_deg,v_col_deg, 2*degree - 2);
            Tmsmall(mtxcol,mtxrow) = coef(eqnum, poly_monom_idx);
          }
        }
      }
    }
    
    // Replace Tmsmall with its QR factorization.
    
    PatchMath::Workspace::RealMatrix hhmat_small1(workspace, tmsmall_col, tmsmall_row);
    PatchMath::Workspace::RealVector betavec_small1(workspace, tmsmall_row);
    PatchMath::Workspace::RealVector select_col_small1(workspace, tmsmall_row);
    
    if (which_factorization == 1) {
      QR_Factor_with_col_pivoting(Tmsmall, select_col_small1, hhmat_small1, betavec_small1);
    }
    else {
      QR_Factor_with_col_pivoting_quadprec(&Tmsmall(0,0), tmsmall_col, tmsmall_row,
        &select_col_small1[0], &hhmat_small1(0,0), &betavec_small1[0],
        &ddworkspace[0]);
    }
    
    // Form Q^T, the transpose of the result of QR factorization
    
    PatchMath::Workspace::RealMatrix Tmsmall2(workspace, tmsmall_row, tmsmall_col);
    
    formQ(hhmat_small1, betavec_small1, Tmsmall2);
    
    // Now perform QR factorization with column pivoting on Tmsmall2 to select
    // monomials.
    //  save select_col.
    
    PatchMath::Workspace::RealMatrix hhmat_small(workspace, tmsmall_row, tmsmall_row);
    PatchMath::Workspace::RealVector betavec_small(workspace, tmsmall_row);
    PatchMath::Workspace::RealVector select_col_small(workspace, tmsmall_col);
    
    QR_Factor_with_col_pivoting(Tmsmall2, select_col_small, hhmat_small, betavec_small);
    
#ifdef DEBUGGING
    if (dump_everything) {
      *debug_str << "Tm = [\n";
      {
        for (int i = 0; i < tmcol; ++i) {
          for (int j = 0; j < tmrow; ++j) {
            *debug_str << Tm(i,j);
            if (j < tmrow - 1)
              *debug_str << ',';
          }
          *debug_str << '\n';
        }
      }
      *debug_str << "];\n% after QR factorization with col pivoting, factors are\nTmp = [\n";
      {
        for (int i = 0; i < tmcol; ++i) {
          for (int j = 0; j < tmrow; ++j) {
            *debug_str << Tmp(i,j) << " ";
          }
          *debug_str << '\n';
        }
      }
      *debug_str << "];\nselect_col_tmp = [";
      {
        for (int i = 0; i < tmrow; ++i)
          *debug_str << static_cast<int>(select_col_tmp[i]) << " ";
      }
      *debug_str << "];\nhhmat_tmp = [";
      {
        for (int i = 0; i < tmcol; ++i) {
          if (i < tmrow) {
            for (int j = 0; j < i + 1; ++j) {
              *debug_str << hhmat_tmp(i,j) << " ";
            }
            for (int j2 = i + 1; j2 < tmrow; ++j2) {
              *debug_str << "0 ";
            }
          }
          else {
            for (int j = 0; j < tmrow; ++j) {
              *debug_str << hhmat_tmp(i,j) << " ";
            }
          }
          *debug_str << '\n';
        }
      }
      *debug_str << "];\nbetavec_tmp = [";
      {
        for (int i = 0; i < tmrow; ++i)
          *debug_str << betavec_tmp[i] << " ";
      }
      *debug_str << "];\n";
      *debug_str << "Tmsmall = [\n";
      {
        for (int i = 0; i < tmsmall_col; ++i) {
          for (int j = 0; j < tmsmall_row; ++j)
            *debug_str << " " << std::setprecision(17) << Tmsmall(i,j);
          *debug_str << "\n";
        }
      }
      *debug_str << "]\n";
      *debug_str << "Tmsmall2 = [\n";
      {
        for (int i = 0; i < tmsmall_row; ++i) {
          for (int j = 0; j < tmsmall_col; ++j)
            *debug_str << " " << std::setprecision(17) << Tmsmall2(i,j);
          *debug_str << "\n";
        }
      }
      *debug_str << "]\n";
      *debug_str << "selectcol = [";
      {
        for (int i = 0; i < tmsmall_row; ++i) 
          *debug_str << static_cast<int>(select_col_small[i]) << " ";
      }
      *debug_str << "]\n";
    }
#endif  
    {
      for (int j = 0; j < monomcount; ++j)
        use_monom[j] = 1.0;
      for (int j2 = 0; j2 < tmsmall_row; ++j2)
        use_monom[static_cast<int>(select_col_small[j2])] = 0.0;
    }
    return rminsv_est;
  }


  // Part 3 of the main algorithm to solve a system of
  // triangular bezier equations.  In this part we
  // take a QR factorization of the top part of the
  // matrix and the list of rows to use, and a line-search
  // direction, and we set up a generalized eigenvalue
  // problem to find all roots.  The roots are returned
  // in act_solution_data.
  // Return -1 if lapack failure, else 0.


  int solve_bez_tri_main2(int degree,
    const PatchMath::Workspace::RealMatrix& coef,
    const PatchMath::Workspace::RealMatrix& hhmat_tmp,
    const PatchMath::Workspace::RealVector& betavec_tmp,
    const PatchMath::Workspace::RealVector& use_monom,
    const Real direc[2],
    PatchMath::Workspace::ComplexMatrix& act_solution_data,
    PatchMath::Workspace& workspace) {
    
    
#ifdef DEBUGGING
    using QMG::Local_to_PatchMath_CPP::debug_str;
    using QMG::Local_to_PatchMath_CPP::dump_everything;
    using std::setprecision;
#endif
    
    using std::abs;

    int tmrow = degree * (degree + 1);
    int tmcol = degree * (2 * degree + 1);
    int evsize = degree * degree; // = tmcol-tmrow
    int monomcount = (2 * degree - 1) * degree;

    
    PatchMath::Workspace::StartPosMarker unnamed_(workspace);
    
    PatchMath::Workspace::RealVector expaArow(workspace, tmcol);
    PatchMath::Workspace::RealVector expaBrow(workspace, tmcol);
    PatchMath::Workspace::ComplexMatrix A(workspace, evsize, evsize);
    PatchMath::Workspace::ComplexMatrix Ao(workspace, evsize, evsize);
    PatchMath::Workspace::ComplexMatrix B(workspace, evsize, evsize);
    PatchMath::Workspace::ComplexMatrix Bo(workspace, evsize, evsize);
    PatchMath::Workspace::ComplexMatrix vr(workspace, evsize, evsize);
    PatchMath::Workspace::ComplexVector transfevec(workspace, tmcol);
    PatchMath::Workspace::ComplexVector lambda(workspace, evsize);
    PatchMath::Workspace::RealVector evec_cond(workspace, evsize);
    int rowcount = 0;
    for (int evrownum = 0; evrownum < monomcount; ++evrownum) {
      if (use_monom[evrownum] == 0.0) continue;
      
      for (int jj = 0; jj < tmcol; ++jj) {
        expaArow[jj] = 0.0;
        expaBrow[jj] = 0.0;
      }
      pair<int,int> rval = colnum_to_deg_pair(evrownum, 2 * degree - 2);
      int ud = rval.first;
      int vd = rval.second;

      int colnum1 = deg_pair_to_colnum(ud + 1, vd, 2*degree - 1);
      expaArow[colnum1] = direc[0];
      colnum1 = deg_pair_to_colnum(ud, vd + 1, 2*degree - 1);
      expaArow[colnum1] = direc[1];
      colnum1 = deg_pair_to_colnum(ud,vd, 2*degree - 1);
      expaBrow[colnum1] = 1.0;
#ifdef DEBUGGING
      if (dump_everything) {
        *debug_str << "expaArow(" << evrownum + 1 << ",:) = [";
        {
          for (int i = 0; i < tmcol; ++i) {
            *debug_str << setprecision(17) << expaArow[i] << " ";
          }
        }
        *debug_str << "];\n expaBrow(" << evrownum + 1 << ",:) = [";
        {
          for (int i = 0; i < tmcol; ++i) {
            *debug_str << setprecision(17) << expaBrow[i] << " ";
          }
        }
        *debug_str << "];\n";
      }
#endif          
      
      // Apply HH transforms 
      
      for (int j = 0; j < tmrow; ++j) {
        Real ip1 = 0.0;
        Real ip2 = 0.0;
        {
          for (int i = j; i < tmcol; ++i) {
            ip1 += hhmat_tmp(i,j) * expaArow[i];
            ip2 += hhmat_tmp(i,j) * expaBrow[i];
          }
        }
        ip1 *= betavec_tmp[j];
        ip2 *= betavec_tmp[j];
        {
          for (int i = j; i < tmcol; ++i) {
            expaArow[i] -= ip1 * hhmat_tmp(i,j);
            expaBrow[i] -= ip2 * hhmat_tmp(i,j);
          }
        }
      }
      
      for (int i = tmrow; i < tmcol; ++i) {
        A(rowcount, i - tmrow) = Complex(expaArow[i], 0.0);
        B(rowcount, i - tmrow) = Complex(expaBrow[i], 0.0);
      }
      ++rowcount;
    }
  
#ifdef DEBUGGING
    if (rowcount != evsize)
      throw_error("Counting problem");
#endif
    
    
    {
      for (int i = 0; i < evsize; ++i)
        for (int j = 0; j < evsize; ++j) {
          Ao(i,j) = A(i,j);
          Bo(i,j) = B(i,j);
        }
    }
    
#ifdef DEBUGGING
    if (dump_everything) {
      *debug_str << "];\nAi = [\n";
      {
        for (int i = 0; i < evsize; ++i) {
          for (int j = 0; j < evsize; ++j) {
            *debug_str << Format_Complex(A(i,j)) << " ";
          }
          *debug_str << '\n';
        }
      }
      *debug_str << "];\nBi = [\n";
      {
        for (int i = 0; i < evsize; ++i) {
          for (int j = 0; j < evsize; ++j) {
            *debug_str << Format_Complex(B(i,j)) << " ";
          }
          *debug_str << '\n';
        }
      }
      *debug_str << "];\n";
    }
#endif
    
    int info = generalized_eigensolver(A, B, vr, evec_cond, lambda, workspace);

    if (info != 0) {
#ifdef DEBUGGING
      *debug_str << "Lapack failure in Solve_BezTri::solve\n";
      *debug_str << "A = ";
      {
        for (int i = 0; i < evsize; ++i) {
          for (int j = 0; j < evsize; ++j) {
            *debug_str << Format_Complex(Ao(i,j)) << " ";
          }
          *debug_str << '\n';
        } 
      }
      *debug_str << "B = ";
      {
        for (int i = 0; i < evsize; ++i) {
          for (int j = 0; j < evsize; ++j) {
            *debug_str << Format_Complex(Bo(i,j)) << " ";
          }
          *debug_str << '\n';
        }
      }  
      *debug_str << "Afinal = ";
      {
        for (int i = 0; i < evsize; ++i) {
          for (int j = 0; j < evsize; ++j) {
            *debug_str << Format_Complex(A(i,j)) << " ";
          }
          *debug_str << '\n';
        } 
      }
      *debug_str << "Bfinal = ";
      {
        for (int i = 0; i < evsize; ++i) {
          for (int j = 0; j < evsize; ++j) {
            *debug_str << Format_Complex(B(i,j)) << " ";
          }
          *debug_str << '\n';
        }
      }
#endif
      return -1;
    }
    
    // Find the polynomial solutions from the eigenvectors.
    for (int vlrow = 0; vlrow < evsize; ++vlrow) {
      
      // Apply the HH transforms to get the transformed eigenvector.          
       {
        for (int j = 0; j < tmrow; ++j)
          transfevec[j] = Complex(0.0,0.0);
      }
      
      {
        for (int j = tmrow; j < tmcol; ++j) 
          transfevec[j] = vr(j - tmrow, vlrow);
      }
      {
        for (int j = tmrow - 1; j >= 0; --j) {
          Complex ip(0.0,0.0);
          {
            for (int i = j; i < tmcol; ++i) {
              ip += Complex(hhmat_tmp(i,j),0.0) * transfevec[i];
            }
          }
          ip *= betavec_tmp[j];
          {
            for (int i = j; i < tmcol; ++i) {
              transfevec[i] -= ip * hhmat_tmp(i,j);
            }
          }
        }
      }
      
      // Find whether u^d, v^d, or w^d is the largest coeff.
      // Then determine u,v,w from powers close to the largest power.
      
      Complex ud = transfevec[deg_pair_to_colnum(2 * degree - 1, 0, 2 * degree - 1)];
      Complex vd = transfevec[deg_pair_to_colnum(0, 2 * degree - 1, 2 * degree - 1)];
      Complex wd = transfevec[deg_pair_to_colnum(0, 0, 2 * degree - 1)];
      Complex u1, v1, w1;
      
      if (abs(ud) >= abs(vd) && abs(ud) >= abs(wd)) {
        u1 = ud;
        v1 = transfevec[deg_pair_to_colnum(2 * degree - 2, 1, 2 * degree - 1)];
        w1 = transfevec[deg_pair_to_colnum(2 * degree - 2, 0, 2 * degree - 1)];
      }
      else if (abs(vd) >= abs(wd)) {
        u1 = transfevec[deg_pair_to_colnum(1, 2 * degree - 2, 2 * degree - 1)];
        v1 = vd;
        w1 = transfevec[deg_pair_to_colnum(0, 2 * degree - 2, 2 * degree - 1)];
      }
      else {
        u1 = transfevec[deg_pair_to_colnum(1, 0, 2 * degree - 1)];
        v1 = transfevec[deg_pair_to_colnum(0, 1, 2 * degree - 1)];
        w1 = wd;
      }
      
      Complex denom0a = (w1 == 0.0)? MACHINE_EPS : w1;
      C2Point sol0(u1/denom0a, v1/denom0a);

      // Now compute the residual.
      
      
      C2Point spfval;      
      spfval.coord[0] = pl_eval(coef, sol0, 0, degree);
      spfval.coord[1] = pl_eval(coef, sol0, 1, degree);
      Real aspfval1 = pow(abs(sol0.coord[0]), degree);
      Real aspfval2 = pow(abs(sol0.coord[1]), degree);
      Real aspfval = (aspfval1 > aspfval2)? aspfval1 : aspfval2;
      aspfval = (aspfval > 1.0)? aspfval : 1.0;
      Real this_resd = spfval.l2norm() / aspfval;
      act_solution_data(0,vlrow) = sol0.coord[0];
      act_solution_data(1,vlrow) = sol0.coord[1];
      act_solution_data(2,vlrow) = Complex(evec_cond[vlrow], 0.0);
      act_solution_data(3,vlrow) = Complex(this_resd, 0.0);
      act_solution_data(4,vlrow) = lambda[vlrow];
#ifdef DEBUGGING
      if (dump_everything) {
        *debug_str << " % vlrow = " << vlrow << "\n";
        *debug_str << " vri = [\n";
        {
          for (int i = 0; i < evsize; ++i) {
            *debug_str << Format_Complex(vr(i,vlrow)) << " ";
          }
        }
        *debug_str << "\n];\n";
        *debug_str << "transfevec = [\n";
        {
          for (int i = 0; i < tmcol; ++i)
            *debug_str << Format_Complex(transfevec[i]) << " ";
        }
        *debug_str << "\n u1i = " 
          << Format_Complex(u1) 
          << ";\n v1i = " << Format_Complex(v1) 
          << ";\n w1i = " << Format_Complex(w1)
          << ";\n actsol = [" << Format_Complex(sol0.coord[0]) << "," << Format_Complex(sol0.coord[1])
          << ";\n eigval = " << Format_Complex(act_solution_data(4,vlrow)) 
          << ";\n resd = " << act_solution_data(3, vlrow).real()
          << ";\n eveccond = " << act_solution_data(2,vlrow).real()
          << ";\n";
      }
#endif
    }
    return 0;
   }


  // Solve a system of two polynomial equations
  // in Bezier triangle form.  Use Macaulay-type
  // result and random u-direction.  Reduce
  // to small system and invoke LAPACK for
  // generalized eigenvalues.
  // Returns true if there was a global degeneracy
  // during the computation.

  
  bool Solve_BezTri::solve(int degree,
    vector<PatchMath::IntersectRecord>& output_stack,
    PatchMath::Workspace& workspace, 
    double scfac, 
    double tol,
    bool newton_improve) const {


#ifdef DEBUGGING

    using QMG::Local_to_PatchMath_CPP::debug_str;
    using QMG::Local_to_PatchMath_CPP::dump_everything;
    using QMG::Local_to_PatchMath_CPP::invocation_count;
    using std::setprecision;
    using std::setiosflags;
    using std::ios;
    using std::cos;

    ++invocation_count;
    bool inner_dump_everything = false;


    // if (invocation_count > 962)
    //  dump_everything = true;
    // dump_everything = (invocation_count == 148);


    if (dump_everything) {
      *debug_str << " In solve beztri degree = " << degree 
        << " scfac = " << scfac << " tol = " << tol << "\n bez  = [\n";
      int ncp = (degree + 1) * (degree + 2) / 2; 
      for (int i = 0; i < 2; ++i) {
        for (int j = 0; j < ncp; ++j) {
          *debug_str << std::setprecision(17) << get_bez_coeff(i,j) << " ";
        }
        *debug_str << '\n';
      }
      *debug_str << "];\n";
      *debug_str << "invocation_count = " << invocation_count << '\n';
    }
#endif

    bool is_degen = false;

    if (degree == 1) {
      Real x[3];
      Real rhs[2];
      Real matrixspace[4];
      Matrix m(2,2,matrixspace);
      for (int eqnum = 0; eqnum < 2; ++eqnum) {
        m(eqnum,0) = get_bez_coeff(eqnum,2) - get_bez_coeff(eqnum,1);
        m(eqnum,1) = get_bez_coeff(eqnum,0) - get_bez_coeff(eqnum,1);
        rhs[eqnum] = -get_bez_coeff(eqnum,1);
      }
      Real invnormest = m.linear_solve(rhs,x);
      if (invnormest >= BIG_REAL) {
        is_degen = true;
        return is_degen;
      }

      x[2] = 1.0 - x[0] - x[1];
      Real condest = scfac * invnormest;
      if (condest > sqrt(tol))
        is_degen = true;
      Real unc = condest * tol;
      if (x[0] < -unc || x[1] < -unc || x[2] < -unc)
        return is_degen;
      PatchMath::IntersectRecord result;
      result.paramcoord[0] = x[0];
      result.paramcoord[1] = x[1];
      result.param_uncertainty = unc;
      result.degen = fabs(x[0]) <= unc || fabs(x[1]) <= unc || fabs(x[2]) <= unc;
      output_stack.push_back(result);
      return is_degen;
    }
    else {    
      using std::abs;
      using std::norm;

      PatchMath::Workspace::StartPosMarker start_pos_marker_(workspace);
      int num_control_pt = (degree + 1) * (degree + 2) / 2;

      
      // Check if (0,0) is in the bounding box; if not then early return.
      {
        for (int eqnum = 0; eqnum < 2; ++eqnum) {
          bool pos_seen = false;
          bool neg_seen = false;
          for (int i = 0; i < num_control_pt; ++i) {
            Real bz = get_bez_coeff(eqnum, i);
            if (bz <= scfac * tol)
              neg_seen = true;
            if (bz >= -scfac * tol)
              pos_seen = true;
          }
          if (!neg_seen || !pos_seen) {
#ifdef DEBUGGING
            if (dump_everything) {
              *debug_str << " early return; equation " << eqnum << " misses 0.\n";
            }
#endif
            return false;
          }
        }
      }
      
      int tmrow = degree * (degree + 1);
      int tmcol = degree * (2 * degree + 1);
      int evsize = degree * degree; // = tmcol-tmrow
      int monomcount = (2 * degree - 1) * degree;


      Real bezmax_poly[2];
      bezmax_poly[0] = 0.0;
      bezmax_poly[1] = 0.0;

      {
        for (int i = 0; i < 2; ++i) {
          for (int j = 0; j < num_control_pt; ++j) {
            Real c1 = fabs(get_bez_coeff(i,j));
            if (c1 > bezmax_poly[i])
              bezmax_poly[i] = c1;
          }
        }
      }
      if (bezmax_poly[0] == 0 || bezmax_poly[1] == 0) {
#ifdef DEBUGGING
        if (dump_everything) {
          *debug_str << " early return; one equation is all zeros\n";
        }
#endif
        is_degen = true;
        return is_degen;
      }

      Real sctol = scfac / bezmax_poly[0];
      if (bezmax_poly[0] > bezmax_poly[1])
        sctol = scfac / bezmax_poly[1];
      if (sctol < 1.0)
        sctol = 1.0;

      
      PatchMath::Workspace::RealMatrix shiftpoly(workspace, 2, num_control_pt); 
      PatchMath::Workspace::RealMatrix acoef(workspace, 2, num_control_pt); 
      PatchMath::Workspace::RealMatrix bcoef(workspace, 2, num_control_pt); 

  
      // Rescale coefficients by max value, and also multiply by binomial
      // coefficients.
      {
        int outer_binom = 1;
        for (int vdeg = 0; vdeg <= degree; ++vdeg) {
          int inner_binom = 1;
          for (int udeg = 0; udeg <= degree - vdeg; ++udeg) {
            int col = deg_pair_to_colnum(udeg, vdeg, degree);
            for (int eqnum = 0; eqnum < 2; ++eqnum) {
              acoef(eqnum,col) = get_bez_coeff(eqnum, col) * 
                inner_binom * outer_binom / bezmax_poly[eqnum];
            }
            inner_binom *= (degree - udeg - vdeg);
            inner_binom /= (udeg + 1);
          }
          outer_binom *= (degree - vdeg);
          outer_binom /= (vdeg + 1);
        }
      }

      Real shiftQ0[9];

      // Compute a shift that changes bezier form to standard form.
      shiftQ0[0] = 1.0;
      shiftQ0[1] = 0.0;
      shiftQ0[2] = -1.0;
      shiftQ0[3] = 0.0;
      shiftQ0[4] = 1.0;
      shiftQ0[5] = -1.0;
      shiftQ0[6] = 0.0;
      shiftQ0[7] = 0.0;
      shiftQ0[8] = 1.0;
      
      tri_poly_frac_lin_trans(degree, shiftQ0, acoef, bcoef, workspace);

      
      Real shiftQ[9];
      
      solve_bez_tri_main0(degree, bcoef,shiftpoly, shiftQ, workspace);
     
#ifdef DEBUGGING
      if (dump_everything) {
        *debug_str << "acoef = [\n";
        {
          for (int i = 0; i < 2; ++i) {
            for (int j = 0; j < num_control_pt; ++j)
              *debug_str << " " << std::setprecision(17)
              << acoef(i,j);
            *debug_str << '\n';
          }
        }
        *debug_str << "];\nbcoef = [\n";
        {
          for (int i = 0; i < 2; ++i) {
            for (int j = 0; j < num_control_pt; ++j)
              *debug_str << " " << std::setprecision(17)
              << bcoef(i,j);
            *debug_str << '\n';
          }
        }
        *debug_str << "];\nshiftQ = [\n";
        {
          for (int i = 0; i < 3; ++i) {
            for (int j = 0; j < 3; ++j) {
              *debug_str << std::setprecision(17) << shiftQ[i + j * 3] << " ";
            }
            *debug_str << "\n";
          }
        }
        *debug_str << "];\n shiftpoly = [\n";
        {
          for (int i = 0; i < 2; ++i) {
            for (int j = 0; j < num_control_pt; ++j)
              *debug_str << " " << std::setprecision(17)
              << shiftpoly(i,j);
            *debug_str << '\n';
          }
        }
        *debug_str << "];\n";
      }
#endif
      
      PatchMath::Workspace::RealMatrix hhmat_tmp(workspace, tmcol, tmrow);
      PatchMath::Workspace::RealVector betavec_tmp(workspace, tmrow);
      PatchMath::Workspace::RealVector use_monom(workspace, monomcount);    
  
      solve_bez_tri_main1(degree, shiftpoly, hhmat_tmp, betavec_tmp, use_monom, workspace);

      PatchMath::Workspace::RealMatrix solution_data(workspace, 4, evsize);
      PatchMath::Workspace::RealMatrix test_solution_data(workspace, 4, evsize);
      PatchMath::Workspace::ComplexMatrix act_solution_data(workspace, 5, evsize);
      PatchMath::Workspace::RealVector tmp_poly(workspace, degree + 1);
      PatchMath::Workspace::RealMatrix zoom_coef(workspace, 2, num_control_pt);
      PatchMath::Workspace::RealMatrix zoom_shift(workspace, 2, num_control_pt);
      PatchMath::Workspace::RealMatrix zoom_hhmat(workspace, tmcol, tmrow);
      PatchMath::Workspace::RealVector zoom_betavec(workspace, tmrow);
      PatchMath::Workspace::RealVector zoom_use_monom(workspace, monomcount);
      PatchMath::Workspace::ComplexMatrix zoom_act_solution_data(workspace, 5, evsize);    
      PatchMath::Workspace::ComplexMatrix jac(workspace, 2, 2);
      

      const int NUM_DIREC_TRY = 3;
      const int LENGTH_DIREC_TRY = NUM_DIREC_TRY * 2;
      double direc_try[LENGTH_DIREC_TRY] = {
      -0.43256481152822,   0.28767642035855,
        -1.66558437823810,  -1.14647135068146,
        0.12533230647483,   1.19091546564300};

      int lapack_failcount = 0;

      int numsol = 0;

      Real best_max_resd2 = BIG_REAL;
      Real resd_cutoff = tol;
      
      for (int direc_try_count = 0; direc_try_count < NUM_DIREC_TRY; 
      ++direc_try_count) {
        
        
        // Now compute the eigenval matrix.
        // Choose some random numbers.

        Real direc[2];
        direc[0] = direc_try[2 * direc_try_count];
        direc[1] = direc_try[2 * direc_try_count + 1];

        // Work out eigenvalue upper and lower bounds.

        // First compute upper and lower bounds on the numerator
        // and denominator.

        Real denom_min = BIG_REAL;
        Real denom_max = -BIG_REAL;
        Real numer_min = BIG_REAL;
        Real numer_max = -BIG_REAL;
        {
          for (int i = 0; i < 3; ++i) {
            Real numeval = direc[0] * shiftQ[i] + direc[1] * shiftQ[i + 3];
            Real denomeval = shiftQ[i + 6];
            if (denomeval < denom_min)
              denom_min = denomeval;
            if (denomeval > denom_max)
              denom_max = denomeval;
            if (numeval < numer_min)
              numer_min = numeval;
            if (numeval > numer_max)
              numer_max = numeval;
          }
        }
        if (denom_min <= 0)
          throw_error("Illegal shift selection");
        Real eigval_ub = numer_max / denom_min;
        Real eigval_lb = (numer_min < 0)? numer_min / denom_min : numer_min / denom_max;

#ifdef DEBUGGING
        if (dump_everything) {
          *debug_str << "direc_try_count = " << direc_try_count << " direc = [";
          for (int i = 0; i < 2; ++i) {
            *debug_str << setprecision(17) << direc[i] << " ";
          }
          *debug_str << "];\neigbounds = [" << eigval_lb << ","
            << eigval_ub << "];\n";
        }
#endif
       
        int returncode = solve_bez_tri_main2(degree, shiftpoly,hhmat_tmp, betavec_tmp,
          use_monom, direc, act_solution_data, workspace);

        if (returncode != 0) {
          ++lapack_failcount;
          if (++lapack_failcount == NUM_DIREC_TRY) {
            throw_error("Repeated lapack failure");
          }
          continue;
        }       
        
        Real max_resd = 0.0;
        Real max_resd2 = 0.0;
        int test_numsol = 0;

        
        for (int vlrow = 0; vlrow < evsize; ++vlrow) {
          Real this_resd = act_solution_data(3,vlrow).real();
          Complex eigval = act_solution_data(4,vlrow);
          bool in_range = eigval.real() > eigval_lb - 0.1 && 
            eigval.real() < eigval_ub + 0.1 &&
            fabs(eigval.imag()) < 0.1;
          Real evec_cond = act_solution_data(2,vlrow).real();
          
          // If the residual is big, try to drive it down with corrective
          // solver iterations.  This loop should be executed infrequently I hope.
          
          int improve_count = 0;
          while (in_range && evec_cond > 1e5 && this_resd > resd_cutoff && ++improve_count < 3) {

            Real dila = sqrt(this_resd);
            Complex transxc = act_solution_data(0,vlrow);
            Complex transyc = act_solution_data(1,vlrow);
            Real transx = transxc.real();
            Real transy = transyc.real();

            // Compute the translated dilated polynomial in zoom_coef, defined
            // to be sum_{j=0}^d sum_{i=0}^{d-j} a_{i,j} (transx+dila*x)^i*(transy+dila*y)^j

            Real shiftQi[9];
            shiftQi[0] = dila;
            shiftQi[1] = 0.0;
            shiftQi[2] = 0.0;
            shiftQi[3] = 0.0;
            shiftQi[4] = dila;
            shiftQi[5] = 0.0;
            shiftQi[6] = transx;
            shiftQi[7] = transy;
            shiftQi[8] = 1.0;            

            tri_poly_frac_lin_trans(degree, shiftQi, shiftpoly, zoom_coef, workspace);

            for (int eqnum = 0; eqnum < 2; ++eqnum) {
              // Rescale so the max is 1
              Real mx = 0.0;
              {
                for (int j = 0; j < num_control_pt; ++j) {
                  Real f = fabs(zoom_coef(eqnum, j));
                  if (f > mx) {
                    mx = f;
                  }
                }
              }
              {
                for (int j = 0; j < num_control_pt; ++j) {
                  zoom_coef(eqnum, j) /= mx;
                }
              }
            }
            
            // Now shift it to make leading coefficients large.
            Real zshiftQ[9];
            
#ifdef DEBUGGING
            bool save_dump_everything = dump_everything;
            if (dump_everything) {
              *debug_str <<" % For vlrow " << vlrow << " improve_count " << improve_count
                << " calling solve_bez_tri_main1\n";
            }
            
            if (!inner_dump_everything)
              dump_everything = false;
            
#endif            
            solve_bez_tri_main0(degree, zoom_coef, zoom_shift, zshiftQ, workspace);
            Real zoom_minsv = solve_bez_tri_main1(degree, zoom_shift, zoom_hhmat,
              zoom_betavec, zoom_use_monom, workspace);
            int returncode = solve_bez_tri_main2(degree, zoom_shift, zoom_hhmat,
              zoom_betavec, zoom_use_monom, direc, zoom_act_solution_data, workspace);
            if (returncode != 0)
              break;
            // Find the correction closest to 0.
            
            Real minac = BIG_REAL;
            C2Point best_actsol;
            int minac_row;
            {
              for (int j = 0; j < evsize; ++j) {
                // Undo the shift               
                
                // Apply the inverse of zshiftQ to (u1,v1,w1) to undo
                // undo effect of shift.
                
                Complex u1 = zoom_act_solution_data(0,j);
                Complex v1 = zoom_act_solution_data(1,j);
                
                Complex u2 = u1 * zshiftQ[0] + v1 * zshiftQ[3] + zshiftQ[6];
                Complex v2 = u1 * zshiftQ[1] + v1 * zshiftQ[4] + zshiftQ[7];
                Complex w2 = u1 * zshiftQ[2] + v1 * zshiftQ[5] + zshiftQ[8];
                
                Complex denom = w2;
                if (abs(denom) == 0.0) {
                  continue;
                }
                
                C2Point test_actsol(u2 / denom, v2 / denom);
                Real f = test_actsol.l2norm();
                if (f < minac) {
                  minac = f;
                  minac_row = j;
                  best_actsol = test_actsol;
                }
              }
            }
            
            // Update the solution.
            
            C2Point newactsol;
            newactsol.coord[0] = best_actsol.coord[0] * dila + transx;
            newactsol.coord[1] = best_actsol.coord[1] * dila + transy;
            C2Point spfval;
            
            // Now compute the residual.
            
            spfval.coord[0] = pl_eval(shiftpoly, newactsol, 0, degree);
            spfval.coord[1] = pl_eval(shiftpoly, newactsol, 1, degree);
            Real aspfval1 =  pow(abs(newactsol.coord[0]),degree);
            Real aspfval2 =  pow(abs(newactsol.coord[1]),degree);
            Real aspfval = (aspfval1 > aspfval2)? aspfval1 : aspfval2;
            aspfval = (aspfval > 1.0)? aspfval : 1.0;
            Real new_resd = spfval.l2norm() / aspfval;    
            Real zdresd =  zoom_act_solution_data(3,minac_row).real();
            

#ifdef DEBUGGING
            dump_everything = save_dump_everything;
            if (dump_everything) { 
              *debug_str << "--- tcorrection #" << improve_count << " vlrow = " << vlrow
                << " direct_try = " << direc_try_count << " invoc count = " << invocation_count
                << "\n eigval = " << Format_Complex(eigval) << "\n  this_resd = " 
                << this_resd << " new_resd = " << new_resd << " cresd = " << zdresd
                << "\n  neweigval = " << Format_Complex(zoom_act_solution_data(4,minac_row))
                << " minac = " << minac << " dila = " << dila << " new_eveccond = " 
                << zoom_act_solution_data(2,minac_row).real()
                << "\n  eveccond = " << evec_cond << " zoom_minsv = " << zoom_minsv
                << "\n best_actsol = [" << Format_Complex(best_actsol.coord[0])
                << "," << Format_Complex(best_actsol.coord[1])
                << "];\n dila*minac = " << dila * minac
                << '\n';
              *debug_str << "  % improvement loop, pass = " << improve_count << '\n';
              *debug_str << " zoom_coef = [\n";
              for (int eqnum = 0; eqnum < 2; ++eqnum) {
                for (int c = 0; c < num_control_pt; ++c) {
                  *debug_str << setprecision(17) << zoom_coef(eqnum,c);
                  if (c < num_control_pt - 1)
                    *debug_str << ",";
                }
                *debug_str << "\n";
              }
              *debug_str << "];\nact_solution_data = [\n";
              for (int vlrow1 = 0; vlrow1 < evsize; ++vlrow1) {
                *debug_str << Format_Complex(zoom_act_solution_data(0,vlrow1)) << "," <<
                  Format_Complex(zoom_act_solution_data(1,vlrow1)) << "," <<
                  zoom_act_solution_data(2,vlrow1).real() << "," << 
                  zoom_act_solution_data(3,vlrow1).real() << "," << 
                  Format_Complex(zoom_act_solution_data(4,vlrow1)) << "\n";
              }
              
              *debug_str << "];\nnewactsol = [" << Format_Complex(newactsol.coord[0]) 
                << "," << Format_Complex(newactsol.coord[1]) 
                << "];\nnew_spfval = " << spfval.l2norm() << ";\nnew_aspfval = " << aspfval
                << ";\nnew_resd = " << new_resd << ";\n";
            }
#endif
            
            // Big corrections not allowed.
            
            if (dila * minac > 0.01) {
              break;              
            }
            if (new_resd > this_resd) {
              break;
            }
            this_resd = new_resd;
            act_solution_data(0,vlrow) = newactsol.coord[0];
            act_solution_data(1,vlrow) = newactsol.coord[1];
          }
          
          max_resd = (max_resd > this_resd)? max_resd : this_resd;
          if (in_range && this_resd > max_resd2)
            max_resd2 = this_resd;
  

          // Apply the inverse of shiftQ to (u1,v1,w1) to undo
          // undo effect of shift.

          Complex u1 = act_solution_data(0,vlrow);
          Complex v1 = act_solution_data(1,vlrow);
          
          Complex u2 = u1 * shiftQ[0] + v1 * shiftQ[3] + shiftQ[6];
          Complex v2 = u1 * shiftQ[1] + v1 * shiftQ[4] + shiftQ[7];
          Complex w2 = u1 * shiftQ[2] + v1 * shiftQ[5] + shiftQ[8];
          
          Complex denom = w2;
          if (abs(denom) == 0.0) {
#ifdef DEBUGGING
            if (dump_everything)
              *debug_str << " discarding, denom = 0\n";
#endif
            continue;
          }
          
          C2Point sol(u2 / denom, v2 / denom);
          
#ifdef DEBUGGING
          if (dump_everything) {
            *debug_str << "\n u1i = " 
              << Format_Complex(u1) 
              << ";\n v1i = " << Format_Complex(v1) 
              << ";\n u2i = " << Format_Complex(u2) 
              << ";\n v2i = " << Format_Complex(v2) 
              << ";\nw2i = " << Format_Complex(w2)
              << ";\n denomi = " << Format_Complex(denom) 
              << ";\n sol = " << Format_Complex(sol.coord[0]) << " " << Format_Complex(sol.coord[1])
              << ";\n";
          }
#endif
          C2Point fval, fvalcopy, direcv, nstep, newsol, newfval;
          
          Real invnormest;
          
          // Try to improve the solution with two steps of Newton's method
          // if newton_improve is true.
          
          
          for (int newtcount = 0; true; ++newtcount) {            
            fval.coord[0] = eval_complex(sol, 0, degree) / bezmax_poly[0];
            fval.coord[1] = eval_complex(sol, 1, degree) / bezmax_poly[1];
            direcv.coord[0] = 1.0;
            direcv.coord[1] = 0.0;
            jac(0,0) = eval_direc_deriv_complex(sol, direcv, 0, degree) / bezmax_poly[0];
            jac(1,0) = eval_direc_deriv_complex(sol, direcv, 1, degree) / bezmax_poly[1];
            direcv.coord[0] = 0.0;
            direcv.coord[1] = 1.0;
            jac(0,1) = eval_direc_deriv_complex(sol, direcv, 0, degree) / bezmax_poly[0];
            jac(1,1) = eval_direc_deriv_complex(sol, direcv, 1, degree) / bezmax_poly[1];
            fvalcopy = fval;
            
            // solve 2-by-2 linear system
            
            invnormest = linear_solve_C2(jac, fval, nstep);

#ifdef DEBUGGING
            if (dump_everything) {
              *debug_str << "soli = [" 
                <<  Format_Complex(sol.coord[0]) << " ," 
                << Format_Complex(sol.coord[1])
                << "];\n fvali = [" << Format_Complex(fval.coord[0])
                << "," << Format_Complex(fval.coord[1])
                << "\n];\n jac1 = [" 
                << Format_Complex(jac(0,0)) 
                << "," << Format_Complex(jac(0,1)) 
                << ";" << Format_Complex(jac(1,0) )
                << " " << Format_Complex(jac(1,1));
              *debug_str << "];\n nstep = [" 
                << Format_Complex(nstep.coord[0])
                << "," << Format_Complex(nstep.coord[1])
                << "];\n invnormest = " << invnormest;
              *debug_str << ";\n";
            }
            
#endif
            if (!newton_improve) break;
            if (newtcount == 2) break;
            if (invnormest >= BIG_REAL) break;
            if (nstep.l2norm() >= 0.01) break;            
            
            newsol.coord[0] = sol.coord[0] - nstep.coord[0];
            newsol.coord[1] = sol.coord[1] - nstep.coord[1];
            newfval.coord[0] = eval_complex(newsol, 0, degree);
            newfval.coord[1] = eval_complex(newsol, 1, degree);
#ifdef DEBUGGING
            if (dump_everything) {
              *debug_str << " newsoli = [" 
                << Format_Complex(newsol.coord[0]) 
                << ", " << Format_Complex(newsol.coord[1])
                << "]; newfvali = [" 
                << Format_Complex(newfval.coord[0]) 
                << ", " << Format_Complex(newfval.coord[1]) 
                << "];\n";
            }
#endif
            if (newfval.l2norm() < fval.l2norm()) {
#ifdef DEBUGGING
              if (dump_everything) {
                *debug_str << " changing\n";
              }
#endif
              sol = newsol;
            }
            else {
              break;
            }
          }

          Real fvalnorm = fval.l2norm();
          Real fvalcutoff = evsize * sctol * tol;       
          Real param_uncertainty = (fvalnorm + fvalcutoff) * invnormest;
          if (param_uncertainty > 1.0)
            param_uncertainty = 1.0;


#ifdef DEBUGGING
          if (dump_everything) {
            *debug_str << " fvalnorm = " << fvalnorm << " fvalcutoff = "
              << fvalcutoff << " param_uncertainty = "
              << param_uncertainty <<  "\n";
          }
#endif
          if (abs(sol.coord[0]) + abs(sol.coord[1]) > 3) {
#ifdef DEBUGGING
            if (dump_everything) {
              *debug_str << " dropping root; norm too large";
            }
#endif
            continue;
          }
          Real delta1 = param_uncertainty;
          if (delta1 > 0.01)
            delta1 = 0.01;

#ifdef DEBUGGING
          if (dump_everything) {
            *debug_str << "testing against delta1 = " << delta1 << '\n';
          }
#endif
          
          if (sol.coord[0].real() >= -delta1 && 
            sol.coord[1].real() >= -delta1 &&
            1.0 - sol.coord[0].real() - sol.coord[1].real() >= -delta1 &&
            fabs(sol.coord[0].imag()) <= delta1 &&
            fabs(sol.coord[1].imag()) <= delta1) {
            
            test_solution_data(0,test_numsol) = sol.coord[0].real();
            test_solution_data(1,test_numsol) = sol.coord[1].real();
            test_solution_data(2,test_numsol) = param_uncertainty;
            
#ifdef DEBUGGING
            if (dump_everything) {
              *debug_str << "test passed, solution #" << test_numsol <<'\n';
              *debug_str << "solution = " << setprecision(17) 
                << test_solution_data(0,test_numsol) << "," << test_solution_data(1,test_numsol)
                << " uncer = " << test_solution_data(2, test_numsol) << '\n';
            }
#endif
            ++test_numsol;

          }
        }
        if (max_resd2 < best_max_resd2) {
#ifdef DEBUGGING
          if (dump_everything) {
            *debug_str << "keeping preceding solutions, test_numsol = "
              << test_numsol << " max_resd2 = " << max_resd2 
              << " previous best = " << best_max_resd2 << '\n';
          }
#endif
          best_max_resd2 = max_resd2;
          numsol = test_numsol;
          for (int j = 0; j < numsol; ++j) {
            for (int i = 0; i < 3; ++i) {
              solution_data(i,j) = test_solution_data(i,j);
            }
          }
        }
      }

      {
        PatchMath::IntersectRecord result;
        for (int j = 0; j < numsol; ++j) {
          result.paramcoord[0] = solution_data(0,j);
          result.paramcoord[1] = solution_data(1,j);
          result.param_uncertainty = solution_data(2,j);
          result.degen = fabs(result.paramcoord[0]) <= result.param_uncertainty ||
            fabs(result.paramcoord[1]) <= result.param_uncertainty ||
            fabs(1.0 - result.paramcoord[0] - result.paramcoord[1]) <= result.param_uncertainty;
          output_stack.push_back(result);
        }
      }
      return is_degen;
    }
  }


  
  // Horner algorithm for evaluating bidegree polynomial.


  Real pl_eval2r(const PatchMath::Workspace::RealMatrix& coef,
    const R2Point& param,
    int eqnum,
    int degree1,
    int degree2) {

    Real u = param.coord[0];
    Real v = param.coord[1];
    Real outer_result = 0.0;
    Real pwrv = 1.0;
    for (int d2 = 0; d2 <= degree2; ++d2) {
      // Compute the inner coefficient.
      Real pwru = 1.0;
      Real inner_result = 0.0;
      for (int d1 = 0; d1 <= degree1; ++d1) {
        inner_result += pwru * coef(eqnum, d2 * (degree1 + 1) + d1);
        pwru *= u;
      }
      outer_result += inner_result * pwrv;
      pwrv *= v;
    }
    return outer_result;
  }
  
  // Horner algorithm for evaluating bidegree polynomial.

  Complex pl_eval2(const PatchMath::Workspace::RealMatrix& coef,
    const C2Point& param,
    int eqnum,
    int degree1,
    int degree2) {

    Complex u = param.coord[0];
    Complex v = param.coord[1];
    Complex outer_result(0.0,0.0);
    Complex pwrv(1.0,0.0);
    for (int d2 = 0; d2 <= degree2; ++d2) {
      // Compute the inner coefficient.
      Complex pwru(1.0,0.0);
      Complex inner_result(0.0,0.0);
      for (int d1 = 0; d1 <= degree1; ++d1) {
        inner_result += pwru * coef(eqnum, d2 * (degree1 + 1) + d1);
        pwru *= u;
      }
      outer_result += inner_result * pwrv;
      pwrv *= v;
    }
    return outer_result;
  }

  // Helper class for math with bezier quad patches.

  class Solve_BezQuad {
  public:
    virtual Real get_bez_coeff(int eqnum, int j) const = 0;
    Solve_BezQuad() { }
    virtual ~Solve_BezQuad() { }
    bool solve(int degree1[2],
      int degree2[2],
      vector<PatchMath::IntersectRecord>& output_stack,
      PatchMath::Workspace& workspace, 
      double scfac, 
      double tol,
      bool newton_improve) const;
    Real eval(const R2Point& param,
      int eqnum,
      int degree1,
      int degree2) const;
    Complex eval_complex(const C2Point& param,
      int eqnum,
      int degree1,
      int degree2) const;
    Real eval_partialu(const R2Point& param,
      int eqnum,
      int degree1,
      int degree2) const;
    Complex eval_partialu_complex(const C2Point& param,
      int eqnum,
      int degree1,
      int degree2) const;
    Real eval_partialv(const R2Point& param,
      int eqnum,
      int degree1,
      int degree2) const;
    Complex eval_partialv_complex(const C2Point& param,
      int eqnum,
      int degree1,
      int degree2) const;
  };

  // Derived class for intersecting a patch with a line.
  
  class Solve_BezQuad1 : public Solve_BezQuad {
  private: 
    const Matrix& lineeq_;
    const Point& linerhs_;
    const PatchMath& pm_;
  public:
    Solve_BezQuad1(const Matrix& lineeq,
      const Point& linerhs,
      const PatchMath& pm) : lineeq_(lineeq), linerhs_(linerhs), pm_(pm) { }
    Real get_bez_coeff(int eqnum, int cpnum) const {
      return lineeq_(eqnum,0) * PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum, 0) +
        lineeq_(eqnum,1) * PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum, 1) +
        lineeq_(eqnum,2) * PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum, 2) -
        linerhs_[eqnum];
    }
  };

  // derived class for directional derivative.

  class Solve_BezQuad2 : public Solve_BezQuad {
  private:
    const Point& direc_;
    int di_;
    const PatchMath& pm_;
    int degree1_;
    int degree2_;
  public:
    Solve_BezQuad2(const Point& direc, 
      int di,
      const PatchMath& pm,
      int degree1,
      int degree2) : direc_(direc), di_(di), pm_(pm), degree1_(degree1), degree2_(degree2)
    { }
    Real get_bez_coeff(int eqnum, int cpnum) const;
  };

  // derived class for evaluating one dimension.
    
  class Solve_BezQuad3 : public Solve_BezQuad {
  private:
    int select_dim_;
    const PatchMath& pm_;
  public:
    Solve_BezQuad3(int select_dim, 
      const PatchMath& pm) : select_dim_(select_dim), pm_(pm) 
    { }
    Real get_bez_coeff(int, int cpnum) const {
      return PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum, select_dim_);
    }
  };

    
  Real Solve_BezQuad2::get_bez_coeff(int eqnum, int cpnum) const {
    Real t = 0.0;
    int ocp1, ocp2, fac;
    if (eqnum == 0) {
      int d2 = cpnum / degree1_;
      int d1 = cpnum % degree1_;
      ocp1 = d1 + 1 + d2 * (degree1_ + 1);
      ocp2 = d1 + d2 * (degree1_ + 1);
      fac = degree1_;
    }
    else {
      int d2 = cpnum / (degree1_ + 1);
      int d1 = cpnum % (degree1_ + 1);
      ocp1 = d1 + (d2 + 1) * (degree1_ + 1);
      ocp2 = d1 + d2 * (degree1_ + 1);
      fac = degree2_;
    }
    for (int j = 0; j < di_; ++j) {
        t += (PatchMath::Bezier_Helper::control_point_coord(pm_, ocp1, j) -
          PatchMath::Bezier_Helper::control_point_coord(pm_, ocp2, j)) * fac 
          * direc_[j];
    }
    return t;
  }

  // Horner algorithm for evaluating bezier quad (parametric
  // to real conversion)


  Real Solve_BezQuad::eval(const R2Point& param,
    int eqnum,
    int degree1,
    int degree2) const {

    Real u = param.coord[0];
    Real oppu = 1.0 - u;
    Real v = param.coord[1];
    Real oppv = 1.0 - v;

    int outer_binom = 1;
    Real outer_result = 0.0;
    Real pwrv = 1.0;
    for (int d2 = 0; d2 <= degree2; ++d2) {
      // Compute the inner coefficient.

      int inner_binom = 1;
      Real pwru = 1.0;
      Real inner_result = 0.0;
      for (int d1 = 0; d1 <= degree1; ++d1) {
        inner_result *= oppu;
        inner_result += get_bez_coeff(eqnum, d2 * (degree1 + 1) + d1) *
          inner_binom * pwru;
        pwru *= u;
        inner_binom *= (degree1 - d1);
        inner_binom /= (d1 + 1);
      }
      outer_result *= oppv;
      outer_result += inner_result * outer_binom * pwrv;
      pwrv *= v;
      outer_binom *= (degree2 - d2);
      outer_binom /= (d2 + 1);
    }
    return outer_result;
  }

  
  // Horner algorithm for evaluating bezier quad (parametric
  // to real conversion)


  Complex Solve_BezQuad::eval_complex(const C2Point& param,
    int eqnum,
    int degree1,
    int degree2) const {

    Complex u = param.coord[0];
    Complex oppu = Complex(1.0,0.0) - u;
    Complex v = param.coord[1];
    Complex oppv = Complex(1.0,0.0) - v;

    int outer_binom = 1;
    Complex outer_result(0.0,0.0);
    Complex pwrv(1.0,0.0);
    for (int d2 = 0; d2 <= degree2; ++d2) {
      // Compute the inner coefficient.

      int inner_binom = 1;
      Complex pwru(1.0,0.0);
      Complex inner_result(0.0,0.0);
      for (int d1 = 0; d1 <= degree1; ++d1) {
        inner_result *= oppu;
        inner_result += Complex(get_bez_coeff(eqnum, d2 * (degree1 + 1) + d1), 0.0) *
          Complex(static_cast<Real>(inner_binom), 0.0) * pwru;
        pwru *= u;
        inner_binom *= (degree1 - d1);
        inner_binom /= (d1 + 1);
      }
      outer_result *= oppv;
      outer_result += inner_result * Complex(static_cast<Real>(outer_binom), 0.0) * pwrv;
      pwrv *= v;
      outer_binom *= (degree2 - d2);
      outer_binom /= (d2 + 1);
    }
    return outer_result;
  }
  
  
  Real Solve_BezQuad::eval_partialu(const R2Point& param,
    int eqnu