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

//.........这里部分代码省略.........
  int nnz2=Askew2->NumMyNonzeros();
  int N=Askew2->NumMyRows();


  // Grab pointers
  IFPACK_CHK_ERR(Askew2->ExtractCrsDataPointers(rowptr_s,colind_s,vals_s));
  IFPACK_CHK_ERR(Aherm2->ExtractCrsDataPointers(rowptr_h,colind_h,vals_h));

  // Sanity checking: Make sure the sparsity patterns are the same
#define SANITY_CHECK
#ifdef SANITY_CHECK
  for(int i=0;i<N;i++)
    if(rowptr_s[i]!=rowptr_h[i]) IFPACK_CHK_ERR(-2);
  for(int i=0;i<nnz2;i++)
    if(colind_s[i]!=colind_h[i]) IFPACK_CHK_ERR(-3);
#endif

  // Dirichlet Detection & Nuking of Aherm2 and Askew2
  // Note: Relies on Aherm2/Askew2 having identical sparsity patterns (see sanity check above)
  if(HaveOAZBoundaries_){
    int numBCRows;
    int* dirRows=FindLocalDiricheltRowsFromOnesAndZeros(*Acrs_,numBCRows);
    Epetra_IntVector* dirCols=FindLocalDirichletColumnsFromRows(dirRows,numBCRows,*Aherm2);
    Apply_BCsToMatrixRowsAndColumns(dirRows,numBCRows,*dirCols,*Aherm2);
    Apply_BCsToMatrixRowsAndColumns(dirRows,numBCRows,*dirCols,*Askew2);
    delete [] dirRows;
    delete dirCols;
  }

  // Grab diagonal of A
  A_->ExtractDiagonalCopy(Adiag);

  // Allocate the diagonal for W
  Epetra_Vector *Wdiag = new Epetra_Vector(*RowMap);

  // Build the W matrix (lower triangle only)
  // Note: Relies on EpetraExt giving me identical sparsity patterns for both Askew2 and Aherm2 (see sanity check above)
  int maxentries=Askew2->MaxNumEntries();
  int_type* gids=new int_type [maxentries];
  double* newvals=new double[maxentries];
  W=new Epetra_CrsMatrix(Copy,*RowMap,0);
  for(int i=0;i<N;i++){
    // Build the - (1+alpha)/2 E - (1-alpha)/2 F part of the W matrix
    int_type rowgid = (int_type) Acrs_->GRID64(i);
    double c_data=0.0;
    double ipdamp=0.0;
    int idx=0;

    for(int j=rowptr_s[i];j<rowptr_s[i+1];j++){
      int_type colgid = (int_type) Askew2->GCID64(colind_s[j]);
      c_data+=fabs(vals_s[j]);
      if(rowgid>colgid){
        // Rely on the fact that off-diagonal entries are always numbered last, dropping the entry entirely.
        if(colind_s[j] < N) {
          gids[idx]=colgid;
          newvals[idx]=vals_h[j]/2 + Alpha_ * vals_s[j]/2;
          idx++;
        }
        else{
          ipdamp+=fabs(vals_h[j]/2 + Alpha_ * vals_s[j]/2);
        }
      }
    }
    if(idx>0)
      IFPACK_CHK_ERR(W->InsertGlobalValues(rowgid,idx,newvals,gids));

    // Do the diagonal
    double w_val= c_data*Alpha_*Gamma_/4 + Adiag[Acrs_->LRID(rowgid)];
    if(UseInterprocDamping_) w_val+=ipdamp;

    W->InsertGlobalValues(rowgid,1,&w_val,&rowgid);
    IFPACK_CHK_ERR(Wdiag->ReplaceGlobalValues(1,&w_val,&rowgid));
  }
  W->FillComplete(A_->OperatorDomainMap(),A_->OperatorRangeMap());
  W_= Teuchos::rcp(W);
  Wdiag_= Teuchos::rcp(Wdiag);

  // Mark as computed
  IsComputed_=true;

  // Global damping, if wanted
  if(UseGlobalDamping_) {
    PowerMethod(10,LambdaMax_);
    if(!A_->Comm().MyPID()) printf("SORa: Global damping parameter = %6.4e (lmax=%6.4e)\n",GetOmega(),LambdaMax_);
  }

  // Cleanup
  delete [] gids;
  delete [] newvals;
  delete Aherm2;
  delete Askew2;
  if(RowMatrixMode) {
    Acrs_=Teuchos::null;
  }

  // Counters
  NumCompute_++;
  ComputeTime_ += Time_.ElapsedTime();
  return 0;
}

//.........这里部分代码省略.........
		    <<  X[testVecIndex][FirstEntryOfLID]
		    << " should = " << newLIDValue << endl;
  newLIDValue = 2.0;
  locerr = X.ReplaceMyValue(testLID, LIDSize-1, testVecIndex, newLIDValue);
  locerr = X.SumIntoMyValue(testLID, LIDSize-1, testVecIndex, newLIDValue);
  if (verbose) cout << "X["<<testVecIndex<<"]["<<FirstEntryOfLID+LIDSize-1<<"] = "
		    <<  X[testVecIndex][FirstEntryOfLID+LIDSize-1]
		    << " should = " << newLIDValue << endl;
  if (X[testVecIndex][FirstEntryOfLID+LIDSize-1]!=(newLIDValue+newLIDValue)) err++;

  ierr += err;

  // ========================================================================
  // Test Post-construction modification of an Epetra_Vector using a vector
  // our multivector X
  // ========================================================================

  if (verbose) cout <<  "\n\nXXXXX Testing Post-construction modification of a vector"
		    << endl << endl;

  Epetra_Vector * x = X(testVecIndex);

  int NumEntries = 2;
  double * VecValues = new double[NumEntries];
  int * VecGIDs = new int[NumEntries];
  VecGIDs[0] = testGID;
  VecGIDs[1] = testGID+1; // Some pathological chance that these GIDs are not valid

  // ========================================================================
  // Test int ReplaceGlobalValues (int NumEntries, double *Values, int *Indices)
  // ========================================================================

  VecValues[0] = 2.0; VecValues[1] = 4.0;
  locerr = x->ReplaceGlobalValues(NumEntries, VecValues, VecGIDs);

  for (i=0; i<NumEntries; i++) {
    testGID = VecGIDs[i];
    if (Map.MyGID(testGID)) {
      LIDOfGID = Map.LID(testGID);
      GIDSize = EPETRA_MIN(GIDSize,Map.ElementSize(LIDOfGID)); // Need this value below
      FirstEntryOfGID = Map.FirstPointInElement(LIDOfGID);
      if ((*x)[FirstEntryOfGID]!=VecValues[i]) err++;
      if (verbose) cout << "x["<<FirstEntryOfGID<<"] = "
			<< (*x)[FirstEntryOfGID]
			<< " should = " << VecValues[i] << endl;
    }
    else
      if (locerr!=1) err++; // Test for GID out of range error (=1)
  }


  // ========================================================================
  // Test int ReplaceGlobalValues (int NumEntries, int BlockOffset, double *Values, int *Indices)
  // ========================================================================

  VecValues[0] = 4.0; VecValues[1] = 8.0;
  locerr = x->ReplaceGlobalValues(NumEntries, GIDSize-1, VecValues, VecGIDs);

  for (i=0; i<NumEntries; i++) {
    testGID = VecGIDs[i];
    if (Map.MyGID(testGID)) {
      LIDOfGID = Map.LID(testGID);
      FirstEntryOfGID = Map.FirstPointInElement(LIDOfGID);
      if ((*x)[FirstEntryOfGID+GIDSize-1]!=VecValues[i]) err++;
      if (verbose) cout << "x["<<FirstEntryOfGID+GIDSize-1<<"] = "
			<< (*x)[FirstEntryOfGID+GIDSize-1]