C++ DMat::resize方法代码示例

//---------------------------------------------------------
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
}

//---------------------------------------------------------
void umAxB(const DMat& A, const DMat& B, DMat& C)
//---------------------------------------------------------
{
  //-------------------------
  // C = A * B
  //-------------------------
  // A = op(A) is (M,K)
  // B = op(B) is (K,N)
  //        C  is (M,N)
  //-------------------------
  int M=A.num_rows(), K=A.num_cols(), N=B.num_cols();
  int LDA=M, LDB=K, LDC=M;
  double one=1.0, zero=0.0;
  if (B.num_rows() != K) { umERROR("umAxB(A,B,C)", "wrong dimensions"); }
  C.resize(M,N);

  GEMM ('N','N',M,N,K, one,A.data(),LDA, 
                           B.data(),LDB, 
                      zero,C.data(),LDC);
}

double ResF4toM2Interface::setDegreeZeroMap(SchreyerFrame& C,
                                       DMat<RingType>& result,
                                       int slanted_degree,
                                       int lev)
// 'result' should be previously initialized, but will be resized.
// return value: -1 means (slanted_degree, lev) is out of range, and the zero matrix was returned.
//   otherwise: the fraction of non-zero elements is returned.
{
  // As above, get the size of the matrix, and 'newcols'
  // Now we loop through the elements of degree 'slanted_degree + lev' at level 'lev'
  const RingType& R = result.ring();
  if (not (lev > 0 and lev <= C.maxLevel()))
    {
      result.resize(0,0);
      return -1;
    }
  assert(lev > 0 and lev <= C.maxLevel());
  int degree = slanted_degree + lev;
  auto& thislevel = C.level(lev);
  int ncols = 0;
  for (auto p=thislevel.begin(); p != thislevel.end(); ++p)
    {
      if (p->mDegree == degree) ncols++;
    }

  auto& prevlevel = C.level(lev-1);
  int* newcomps = new int[prevlevel.size()];
  int nrows = 0;
  for (int i=0; i<prevlevel.size(); i++)
    if (prevlevel[i].mDegree == degree)
      newcomps[i] = nrows++;
    else
      newcomps[i] = -1;

  result.resize(nrows, ncols);

  int col = 0;
  long nnonzeros = 0;  
  for (auto p=thislevel.begin(); p != thislevel.end(); ++p)
    {
      if (p->mDegree != degree) continue;
      auto& f = p->mSyzygy;
      auto end = poly_iter(C.ring(), f, 1);
      auto i = poly_iter(C.ring(), f);
      for ( ; i != end; ++i)
        {
          long comp = C.monoid().get_component(i.monomial());
          if (newcomps[comp] >= 0)
            {
              R.set_from_long(result.entry(newcomps[comp], col), C.gausser().coeff_to_int(i.coefficient()));
              nnonzeros++;
            }
        }
      ++col;
    }
  double frac_nonzero = (nrows*ncols);
  frac_nonzero = static_cast<double>(nnonzeros) / frac_nonzero;

  delete[] newcomps;

  return frac_nonzero;
}

// DGELSS computes minimum norm solution to a real linear 
// least squares problem:   Minimize 2-norm(| b - A*x |).   
// using the singular value decomposition (SVD) of A. 
// A is an M-by-N matrix which may be rank-deficient.   
//---------------------------------------------------------
void umSOLVE_LS(const DMat& mat, const DMat& B, DMat& X)
//---------------------------------------------------------
{
  if (!mat.ok()) {umWARNING("umSOLVE_LS()", "system is empty"); return;}

  DMat A(mat);    // work with copy of input.

  int rows=A.num_rows(), cols=A.num_cols(), mmn=A.min_mn();
  int LDB=A.max_mn(), NRHS=B.num_cols();
  if (rows!=B.num_rows()) {umERROR("umSOLVE_LS(A,B)", "Inconsistant matrix sizes.");}

  DVec s(mmn);    // allocate array for singular values

  // X must be big enough to store various results.
  // Resize X so that its leading dimension = max(M,N), 
  // then load the set of right hand sides.

  X.resize(LDB,NRHS, true, 0.0);

  for (int j=1; j<=NRHS; ++j)     // loop across colums
    for (int i=1; i<=rows; ++i)   // loop down rows
      X(i,j) = B(i,j);

  // RCOND is used to determine the effective rank of A.   
  // Singular values S(i) <= RCOND*S(1) are treated as zero.   
  // If RCOND < 0, machine precision is used instead.   

//double rcond =  1.0 / 1.0e16;
  double rcond = -1.0;

  // NBN: ACML does not use the work vector.
  int mnLo=A.min_mn(), mnHi=A.max_mn(), rank=1, info=1;
  int lwork = 10*mnLo + std::max(2*mnLo, std::max(mnHi, NRHS));
  DVec work(lwork); 

  // Solve the system
  GELSS (rows, cols, NRHS, A.data(), rows, X.data(), LDB, s.data(), rcond, rank, work.data(), lwork, info);

  //---------------------------------------------
  // Report:
  //---------------------------------------------

  if (info == 0) {
    umLOG(1, "umSOLVE_LS reports successful LS-solution."
             "\nRCOND = %0.6e, "
             "\nOptimal length of work array was %d\n", rcond, lwork);
  } 
  else 
  {
    if (info < 0) { 
      X = 0.0;
      umERROR("umSOLVE_LS(DMat&, DMat&)", 
              "Error in input argument (%d)\nNo solution or error bounds computed.", -info);

    } else if (info > 0) {
      X = 0.0;
      umERROR("umSOLVE_LS(DMat&, DMat&)", 
          "\nThe algorithm for computing the SVD failed to converge.\n"
          "\n%d off-diagonal elements of an intermediate "
          "\nbidiagonal form did not converge to zero.\n "
          "\nRCOND = %0.6e, "
          "\nOptimal length of work array was %d.\n", info, rcond, lwork);
    }
  }
}

//---------------------------------------------------------
void NDG2D::OutputSampleXYZ
(
        int sample_N,
        DMat &newX, 
        DMat &newY, 
        DMat &newZ,     // e.g. triangles on a sphere
  const DMat &FData,    // old field data
        DMat &newFData, // new field data
        int zfield      // if>0, use as z-elevation
)
//---------------------------------------------------------
{
  DVec newR, newS, newT;
  DMat newVDM;
  int newNpts = 0;

  // Triangles
  OutputSampleNodes2D(sample_N, newR, newS);
  newNpts = newR.size();
  newVDM = Vandermonde2D(this->N, newR, newS);

  const DMat& oldV = this->V;
  DMat oldtonew(newNpts, this->Np, "OldToNew");
  oldtonew = trans(trans(oldV) | trans(newVDM));

  //-----------------------------------
  // interpolate the field data
  //-----------------------------------
  int Nfields = FData.num_cols();
  newFData.resize(newNpts*this->K, Nfields);
  //DVec scales(Nfields);

  // For each field, use tOldF to wrap field i.
  // Use tNewF to load the interpolated field
  // directly into column i of the output array.
  DMat tOldF, tNewF;
  for (int i=1; i<=Nfields; ++i) {
    tOldF.borrow(this->Np, this->K, (double*)   FData.pCol(i));
    tNewF.borrow(newNpts,  this->K, (double*)newFData.pCol(i));
    tNewF = oldtonew * tOldF;
  //scales(i) = tNewF.max_col_val_abs(i);
  }

  //-----------------------------------
  // interpolate the vertices
  //-----------------------------------
  newX = oldtonew * this->x;
  newY = oldtonew * this->y;

  if (this->bCoord3D) {
    newZ = oldtonew * this->z;
  } 
  else 
  {
    if (zfield>=1 && zfield<=Nfields) {
      // use field data for z-height
      newZ.load(newNpts, K, newFData.pCol(Nfields));
    } else {
      // set z-data to 0.0
      newZ.resize(newNpts, K, true, 0.0);
    }
  }
}