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;
}

void  ZeroLengthContact3D::formResidAndTangent( int tang_flag ) 

{



	// trial displacement vectors

 	Vector DispTrialS(3); // trial disp for slave node

	Vector DispTrialM(3); // trial disp for master node

	// trial frictional force vectors (in local coordinate)

    Vector t_trial(2);

    double TtrNorm;



    // Coulomb friction law surface

	double Phi;     



    int i, j;



    //zero stiffness and residual 

    stiff.Zero( ) ;

    resid.Zero( ) ;

   

	// detect contact and set flag

    ContactFlag = contactDetect(); 



	//opserr<<this->getTag()<< " ZeroLengthContact3D::ContactFlag=" << ContactFlag<<endln;





	if (ContactFlag == 1) // contacted

	{  

       // contact presure;

		pressure = Kn*gap;   // Kn : normal penalty

  

		DispTrialS=nodePointers[0]->getTrialDisp();

        DispTrialM=nodePointers[1]->getTrialDisp();



       //nodal displacements 

        double ul[6];



		ul[0]=DispTrialS(0);

		ul[1]=DispTrialS(1);

		ul[2]=DispTrialS(2);

		ul[3]=DispTrialM(0);

		ul[4]=DispTrialM(1);

		ul[5]=DispTrialM(2);



		t_trial.Zero();

        xi.Zero();	



		for (i=0; i<6; i++){

			xi(0) += T1(i)*ul[i];

			xi(1) += T2(i)*ul[i];

		}

//.........这里部分代码省略.........

//.........这里部分代码省略.........

  Epetra_CrsMatrix A(Copy, Map, &NumNz[0]);
  EPETRA_TEST_ERR(A.IndicesAreGlobal(),ierr);
  EPETRA_TEST_ERR(A.IndicesAreLocal(),ierr);

  // Add  rows one-at-a-time
  // Need some vectors to help
  // Off diagonal Values will always be -1


  vector<double> Values(2);
  Values[0] = -1.0;
	Values[1] = -1.0;
	vector<long long> Indices(2);
  double two = 2.0;
  int NumEntries;

  forierr = 0;
  for(i = 0; i < NumMyEquations; i++) {
    if(MyGlobalElements[i] == 0) {
			Indices[0] = 1;
			NumEntries = 1;
		}
    else if (MyGlobalElements[i] == NumGlobalEquations-1) {
			Indices[0] = NumGlobalEquations-2;
			NumEntries = 1;
		}
    else {
			Indices[0] = MyGlobalElements[i]-1;
			Indices[1] = MyGlobalElements[i]+1;
			NumEntries = 2;
		}
		forierr += !(A.InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0])==0);
		forierr += !(A.InsertGlobalValues(MyGlobalElements[i], 1, &two, &MyGlobalElements[i])>0); // Put in the diagonal entry
  }
  EPETRA_TEST_ERR(forierr,ierr);

  // Finish up
  A.FillComplete();
  A.OptimizeStorage();

  Epetra_JadMatrix JadA(A);
  Epetra_JadMatrix JadA1(A);
  Epetra_JadMatrix JadA2(A);

  // Create vectors for Power method

  Epetra_Vector q(Map);
  Epetra_Vector z(Map); z.Random();
  Epetra_Vector resid(Map);

  Epetra_Flops flopcounter;
  A.SetFlopCounter(flopcounter);
  q.SetFlopCounter(A);
  z.SetFlopCounter(A);
  resid.SetFlopCounter(A);
  JadA.SetFlopCounter(A);
  JadA1.SetFlopCounter(A);
  JadA2.SetFlopCounter(A);


  if (verbose) cout << "=======================================" << endl
		    << "Testing Jad using CrsMatrix as input..." << endl
		    << "=======================================" << endl;

  A.ResetFlops();
  powerMethodTests(A, JadA, Map, q, z, resid, verbose);

  // Increase diagonal dominance

  if (verbose) cout << "\n\nIncreasing the magnitude of first diagonal term and solving again\n\n"
		    << endl;


  if (A.MyGlobalRow(0)) {
    int numvals = A.NumGlobalEntries(0);
    vector<double> Rowvals(numvals);
    vector<long long> Rowinds(numvals);
    A.ExtractGlobalRowCopy(0, numvals, numvals, &Rowvals[0], &Rowinds[0]); // Get A[0,0]

    for (i=0; i<numvals; i++) if (Rowinds[i] == 0) Rowvals[i] *= 10.0;

    A.ReplaceGlobalValues(0, numvals, &Rowvals[0], &Rowinds[0]);
  }
  JadA.UpdateValues(A);
  A.ResetFlops();
  powerMethodTests(A, JadA, Map, q, z, resid, verbose);

  if (verbose) cout << "================================================================" << endl
		          << "Testing Jad using Jad matrix as input matrix for construction..." << endl
		          << "================================================================" << endl;
  JadA1.ResetFlops();
  powerMethodTests(JadA1, JadA2, Map, q, z, resid, verbose);

#ifdef EPETRA_MPI
  MPI_Finalize() ;
#endif

return ierr ;
}

//.........这里部分代码省略.........
    serialA = &transposeA ; 
  } else {
    serialA = readA ; 
  }

  // Create uniform distributed map
  Epetra_Map map(readMap->NumGlobalElements(), 0, Comm);
  Epetra_Map* map_;

  if( NonContiguousMap ) {
    //
    //  map gives us NumMyElements and MyFirstElement;
    //
    int NumGlobalElements =  readMap->NumGlobalElements();
    int NumMyElements = map.NumMyElements();
    int MyFirstElement = map.MinMyGID();
    std::vector<int> MapMap_( NumGlobalElements );
    readMap->MyGlobalElements( &MapMap_[0] ) ;
    Comm.Broadcast( &MapMap_[0], NumGlobalElements, 0 ) ; 
    map_ = new Epetra_Map( NumGlobalElements, NumMyElements, &MapMap_[MyFirstElement], 0, Comm);
  } else {
    map_ = new Epetra_Map( map ) ; 
  }


  // Create Exporter to distribute read-in matrix and vectors
  Epetra_Export exporter(*readMap, *map_);
  Epetra_CrsMatrix A(Copy, *map_, 0);

  Epetra_RowMatrix * passA = 0; 
  Epetra_MultiVector * passx = 0; 
  Epetra_MultiVector * passb = 0;
  Epetra_MultiVector * passxexact = 0;
  Epetra_MultiVector * passresid = 0;
  Epetra_MultiVector * passtmp = 0;

  Epetra_MultiVector x(*map_,numsolves);
  Epetra_MultiVector b(*map_,numsolves);
  Epetra_MultiVector xexact(*map_,numsolves);
  Epetra_MultiVector resid(*map_,numsolves);
  Epetra_MultiVector tmp(*map_,numsolves);

  Epetra_MultiVector serialx(*readMap,numsolves);
  Epetra_MultiVector serialb(*readMap,numsolves);
  Epetra_MultiVector serialxexact(*readMap,numsolves);
  Epetra_MultiVector serialresid(*readMap,numsolves);
  Epetra_MultiVector serialtmp(*readMap,numsolves);

  bool distribute_matrix = ( matrix_type == AMESOS_Distributed ) ; 
  if ( distribute_matrix ) { 
    //
    //  Initialize x, b and xexact to the values read in from the file
    //
    
    A.Export(*serialA, exporter, Add);
    Comm.Barrier();

    assert(A.FillComplete()==0);    
    Comm.Barrier();

    passA = &A; 
    passx = &x; 
    passb = &b;
    passxexact = &xexact;
    passresid = &resid;
    passtmp = &tmp;

void mglin(float ***u, int n, int ncycle){
/*
  Full Multigrid Algorithm for solution of the steady state heat
  equation with forcing.  On input u[1..n][1..n] contains the
  right-hand side ρ, while on output it returns the solution.  The
  dimension n must be of the form 2j + 1 for some integer j. (j is
  actually the number of grid levels used in the solution, called ng
  below.) ncycle is the number of V-cycles to be used at each level.
*/
  unsigned int j,jcycle,jj,jpost,jpre,nf,ng=0,ngrid,nn;
  /*** setup multigrid jagged arrays ***/
  float ***iu[NGMAX+1];   /* stores solution at each grid level */
  float ***irhs[NGMAX+1]; /* stores rhs at each grid level */
  float ***ires[NGMAX+1]; /* stores residual at each grid level */
  float ***irho[NGMAX+1]; /* stores rhs during intial solution of FMG */
  
  /*** use bitshift to find the number of grid levels, stored in ng ***/
  nn=n;                   
  while (nn >>= 1) ng++;     
  
  /*** some simple input checks ***/
  if (n != 1+(1L << ng)) nrerror("n-1 must be a power of 2 in mglin.");
  if (ng > NGMAX) nrerror("increase NGMAX in mglin.");
  
  /***restrict solution to next coarsest grid (irho[ng-1])***/
  nn=n/2+1;
  ngrid=ng-1;
  irho[ngrid]=f3tensor(1,nn,1,nn,1,nn); 
  rstrct(irho[ngrid],u,nn);/* coarsens rhs (u at this point) to irho on mesh size nn */
  
  /***continue setting up coarser grids down to coarsest level***/
  while (nn > 3) { 
    nn=nn/2+1; 
    irho[--ngrid]=f3tensor(1,nn,1,nn,1,nn);
    rstrct(irho[ngrid],irho[ngrid+1],nn); 
  }
  
  /***now setup and solve coarsest level iu[1],irhs[1] ***/
  nn=3;
  iu[1]=f3tensor(1,nn,1,nn,1,nn);
  irhs[1]=f3tensor(1,nn,1,nn,1,nn);
  slvsml(iu[1],irho[1]);          /* solve the small system directly */
  free_f3tensor(irho[1],1,nn,1,nn,1,nn);
  ngrid=ng;                       /* reset ngrid to original size */

  for (j=2;j<=ngrid;j++) {        /* loop over coarse to fine, starting at level 2 */
    printf("at grid level %d\n",j);
    nn=2*nn-1;                     
    iu[j]=f3tensor(1,nn,1,nn,1,nn);     /* setup grids for lhs,rhs, and residual */
    irhs[j]=f3tensor(1,nn,1,nn,1,nn);
    ires[j]=f3tensor(1,nn,1,nn,1,nn);
    interp(iu[j],iu[j-1],nn);
    /* irho contains rhs except on fine grid where it is in u */
    copy(irhs[j],(j != ngrid ? irho[j] : u),nn); 
    /* v-cycle at current grid level */
    for (jcycle=1;jcycle<=ncycle;jcycle++) {  
      /* nf is # points on finest grid for current v-sweep */
      nf=nn;                                  
      for (jj=j;jj>=2;jj--) {                 
	for (jpre=1;jpre<=NPRE;jpre++)  /* NPRE g-s sweeps on the finest (relatively) scale */
	  relax(iu[jj],iu[jj-1],irhs[jj],nf); //need iu[jj-1] for jacobi
	resid(ires[jj],iu[jj],irhs[jj],nf); /* compute res on finest scale, store in ires */
	nf=nf/2+1;                        /* next coarsest scale */
	rstrct(irhs[jj-1],ires[jj],nf);  /* restrict residuals as rhs of next coarsest scale */
	fill0(iu[jj-1],nf);              /* set the initial solution guess to zero */
      } 
      slvsml(iu[1],irhs[1]);                  /* solve the small problem exactly */
      nf=3;                                   /* fine scale now n=3 */
      for (jj=2;jj<=j;jj++) {                 /* work way back up to current finest grid */
	nf=2*nf-1;                            /* next finest scale */
	addint(iu[jj],iu[jj-1],ires[jj],nf);  /* inter error and add to previous solution guess */
	for (jpost=1;jpost<=NPOST;jpost++)    /* do NPOST g-s sweeps */
	  relax(iu[jj],iu[jj-1],irhs[jj],nf);
      }
    }
  }

  copy(u,iu[ngrid],n);              /* copy solution into input array (implicitly returned) */

  /*** clean up memory ***/
  for (nn=n,j=ng;j>=2;j--,nn=nn/2+1) {       
    free_f3tensor(ires[j],1,nn,1,nn,1,nn);      
    free_f3tensor(irhs[j],1,nn,1,nn,1,nn);
    free_f3tensor(iu[j],1,nn,1,nn,1,nn);
    if (j != ng) free_f3tensor(irho[j],1,nn,1,nn,1,nn);
  }
  free_f3tensor(irhs[1],1,3,1,3,1,3);
  free_f3tensor(iu[1],1,3,1,3,1,3);
}

//get residual with inertia terms
const XC::Vector &XC::Twenty_Node_Brick::getResistingForceIncInertia(void) const
  {
        static XC::Vector res(60);

//        printf("getResistingForceIncInertia()\n");

        int i, j;
        static double a[60];

        for(i=0; i<nenu; i++) {
                const XC::Vector &accel = theNodes[i]->getTrialAccel();
                if( 3 != accel.Size() ) {
                        std::cerr << "XC::Twenty_Node_Brick::getResistingForceIncInertia matrix and vector sizes are incompatable\n";
                        exit(-1);
                }

                a[i*3] = accel(0);
                a[i*3+1] = accel(1);
                a[i*3+2] = accel(2);
        }
        // Compute the current resisting force
        this->getResistingForce();
//        std::cerr<<"K "<<resid<<std::endl;

        // Compute the mass matrix
        this->getMass();

        for(i = 0; i < 60; i++) {
                for(j = 0; j < 60; j++){
                        resid(i) += mass(i,j)*a[j];
                }
        }
//        printf("\n");
        //std::cerr<<"K+M "<<P<<std::endl;


        for(i=0; i<nenu; i++) {
                const XC::Vector &vel = theNodes[i]->getTrialVel();
                if( 3!= vel.Size() ) {
                        std::cerr << "XC::Twenty_Node_Brick::getResistingForceIncInertia matrix and vector sizes are incompatable\n";
                        exit(-1);
                }
                a[i*3] = vel(0);
                a[i*3+1] = vel(1);
                a[i*3+2] = vel(2);
        }

        this->getDamp();

        for(i = 0; i < 60; i++) {
                for(j = 0; j < 60; j++) {
                        resid(i) += damp(i,j)*a[j];
                }
        }
//        std::cerr<<"Pd"<<Pd<<std::endl;

        res = resid;
//        std::cerr<<"res "<<res<<std::endl;

//        exit(-1);
    if(isDead())
      res*=dead_srf;
    return res;
  }