C++ Dim函数代码示例

  void DenseSymMatrix::HighRankUpdateTranspose(Number alpha,
      const MultiVectorMatrix& V1,
      const MultiVectorMatrix& V2,
      Number beta)
  {
    DBG_ASSERT(Dim()==V1.NCols());
    DBG_ASSERT(Dim()==V2.NCols());
    DBG_ASSERT(beta==0. || initialized_);

    const Index dim = Dim();
    if (beta==0.) {
      for (Index j=0; j<dim; j++) {
        for (Index i=j; i<dim; i++) {
          values_[i+j*dim] = alpha*V1.GetVector(i)->Dot(*V2.GetVector(j));
        }
      }
    }
    else {
      for (Index j=0; j<dim; j++) {
        for (Index i=j; i<dim; i++) {
          values_[i+j*dim] = alpha*V1.GetVector(i)->Dot(*V2.GetVector(j))
                             + beta*values_[i+j*dim];
        }
      }
    }
    initialized_ = true;
    ObjectChanged();
  }

void Matrix::Set(int row, int col, double val)
{
	vector<int> indices(2,0);
	indices[0] = row;
	indices[1] = col;
	assert(row < Dim(0));
	assert(col < Dim(1));
	Tensor::Set(indices, val);
}

  void SymTMatrix::MultVectorImpl(Number alpha, const Vector &x, Number beta,
                                  Vector &y) const
  {
    //  A few sanity checks
    DBG_ASSERT(Dim()==x.Dim());
    DBG_ASSERT(Dim()==y.Dim());

    // Take care of the y part of the addition
    DBG_ASSERT(initialized_);
    if( beta!=0.0 ) {
      y.Scal(beta);
    }
    else {
      y.Set(0.0);  // In case y hasn't been initialized yet
    }

    // See if we can understand the data
    const DenseVector* dense_x = dynamic_cast<const DenseVector*>(&x);
    DBG_ASSERT(dense_x); /* ToDo: Implement others */
    DenseVector* dense_y = dynamic_cast<DenseVector*>(&y);
    DBG_ASSERT(dense_y); /* ToDo: Implement others */

    if (dense_x && dense_y) {
      const Index*  irn=Irows();
      const Index*  jcn=Jcols();
      const Number* val=values_;
      Number* yvals=dense_y->Values();

      if (dense_x->IsHomogeneous()) {
        Number as = alpha *  dense_x->Scalar();
        for(Index i=0; i<Nonzeros(); i++) {
          yvals[*irn-1] += as * (*val);
          if (*irn!=*jcn) {
            // this is not a diagonal element
            yvals[*jcn-1] += as * (*val);
          }
          val++;
          irn++;
          jcn++;
        }
      }
      else {
        const Number* xvals=dense_x->Values();
        for(Index i=0; i<Nonzeros(); i++) {
          yvals[*irn-1] += alpha* (*val) * xvals[*jcn-1];
          if (*irn!=*jcn) {
            // this is not a diagonal element
            yvals[*jcn-1] += alpha* (*val) * xvals[*irn-1];
          }
          val++;
          irn++;
          jcn++;
        }
      }
    }
  }

 inline
 Number* DenseVectorSpace::AllocateInternalStorage() const
 {
   if (Dim()>0) {
     return new Number[Dim()];
   }
   else {
     return NULL;
   }
 }

void ANCFBeamBE2D::GetdPosdqT(const Vector2D& p_loc, Matrix& dpdqi)
{
	//p = S(p.x,p.y,p.z)*q; d(p)/dq
	dpdqi.SetSize(SOS(),Dim());
	dpdqi.FillWithZeros();
	//d = S + ...
	for (int i = 1; i <= NS(); i++)
	{
		double s = GetS0(p_loc.X(), i);
		dpdqi((i-1)*Dim()+1,1) = s;
		dpdqi((i-1)*Dim()+2,2) = s;
	}
	if (p_loc.Y() != 0)
	{
		double y = p_loc.Y();
		Vector2D rx = GetPosx2D(p_loc.X());
		Vector2D n(-rx.X(), rx.Y());
		n /= rx.Norm();

		for (int i = 1; i <= NS(); i++)
		{
			double sx = GetS0x(p_loc.X(), i) * 2./GetLx();
			//y/|n|*dn/dq
			dpdqi((i-1)*Dim()+1,2) +=  y*sx;
			dpdqi((i-1)*Dim()+2,1) += -y*sx;

			//y*n/|n|*(r_x1*S_x1 + r_x2*S_x2)
			dpdqi((i-1)*Dim()+1,1) +=  y*n.X()*(rx.X()*sx);
			dpdqi((i-1)*Dim()+1,2) +=  y*n.Y()*(rx.X()*sx);
			dpdqi((i-1)*Dim()+2,1) +=  y*n.X()*(rx.Y()*sx);
			dpdqi((i-1)*Dim()+2,2) +=  y*n.Y()*(rx.Y()*sx);
		}
	}
};

void TestExprSplitOcc::test07() {
	const ExprSymbol& x=ExprSymbol::new_("x",Dim(1,3,1));

	// several occurrences of an index (but different nodes)
	const ExprNode& e1=x[0];
	const ExprNode& e2=x[1];
	Function f1(x,(x[0]+(-x)[1])+x[1]+e1+e2);
	ExprSplitOcc eso(f1.args(),f1.expr());

	const Array<const ExprSymbol>& x2=eso.get_x();
	const ExprNode& y2=eso.get_y();
	TEST_ASSERT(x2.size()==4);
	TEST_ASSERT(sameExpr(y2,"((((x[0]+(-x_0_)[1])+x[1])+x[0]_1_)+x[1]_1_)"));
	TEST_ASSERT(&eso.node(x2[0])==&x); // the "special node" (which is inserted first)
	TEST_ASSERT(&eso.node(x2[1])==&x);
	TEST_ASSERT(&eso.node(x2[2])==&e1);
	TEST_ASSERT(&eso.node(x2[3])==&e2);
	int* var;
	int n=eso.var_map(var);
	TEST_ASSERT(n==8);
	TEST_ASSERT(var[0]==0);
	TEST_ASSERT(var[1]==1);
	TEST_ASSERT(var[2]==2);
	TEST_ASSERT(var[3]==0);
	TEST_ASSERT(var[4]==1);
	TEST_ASSERT(var[5]==2);
	TEST_ASSERT(var[6]==0);
	TEST_ASSERT(var[7]==1);
}

static inline void SolveEqs(Cut *cut, count ncut,
  creal *delta, creal diff)
{
  real last = 0;
  real r = 1;
  Cut *c;

  for( c = cut; ; ++c ) {
    ccount dim = Dim(c->i);
    c->row = r -=
      Div(diff, (delta[Lower(dim)] + delta[Upper(dim)])*c->df);
    if( --ncut == 0 ) break;
    last += r*c->lhs;
  }

  last = Div(c->lhs - last, r);

  for( ; c >= cut; last += (--c)->lhs ) {
    creal delmin = -(c->delta = delta[c->i]);
    creal delmax = FRACT*(delmin + c->save);
    c->sol = Div(last, c->df);
    if( c->sol > delmax ) c->sol = .75*delmax;
    if( c->sol < delmin ) c->sol = .75*delmin;
  }
}

  CompoundVector::CompoundVector(const CompoundVectorSpace* owner_space, bool create_new)
      :
      Vector(owner_space),
      comps_(owner_space->NCompSpaces()),
      const_comps_(owner_space->NCompSpaces()),
      owner_space_(owner_space),
      vectors_valid_(false)
  {
    Index dim_check = 0;
    for (Index i=0; i<NComps(); i++) {
      SmartPtr<const VectorSpace> space = owner_space_->GetCompSpace(i);
      DBG_ASSERT(IsValid(space));
      dim_check += space->Dim();

      if (create_new) {
        comps_[i] = space->MakeNew();
      }
    }

    DBG_ASSERT(dim_check == Dim());

    if (create_new) {
      vectors_valid_ = VectorsValid();
    }
  }

void ANCFBeamBE2D::GetH(Matrix& H)
{
	if (Hmatrix.Getrows() == SOS())
	{
		H = Hmatrix;
		return;
	}
	else
	{
		double A = this->GetMaterial().BeamRhoA() / this->GetMaterial().Density();

		H.SetSize(SOS(), Dim());
		H.SetAll(0);

		for (IntegrationPointsIterator ip(integrationRuleLoad); !ip.IsEnd(); ++ip)
		{
			double x = ip.Point2D().X();

			// jacobi determinant
			//Vector2D rx0 = GetRefPosx2D(x);
			double det = 0.5*GetLx(); //*rx0.Norm();

			double d = A * det * ip.Weight();

			for (int i = 1; i <= NS(); i++)
			{
				double s = GetS0(x, i);
				H(2*i-1, 1) += d * s;
				H(2*i,   2) += d * s; //H(2*i-1, 1);
			}	
		}
		Hmatrix = H;
	}
};

bool SemiglobalLabMatcher::update()
{
    WriteGuard<ReadWritePipe<FloatImage, FloatImage> > wguard(m_wpipe);
    FloatImage leftImg, rightImg;   
    if ( m_lrpipe->read(&leftImg) && m_rrpipe->read(&rightImg) )
    {
        Dim dsiDim(leftImg.dim().width(), leftImg.dim().height(), 
                   m_maxDisparity);
        
        float *leftImg_d = leftImg.devMem();
        float *rightImg_d = rightImg.devMem();   
        FloatImage dispImage = FloatImage::CreateDev(
            Dim(dsiDim.width(), dsiDim.height()));
     
        cudaPitchedPtr aggregDSI = m_aggregDSI.mem(dsiDim);
        SGPath *paths = m_sgPaths.getDescDev(dispImage.dim());
                   
        
        SemiGlobalLabDevRun(dsiDim, paths, m_sgPaths.pathCount(),
                            leftImg_d, rightImg_d,
                            aggregDSI, dispImage.devMem(), m_zeroAggregDSI);
        m_zeroAggregDSI = false;
        
        dispImage.cpuMem();
        wguard.write(dispImage);        
    }
    
    return wguard.wasWrite();
}

  void DenseSymMatrix::PrintImpl(const Journalist& jnlst,
                                 EJournalLevel level,
                                 EJournalCategory category,
                                 const std::string& name,
                                 Index indent,
                                 const std::string& prefix) const
  {
    jnlst.Printf(level, category, "\n");
    jnlst.PrintfIndented(level, category, indent,
                         "%sDenseSymMatrix \"%s\" of dimension %d (only lower triangular part printed):\n",
                         prefix.c_str(), name.c_str(), Dim());

    if (initialized_) {
      for (Index j=0; j<NCols(); j++) {
        for (Index i=j; i<NRows(); i++) {
          jnlst.PrintfIndented(level, category, indent,
                               "%s%s[%5d,%5d]=%23.16e\n",
                               prefix.c_str(), name.c_str(), i, j, values_[i+NRows()*j]);
        }
      }
    }
    else {
      jnlst.PrintfIndented(level, category, indent,
                           "The matrix has not yet been initialized!\n");
    }
  }

  typename ImageFactory<T>::view_type* difference_of_exponential_crack_edge_image(const T& src, double scale, double gradient_threshold, unsigned int min_edge_length, unsigned int close_gaps, unsigned int beautify) {
    if ((scale < 0) || (gradient_threshold < 0))
      throw std::runtime_error("The scale and gradient threshold must be greater than 0");

    typename ImageFactory<T>::data_type* dest_data =
      new typename ImageFactory<T>::data_type(Dim(src.ncols() * 2, src.nrows() * 2), src.origin());

    typename ImageFactory<T>::view_type* dest =
      new typename ImageFactory<T>::view_type(*dest_data);

    try {
      vigra::differenceOfExponentialCrackEdgeImage(src_image_range(src), dest_image(*dest), scale, gradient_threshold, NumericTraits<typename T::value_type>::one());
    
      if (min_edge_length > 0)
        vigra::removeShortEdges(dest_image_range(*dest), min_edge_length, NumericTraits<typename T::value_type>::one());
    
      if (close_gaps)
        vigra::closeGapsInCrackEdgeImage(dest_image_range(*dest), NumericTraits<typename T::value_type>::one());
    
      if (beautify)
        vigra::beautifyCrackEdgeImage(dest_image_range(*dest), NumericTraits<typename T::value_type>::one(), NumericTraits<typename T::value_type>::zero());
    } catch (std::exception e) {
      delete dest;
      delete dest_data;
      throw;
    }
    return dest;
  }

void ANCFBeamBE2D::EvalM(Matrix& m, double t)
{
	if (this->massmatrix.Getrows() != 0)	// mass matrix already computed
	{
		m = massmatrix;
		return;
	}
	
	Matrix H(Dim(), SOS());
	massmatrix.SetSize(SOS(), SOS());
	massmatrix.SetAll(0.);

	double rhoA = GetBeamRhoA();
	
	for (IntegrationPointsIterator ip(integrationRuleMass); !ip.IsEnd(); ++ip)
	{
		for (int i = 1; i <= NS(); i++)
		{
			H(1, 2*i-1) = GetS0(ip.Point2D().X(), i);
			H(2, 2*i) = H(1, 2*i-1);
		}
		H = (0.5*GetLx()*rhoA*ip.Weight()) * (H.GetTp() * H);
		massmatrix += H;
	}
	m = massmatrix;
};

Dim Dim::index_dim() const {

  const Dim& dim=*this;

  if (dim.dim2==1 && dim.dim3==1) {
	  return Dim::scalar();
	  // we allow x[1] for x scalar
	  // old error message -> "Too many subscripts (e.g., a vector symbol cannot be indexed twice)";
  }
  if (dim.dim1>1) // array of matrices
	  return Dim(1,dim.dim2,dim.dim3);
  else
	  if (dim.dim2==1 || dim.dim3==1) // vector
		  return Dim(1,1,1);
	  else // matrix
		  return Dim(1,1,dim.dim3); // return a row vector
}

Image* scale(T& image, double scaling, int resize_quality) {
    // nrows, ncols are cast to a double so that the multiplication happens
    // exactly as it does in Python
    return resize(image,
                  Dim(size_t(double(image.ncols()) * scaling),
                      size_t(double(image.nrows()) * scaling)),
                  resize_quality);
}