C++ Epetra_Vector::Import方法代码示例

/*----------------------------------------------------------------------*
 |  Constructor (public)                                     m.gee 12/04|
 |  IMPORTANT:                                                          |
 |  No matter on which level we are here, the vector xfine is ALWAYS    |
 |  a fine grid vector here!                                            |
 *----------------------------------------------------------------------*/
ML_NOX::ML_Nox_MatrixfreeLevel::ML_Nox_MatrixfreeLevel(int level, int nlevel, int plevel, 
                                 ML* ml, ML_Aggregate* ag, Epetra_CrsMatrix** P,
                                 ML_NOX::Ml_Nox_Fineinterface& interface, const Epetra_Comm& comm,
                                 const Epetra_Vector& xfine, double fd_alpha, double fd_beta,
                                 bool fd_centered, bool isDiagonalOnly, int bsize) 
: fineinterface_(interface),
comm_(comm)                                                          
{
   level_          = level;
   nlevel_         = nlevel;
   ml_printlevel_  = plevel;
   ml_             = ml;
   ag_             = ag;
   fd_alpha_       = fd_alpha;
   fd_beta_        = fd_beta;
   fd_centered_    = fd_centered;
   isDiagonalOnly_ = isDiagonalOnly;
   A_              = 0;
   coarseinterface_= 0;
   bsize_          = bsize;
   

   // we need the graph of the operator on this level. On the fine grid we can just ask the
   // fineinterface for it, on the coarser levels it has to be extracted from the ML-hierachy
   if (level_==0)
   {
      // the Epetra_CrsGraph-copy-constructor shares data with the original one.
      // We want a really deep copy here so we cannot use it
      // graph_ will be given to the FiniteDifferencing class and will be destroyed by it
      graph_ = ML_NOX::deepcopy_graph(interface.getGraph());
   }
   else
   {
      // Note that ML has no understanding of global indices, so it makes up new GIDs
      // (This also holds for the Prolongators P)
      Epetra_CrsMatrix* tmpMat  = 0;
      int               maxnnz  = 0;
      double            cputime = 0.0;
      ML_Operator2EpetraCrsMatrix(&(ml_->Amat[level_]), tmpMat, maxnnz, false, cputime);
      // copy the graph
      double t0 = GetClock();
      graph_ = ML_NOX::deepcopy_graph(&(tmpMat->Graph()));
      // delete the copy of the Epetra_CrsMatrix
      if (tmpMat) delete tmpMat; tmpMat = 0;
      double t1 = GetClock();
      if (ml_printlevel_ > 0 && 0 == comm_.MyPID())
         cout << "matrixfreeML (level " << level_ << "): extraction/copy of Graph in " << cputime+t1-t0 << " sec\n"
              << "                        max-nonzeros in Graph: " << maxnnz << "\n";
   }
   
   // create this levels coarse interface
   coarseinterface_ = new ML_NOX::Nox_CoarseProblem_Interface(fineinterface_,level_,ml_printlevel_,
                                                              P,&(graph_->RowMap()),nlevel_);

   // restrict the xfine-vector to this level 
   Epetra_Vector* xthis = coarseinterface_->restrict_fine_to_this(xfine);
   if (!xthis)
   {
     cout << "**ERR**: ML_Epetra::ML_Nox_MatrixfreeLevel::ML_Nox_MatrixfreeLevel:\n"
          << "**ERR**: ML_Epetra::Nox_CoarseProblem_Interface::restrict_fine_to_this returned NULL on level " 
          << level_ << "\n"
          << "**ERR**: file/line: " << __FILE__ << "/" << __LINE__ << "\n"; throw -1;
   }

   Epetra_Vector* xc = new Epetra_Vector(graph_->RowMap(),false);
   // FIXME: after intesive testing, this test might be obsolet
#if 0
   bool samemap = xc->Map().PointSameAs(xthis->Map());
   if (samemap)
   {
#endif
      xc->Update(1.0,*xthis,0.0);
#if 0
   }
   else
   {
      cout << "**WRN** Maps are not equal in\n"
           << "**WRN** file/line: " << __FILE__ << "/" << __LINE__ << "\n";
      // import the xthis vector in the Map that ML produced for graph_
      Epetra_Import* importer = new Epetra_Import(graph_->RowMap(),xthis->Map());  
      int ierr = xc->Import(*xthis,*importer,Insert);
      if (ierr)
      {
         cout << "**ERR**: ML_Epetra::ML_Nox_MatrixfreeLevel::ML_Nox_MatrixfreeLevel:\n"
              << "**ERR**: export from xthis to xc returned err=" << ierr <<"\n"
              << "**ERR**: file/line: " << __FILE__ << "/" << __LINE__ << "\n"; throw -1;
      }
      if (importer) delete importer; importer = 0;
   }
#endif   
   if (xthis) delete xthis; xthis = 0;

   // create the coloring of the graph
   if (ml_printlevel_>0 && comm_.MyPID()==0) 
//.........这里部分代码省略.........

//.........这里部分代码省略.........
      coarseinterface_->recreate(ml_printlevel_,P,&(graph_->RowMap()));
   else
   {
      delete coarseinterface_;
      coarseinterface_ = new ML_NOX::Nox_CoarseProblem_Interface(fineinterface_,level_,ml_printlevel_,
                                                                 P,&(graph_->RowMap()),nlevel_);
   }
   
   // restrict the xfine-vector to this level 
   Epetra_Vector* xthis = coarseinterface_->restrict_fine_to_this(xfine);
   if (!xthis)
   {
     cout << "**ERR**: ML_Epetra::ML_Nox_MatrixfreeLevel::ML_Nox_MatrixfreeLevel:\n"
          << "**ERR**: ML_Epetra::Nox_CoarseProblem_Interface::restrict_fine_to_this returned NULL on level " 
          << level_ << "\n"
          << "**ERR**: file/line: " << __FILE__ << "/" << __LINE__ << "\n"; throw -1;
   }

   Epetra_Vector* xc = new Epetra_Vector(graph_->RowMap(),false);
   // FIXME: after intesive testing, this test might be obsolet
#if 0
   bool samemap = xc->Map().PointSameAs(xthis->Map());
   if (samemap)
   {
#endif
      xc->Update(1.0,*xthis,0.0);
#if 0
   }
   else
   {
      cout << "**WRN** Maps are not equal in\n"
           << "**WRN** file/line: " << __FILE__ << "/" << __LINE__ << "\n";
      // import the xthis vector in the Map that ML produced for graph_
      Epetra_Import* importer = new Epetra_Import(graph_->RowMap(),xthis->Map());  
      int ierr = xc->Import(*xthis,*importer,Insert);
      if (ierr)
      {
         cout << "**ERR**: ML_Epetra::ML_Nox_MatrixfreeLevel::ML_Nox_MatrixfreeLevel:\n"
              << "**ERR**: export from xthis to xc returned err=" << ierr <<"\n"
              << "**ERR**: file/line: " << __FILE__ << "/" << __LINE__ << "\n"; throw -1;
      }
      if (importer) delete importer; importer = 0;
   }
#endif   
   if (xthis) delete xthis; xthis = 0;

   // create the coloring of the graph
   if (ml_printlevel_>0 && comm_.MyPID()==0)
   {
      cout << "matrixfreeML (level " << level_ << "): Entering Recoloring on level " << level_ << "\n";
      fflush(stdout);
   }
   double t0 = GetClock();
   if (!same) // te graph has obviously changed, so we need to recolor
   {
      if (colorMap_) delete colorMap_; colorMap_ = 0;
      if (colorMapIndex_) delete colorMapIndex_; colorMapIndex_ = 0;
      if (colorcolumns_) delete colorcolumns_; colorcolumns_ = 0;

      colorMap_ = ML_NOX::ML_Nox_collapsedcoloring(graph_,bsize_,isDiagonalOnly_,ml_printlevel_);
      if (!colorMap_) colorMap_ = ML_NOX::ML_Nox_standardcoloring(graph_,isDiagonalOnly_);
      colorMapIndex_ = new EpetraExt::CrsGraph_MapColoringIndex(*colorMap_);
      colorcolumns_  = &(*colorMapIndex_)(*graph_);
    }
    else if (ml_printlevel_>0 && comm_.MyPID()==0)
      cout << "matrixfreeML (level " << level_ << "): Reusing existing Coloring on level " << level_ << "\n";

bool FiniteDifferenceColoringWithUpdate::differenceProbe(const Epetra_Vector& x, Epetra_CrsMatrix& jac,const Epetra_MapColoring& colors){

  // Allocate space for perturbation, get column version of x for scaling
  Epetra_Vector xp(x);
  Epetra_Vector *xcol;
  int N=jac.NumMyRows();

  if(jac.ColMap().SameAs(x.Map()))
     xcol=const_cast<Epetra_Vector*>(&x);
  else{
    xcol=new Epetra_Vector(jac.ColMap(),true);//zeros out by default
    xcol->Import(x,*jac.Importer(),InsertAdd);
  }

  // Counters for probing diagnostics
  double tmp,probing_error_lower_bound=0.0,jc_norm=0.0;

  // Grab coloring info (being very careful to ignore color 0)
  int Ncolors=colors.MaxNumColors()+1;
  int num_c0_global,num_c0_local=colors.NumElementsWithColor(0);
  colors.Comm().MaxAll(&num_c0_local,&num_c0_global,1);
  if(num_c0_global>0) Ncolors--;

  if(Ncolors==0) return false;

  // Pointers for Matrix Info
  int entries, *indices;
  double *values;

  // NTS: Fix me
  if ( diffType == Centered ) exit(1);

  double scaleFactor = 1.0;
  if ( diffType == Backward )
    scaleFactor = -1.0;

  // Compute RHS at initial solution
  computeF(x,fo,NOX::Epetra::Interface::Required::FD_Res);

  /* Probing, vector by vector since computeF does not have a MultiVector interface */
  // Assume that anything with Color 0 gets ignored.
  for(int j=1;j<Ncolors;j++){
    xp=x;
    for(int i=0;i<N;i++){
      if(colors[i]==j)
    xp[i] += scaleFactor*(alpha*abs(x[i])+beta);
    }

    computeF(xp, fp, NOX::Epetra::Interface::Required::FD_Res);

    // Do the subtraction to estimate the Jacobian (w/o including step length)
    Jc.Update(1.0, fp, -1.0, fo, 0.0);

    // Relative error in probing
     if(use_probing_diags){
       Jc.Norm2(&tmp);
       jc_norm+=tmp*tmp;
     }

    for(int i=0;i<N;i++){
      // Skip for uncolored row/columns, else update entries
      if(colors[i]==0) continue;

      jac.ExtractMyRowView(i,entries,values,indices);
      for(int k=0;k<jac.NumMyEntries(i);k++){
    if(colors[indices[k]]==j){
      values[k]=Jc[i] / (scaleFactor*(alpha*abs((*xcol)[indices[k]])+beta));
      // If probing diagnostics are on, zero out the entries as they are used
      if(use_probing_diags) Jc[i]=0.0;
      break;// Only one value per row...
    }
      }
    }
    if(use_probing_diags){
      Jc.Norm2(&tmp);
      probing_error_lower_bound+=tmp*tmp;
    }
  }

  // If diagnostics are requested, output Frobenius norm lower bound
  if(use_probing_diags && !x.Comm().MyPID()) printf("Probing Error Lower Bound (Frobenius) abs = %6.4e rel = %6.4e\n",sqrt(probing_error_lower_bound),sqrt(probing_error_lower_bound)/sqrt(jc_norm));

  // Cleanup
  if(!jac.ColMap().SameAs(x.Map()))
    delete xcol;

  return true;
}

// ====================================================================== 
int Amesos_Lapack::GEEV(Epetra_Vector& Er, Epetra_Vector& Ei)
{
  if (IsSymbolicFactorizationOK_ == false)
    AMESOS_CHK_ERR(SymbolicFactorization());

  if (MyPID_ == 0)
    AMESOS_CHK_ERR(DenseMatrix_.Shape(static_cast<int>(NumGlobalRows64()),static_cast<int>(NumGlobalRows64())));

  AMESOS_CHK_ERR(DistributedToSerial());
  AMESOS_CHK_ERR(SerialToDense());

  Teuchos::RCP<Epetra_Vector> LocalEr;
  Teuchos::RCP<Epetra_Vector> LocalEi;

  if (NumProcs_ == 1)
  {
    LocalEr = Teuchos::rcp(&Er, false);
    LocalEi = Teuchos::rcp(&Ei, false);
  }
  else
  {
    LocalEr = Teuchos::rcp(new Epetra_Vector(*SerialMap_));
    LocalEi = Teuchos::rcp(new Epetra_Vector(*SerialMap_));
  }

  if (MyPID_ == 0) 
  {
    int n = static_cast<int>(NumGlobalRows64());
    char jobvl = 'N'; /* V/N to calculate/not calculate left eigenvectors
                         of matrix H.*/
    char jobvr = 'N'; /* As above, but for right eigenvectors. */
    int info = 0;
    int ldvl = n;
    int ldvr = n;

    Er.PutScalar(0.0);
    Ei.PutScalar(0.0);

    Epetra_LAPACK LAPACK;

    std::vector<double> work(1);
    int lwork = -1;

    LAPACK.GEEV(jobvl, jobvr, n, DenseMatrix_.A(), n, 
                LocalEr->Values(), LocalEi->Values(), NULL,
                ldvl, NULL, 
                ldvr, &work[0], 
                lwork, &info);

    lwork = (int)work[0];
    work.resize(lwork);
    LAPACK.GEEV(jobvl, jobvr, n, DenseMatrix_.A(), n, 
                LocalEr->Values(), LocalEi->Values(), NULL,
                ldvl, NULL, 
                ldvr, &work[0], 
                lwork, &info);

    if (info)
      AMESOS_CHK_ERR(info);
  }

  if (NumProcs_ != 1)
  {
    // I am not really sure that exporting the results make sense... 
    // It is just to be coherent with the other parts of the code.
    Er.Import(*LocalEr, Epetra_Import(Er.Map(), SerialMap()), Insert);
    Ei.Import(*LocalEi, Epetra_Import(Ei.Map(), SerialMap()), Insert);
  }

  return(0);
}