C++ MultiVector类代码示例

// ====================================================================== 
void Eig(const Operator& Op, MultiVector& ER, MultiVector& EI)
{
  int ierr;
  if (Op.GetDomainSpace() != Op.GetRangeSpace())
    ML_THROW("Matrix is not square", -1);

  ER.Reshape(Op.GetDomainSpace());
  EI.Reshape(Op.GetDomainSpace());

  Epetra_LinearProblem Problem;
  Problem.SetOperator(const_cast<Epetra_RowMatrix*>(Op.GetRowMatrix()));
  Amesos_Lapack Lapack(Problem);

  Epetra_Vector ER_Epetra(Op.GetRowMatrix()->RowMatrixRowMap());
  Epetra_Vector EI_Epetra(Op.GetRowMatrix()->RowMatrixRowMap());

  ierr = Lapack.GEEV(ER_Epetra, EI_Epetra);

  if (ierr)
    ML_THROW("GEEV returned error code = " + GetString(ierr), -1);
  
  for (int i = 0 ; i < ER.GetMyLength() ; ++i) {
    ER(i) = ER_Epetra[i];
    EI(i) = EI_Epetra[i];
  }
}

//---------------------------------------------------------------------------
const std::vector<int> ManyDigitIndexTable::FindInternal(
  const std::vector<int>& indices,
  const MultiVector<int>& v,
  const int value)
{
  //Check indices
  {
    const auto begin = v.PeekIndices().begin();
    const auto end   = v.PeekIndices().end();
    const auto x = std::find(begin,end,value);
    if (x!=end)
    {
      const int index_found = std::distance(begin,x);
      assert(index_found >= 0);
      assert(index_found < boost::numeric_cast<int>(v.PeekIndices().size()));
      std::vector<int> result = indices;
      result.push_back(index_found);
      return result;
    }
  }

  //Check MultiVector
  const auto& mvs = v.PeekMultiVectors();
  const int size = mvs.size();
  for (int i=0; i!=size; ++i)
  {
    std::vector<int> indices_deeper = indices;
    indices_deeper.push_back(i);
    const auto result
      = FindInternal(
        indices_deeper,mvs[i],value);
    if (!result.empty()) return result;
  }
  return std::vector<int>();
}

/*----------------------------------------------------------------------*
 |  apply multigrid linear preconditioner (private)          m.gee 03/06|
 *----------------------------------------------------------------------*/
int MOERTEL::Mortar_ML_Preconditioner::MultiLevelSA(
    const MultiVector& b1_f,
    const MultiVector& b2_f,
    MultiVector& x1_f,
    MultiVector& x2_f,
    int level) const
{
    MultiVector r1_f(b1_f.GetVectorSpace(),1,false);
    MultiVector z1_f(b1_f.GetVectorSpace(),1,false);

    // presmoothing
    x1_f = 0;
    G_.Apply(b1_f,x1_f);
    x2_f = mlapiMT_ * x1_f;
    x2_f.Scale(-1.0);

    // compute residual
    r1_f = b1_f - mlapiAhat11_ * x1_f;

    // postsmoothing
    z1_f = 0;
    G_.Apply(r1_f,z1_f);
    x1_f = x1_f + z1_f;
    x2_f = mlapiMT_ * x1_f;
    x2_f.Scale(-1.0);



    return 0;
}

  void IfpackSmoother::Apply(MultiVector& X, const MultiVector& B, bool InitialGuessIsZero) const {
    TEUCHOS_TEST_FOR_EXCEPTION(SmootherPrototype::IsSetup() == false, Exceptions::RuntimeError, "MueLu::IfpackSmoother::Apply(): Setup() has not been called");

    // Forward the InitialGuessIsZero option to Ifpack
    Teuchos::ParameterList  paramList;
    if (type_ == "Chebyshev")
      paramList.set("chebyshev: zero starting solution",  InitialGuessIsZero);

    else if (type_ == "point relaxation stand-alone")
      paramList.set("relaxation: zero starting solution", InitialGuessIsZero);

    SetPrecParameters(paramList);

    // Apply
    if (InitialGuessIsZero) {
      Epetra_MultiVector&       epX = Utils::MV2NonConstEpetraMV(X);
      const Epetra_MultiVector& epB = Utils::MV2EpetraMV(B);
      prec_->ApplyInverse(epB, epX);

    } else {
      RCP<MultiVector> Residual = Utils::Residual(*A_,X,B);
      RCP<MultiVector> Correction = MultiVectorFactory::Build(A_->getDomainMap(), X.getNumVectors());
      Epetra_MultiVector&       epX = Utils::MV2NonConstEpetraMV(*Correction);
      const Epetra_MultiVector& epB = Utils::MV2EpetraMV(*Residual);
      prec_->ApplyInverse(epB, epX);
      X.update(1.0, *Correction, 1.0);
    }
  }

  void ZoltanInterface<LocalOrdinal, GlobalOrdinal, Node, LocalMatOps>::
  GetProblemGeometry(void *data, int numGIDEntries, int numLIDEntries, int numObjectIDs,
                     ZOLTAN_ID_PTR gids, ZOLTAN_ID_PTR lids, int dim, double *coordinates, int *ierr)
  {
    if (data == NULL) {
      *ierr = ZOLTAN_FATAL;
      return;
    }

    MultiVector *Coords = (MultiVector*) data;

    if (dim != Teuchos::as<int>(Coords->getNumVectors())) {
      //FIXME I'm assuming dim should be 1, 2, or 3 coming in?!
      *ierr = ZOLTAN_FATAL;
      return;
    }

    TEUCHOS_TEST_FOR_EXCEPTION(numObjectIDs != Teuchos::as<int>(Coords->getLocalLength()), Exceptions::Incompatible, "Length of coordinates must be the same as the number of objects");

    ArrayRCP<ArrayRCP<const SC> > CoordsData(dim);
    for (int j = 0; j < dim; ++j)
      CoordsData[j] = Coords->getData(j);

    size_t numElements = Coords->getLocalLength();
    for (size_t i = 0; i < numElements; ++i)
      for (int j = 0; j < dim; ++j)
        coordinates[i*dim+j] = (double) CoordsData[j][i];

    *ierr = ZOLTAN_OK;

  } //GetProblemGeometry

void test_rotate()
{
    QVector<Vector> axis(MultiVector::count()), vec(MultiVector::count());
    QVector<double> angle(MultiVector::count());
    QVector<Vector> result(MultiVector::count());

    for (int i=0; i<MultiVector::count(); ++i)
    {
        axis[i] = rangen.vectorOnSphere(1);
        vec[i] = rangen.vectorOnSphere(5);
        angle[i] = rangen.rand(-2.0*SireMaths::pi, 2.0*SireMaths::pi);
        result[i] = Quaternion(Angle(angle[i]),axis[i]).rotate(vec[i]);
    }

    MultiVector maxis(axis);
    MultiVector mvec(vec);
    MultiDouble mangle(angle);

    MultiVector mresult = MultiQuaternion(mangle, maxis).rotate(mvec);

    for (int i=0; i<MultiVector::count(); ++i)
    {
        for (int j=0; j<3; ++j)
        {
            assert_nearly_equal( mresult.at(i)[j], result[i][j], 1e-5 );
        }
    }
}

// ====================================================================== 
// FIXME: Add List
void Eigs(const Operator& A, int NumEigenvalues, 
          MultiVector& ER, MultiVector& EI)
{

  if (A.GetDomainSpace() != A.GetRangeSpace())
    ML_THROW("Input Operator is not square", -1);

  double time;

  time = GetClock();

  int length = NumEigenvalues;
  double tol = 1e-3;
  int restarts = 1;
  int output = 10;
  bool PrintStatus = true;

  // 1.- set parameters for Anasazi
  Teuchos::ParameterList AnasaziList;
  // MatVec should be either "A" or "ML^{-1}A"
  AnasaziList.set("eigen-analysis: matrix operation", "A");
  AnasaziList.set("eigen-analysis: use diagonal scaling", false);
  AnasaziList.set("eigen-analysis: symmetric problem", false);
  AnasaziList.set("eigen-analysis: length", length);
  AnasaziList.set("eigen-analysis: block-size",1);
  AnasaziList.set("eigen-analysis: tolerance", tol);
  AnasaziList.set("eigen-analysis: restart", restarts);
  AnasaziList.set("eigen-analysis: output", output);
  AnasaziList.get("eigen-analysis: print current status",PrintStatus);

  // data to hold real and imag for eigenvalues and eigenvectors
  Space ESpace(-1, NumEigenvalues);
  ER.Reshape(ESpace);
  EI.Reshape(ESpace);

  // this is the starting value -- random
  Epetra_MultiVector EigenVectors(A.GetRowMatrix()->OperatorDomainMap(),
                                  NumEigenvalues);
  EigenVectors.Random();
#ifdef HAVE_ML_ANASAxI
  //int NumRealEigenvectors, NumImagEigenvectors;
#endif

  AnasaziList.set("eigen-analysis: action", "LM");

#ifdef HAVE_ML_ANASAxI
  ML_THROW("fixme...", -1);
  /* FIXME
  ML_Anasazi::Interface(A.GetRowMatrix(),EigenVectors,ER.GetValues(),
			EI.GetValues(), AnasaziList, 0, 0,
			&NumRealEigenvectors, &NumImagEigenvectors, 0);
                        */
#else
  ML_THROW("Anasazi is no longer supported", -1);
#endif

  return;
}

  void IndefBlockedDiagonalSmoother<Scalar, LocalOrdinal, GlobalOrdinal, Node, LocalMatOps>::Apply(MultiVector& X, const MultiVector& B, bool InitialGuessIsZero) const
  {
    TEUCHOS_TEST_FOR_EXCEPTION(SmootherPrototype::IsSetup() == false, Exceptions::RuntimeError, "MueLu::IndefBlockedDiagonalSmoother::Apply(): Setup() has not been called");

    Teuchos::RCP<Teuchos::FancyOStream> fos = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));

    SC zero = Teuchos::ScalarTraits<SC>::zero(), one = Teuchos::ScalarTraits<SC>::one();

    // extract parameters from internal parameter list
    const ParameterList & pL = Factory::GetParameterList();
    LocalOrdinal nSweeps = pL.get<LocalOrdinal>("Sweeps");
    Scalar omega = pL.get<Scalar>("Damping factor");

    // wrap current solution vector in RCP
    RCP<MultiVector> rcpX = Teuchos::rcpFromRef(X);

    // create residual vector
    // contains current residual of current solution X with rhs B
    RCP<MultiVector> residual = MultiVectorFactory::Build(B.getMap(), B.getNumVectors());

    // incrementally improve solution vector X
    for (LocalOrdinal run = 0; run < nSweeps; ++run) {
      // 1) calculate current residual
      residual->update(one,B,zero); // residual = B
      A_->apply(*rcpX, *residual, Teuchos::NO_TRANS, -one, one);

      // split residual vector
      Teuchos::RCP<MultiVector> r1 = rangeMapExtractor_->ExtractVector(residual, 0);
      Teuchos::RCP<MultiVector> r2 = rangeMapExtractor_->ExtractVector(residual, 1);

      // 2) solve F * \Delta \tilde{x}_1 = r_1
      //    start with zero guess \Delta \tilde{x}_1
      RCP<MultiVector> xtilde1 = MultiVectorFactory::Build(F_->getRowMap(),1);
      xtilde1->putScalar(zero);
      velPredictSmoo_->Apply(*xtilde1,*r1);

      // 3) solve SchurComp equation
      //    start with zero guess \Delta \tilde{x}_2
      RCP<MultiVector> xtilde2 = MultiVectorFactory::Build(Z_->getRowMap(),1);
      xtilde2->putScalar(zero);
      schurCompSmoo_->Apply(*xtilde2,*r2);

      // 4) extract parts of solution vector X
      Teuchos::RCP<MultiVector> x1 = domainMapExtractor_->ExtractVector(rcpX, 0);
      Teuchos::RCP<MultiVector> x2 = domainMapExtractor_->ExtractVector(rcpX, 1);

      // 5) update solution vector with increments xhat1 and xhat2
      //    rescale increment for x2 with omega_
      x1->update(omega,*xtilde1,one); // x1 = x1_old + omega xtilde1
      x2->update(omega,*xtilde2,one); // x2 = x2_old + omega xtilde2

      // write back solution in global vector X
      domainMapExtractor_->InsertVector(x1, 0, rcpX);
      domainMapExtractor_->InsertVector(x2, 1, rcpX);

    }

  }

// ======================================================================
MultiVector Extract(const MultiVector& y, const int v)
{
  if ((v < 0) || v >= y.GetNumVectors())
    ML_THROW("Wrong input parameter v (" +
             GetString(v) + ")", -1);

  MultiVector x(y.GetVectorSpace(), y.GetRCPValues(v));

  return(x);
}

/** Use this quaternion to rotate 'p' */
MultiVector MultiQuaternion::rotate(const MultiVector &p) const
{
    const MultiDouble sx2 = sc[0]*sc[0];
    const MultiDouble sy2 = sc[1]*sc[1];
    const MultiDouble sz2 = sc[2]*sc[2];

    const MultiDouble sxy = sc[0]*sc[1];
    const MultiDouble sxz = sc[0]*sc[2];
    const MultiDouble syz = sc[1]*sc[2];

    const MultiDouble swx = sc[0]*sc[3];
    const MultiDouble swy = sc[1]*sc[3];
    const MultiDouble swz = sc[2]*sc[3];

    const MultiDouble two(2.0);
    const MultiDouble half(0.5);

    return MultiVector( two*( ( half - sy2 - sz2 )*p.x()
                           + ( sxy - swz )       *p.y()
                           + ( sxz + swy )       *p.z()),

                        two*( ( sxy + swz )      *p.x()
                           + ( half - sx2 - sz2 ) *p.y()
                           + ( syz - swx )       *p.z()),

                        two*( ( sxz - swy )      *p.x()
                           + ( syz + swx )       *p.y()
                           + ( half - sx2 - sy2 ) *p.z()) );
}

// ======================================================================
MultiVector Duplicate(const MultiVector& y, const int v)
{
  if ((v < 0) || v >= y.GetNumVectors())
    ML_THROW("Wrong input parameter v (" +
             GetString(v) + ")", -1);

  // FIXME: use Extract
  MultiVector x(y.GetVectorSpace(), 1);
  for (int i = 0 ; i < x.GetMyLength() ; ++i)
    x(i) = y(i,v);

  return(x);
}

void LinearCheckNodeUpdater::updateNode(
										unsigned int nodeInd,
										MultiVector<float>& messages)

{
	uint32_t begin = messages.begin(nodeInd);
	uint32_t end = messages.end(nodeInd);
	uint32_t N = end - begin;

	if(N == 1) {
		messages[begin] = m_checkNodesPriorLogPmQ[nodeInd].getFloat();
		return;
	}

	// Convert all message values into log(p-q) values
	for (unsigned int j = begin; j < end; j++) {
		m_messagesPmQ.push_back(LogPmQ(messages[j]));
	}

	// calculate the multiplication from the right of p-q values:
	// rightSum[j] = encodedSoftBit[i] * messages[-1].p-q * messages[-2].p-q * ... * messages[-j].p-q
	LogPmQ curSum = m_checkNodesPriorLogPmQ[nodeInd]; // get log(p-q)

	m_rightSum.push_back(curSum);
	// go over all messages except the first one (which will never be
	// needed anyway)
	for(unsigned int rightInd = N - 1; rightInd > 0; rightInd--) {
		curSum += m_messagesPmQ[rightInd];

		m_rightSum.push_back(curSum);
	}


	// special case for j = 0
	messages[begin] = m_rightSum[N - 1].getFloat();

	// now we use curSum as the cumulative multiplication from the
	// left rather than right
	curSum = m_messagesPmQ[0];

	for (unsigned int j = 1; j < N; j++) {
		messages[begin + j] = (curSum + m_rightSum[N - j - 1]).getFloat();

		//assert(isfinite(messages[begin + j]));

		curSum += m_messagesPmQ[j];
	}

	m_rightSum.clear();
	m_messagesPmQ.clear();
}

      void
      apply (const MultiVector<S,LO,GO,Node> &X, 
	     MultiVector<S,LO,GO,Node> &Y,
	     Teuchos::ETransp mode = Teuchos::NO_TRANS, 
	     S alpha = Teuchos::ScalarTraits<S>::one (), 
	     S beta = Teuchos::ScalarTraits<S>::zero ()) const
      {
	const size_t numVectors = X.getNumVectors ();
	RCP<MultiVector<S,LO,GO,Node> > mvec_inout;
	RCP<const MultiVector<S,LO,GO,Node> > mvec_in2;

	if (_importer != null) {
	  if (_importMV != null && _importMV->getNumVectors () != numVectors) {
	    _importMV = null;
	  }
	  if (_importMV == null) {
	    _importMV = createMultiVector<S> (_importer->getTargetMap (), numVectors);
	  }
	  _importMV->doImport (X, *_importer, INSERT);
	  mvec_in2 = _importMV;
	}
	else {
	  mvec_in2 = rcpFromRef(X);
	}

	if (_exporter != null) {
	  if (_exportMV != null && _exportMV->getNumVectors () != numVectors) {
	    _exportMV = null;
	  }
	  if (_exportMV == null) {
	    _exportMV = createMultiVector<S> (_exporter->getSourceMap (), numVectors);
	  }
	  mvec_inout = _exportMV;
	}
	else {
	  mvec_inout = rcpFromRef (Y);
	}
	_kernel.setAlphaBeta (alpha, beta);
	//
	for (size_t j=0; j < numVectors; ++j) {
	  RCP<       Vector<S,LO,GO,Node> > vec_inout = mvec_inout->getVectorNonConst(j);
	  RCP< const Vector<S,LO,GO,Node> > vec_in2   = mvec_in2->getVector(j);
	  Tpetra::RTI::detail::binary_transform( *vec_inout, *vec_in2, _kernel );
	}
	// export
	if (_exporter != null) {
	  Y.doExport (*_exportMV, *_exporter, ADD);
	}
      }

/** \brief Copy the contents of a multivector to a destination vector.
  *
  * Copy the contents of a multivector to a new vector. If the destination
  * vector is null, a deep copy of the source multivector is made to a newly allocated
  * vector. Also, if the destination and the source do not match, a new destination
  * object is allocated and returned to the user.
  *
  * \param[in] src Source multivector to be copied.
  * \param[in] dst Destination multivector.  If null a new multivector will be allocated.
  *
  * \returns A copy of the source multivector. If dst is not null a pointer to this object
  *          is returned. Otherwise a new multivector is returned.
  */
inline MultiVector datacopy(const MultiVector & src,MultiVector & dst)
{ 
   if(dst==Teuchos::null)
      return deepcopy(src);

   bool rangeCompat = src->range()->isCompatible(*dst->range());
   bool domainCompat = src->domain()->isCompatible(*dst->domain());
 
   if(not (rangeCompat && domainCompat))
      return deepcopy(src);

   // perform data copy
   Thyra::assign<double>(dst.ptr(),*src);
   return dst;
}

// ====================================================================== 
Operator GetPtent1D(const MultiVector& D, const int offset = 0)
{
  if (D.GetNumVectors() != 1)
    ML_THROW("D.GetNumVectors() != 1", -1);

  int size = D.GetMyLength();
  if (size == 0)
    ML_THROW("empty diagonal vector in input", -1);

  double* diag = new double[size];
  for (int i = 0 ; i < size ; ++i)
    diag[i] = D(i);

  // creates the ML operator and store the diag pointer,
  // as well as the function pointers
  ML_Operator* MLDiag = ML_Operator_Create(GetML_Comm());

  int invec_leng = size / 3 + size % 3;
  int outvec_leng = size;

  MLDiag->invec_leng = invec_leng;
  MLDiag->outvec_leng = outvec_leng;
  MLDiag->data = (void*)diag;
  MLDiag->data_destroy = Ptent1D_destroy;
  MLDiag->matvec->func_ptr = Ptent1D_matvec;

  MLDiag->matvec->ML_id = ML_NONEMPTY;
  MLDiag->matvec->Nrows = outvec_leng;
  MLDiag->from_an_ml_operator = 0;

  MLDiag->getrow->func_ptr = Ptent1D_getrows;

  MLDiag->getrow->ML_id = ML_NONEMPTY;
  MLDiag->getrow->Nrows = outvec_leng;

  // creates the domain space
  vector<int> MyGlobalElements(invec_leng);
  for (int i = 0 ; i < invec_leng ; ++i) 
    MyGlobalElements[i] = D.GetVectorSpace()(i * 3) / 3;
  Space DomainSpace(invec_leng, -1, &MyGlobalElements[0]);
  Space RangeSpace = D.GetVectorSpace();

  // creates the MLAPI wrapper
  Operator Diag(DomainSpace,RangeSpace,MLDiag,true);
  return(Diag);
}