C++ resid函数代码示例

static void mg3P(double **u, double *v, double **r, double a[4],
		 double c[4], int n1, int n2, int n3, int k) {

/*--------------------------------------------------------------------
c-------------------------------------------------------------------*/

/*--------------------------------------------------------------------
c     multigrid V-cycle routine
c-------------------------------------------------------------------*/

    int j;

/*--------------------------------------------------------------------
c     down cycle.
c     restrict the residual from the find grid to the coarse
c-------------------------------------------------------------------*/

    for (k = lt; k >= lb+1; k--) {
	j = k-1;
	rprj3(r[k], m1[k], m2[k], m3[k],
	      r[j], m1[j], m2[j], m3[j], k);
    }

    k = lb;
/*--------------------------------------------------------------------
c     compute an approximate solution on the coarsest grid
c-------------------------------------------------------------------*/
    zero3(u[k], m1[k], m2[k], m3[k]);
    psinv(r[k], u[k], m1[k], m2[k], m3[k], c, k);

    for (k = lb+1; k <= lt-1; k++) {
	j = k-1;
/*--------------------------------------------------------------------
c        prolongate from level k-1  to k
c-------------------------------------------------------------------*/
	zero3(u[k], m1[k], m2[k], m3[k]);
	interp(u[j], m1[j], m2[j], m3[j],
	       u[k], m1[k], m2[k], m3[k], k);
/*--------------------------------------------------------------------
c        compute residual for level k
c-------------------------------------------------------------------*/
	resid(u[k], r[k], r[k], m1[k], m2[k], m3[k], a, k);
/*--------------------------------------------------------------------
c        apply smoother
c-------------------------------------------------------------------*/
	psinv(r[k], u[k], m1[k], m2[k], m3[k], c, k);
    }

    j = lt - 1;
    k = lt;
    interp(u[j], m1[j], m2[j], m3[j], u[lt], n1, n2, n3, k);
    resid(u[lt], v, r[lt], n1, n2, n3, a, k);
    psinv(r[lt], u[lt], n1, n2, n3, c, k);
}

void EBPoissonOp::
levelJacobi(LevelData<EBCellFAB>&       a_phi,
            const LevelData<EBCellFAB>& a_rhs)
{
  CH_TIME("EBPoissonOp::levelJacobi");

  // Overhauled from previous in-place Gauss-Seidel-like method
  // This implementation is very inefficient in terms of memory usage,
  // and could be greatly improved.  It has the advantage, though,
  // of being very simple and easy to verify as correct.

  EBCellFactory factory(m_eblg.getEBISL());
  // Note: this is hardcoded to a single variable (component),
  // like some other code in this class
  LevelData<EBCellFAB> resid(m_eblg.getDBL(), 1, m_ghostCellsRHS, factory);

  residual(resid, a_phi, a_rhs, true);

  LevelData<EBCellFAB> relaxationCoeff(m_eblg.getDBL(), 1, m_ghostCellsRHS, factory);
  getJacobiRelaxCoeff(relaxationCoeff);

  EBLevelDataOps::scale(resid, relaxationCoeff);

  EBLevelDataOps::incr(a_phi, resid, 1.0);
}

void TrdHSegmentR::calChi2(){
  Chi2=0.;
  for(int i=0;i<nTrdRawHit();i++){
    TRDHitRZD rzd=TRDHitRZD(*pTrdRawHit(i));
    Chi2+=pow(resid(rzd.r,rzd.z,rzd.d)/ (0.62/sqrt(12.)),2);
  }
}

size_t
suiolite::capacity () const
{
  int inuse = resid ();
  assert (inuse <= len);
  return len - inuse;
}

int residuals_b_predicts_a(double ax[], double ay[], double bx[], double by[],
			   int use[], int n, double residuals[], double *rms)
{
    resid(ax, ay, bx, by, use, n, residuals, rms, 0);

    return 0;
}

bool
AmesosBTFGlobal_LinearProblem::
rvs()
{
  //  cout << "AmesosBTFGlobal_LinearProblem: NewLHS_" << endl;
  //  cout << *NewLHS_ << endl;

  OldLHS_->Import( *NewLHS_, *Importer2_, Insert );
  int numrhs = OldLHS_->NumVectors();
  std::vector<double> actual_resids( numrhs ), rhs_norm( numrhs );
  Epetra_MultiVector resid( OldLHS_->Map(), numrhs );
  OldMatrix_->Apply( *OldLHS_, resid );
  resid.Update( -1.0, *OldRHS_, 1.0 );
  resid.Norm2( &actual_resids[0] );
  OldRHS_->Norm2( &rhs_norm[0] );
  if (OldLHS_->Comm().MyPID() == 0 ) {
    for (int i=0; i<numrhs; i++ ) {
      std::cout << "Problem " << i << " (in AmesosBTFGlobal): \t" << actual_resids[i]/rhs_norm[i] << std::endl;
    }
  }
  //cout << "AmesosBTFGlobal_LinearProblem: OldLHS_" << endl;
  //cout << *OldLHS_ << endl;

  return true;
}

// Simple Power method algorithm
double power_method(const Epetra_CrsMatrix& A) {  
  // variable needed for iteration
  double lambda = 0.0;
  int niters = A.RowMap().NumGlobalElements()*10;
  double tolerance = 1.0e-10;
  // Create vectors
  Epetra_Vector q(A.RowMap());
  Epetra_Vector z(A.RowMap());
  Epetra_Vector resid(A.RowMap());
  // Fill z with random Numbers
  z.Random();
  // variable needed for iteration
  double normz;
  double residual = 0;
  int iter = 0;
  while (iter==0 || (iter < niters && residual > tolerance)) {
    z.Norm2(&normz); // Compute 2-norm of z
    q.Scale(1.0/normz, z);
    A.Multiply(false, q, z); // Compute z = A*q
    q.Dot(z, &lambda); // Approximate maximum eigenvalue
    if (iter%10==0 || iter+1==niters) {
      // Compute A*q - lambda*q every 10 iterations
      resid.Update(1.0, z, -lambda, q, 0.0);
      resid.Norm2(&residual);
      if (q.Map().Comm().MyPID()==0)
	std::cout << "Iter = " << iter << "  Lambda = " << lambda 
	     << "  Two-norm of A*q - lambda*q = " 
	     << residual << std::endl;
    } 
    iter++;
  }
  return(lambda);
}

int TpetraLinearSolver::residual_norm(int whichNorm, Teuchos::RCP<LinSys::Vector> sln, double& norm)
{
  LinSys::Vector resid(rhs_->getMap());
  ThrowRequire(! (sln.is_null()  || rhs_.is_null() ) );

  if (matrix_->isFillActive() )
  {
    // FIXME
    //!matrix_->fillComplete(map_, map_);
    throw std::runtime_error("residual_norm");
  }
  matrix_->apply(*sln, resid);

  LinSys::OneDVector rhs = rhs_->get1dViewNonConst ();
  LinSys::OneDVector res = resid.get1dViewNonConst ();
  for (int i=0; i<rhs.size(); ++i)
    res[i] -= rhs[i];
  if ( whichNorm == 0 )
    norm = resid.normInf();
  else if ( whichNorm == 1 )
    norm = resid.norm1();
  else if ( whichNorm == 2 )
    norm = resid.norm2();
  else
    return 1;

  return 0;
}

/* -----------------------------------------------------------------------------------
 * Solve: M x = b
 *
 */
int
gsl_sparse_matrix_complex_LU_solve(gsl_sparse_matrix_complex *spmat, double *b_real, double *b_imag, double *x_real, double *x_imag)
{
    void *Symbolic, *Numeric;
    int status;

    //gsl_sparse_matrix_complex_print_col(spmat);

    /* --- symbolic factorization --- */

    status = umfpack_zl_symbolic(spmat->n, spmat->n, spmat->Ap, spmat->Ai, spmat->Ax, spmat->Az, &Symbolic, spmat->Control, spmat->Info);
   
    if (status < 0)
        {
            umfpack_zl_report_info(spmat->Control, spmat->Info);
            //umfpack_zl_report_status(spmat->Control, status);
            fprintf(stderr, "%s: umfpack_zl_symbolic failed\n", __PRETTY_FUNCTION__);
        return -1;
        }

        //printf("%s: Symbolic factorization of A: ", __PRETTY_FUNCTION__);
    //umfpack_zl_report_symbolic(Symbolic, spmat->Control);

    /* --- numeric factorization --- */

    status = umfpack_zl_numeric(spmat->Ap, spmat->Ai, spmat->Ax, spmat->Az, Symbolic, &Numeric, spmat->Control, spmat->Info);

    if (status < 0)
    {
            umfpack_zl_report_info(spmat->Control, spmat->Info);
            umfpack_zl_report_status(spmat->Control, status);
            fprintf(stderr, "%s: umfpack_zl_numeric failed", __PRETTY_FUNCTION__);
        return -1;
        }

        //printf("%s: Numeric factorization of A: ", __PRETTY_FUNCTION__);
    //umfpack_zl_report_numeric(Numeric, spmat->Control);    

    /* --- Solve M x = b --- */
    
    status = umfpack_zl_solve(UMFPACK_A, spmat->Ap, spmat->Ai, spmat->Ax, spmat->Az, x_real, x_imag, b_real, b_imag, Numeric, spmat->Control, spmat->Info);

    //umfpack_zl_report_info(spmat->Control, spmat->Info);
        //umfpack_zl_report_status(spmat->Control, status);

        if (status < 0)
        {
            fprintf(stderr, "%s: umfpack_zl_solve failed\n", __PRETTY_FUNCTION__);
        }
       //printf("%s: x (solution of Ax=b): ") ;
    //umfpack_zl_report_vector(spmat->n, x_real, x_imag, spmat->Control);
    
    {
        double rnorm = resid(FALSE, spmat->Ap, spmat->Ai, spmat->Ax, spmat->Az, x_real, x_imag, b_real, b_imag, spmat->n);
           printf ("maxnorm of residual: %g\n\n", rnorm) ;
    }
    
    return 0;
}

bool IterativeSolvers::pcg(const IRCMatrix &A,
                           Vector &x,
                           const Vector &b,
                           const Preconditioner &M) {
    /*!
      Solves Ax=b using the preconditioned conjugate gradient method.
      */
    const idx N = x.getLength();
    real resid(100.0);
    Vector p(N), z(N), q(N);
    real alpha;
    real normr(0);
    real normb = norm(b);
    real rho(0), rho_1(0), beta(0);
    Vector r = b - A * x;
    if (normb == 0.0)
        normb = 1;
    resid = norm(r) / normb;
    if (resid <= IterativeSolvers::toler) {
        IterativeSolvers::toler = resid;
        IterativeSolvers::maxIter = 0;
        return true;
    }
    // MAIN LOOP
    idx i = 1;
    for (; i <= IterativeSolvers::maxIter; i++) {
        M.solveMxb(z, r);
        rho = dot(r, z);
        if (i == 1)
            p = z;
        else {
            beta = rho / rho_1;
            aypx(beta, p, z); // p = beta*p + z;
        }
        // CALCULATES q = A*p AND dp = dot(q,p)
        real dp = multiply_dot(A, p, q);
        alpha = rho / dp;
        normr = 0;
#ifdef USES_OPENMP
        #pragma omp parallel for reduction(+:normr)
#endif
        for (idx j = 0 ; j < N ; ++j) {
            x[j] += alpha * p[j]; // x + alpha(0) * p;
            r[j] -= alpha * q[j]; // r - alpha(0) * q;
            normr += r[j] * r[j];
        }
        normr = sqrt(normr);
        resid = normr / normb;
        if (resid <= IterativeSolvers::toler) {
            IterativeSolvers::toler = resid;
            IterativeSolvers::maxIter = i;
            return true;
        }
        rho_1 = rho;
    }
    IterativeSolvers::toler = resid;
    return false;
}

void
suiolite::grow (size_t ns)
{
  assert (resid () == 0);
  assert (bytes_read == 0);
  xfree (buf);
  len = min<int> (ns, SUIOLITE_MAX_BUFLEN);
  buf = static_cast<char *> (xmalloc (len));
  clear ();
}

//---------------------------------------------------------
void EulerShock2D::Report(bool bForce)
//---------------------------------------------------------
{
  CurvedEuler2D::Report(bForce);

  if (tstep>=1 && tstep <= resid.size()) 
  {
    // calculate residual
    // resid(tstep) = sqrt(sum(sum(sum((Q-oldQ).^2)))/(4*K*Np));
    // resid(tstep) = sqrt(sum(sum(sum((Q-oldQ).^2)))/(4*K*Np))/dt;


    DMat Qresid = sqr(Q-oldQ); double d4KNp = double(4*K*Np);
    resid(tstep) = sqrt(Qresid.sum()/d4KNp);

    if (eScramInlet == sim_type) {
      // scale residual
      resid(tstep) /= dt;
    }
  }
}

DoubleVector RowDoubleMatrix::conjugateGradient(const DoubleVector &B, double epsilon, unsigned int niter, bool printMessages, unsigned int messageStep) const
{
    DoubleVector X(size_, 0.0); // начальное приближение - вектор нулей
    DoubleVector resid(size_); // невязка
    DoubleVector direction; // направление поиска
    DoubleVector temp(size_); // ременное хранилище для обмена данными
    double resid_norm; // норма невязки
    double alpha;
    double beta;

    double resid_resid, resid_resid_new;

    residual(X, B, resid);

    direction = resid;

    resid_norm = resid.norm_2();

    if (printMessages) std::cout << "Начальная невязка: " << resid_norm << std::endl;
    if (resid_norm > epsilon)
    {
        resid_resid = resid * resid;
        for (unsigned int i = 0; i < niter; i++)
        {
            product(direction, temp);
//            std::cout << direction.norm_2() << "    " << temp.norm_2() << std::endl;
            alpha = (resid_resid) / (direction * temp);
            X += alpha * direction;
            resid -= alpha * temp;
            resid_resid_new = resid * resid;
            resid_norm = sqrt(resid_resid_new);
            if (resid_norm <= epsilon)
            {
                if (printMessages)
                    std::cout << "Решение найдено. Итераций: " << i << ", невязка: " << resid_norm << std::endl;
                break;
            }
            if (printMessages && (i % messageStep == 0))
                std::cout << i << ", невязка: " << resid_norm << std::endl;

            beta = (resid_resid_new) / (resid_resid);
            // d = r + d*beta
            direction.scale(beta);
            direction += resid;
            //
            resid_resid = resid_resid_new;
        }
    }
    return X;
}

int checkResults(Epetra_RowMatrix * A, Epetra_CrsMatrix * transA, 
                 Epetra_Vector * xexact, bool verbose) {

  int n = A->NumGlobalRows();

  if (n<100) cout << "A transpose = " << endl << *transA << endl;

  Epetra_Vector x1(View,A->OperatorDomainMap(), &((*xexact)[0]));
  Epetra_Vector b1(A->OperatorRangeMap());

  A->SetUseTranspose(true);

  Epetra_Time timer(A->Comm());
  double start = timer.ElapsedTime();
  A->Apply(x1, b1);
  if (verbose) cout << "\nTime to compute b1: matvec with original matrix using transpose flag  = " << timer.ElapsedTime() - start << endl;

  if (n<100) cout << "b1 = " << endl << b1 << endl;
  Epetra_Vector x2(View,transA->OperatorRangeMap(), &((*xexact)[0]));
  Epetra_Vector b2(transA->OperatorDomainMap());
  start = timer.ElapsedTime();
  transA->Multiply(false, x2, b2);
  if (verbose) cout << "\nTime to compute b2: matvec with transpose matrix                      = " << timer.ElapsedTime() - start << endl;

  if (n<100) cout << "b1 = " << endl << b1 << endl;

  double residual;
  Epetra_Vector resid(A->OperatorDomainMap());

  resid.Update(1.0, b1, -1.0, b2, 0.0);
  resid.Norm2(&residual);
  if (verbose) cout << "Norm of b1 - b2 = " << residual << endl;

  int ierr = 0;

  if (residual > 1.0e-10) ierr++;

  if (ierr!=0 && verbose) cerr << "Status: Test failed" << endl;
  else if (verbose) cerr << "Status: Test passed" << endl;

  return(ierr);
}

int
NineNodeMixedQuad::addInertiaLoadToUnbalance(const Vector &accel)
{
  static const int numberGauss = 9 ;
  static const int numberNodes = 9 ;
  static const int ndf = 2 ; 

  int i;

  // check to see if have mass
  int haveRho = 0;
  for (i = 0; i < numberGauss; i++) {
    if (materialPointers[i]->getRho() != 0.0)
      haveRho = 1;
  }

  if (haveRho == 0)
    return 0;

  // Compute mass matrix
  int tangFlag = 1 ;
  formInertiaTerms( tangFlag ) ;

  // store computed RV fro nodes in resid vector
  int count = 0;

  for (i=0; i<numberNodes; i++) {
    const Vector &Raccel = nodePointers[i]->getRV(accel);
    for (int j=0; j<ndf; j++)
      resid(count++) = Raccel(i);
  }

  // create the load vector if one does not exist
  if (load == 0) 
    load = new Vector(numberNodes*ndf);

  // add -M * RV(accel) to the load vector
  load->addMatrixVector(1.0, mass, resid, -1.0);
  
  return 0;
}