C++ DVec::data方法代码示例

// compute eigensystem of a real symmetric matrix
//---------------------------------------------------------
void eig_sym(const DMat& A, DVec& ev, DMat& Q, bool bDoEVecs)
//---------------------------------------------------------
{
  if (!A.is_square()) { umERROR("eig_sym(A)", "matrix is not square."); }

  int N = A.num_rows();
  int LDA=N, LDVL=N, LDVR=N, ldwork=10*N, info=0;
  DVec work(ldwork, 0.0, OBJ_temp, "work_TMP");

  Q = A;          // Calculate eigenvectors in Q (optional)
  ev.resize(N);   // Calculate eigenvalues in ev

  char jobV = bDoEVecs ? 'V' : 'N';

  SYEV (jobV,'U', N, Q.data(), LDA, ev.data(), work.data(), ldwork, info);  

  if (info < 0) { 
    umERROR("eig_sym(A, Re,Im)", "Error in input argument (%d)\nNo solution computed.", -info);
  } else if (info > 0) {
    umLOG(1, "eig_sym(A, W): ...\n"
             "\nthe algorithm failed to converge;"
             "\n%d off-diagonal elements of an intermediate"
             "\ntridiagonal form did not converge to zero.\n", info);
  }
}

// DPOSV uses Cholesky factorization A=U^T*U, A=L*L^T 
// to compute the solution to a real system of linear 
// equations A*X=B, where A is a square, (N,N) symmetric 
// positive definite matrix and X and B are (N,NRHS).
//---------------------------------------------------------
void umSOLVE_CH(const DMat& mat, const DVec& b, DVec& x)
//---------------------------------------------------------
{
  // check args
  assert(mat.is_square());            // symmetric
  assert(b.size() >= mat.num_rows()); // is b consistent?
  assert(b.size() <= x.size());       // can x store solution?
  
  DMat A(mat);    // work with copy of input
  x = b;          // allocate solution vector

  int rows=A.num_rows(), LDA=A.num_rows(), cols=A.num_cols();
  int  LDB=b.size(), NRHS=1, info=0;
  if (rows<1) {umWARNING("umSOLVE_CH()", "system is empty"); return;}

  // Solve the system.
  POSV('U', rows, NRHS, A.data(), LDA, x.data(), LDB, info);

  if (info < 0) { 
    x = 0.0;
    umERROR("umSOLVE_CH(A,b, x)", 
            "Error in input argument (%d)\nNo solution computed.", -info);
  } else if (info > 0) {
    x = 0.0;
    umERROR("umSOLVE_CH(A,b, x)", 
            "\nINFO = %d.  The leading minor of order %d of A"
            "\nis not positive definite, so the factorization" 
            "\ncould not be completed. No solution computed.", 
              info, info);
  }
}

//---------------------------------------------------------
void Poly3D::AddPoint(const DVec& point) 
//---------------------------------------------------------
{
  if (!HavePoint(point)) {
    // append this point
    m_xyz.append_col(3, (double*)point.data());
    ++m_N;
  }
}

//---------------------------------------------------------
void eig(const DMat& A, DVec& Re, DMat& VL, DMat& VR, bool bL, bool bR)
//---------------------------------------------------------
{
  // Compute eigensystem of a real general matrix
  // Currently NOT returning imaginary components

  static DMat B;

  if (!A.is_square()) { umERROR("eig(A)", "matrix is not square."); }

  int N = A.num_rows();
  int LDA=N, LDVL=N, LDVR=N, ldwork=10*N, info=0;

  Re.resize(N);     // store REAL components of eigenvalues in Re
  VL.resize(N,N);   // storage for LEFT eigenvectors
  VR.resize(N,N);   // storage for RIGHT eigenvectors
  DVec Im(N);     // NOT returning imaginary components
  DVec work(ldwork, 0.0);

  // Work on a copy of A
  B = A;

  char jobL = bL ? 'V' : 'N';   // calc LEFT eigenvectors?
  char jobR = bR ? 'V' : 'N';   // calc RIGHT eigenvectors?

  GEEV (jobL,jobR, N, B.data(), LDA, Re.data(), Im.data(), 
        VL.data(), LDVL, VR.data(), LDVR, work.data(), ldwork, info);

  if (info < 0) { 
    umERROR("eig(A, Re,Im)", "Error in input argument (%d)\nNo solution computed.", -info);
  } else if (info > 0) {
    umLOG(1, "eig(A, Re,Im): ...\n"
             "\nThe QR algorithm failed to compute all the"
             "\neigenvalues, and no eigenvectors have been" 
             "\ncomputed;  elements %d+1:N of WR and WI contain"
             "\neigenvalues which have converged.\n", info);
  }

#if (0)
  // Return (Re,Imag) parts of eigenvalues as columns of Ev
  Ev.resize(N,2);
  Ev.set_col(1, Re);
  Ev.set_col(2, Im);
#endif

#ifdef _DEBUG
    //#####################################################
    // Check for imaginary components in eigenvalues
    //#####################################################
    double im_max = Im.max_val_abs();
    if (im_max > 1e-6) {
      umERROR("eig(A)", "imaginary components in eigenvalues.");
    }
    //#####################################################
#endif
}

//---------------------------------------------------------
bool Poly3D::HavePoint(const DVec& point) 
//---------------------------------------------------------
{
  const double* p1 = point.data();

  double tol = 1e-6, normi=0.0; DVec pnti,tv;
  for (int i=1; i<=m_N; ++i) 
  {
    const double* p2 = m_xyz.pCol(i);
    normi = sqrt( SQ(p1[0]-p2[0]) + 
                  SQ(p1[1]-p2[1]) + 
                  SQ(p1[2]-p2[2]) );
    if (normi < tol) { return true; }
  }
  return false;
}

//---------------------------------------------------------
DVec& chol_solve(const DMat& ch, const DVec& b)
//---------------------------------------------------------
{
  // Solves a linear system using Cholesky-factored 
  // symmetric positive-definite matrix, A = U^T U.

  if (FACT_CHOL != ch.get_factmode()) {umERROR("chol_solve(ch,b)", "matrix is not factored.");}
  int M=ch.num_rows(), lda=ch.num_rows(); 
  int nrhs=1, ldb=b.size();   assert(ldb == M);
  char uplo = 'U';  int info=0; 
  double* ch_data = const_cast<double*>(ch.data());

  // copy RHS into x, then overwrite x with solution
  DVec* x = new DVec(b, OBJ_temp);
  POTRS (uplo, M, nrhs, ch_data, lda, x->data(), ldb, info);
  if (info) { umERROR("chol_solve(ch,b)", "dpotrs reports: info = %d", info); }
  return (*x);
}

// Computes an SVD factorization of a real MxN matrix.
// Returns the vector of singular values.
// Also, factors U, VT, where A = U * D * VT.
//---------------------------------------------------------
DVec& svd
(
  const DMat& mat,  // [in]
        DMat& U,    // [out: left singular vectors]
        DMat& VT,   // [out: right singular vectors]
        char ju,    // [in: want U?]
        char jvt    // [in: want VT?]
)
//---------------------------------------------------------
{
  // Work with a copy of the input matrix.
  DMat A(mat, OBJ_temp, "svd.TMP");

  // A(MxN)
  int m=A.num_rows(), n=A.num_cols();
  int mmn=A.min_mn(), xmn=A.max_mn();

  // resize parameters
  U.resize (m,m, true, 0.0);
  VT.resize(n,n, true, 0.0);
  DVec* s = new DVec(mmn, 0.0, OBJ_temp, "s.TMP");
  char jobu  = ju;
  char jobvt = jvt;
  int info = 0;

  // NBN: ACML does not use the work vector.
  int lwork = 2 * std::max(3*mmn+xmn, 5*mmn);
  DVec work(lwork, 0.0, OBJ_temp, "work.TMP");
  GESVD (jobu, jobvt, m, n, A.data(), m, s->data(), U.data(), m, VT.data(), n, work.data(), lwork, info);

  if (info < 0) { 
    (*s) = 0.0;
    umERROR("SVD", "Error in input argument (%d)\nNo solution computed.", -info);
  } else if (info > 0) {
    (*s) = 0.0;
    umLOG(1, "DBDSQR did not converge."
             "\n%d superdiagonals of an intermediate bidiagonal"
             "\nform B did not converge to zero.\n", info);
  }

  return (*s);
}

// DGESV computes the solution to a real system of linear 
// equations, A*x = b, where A is an N-by-N matrix, and 
// x and b are N-by-1 vectors.  The LU decomposition 
// with partial pivoting and row interchanges is used to 
// factor A as A = P*L*U, where P is a permutation matrix,
// L is unit lower triangular, and U is upper triangular.
// The system is solved using this factored form of A.
//---------------------------------------------------------
void umSOLVE(const DMat& mat, const DVec& b, DVec& x)
//---------------------------------------------------------
{
  // Work with copies of input arrays.
  DMat A(mat);
  x = b;

  int NRHS = 1;
  int LDA  = A.num_rows();
  int rows = A.num_rows();
  int cols = A.num_cols();
  int info = 0;
  
  if (rows != cols) {
    umERROR("umSOLVE(DMat, DVec)", 
            "Matrix A (%d,%d) is not square.\n"
            "For a Least-Squares solution, see umSOLVE_LS(A,B).", 
            rows, cols);
  }

  if (rows < 1) {
    umLOG(1, "Empty system passed into umSOLVE().\n");
    return;
  }

  IVec ipiv(rows, 0);

  GESV (rows, NRHS, A.data(), LDA, ipiv.data(), x.data(), rows, info);

  if (info < 0) { 
    x = 0.0;
    umERROR("umSOLVE(DMat&, DVec&)", 
            "Error in input argument (%d)\nNo solution computed.", -info);
  } else if (info > 0) {
    x = 0.0;
    umERROR("umSOLVE(DMat&, DVec&)", 
            "\nINFO = %d.  U(%d,%d) was exactly zero."
            "\nThe factorization has been completed, but the factor U is "
            "\nexactly singular, so the solution could not be computed.", 
              info, info, info);
  }
}

//---------------------------------------------------------
void umPOLISH(DVec& V, double eps)
//---------------------------------------------------------
{
  // round elements close to certain values

  int N = V.size();
  double *p = V.data();

  for (int i=0; i<N; ++i) 
  {
    if (fabs(p[i]) < eps) 
    {
      p[i] = 0.0;
    }
    else
    {
      if (p[i] > 0.0) 
      {
        // check for proximity to certain positive values
        if      (fabs (p[i] - 0.10) < eps) { p[i] = 0.10; }
        else if (fabs (p[i] - 0.20) < eps) { p[i] = 0.20; }
        else if (fabs (p[i] - 0.25) < eps) { p[i] = 0.25; }
        else if (fabs (p[i] - 0.50) < eps) { p[i] = 0.50; }
        else if (fabs (p[i] - 0.75) < eps) { p[i] = 0.75; }
        else if (fabs (p[i] - 0.80) < eps) { p[i] = 0.80; }
        else if (fabs (p[i] - 0.90) < eps) { p[i] = 0.90; }
        else if (fabs (p[i] - 1.00) < eps) { p[i] = 1.00; }
        else if (fabs (p[i] - 2.00) < eps) { p[i] = 2.00; }
        else if (fabs (p[i] - 4.00) < eps) { p[i] = 4.00; }
        else if (fabs (p[i] - 4.50) < eps) { p[i] = 4.50; }
        else if (fabs (p[i] - 5.00) < eps) { p[i] = 5.00; }

        else if (fabs (p[i] - M_PI  ) < eps) { p[i] = M_PI  ; }
        else if (fabs (p[i] - M_PI_2) < eps) { p[i] = M_PI_2; }
        else if (fabs (p[i] - M_PI_4) < eps) { p[i] = M_PI_4; }
        else if (fabs (p[i] - M_E   ) < eps) { p[i] = M_E   ; }
      }
      else
      {
        // check for proximity to certain negative values
        if      (fabs (p[i] + 0.10) < eps) { p[i] = -0.10; }
        else if (fabs (p[i] + 0.20) < eps) { p[i] = -0.20; }
        else if (fabs (p[i] + 0.25) < eps) { p[i] = -0.25; }
        else if (fabs (p[i] + 0.50) < eps) { p[i] = -0.50; }
        else if (fabs (p[i] + 0.75) < eps) { p[i] = -0.75; }
        else if (fabs (p[i] + 0.80) < eps) { p[i] = -0.80; }
        else if (fabs (p[i] + 0.90) < eps) { p[i] = -0.90; }
        else if (fabs (p[i] + 1.00) < eps) { p[i] = -1.00; }
        else if (fabs (p[i] + 2.00) < eps) { p[i] = -2.00; }
        else if (fabs (p[i] + 4.00) < eps) { p[i] = -4.00; }
        else if (fabs (p[i] + 4.50) < eps) { p[i] = -4.50; }
        else if (fabs (p[i] + 5.00) < eps) { p[i] = -5.00; }

        else if (fabs (p[i] + M_PI  ) < eps) { p[i] = -M_PI  ; }
        else if (fabs (p[i] + M_PI_2) < eps) { p[i] = -M_PI_2; }
        else if (fabs (p[i] + M_PI_4) < eps) { p[i] = -M_PI_4; }
        else if (fabs (p[i] + M_E   ) < eps) { p[i] = -M_E   ; }
      }
    }
  }
}

// function call ///////////////////////////////////////////////////////////////
double DoubleModel::operator()(const DVec &arg, const DVec &par) const
{
    checkSize(arg.size(), par.size());
    
    return (*this)(arg.data(), par.data());
}