C++ DVec类代码示例

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

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

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

// access //////////////////////////////////////////////////////////////////////
void TabFunction::setData(const DVec &x, const DVec &y)
{
    if (x.size() != y.size())
    {
        LATAN_ERROR(Size, "tabulated function x/y data size mismatch");
    }
    FOR_VEC(x, i)
    {
        value_[x(i)] = y(i);
    }

int main(int argc, char *argv[]) {


  // preprocess data
  NodeFiles node_files;
  std::vector<Block> col_blocks;
  // std::vector<string> in_files;
  AssignData(&node_files, &col_blocks);
  // PreprocessData(node_files, col_blocks, &in_files);

  std::vector<string> in_files = node_files[FLAGS_my_row_rank];
  int nf = in_files.size();

  Block col = col_blocks[FLAGS_my_col_rank];

  RSpMat<size_t> tmp;
  RSpMat<uint32_t> *adjs = new RSpMat<uint32_t>[nf];
  for (int i = 0; i < nf; ++i) {
    // load data
    tmp.Load(in_files[i]);
    tmp.VSlice(col, adjs+i);
  }


  // do actual computing
  // w = alpha * X * w + (1-alpha)*1/n
  // TODO
  CHECK_EQ(nf, 1);
  CHECK(adjs[0].square());

    int f = 0;

    uint32_t* index = adjs[f].index();
    size_t* offset = adjs[f].offset();
    double penalty = (1 - FLAGS_alpha) / (double) adjs[f].rows();

    DVec w = DVec::Ones(col.size()) * penalty;
    DVec u(adjs[f].rows());

  for (int it = 0; it < 40; ++it) {
    for (size_t i = 0; i < adjs[f].rows(); ++i) {
      double v = 0;
      double degree = offset[i+1] - offset[i];
      for (uint32_t j = offset[i]; j < offset[i+1]; ++j) {
        v += w[index[j]] / degree;
      }
      u[i] = v*FLAGS_alpha + penalty;
    }
    LL << "iter " << it << " err " << (u-w).norm() / w.norm()
       << " 1-norm " << w.cwiseAbs().sum();
    w = u;
  }

  return 0;
}

//---------------------------------------------------------
void CurvedINS2D::INScylinderBC2D
(
  const DVec&   xin,    // [in]
  const DVec&   yin,    // [in]
  const DVec&   nxi,    // [in]
  const DVec&   nyi,    // [in]
  const IVec&   MAPI,   // [in]
  const IVec&   MAPO,   // [in]
  const IVec&   MAPW,   // [in]
  const IVec&   MAPC,   // [in]
        double  ti,     // [in]
        double  nu,     // [in]
        DVec&   BCUX,   // [out]
        DVec&   BCUY,   // [out]
        DVec&   BCPR,   // [out]
        DVec&   BCDUNDT // [out]
)
//---------------------------------------------------------
{
  // function [bcUx, bcUy, bcPR, bcdUndt] = INScylinderBC2D(x, y, nx, ny, mapI, mapO, mapW, mapC, time, nu)
  // Purpose: evaluate boundary conditions for channel bounded cylinder flow with walls at y=+/- .15

  // TEST CASE: from 
  // V. John "Reference values for drag and lift of a two-dimensional time dependent flow around a cylinder", 
  // Int. J. Numer. Meth. Fluids 44, 777 - 788, 2004

  DVec yI("yI");  int len = xin.size();
  BCUX.resize(len); BCUY.resize(len);     // resize result arrays
  BCPR.resize(len); BCDUNDT.resize(len);  // and set to zero

  // inflow
#ifdef _MSC_VER
  yI = yin(MAPI);  yI += 0.20;
#else
  int Ni=MAPI.size(); yI.resize(Ni);
  for(int n=1;n<=Ni;++n) {yI(n)=yin(MAPI(n))+0.2;}
#endif

  BCUX(MAPI)    =  SQ(1.0/0.41)*6.0 * yI.dm(0.41 - yI);
  BCUY(MAPI)    =  0.0;
  BCDUNDT(MAPI) = -SQ(1.0/0.41)*6.0 * yI.dm(0.41 - yI);

  // wall
  BCUX(MAPW) = 0.0;
  BCUY(MAPW) = 0.0;

  // cylinder
  BCUX(MAPC) = 0.0;
  BCUY(MAPC) = 0.0;

  // outflow
  BCUX(MAPO)    = 0.0;
  BCUY(MAPO)    = 0.0;
  BCDUNDT(MAPO) = 0.0;
}

//---------------------------------------------------------
DVec& lu_solve(DMat& LU, const DVec& b)
//---------------------------------------------------------
{
  // Solve a linear system using lu-factored square matrix.

  DVec *x = new DVec("x", OBJ_temp);
  try {
    LU.solve_LU(b, (*x), false, false);
  } catch(...) { x->Free(); }
  return (*x);
}

DVec OCCFace::inertia() {
    DVec ret;
    GProp_GProps prop;
    BRepGProp::SurfaceProperties(this->getShape(), prop);
    gp_Mat mat = prop.MatrixOfInertia();
    ret.push_back(mat(1,1)); // Ixx
    ret.push_back(mat(2,2)); // Iyy
    ret.push_back(mat(3,3)); // Izz
    ret.push_back(mat(1,2)); // Ixy
    ret.push_back(mat(1,3)); // Ixz
    ret.push_back(mat(2,3)); // Iyz
    return ret;
}

// load {R,S} output nodes
void OutputSampleNodes2D(int sample_N, DVec &R, DVec &S)
{
  int Npts = OutputSampleNpts2D(sample_N);
  R.resize(Npts); S.resize(Npts);
  double denom = (double)(sample_N);
  int sampleNq = sample_N+1;
  for (int sk=0, i=0; i<sampleNq; ++i) {
    for (int j=0; j<sampleNq-i; ++j, ++sk) {
      R[sk] = -1. + (2.*j)/denom;
      S[sk] = -1. + (2.*i)/denom;
    }
  }
}

//---------------------------------------------------------
void GradSimplex3DP
(
  const DVec& a,      // [in]
  const DVec& b,      // [in]
  const DVec& c,      // [in]
        int   id,     // [in]
        int   jd,     // [in]
        int   kd,     // [in]
        DVec& V3Dr,   // [out]
        DVec& V3Ds,   // [out]
        DVec& V3Dt    // [out]
)
//---------------------------------------------------------
{
  // function [V3Dr, V3Ds, V3Dt] = GradSimplex3DP(a,b,c,id,jd,kd)
  // Purpose: Return the derivatives of the modal basis (id,jd,kd) 
  //          on the 3D simplex at (a,b,c)

  DVec fa, dfa, gb, dgb, hc, dhc, tmp;


  fa = JacobiP(a,0,0,id);           dfa = GradJacobiP(a,0,0,id);
  gb = JacobiP(b,2*id+1,0,jd);      dgb = GradJacobiP(b,2*id+1,0,jd);
  hc = JacobiP(c,2*(id+jd)+2,0,kd); dhc = GradJacobiP(c,2*(id+jd)+2,0,kd);

  // r-derivative
  V3Dr = dfa.dm(gb.dm(hc));
  if (id>0)    { V3Dr *= pow(0.5*(1.0-b), double(id-1)); }
  if (id+jd>0) { V3Dr *= pow(0.5*(1.0-c), double(id+jd-1)); }

  // s-derivative 
  V3Ds = V3Dr.dm(0.5*(1.0+a));
  tmp = dgb.dm(pow(0.5*(1.0-b), double(id)));

  if (id>0)    { tmp -= (0.5*id) * gb.dm(pow(0.5*(1.0-b), (id-1.0))); }
  if (id+jd>0) { tmp *= pow(0.5*(1.0-c),(id+jd-1.0)); }
  tmp = fa.dm(tmp.dm(hc));
  V3Ds += tmp;

  // t-derivative 
  V3Dt = V3Dr.dm(0.5*(1.0+a)) + tmp.dm(0.5*(1.0+b));
  tmp = dhc.dm(pow(0.5*(1.0-c), double(id+jd)));
  if (id+jd>0) { tmp -= (0.5*(id+jd)) * hc.dm(pow(0.5*(1.0-c), (id+jd-1.0))); }
  tmp = fa.dm(gb.dm(tmp));
  tmp *= pow(0.5*(1.0-b), double(id));
  V3Dt += tmp;

  // normalize
  double fac = pow(2.0, (2.0*id + jd+1.5));
  V3Dr *= fac;  V3Ds *= fac;  V3Dt *= fac;
}

//---------------------------------------------------------
void xyztorst
(
  const DVec& X,  // [in]
  const DVec& Y,  // [in]
  const DVec& Z,  // [in]
        DVec& r,  // [out]
        DVec& s,  // [out]
        DVec& t   // [out]
)
//---------------------------------------------------------
{
  // function [r,s,t] = xyztorst(x, y, z)
  // Purpose : Transfer from (x,y,z) in equilateral tetrahedron
  //           to (r,s,t) coordinates in standard tetrahedron

  double sqrt3=sqrt(3.0), sqrt6=sqrt(6.0); int Nc=X.size();
  DVec v1(3),v2(3),v3(3),v4(3);
  DMat tmat1(3,Nc), A(3,3), rhs; 
  
  v1(1)=(-1.0);  v1(2)=(-1.0/sqrt3);  v1(3)=(-1.0/sqrt6);
  v2(1)=( 1.0);  v2(2)=(-1.0/sqrt3);  v2(3)=(-1.0/sqrt6);
  v3(1)=( 0.0);  v3(2)=( 2.0/sqrt3);  v3(3)=(-1.0/sqrt6);
  v4(1)=( 0.0);  v4(2)=( 0.0      );  v4(3)=( 3.0/sqrt6);

  // back out right tet nodes
  tmat1.set_row(1,X); tmat1.set_row(2,Y); tmat1.set_row(3,Z);
  rhs = tmat1 - 0.5*outer(v2+v3+v4-v1, ones(Nc));
  A.set_col(1,0.5*(v2-v1)); A.set_col(2,0.5*(v3-v1)); A.set_col(3,0.5*(v4-v1));

  DMat RST = A|rhs;

  r=RST.get_row(1); s=RST.get_row(2); t=RST.get_row(3);
}

//==================================================================
bool ParamsFindP(	ParamList &params,
					const SymbolList &globalSymbols,
					DVec<Float3> &out_vectorP,
					int fromIdx )
{
	bool	gotP = false;

	for (size_t i=fromIdx; i < params.size(); i += 2)
	{
		DASSERT( params[i].type == Param::STR );

		const Symbol* pSymbol = globalSymbols.FindSymbol( params[i] );
		if ( pSymbol && pSymbol->IsName( "P" ) )
		{
			DASSTHROW( (i+1) < params.size(), "Invalid number of arguments" );
			
			const FltVec	&fltVec = params[ i+1 ].NumVec();
			
			DASSTHROW( (fltVec.size() % 3) == 0, "Invalid number of arguments" );
			
			out_vectorP.resize( fltVec.size() / 3 );

			for (size_t iv=0, id=0; iv < fltVec.size(); iv += 3)
				out_vectorP[id++] = Float3( &fltVec[ iv ] );

			return true;
		}
	}
	
	return false;
}

//---------------------------------------------------------
void Maxwell2D::SetIC()
//---------------------------------------------------------
{
  // Set initial conditions for simulation

  mmode = 1.0; nmode = 1.0;

  //Ez = sin(mmode*pi*x) .* sin(nmode*pi*y);
  DVec tsinx = apply(sin, (mmode*pi*x)),
       tsiny = apply(sin, (nmode*pi*y));
  
  Ezinit = tsinx.dm(tsiny);

  Ez = Ezinit;
  Hx = 0.0;
  Hy = 0.0;
}

//==================================================================
void MemFile::InitExclusiveOwenership( DVec<U8> &fromData )
{
	mDataSize = fromData.size();
	mOwnData.get_ownership( fromData );

	mpData		= &mOwnData[0];
	mReadPos	= 0;
	mIsReadOnly	= true;
}

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