C++ MultiVector::Reshape方法代码示例

// ====================================================================== 
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];
  }
}

// ====================================================================== 
// 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 GetPtent(const Operator& A, Teuchos::ParameterList& List,
              const MultiVector& ThisNS,
              Operator& Ptent, MultiVector& NextNS)
{
  std::string CoarsenType     = List.get("aggregation: type", "Uncoupled");
  /* old version
  int    NodesPerAggr    = List.get("aggregation: per aggregate", 64);
  */
  double Threshold       = List.get("aggregation: threshold", 0.0);
  int    NumPDEEquations = List.get("PDE equations", 1);

  ML_Aggregate* agg_object;
  ML_Aggregate_Create(&agg_object);
  ML_Aggregate_Set_MaxLevels(agg_object,2);
  ML_Aggregate_Set_StartLevel(agg_object,0);
  ML_Aggregate_Set_Threshold(agg_object,Threshold);
  //agg_object->curr_threshold = 0.0;

  ML_Operator* ML_Ptent = 0;
  ML_Ptent = ML_Operator_Create(GetML_Comm());

  if (ThisNS.GetNumVectors() == 0)
    ML_THROW("zero-dimension null space", -1);

  int size = ThisNS.GetMyLength();

  double* null_vect = 0;
  ML_memory_alloc((void **)(&null_vect), sizeof(double) * size * ThisNS.GetNumVectors(), "ns");

  int incr = 1;
  for (int v = 0 ; v < ThisNS.GetNumVectors() ; ++v)
    DCOPY_F77(&size, (double*)ThisNS.GetValues(v), &incr,
              null_vect + v * ThisNS.GetMyLength(), &incr);


  ML_Aggregate_Set_NullSpace(agg_object, NumPDEEquations,
                             ThisNS.GetNumVectors(), null_vect,
                             ThisNS.GetMyLength());

  if (CoarsenType == "Uncoupled")
    agg_object->coarsen_scheme = ML_AGGR_UNCOUPLED;
  else if (CoarsenType == "Uncoupled-MIS")
    agg_object->coarsen_scheme = ML_AGGR_HYBRIDUM;
  else if (CoarsenType == "MIS") {
   /* needed for MIS, otherwise it sets the number of equations to
    * the null space dimension */
    agg_object->max_levels  = -7;
    agg_object->coarsen_scheme = ML_AGGR_MIS;
  }
  else if (CoarsenType == "METIS")
    agg_object->coarsen_scheme = ML_AGGR_METIS;
  else {
    ML_THROW("Requested aggregation scheme (" + CoarsenType +
             ") not recognized", -1);
  }

  int NextSize = ML_Aggregate_Coarsen(agg_object, A.GetML_Operator(),
                                      &ML_Ptent, GetML_Comm());

  /* This is the old version
  int NextSize;

  if (CoarsenType == "Uncoupled") {
    NextSize = ML_Aggregate_CoarsenUncoupled(agg_object, A.GetML_Operator(),
  }
  else if (CoarsenType == "MIS") {
    NextSize = ML_Aggregate_CoarsenMIS(agg_object, A.GetML_Operator(),
                                       &ML_Ptent, GetML_Comm());
  }
  else if (CoarsenType == "METIS") {
    ML ml_object;
    ml_object.ML_num_levels = 1; // crap for line below
    ML_Aggregate_Set_NodesPerAggr(&ml_object,agg_object,0,NodesPerAggr);
    NextSize = ML_Aggregate_CoarsenMETIS(agg_object, A.GetML_Operator(),
                                         &ML_Ptent, GetML_Comm());
  }
  else {
    ML_THROW("Requested aggregation scheme (" + CoarsenType +
             ") not recognized", -1);
  }
  */

  ML_Operator_ChangeToSinglePrecision(ML_Ptent);

  int NumMyElements = NextSize;
  Space CoarseSpace(-1,NumMyElements);
  Ptent.Reshape(CoarseSpace,A.GetRangeSpace(),ML_Ptent,true);

  assert (NextSize * ThisNS.GetNumVectors() != 0);

  NextNS.Reshape(CoarseSpace, ThisNS.GetNumVectors());

  size = NextNS.GetMyLength();
  for (int v = 0 ; v < NextNS.GetNumVectors() ; ++v)
    DCOPY_F77(&size, agg_object->nullspace_vect + v * size, &incr,
              NextNS.GetValues(v), &incr);

  ML_Aggregate_Destroy(&agg_object);
  ML_memory_free((void**)(&null_vect));
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    Operator mlapiAtildesplit(space,space,Atildesplit_.get(),false);
    Operator mlapiAhat(space,space,Ahat_.get(),false);
    Operator mlapiBWT(space,space,BWT.get(),false);
    Operator mlapiBWTcoarse;
    Operator ImBWTfine = GetIdentity(space,space) - mlapiBWT;
    Operator ImBWTcoarse;
    Operator Ptent;
    Operator P;
    Operator Pmod;
    Operator Rtent;
    Operator R;
    Operator Rmod;
    Operator IminusA;
    Operator C;
    InverseOperator S;

    mlapiAtildesplit_.resize(maxlevels);
    mlapiAhat_.resize(maxlevels);
    mlapiImBWT_.resize(maxlevels);
    mlapiImWBT_.resize(maxlevels);
    mlapiRmod_.resize(maxlevels);
    mlapiPmod_.resize(maxlevels);
    mlapiS_.resize(maxlevels);

    mlapiAtildesplit_[0] = mlapiAtildesplit;
    mlapiAhat_[0]        = mlapiAhat;
    mlapiImBWT_[0]       = ImBWTfine;
    mlapiImWBT_[0]       = GetTranspose(ImBWTfine);


    // build nullspace
    MultiVector NS;
    MultiVector NextNS;
    NS.Reshape(mlapiAsplit.GetRangeSpace(),nsdim);
    if (nullspace)
    {
        for (int i=0; i<nsdim; ++i)
            for (int j=0; j<NS.GetMyLength(); ++j)
                NS(j,i) = nullspace[i*NS.GetMyLength()+j];
    }
    else
    {
        if (numpde==1) NS = 1.0;
        else
        {
            NS = 0.0;
            for (int i=0; i<NS.GetMyLength(); ++i)
                for (int j=0; j<numpde; ++j)
                    if ( i % numpde == j)
                        NS(i,j) = 1.0;
        }
    }

    double lambdamax;

    // construct the hierarchy
    int level=0;
    for (level=0; level<maxlevels-1; ++level)
    {
        // this level's smoothing operator
        mlapiAtildesplit = mlapiAtildesplit_[level];

        // build smoother
        if (Comm().MyPID()==0)
        {
            ML_print_line("-", 78);