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

void NOX::Epetra::Scaling::computeScaling(const Epetra_LinearProblem& problem)
{

  Epetra_Vector* diagonal = 0;
  for (unsigned int i = 0; i < scaleVector.size(); i ++) {

    if (sourceType[i] == RowSum) {

      diagonal = scaleVector[i].get();

      // Make sure the Jacobian is an Epetra_RowMatrix, otherwise we can't
      // perform a row sum scale!
      const Epetra_RowMatrix* test = 0;
      test = dynamic_cast<const Epetra_RowMatrix*>(problem.GetOperator());
      if (test == 0) {
    std::cout << "ERROR: NOX::Epetra::Scaling::scaleLinearSystem() - "
         << "For \"Row Sum\" scaling, the Matrix must be an "
         << "Epetra_RowMatrix derived object!" << std::endl;
    throw "NOX Error";
      }

      test->InvRowSums(*diagonal);
      diagonal->Reciprocal(*diagonal);

    }

    else if (sourceType[i] == ColSum) {

      diagonal = scaleVector[i].get();

      // Make sure the Jacobian is an Epetra_RowMatrix, otherwise we can't
      // perform a row sum scale!
      const Epetra_RowMatrix* test = 0;
      test = dynamic_cast<const Epetra_RowMatrix*>(problem.GetOperator());
      if (test == 0) {
    std::cout << "ERROR: NOX::Epetra::Scaling::scaleLinearSystem() - "
         << "For \"Column Sum\" scaling, the Matrix must be an "
         << "Epetra_RowMatrix derived object!" << std::endl;
    throw "NOX Error";
      }

      test->InvColSums(*diagonal);
      diagonal->Reciprocal(*diagonal);

    }

  }

}

void AZOO_iterate(double * xsolve, double * b, 
		  int * options, double * params, 
		  double * status, int *proc_config,
		  AZ_MATRIX * Amat,
		  AZ_PRECOND *precond, struct AZ_SCALING *scaling)
{
  (void)precond;
  (void)scaling;
  bool verbose = (options[AZ_output]!=AZ_none); // Print info unless all output is turned off

  Epetra_Comm * comm;
  Epetra_BlockMap * map;
  Epetra_RowMatrix * A;
  Epetra_Vector * px;
  Epetra_Vector * pb;
  int * global_indices;

  int ierr = Aztec2Petra(proc_config, Amat, xsolve, b, comm, map, A, px, pb, &global_indices);
  if (ierr!=0) {
    cerr << "Error detected in Aztec2Petra. Value = " << ierr << endl;
    exit(1);
  }


  Epetra_LinearProblem problem(A, px, pb);

  Epetra_Vector * leftScaleVec = 0;
  Epetra_Vector * rightScaleVec = 0;
  bool doRowScaling = false;
  bool doColScaling = false;
  
  if ((options[AZ_scaling]==AZ_Jacobi) || options[AZ_scaling]==AZ_BJacobi) {
    doRowScaling = true;
    leftScaleVec = new Epetra_Vector(*map);
    A->ExtractDiagonalCopy(*leftScaleVec); // Extract diagonal of matrix
    leftScaleVec->Reciprocal(*leftScaleVec); // invert it
  }

  else if (options[AZ_scaling]==AZ_row_sum) {
    doRowScaling = true;
    leftScaleVec = new Epetra_Vector(*map);
    A->InvRowSums(*leftScaleVec);
  }
  else if (options[AZ_scaling]==AZ_sym_diag) {
    doRowScaling = true;
    doColScaling = true;
    leftScaleVec = new Epetra_Vector(*map);
    A->ExtractDiagonalCopy(*leftScaleVec); // Extract diagonal of matrix

    int length = leftScaleVec->MyLength();
    for (int i=0; i<length; i++) (*leftScaleVec)[i] = sqrt(fabs((*leftScaleVec)[i])); // Take its sqrt

    rightScaleVec = leftScaleVec; // symmetric, so left and right the same
    leftScaleVec->Reciprocal(*leftScaleVec); // invert it
  }
  else if (options[AZ_scaling]==AZ_sym_row_sum) {
    doRowScaling = true;
    doColScaling = true;
    leftScaleVec = new Epetra_Vector(*map);
    A->InvRowSums(*leftScaleVec);
    int length = leftScaleVec->MyLength();
    for (int i=0; i<length; i++) (*leftScaleVec)[i] = sqrt(fabs((*leftScaleVec)[i])); // Take its sqrt

    rightScaleVec = leftScaleVec; // symmetric, so left and right the same
  }
  if ((doRowScaling || doColScaling) && verbose) {
    double norminf = A->NormInf();
    double normone = A->NormOne();
    if (comm->MyPID()==0) 
      cout << "\n Inf-norm of A before scaling = " << norminf 
	   << "\n One-norm of A before scaling = " << normone<< endl << endl;
  }
  if (doRowScaling) problem.LeftScale(*leftScaleVec);
  if (doColScaling) problem.RightScale(*rightScaleVec);

  if ((doRowScaling || doColScaling) && verbose) {
    double norminf = A->NormInf();
    double normone = A->NormOne();
    if (comm->MyPID()==0) 
      cout << "\n Inf-norm of A after  scaling = " << norminf  
	   << "\n One-norm of A after  scaling = " << normone << endl << endl;
  }



  AztecOO solver(problem);

  solver.SetAllAztecParams(params); // set all AztecOO params with user-provided params
  solver.SetAllAztecOptions(options); // set all AztecOO options with user-provided options

  solver.CheckInput();
  solver.SetAztecOption(AZ_scaling, AZ_none); // Always must have scaling off
  solver.Iterate(options[AZ_max_iter], params[AZ_tol]);
  solver.GetAllAztecStatus(status);
  
  if (doColScaling) {
    rightScaleVec->Reciprocal(*rightScaleVec);
    problem.RightScale(*rightScaleVec);
  }
  if (doRowScaling) {
//.........这里部分代码省略.........