qpatchmath.cpp

来自「算断裂的」· C++ 代码 · 共 2,107 行 · 第 1/5 页

 // ------------------------------------------------------------------
// 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 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, November 10, 1999.  Steve Vavasis
// ------------------------------------------------------------------
#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;
  }
}
#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);


// 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::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];
  }
  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_];
  }
  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];
  }
  Lapack_doublecomplex* make_fortran_complex() {
    return reinterpret_cast<Lapack_doublecomplex*>(
      static_cast<Complex*>(&(wkspa_.cwkspa_[base_])));
  }

};

class QMG::PatchMath::Workspace::RealMatrix {
private:
  Workspace& wkspa_;
  int base_;
  int nr_;
  int nc_;
  // 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) {
    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_];
  }
  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_;
  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()) {
    ++wkspa.marker_active_;
  }
  ~StartPosMarker() {
    wkspa_.rwkspa_.resize(rpos_);
    wkspa_.cwkspa_.resize(cpos_);
    --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;
  }





  // ------------------------------------------------------------------
  // This routine carries out QR factorization on a matrix.
  // If sing_cutoff is nonzero, it also tries to find
  // dependent columns and quits as soon as a dependent column is found.
  // If it quits early then it returns the number of columns processed.


  void QR_Factor(PatchMath::Workspace::RealMatrix& A, //matrix to factor
    PatchMath::Workspace::RealMatrix& hhmat,
    PatchMath::Workspace::RealVector& betavec) {

    int m = A.numrow();
    int n = A.numcol();
    int nstep = (m < n)? m : n;
    for (int k = 0; k < nstep; ++k) {
      Real* hhmat_col = &hhmat(0,k);

      Real beta = 
        Matrix::compute_hh_transform(&hhmat_col[k], &A(k,k), 1, m - k);
      Real mu = (beta == 0.0)? 0 : (1.0 / fabs(beta * hhmat_col[k]));
      betavec[k] = beta;
      {
        for (int j = k; j < n; ++j) {
          Real ip1 = 0.0;
          {
            for (int i = k; i < m; ++i)
              ip1 += A(i,j) * hhmat_col[i];
          }
          ip1 *= beta;
          {
            for (int i = k; i < m; ++i)
              A(i,j) -= ip1 * hhmat_col[i];
          }
        }
      }
    }
  }


  void QR_Factor_with_col_pivoting(PatchMath::Workspace::RealMatrix& A,