C++ teuchos::LAPACK类代码示例

int main(int argc, char* argv[])
{
#ifdef HAVE_MPI
  MPI_Init(&argc,&argv);
#endif

  // Creating an instance of the LAPACK class for double-precision routines looks like:
  Teuchos::LAPACK<int, double> lapack;

  // This instance provides the access to all the LAPACK routines.
  Teuchos::SerialDenseMatrix<int, double> My_Matrix(4,4);
  Teuchos::SerialDenseVector<int, double> My_Vector(4);
  My_Matrix.random();
  My_Vector.random();

  // Perform an LU factorization of this matrix. 
  int ipiv[4], info;
  char TRANS = 'N';
  lapack.GETRF( 4, 4, My_Matrix.values(), My_Matrix.stride(), ipiv, &info ); 
  
  // Solve the linear system.
  lapack.GETRS( TRANS, 4, 1, My_Matrix.values(), My_Matrix.stride(), 
		ipiv, My_Vector.values(), My_Vector.stride(), &info );  

  // Print out the solution.
  cout << My_Vector << endl;

#ifdef HAVE_MPI
  MPI_Finalize();
#endif
  return 0;
}

 void solve_state(std::vector<std::vector<Real> > &U, const std::vector<Real> &z) {
   // Get Diagonal and Off-Diagonal Entries of PDE Jacobian
   std::vector<Real> d(nx_,4.0*dx_/6.0 + dt_*2.0/dx_);
   d[0]           = dx_/3.0 + dt_/dx_;
   d[nx_-1] = dx_/3.0 + dt_/dx_;
   std::vector<Real> o(nx_-1,dx_/6.0 - dt_/dx_);
   // Perform LDL factorization
   Teuchos::LAPACK<int,Real> lp;
   int info;
   int ldb  = nx_;
   int nhrs = 1;
   lp.PTTRF(nx_,&d[0],&o[0],&info);
   // Initialize State Storage
   U.clear();
   U.resize(nt_+1);
   (U[0]).assign(u0_.begin(),u0_.end());
   // Time Step Using Implicit Euler
   std::vector<Real> b(nx_,0.0);
   for ( uint t = 0; t < nt_; t++ ) {
     // Get Right Hand Side
     apply_mass(b,U[t]);
     b[nx_-1] += dt_*z[t];
     // Solve Tridiagonal System Using LAPACK's SPD Tridiagonal Solver
     lp.PTTRS(nx_,nhrs,&d[0],&o[0],&b[0],ldb,&info);
     // Update State Storage
     (U[t+1]).assign(b.begin(),b.end());
   }
 }

 void solve_adjoint_sensitivity(std::vector<std::vector<Real> > &P, const std::vector<std::vector<Real> > &U) {
   // Get Diagonal and Off-Diagonal Entries of PDE Jacobian
   std::vector<Real> d(nx_,4.0*dx_/6.0 + dt_*2.0/dx_);
   d[0]           = dx_/3.0 + dt_/dx_;
   d[nx_-1] = dx_/3.0 + dt_/dx_;
   std::vector<Real> o(nx_-1,dx_/6.0 - dt_/dx_);
   // Perform LDL factorization
   Teuchos::LAPACK<int,Real> lp;
   int info;
   int ldb  = nx_;
   int nhrs = 1;
   lp.PTTRF(nx_,&d[0],&o[0],&info);
   // Initialize State Storage
   P.clear();
   P.resize(nt_);
   // Time Step Using Implicit Euler
   std::vector<Real> b(nx_,0.0);
   for ( uint t = nt_; t > 0; t-- ) {
     // Get Right Hand Side
     if ( t == nt_ ) {
       apply_mass(b,U[t-1]);
     }
     else {
       apply_mass(b,P[t]);
     }
     // Solve Tridiagonal System Using LAPACK's SPD Tridiagonal Solver
     lp.PTTRS(nx_,nhrs,&d[0],&o[0],&b[0],ldb,&info);
     // Update State Storage
     (P[t-1]).assign(b.begin(),b.end());
   }
 }

  KOKKOS_INLINE_FUNCTION
  int
  Chol<Uplo::Upper,
       AlgoChol::ExternalLapack,Variant::One>
  ::invoke(PolicyType &policy,
           const MemberType &member,
           DenseExecViewTypeA &A) {
    // static_assert( Kokkos::Impl::is_same<
    //                typename DenseMatrixTypeA::space_type,
    //                Kokkos::Cuda
    //                >::value,
    //                "Cuda space is not available for calling external BLAS" );

    //typedef typename DenseExecViewTypeA::space_type   space_type;
    typedef typename DenseExecViewTypeA::ordinal_type ordinal_type;
    typedef typename DenseExecViewTypeA::value_type   value_type;

    int r_val = 0;      
    if (member.team_rank() == 0) {
#ifdef HAVE_SHYLUTACHO_TEUCHOS
      Teuchos::LAPACK<ordinal_type,value_type> lapack;

      lapack.POTRF('U',
                   A.NumRows(),
                   A.ValuePtr(), A.BaseObject().ColStride(),
                   &r_val);
#else
    TACHO_TEST_FOR_ABORT( true, MSG_NOT_HAVE_PACKAGE("Teuchos") );
#endif
    }

    return r_val;
  }

MxDimMatrix<T, DIM> MxDimMatrix<T, DIM>::inv() const {

  MxDimMatrix<T, DIM> thisT = this->transpose(); 
  MxDimMatrix<T, DIM> solnT(MxDimVector<T, DIM>(1));

  int info;
  int ipiv[DIM];
  Teuchos::LAPACK<int, T> lapack;
  lapack.GESV(DIM, DIM, &thisT(0, 0), DIM, ipiv, &solnT(0, 0), DIM, &info);
  
  if (info != 0)
    std::cout << "MxDimMatrix::inv(): error in lapack routine. Return code: " << info << ".\n";

  return solnT.transpose();
}

int MxDimMatrix<T, DIM>::eig(MxDimVector<std::complex<T>, DIM> & evals,
MxDimMatrix<std::complex<T>, DIM> & levecs,
MxDimMatrix<std::complex<T>, DIM> & revecs) const {

  MxDimMatrix<T, DIM> thisT = this->transpose(); 

  MxDimVector<T, DIM> rva, iva;
  MxDimMatrix<T, DIM> rve, lve;

  int info;
  int lwork = 5 * DIM;
  T work[lwork];
  Teuchos::LAPACK<int, T> lapack;
  lapack.GEEV('V', 'V', DIM, &thisT(0, 0), DIM, &rva[0], &iva[0], &lve(0, 0), DIM, &rve(0, 0), DIM, work, lwork, &info);
  
  if (info != 0)
    std::cout << "MxDimMatrix::eig(): error in lapack routine. Return code: " << info << ".\n";

  // gather results
  for (size_t i = 0; i < DIM; ++i) {
    if (iva[i] != 0) {
      for (size_t k = 0; k < 2; ++k) {
        evals[i + k].real() = rva[i + k];
        evals[i + k].imag() = iva[i + k];
        for (size_t j = 0; j < DIM; ++j) {
          levecs(i + k, j).real() = lve(j, i);
          levecs(i + k, j).imag() = (k == 0 ? 1.0 : -1.0) * lve(j, i + 1);
          revecs(i + k, j).real() = rve(j, i);
          revecs(i + k, j).imag() = (k == 0 ? 1.0 : -1.0) * rve(j, i + 1);
        }
      }
      i++;
    }
    else {
      evals[i].real() = rva[i];
      evals[i].imag() = iva[i];
      for (size_t j = 0; j < DIM; ++j) {
        levecs(i, j).real() = lve(j, i);
        revecs(i, j).real() = rve(j, i);
      }
    }
  }

  return info;
}

  void Constraint<Scalar, LocalOrdinal, GlobalOrdinal, Node, LocalMatOps>::Setup(const MultiVector& B, const MultiVector& Bc, RCP<const CrsGraph> Ppattern) {
    Ppattern_ = Ppattern;

    const RCP<const Map> uniqueMap    = Ppattern_->getDomainMap();
    const RCP<const Map> nonUniqueMap = Ppattern_->getColMap();
    RCP<const Import> importer = ImportFactory::Build(uniqueMap, nonUniqueMap);

    const size_t NSDim = Bc.getNumVectors();
    X_ = MultiVectorFactory::Build(nonUniqueMap, NSDim);
    X_->doImport(Bc, *importer, Xpetra::INSERT);

    size_t numRows = Ppattern_->getNodeNumRows();
    XXtInv_.resize(numRows);
    Teuchos::SerialDenseVector<LO,SC> BcRow(NSDim, false);
    for (size_t i = 0; i < numRows; i++) {
      Teuchos::ArrayView<const LO> indices;
      Ppattern_->getLocalRowView(i, indices);

      size_t nnz = indices.size();

      Teuchos::SerialDenseMatrix<LO,SC> locX(NSDim, nnz, false);
      for (size_t j = 0; j < nnz; j++) {
        for (size_t k = 0; k < NSDim; k++)
          BcRow[k] = X_->getData(k)[indices[j]];

        Teuchos::setCol(BcRow, (LO)j, locX);
      }

      XXtInv_[i] = Teuchos::SerialDenseMatrix<LO,SC>(NSDim, NSDim, false);

      Teuchos::BLAS<LO,SC> blas;
      blas.GEMM(Teuchos::NO_TRANS, Teuchos::CONJ_TRANS, NSDim, NSDim, nnz,
                Teuchos::ScalarTraits<SC>::one(), locX.values(), locX.stride(),
                locX.values(), locX.stride(), Teuchos::ScalarTraits<SC>::zero(),
                XXtInv_[i].values(), XXtInv_[i].stride());

      Teuchos::LAPACK<LO,SC> lapack;
      LO info, lwork = 3*NSDim;
      ArrayRCP<LO> IPIV(NSDim);
      ArrayRCP<SC> WORK(lwork);
      lapack.GETRF(NSDim, NSDim, XXtInv_[i].values(), XXtInv_[i].stride(), IPIV.get(), &info);
      lapack.GETRI(NSDim, XXtInv_[i].values(), XXtInv_[i].stride(), IPIV.get(), WORK.get(), lwork, &info);
    }
  }

int MxDimMatrix<T, DIM>::solve(MxDimVector<T, DIM> & x, const MxDimVector<T, DIM> & b) const {
  x = b;

  MxDimMatrix<T, DIM> copy(*this);

  Teuchos::LAPACK<int, T> lapack;
  MxDimVector<int, DIM> ipiv;
  //int ipiv[DIM];
  int info;

  lapack.GETRF(DIM, DIM, &copy(0, 0), DIM, &ipiv[0], &info);
  if (info != 0)
    std::cout << "MxDimMatrix::solve(...): error in lapack routine getrf. Return code: " << info << ".\n";

  lapack.GETRS('T', DIM, 1, &copy(0, 0), DIM, &ipiv[0], &x[0], DIM, &info);
  if (info != 0)
    std::cout << "MxDimMatrix::solve(...): error in lapack routine getrs. Return code: " << info << ".\n";

  return info;
}

void updateGuess(Teuchos::SerialDenseVector<int, std::complex<double> >& myCurrentGuess,
		Teuchos::SerialDenseVector<int, std::complex<double> >& myTargetsCalculated,
		Teuchos::SerialDenseMatrix<int, std::complex<double> >& myJacobian, 
		Teuchos::LAPACK<int, std::complex<double> >& myLAPACK
		 )
{
	//v = J(inverse) * (-F(x))
	//new guess = v + old guess
	myTargetsCalculated *= -1.0;

	//Perform an LU factorization of this matrix. 
	int ipiv[NUMDIMENSIONS], info;
	char TRANS = 'N';
	myLAPACK.GETRF( NUMDIMENSIONS, NUMDIMENSIONS, myJacobian.values(), myJacobian.stride(), ipiv, &info ); 

	// Solve the linear system.
	myLAPACK.GETRS( TRANS, NUMDIMENSIONS, 1, myJacobian.values(), myJacobian.stride(),
		       	ipiv, myTargetsCalculated.values(), myTargetsCalculated.stride(), &info );  

	//We have overwritten myTargetsCalculated with guess update values
	//myBLAS.AXPY(NUMDIMENSIONS, 1.0, myGuessAdjustment.values(), 1, myCurrentGuess.values(), 1);
	myCurrentGuess += myTargetsCalculated;
}

void DenseContainer<MatrixType, LocalScalarType>::factor ()
{
    Teuchos::LAPACK<int, local_scalar_type> lapack;
    int INFO = 0;
    lapack.GETRF (diagBlock_.numRows (), diagBlock_.numCols (),
                  diagBlock_.values (), diagBlock_.stride (),
                  ipiv_.getRawPtr (), &INFO);
    // INFO < 0 is a bug.
    TEUCHOS_TEST_FOR_EXCEPTION(
        INFO < 0, std::logic_error, "Ifpack2::DenseContainer::factor: "
        "LAPACK's _GETRF (LU factorization with partial pivoting) was called "
        "incorrectly.  INFO = " << INFO << " < 0.  "
        "Please report this bug to the Ifpack2 developers.");
    // INFO > 0 means the matrix is singular.  This is probably an issue
    // either with the choice of rows the rows we extracted, or with the
    // input matrix itself.
    TEUCHOS_TEST_FOR_EXCEPTION(
        INFO > 0, std::runtime_error, "Ifpack2::DenseContainer::factor: "
        "LAPACK's _GETRF (LU factorization with partial pivoting) reports that the "
        "computed U factor is exactly singular.  U(" << INFO << "," << INFO << ") "
        "(one-based index i) is exactly zero.  This probably means that the input "
        "matrix has a singular diagonal block.");
}

NOX::Abstract::Group::ReturnType
LOCA::EigenvalueSort::LargestMagnitude::sort(int n, double* r_evals,
                         double* i_evals,
                         std::vector<int>* perm) const
{
  int i, j;
  int tempord = 0;
  double temp, tempr, tempi;
  Teuchos::LAPACK<int,double> lapack;

  //
  // Reset the index
  if (perm) {
    for (i=0; i < n; i++) {
      (*perm)[i] = i;
    }
  }

  //---------------------------------------------------------------
  // Sort eigenvalues in decreasing order of magnitude
  //---------------------------------------------------------------
  for (j=1; j < n; ++j) {
    tempr = r_evals[j]; tempi = i_evals[j];
    if (perm)
      tempord = (*perm)[j];
    temp=lapack.LAPY2(r_evals[j],i_evals[j]);
    for (i=j-1; i>=0 && lapack.LAPY2(r_evals[i],i_evals[i])<temp; --i) {
      r_evals[i+1]=r_evals[i]; i_evals[i+1]=i_evals[i];
      if (perm)
    (*perm)[i+1]=(*perm)[i];
    }
    r_evals[i+1] = tempr; i_evals[i+1] = tempi;
    if (perm)
      (*perm)[i+1] = tempord;
  }
  return NOX::Abstract::Group::Ok;
}

void gauss( const Teuchos::LAPACK<int,Real> &lapack,
            const ROL::Vector<Real> &a,
            const ROL::Vector<Real> &b,
            ROL::Vector<Real> &x,
            ROL::Vector<Real> &w ) {
    int INFO;

    Teuchos::RCP<const std::vector<Real> > ap = 
        (Teuchos::dyn_cast<ROL::StdVector<Real> >(const_cast<ROL::Vector<Real> &>(a))).getVector();
    Teuchos::RCP<const std::vector<Real> > bp = 
        (Teuchos::dyn_cast<ROL::StdVector<Real> >(const_cast<ROL::Vector<Real> &>(b))).getVector();
    Teuchos::RCP<std::vector<Real> > xp = 
        Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<ROL::StdVector<Real> >(x)).getVector()); 
    Teuchos::RCP<std::vector<Real> > wp = 
        Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<ROL::StdVector<Real> >(w)).getVector()); 

    const int N = ap->size();  
    const int LDZ = N;
    const char COMPZ = 'I';

    Teuchos::RCP<std::vector<Real> > Dp = Teuchos::rcp(new std::vector<Real>(N,0.0));
    Teuchos::RCP<std::vector<Real> > Ep = Teuchos::rcp(new std::vector<Real>(N,0.0));
    Teuchos::RCP<std::vector<Real> > WORKp = Teuchos::rcp(new std::vector<Real>(4*N,0.0));

    // Column-stacked matrix of eigenvectors needed for weights
    Teuchos::RCP<std::vector<Real> > Zp = Teuchos::rcp(new std::vector<Real>(N*N,0));

    // D = a
    std::copy(ap->begin(),ap->end(),Dp->begin());
     
    for(int i=0;i<N-1;++i) {
        (*Ep)[i] = sqrt((*bp)[i+1]);  
    }

    // Eigenvalue Decomposition 
    lapack.STEQR(COMPZ,N,&(*Dp)[0],&(*Ep)[0],&(*Zp)[0],LDZ,&(*WORKp)[0],&INFO);

    for(int i=0;i<N;++i){
        (*xp)[i] = (*Dp)[i];
        (*wp)[i] = (*bp)[0]*pow((*Zp)[N*i],2);
    } 
}

Stokhos::MonoProjPCEBasis<ordinal_type, value_type>::
MonoProjPCEBasis(
   ordinal_type p,
   const Stokhos::OrthogPolyApprox<ordinal_type, value_type>& pce,
   const Stokhos::Quadrature<ordinal_type, value_type>& quad,
   const Stokhos::Sparse3Tensor<ordinal_type, value_type>& Cijk,
   bool limit_integration_order_) :
  RecurrenceBasis<ordinal_type, value_type>("Monomial Projection", p, true),
  limit_integration_order(limit_integration_order_),
  pce_sz(pce.basis()->size()),
  pce_norms(pce.basis()->norm_squared()),
  a(pce_sz), 
  b(pce_sz),
  basis_vecs(pce_sz, p+1),
  new_pce(p+1)
{
  // If the original basis is normalized, we can use the standard QR
  // factorization.  For simplicity, we renormalize the PCE coefficients
  // for a normalized basis
  Stokhos::OrthogPolyApprox<ordinal_type, value_type> normalized_pce(pce);
  for (ordinal_type i=0; i<pce_sz; i++) {
    pce_norms[i] = std::sqrt(pce_norms[i]);
    normalized_pce[i] *= pce_norms[i];
  }

  // Evaluate PCE at quad points
  ordinal_type nqp = quad.size();
  Teuchos::Array<value_type> pce_vals(nqp);
  const Teuchos::Array<value_type>& weights = quad.getQuadWeights();
  const Teuchos::Array< Teuchos::Array<value_type> >& quad_points =
    quad.getQuadPoints();
  const Teuchos::Array< Teuchos::Array<value_type> >& basis_values =
    quad.getBasisAtQuadPoints();
  for (ordinal_type i=0; i<nqp; i++) {
    pce_vals[i] = normalized_pce.evaluate(quad_points[i], basis_values[i]);
  }

  // Form Kylov matrix up to order pce_sz
  matrix_type K(pce_sz, pce_sz);

  // Compute matrix
  matrix_type A(pce_sz, pce_sz);
  typedef Stokhos::Sparse3Tensor<ordinal_type, value_type> Cijk_type;
  for (typename Cijk_type::k_iterator k_it = Cijk.k_begin();
       k_it != Cijk.k_end(); ++k_it) {
    ordinal_type k = index(k_it);
    for (typename Cijk_type::kj_iterator j_it = Cijk.j_begin(k_it); 
	 j_it != Cijk.j_end(k_it); ++j_it) {
      ordinal_type j = index(j_it);
      value_type val = 0;
      for (typename Cijk_type::kji_iterator i_it = Cijk.i_begin(j_it);
	   i_it != Cijk.i_end(j_it); ++i_it) {
	ordinal_type i = index(i_it);
	value_type c = value(i_it) / (pce_norms[j]*pce_norms[k]);
	val += pce[i]*c;
      }
      A(k,j) = val;
    }
  }

  // Each column i is given by projection of the i-th order monomial 
  // onto original basis
  vector_type u0 = Teuchos::getCol(Teuchos::View, K, 0);
  u0(0) = 1.0;
  for (ordinal_type i=1; i<pce_sz; i++)
    u0(i) = 0.0;
  for (ordinal_type k=1; k<pce_sz; k++) {
    vector_type u = Teuchos::getCol(Teuchos::View, K, k);
    vector_type up = Teuchos::getCol(Teuchos::View, K, k-1);
    u.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, A, up, 0.0);
  }
  /*
  for (ordinal_type j=0; j<pce_sz; j++) {
    for (ordinal_type i=0; i<pce_sz; i++) {
      value_type val = 0.0;
      for (ordinal_type k=0; k<nqp; k++)
	val += weights[k]*std::pow(pce_vals[k],j)*basis_values[k][i];
      K(i,j) = val;
    }
  }
  */

  std::cout << K << std::endl << std::endl;

  // Compute QR factorization of K
  ordinal_type ws_size, info;
  value_type ws_size_query;
  Teuchos::Array<value_type> tau(pce_sz);
  Teuchos::LAPACK<ordinal_type,value_type> lapack;
  lapack.GEQRF(pce_sz, pce_sz, K.values(), K.stride(), &tau[0], 
	       &ws_size_query, -1, &info);
  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error, 
		     "GEQRF returned value " << info);
  ws_size = static_cast<ordinal_type>(ws_size_query);
  Teuchos::Array<value_type> work(ws_size);
  lapack.GEQRF(pce_sz, pce_sz, K.values(), K.stride(), &tau[0], 
	       &work[0], ws_size, &info);
  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error, 
		     "GEQRF returned value " << info);
  
//.........这里部分代码省略.........

//.........这里部分代码省略.........
              
              // get normal direction, this has the edge weight factored into it already
              CellTools<double>::getReferenceSideNormal(side_normal , 
                                                        (int)side_cur,cell );

              // v.n at cub points on side
              vectorBasis->getValues(value_of_v_basis_at_cub_points_side ,
                                    cub_points_side_refcell ,
                                    OPERATOR_VALUE );

              for (int i=0;i<numVectorFields;i++) {
                for (int j=0;j<numCubPointsSide;j++) {
                  n_of_v_basis_at_cub_points_side(i,j) = 0.0;
                  for (int k=0;k<cellDim;k++) {
                    n_of_v_basis_at_cub_points_side(i,j) += side_normal(k) * 
                      value_of_v_basis_at_cub_points_side(i,j,k);
                  }
                } 
              }
              
              cub_weights_side.resize(1,numCubPointsSide);
              FunctionSpaceTools::multiplyMeasure<double>(w_n_of_v_basis_at_cub_points_side,
                                                          cub_weights_side,
                                                          n_of_v_basis_at_cub_points_side);
              cub_weights_side.resize(numCubPointsSide);
              
              FunctionSpaceTools::integrate<double>(rhs_vector_vec,
                                                    diri_data_at_cub_points_side,
                                                    w_n_of_v_basis_at_cub_points_side,
                                                    COMP_BLAS,
                                                    false);

              for (int i=0;i<numVectorFields;i++) {
                rhs_and_soln_vec(0,i) -= rhs_vector_vec(0,i);
              }
              
            }
            
            // solve linear system
            int info = 0;
            Teuchos::LAPACK<int, double> solver;
            solver.GESV(numTotalFields, 1, &fe_matrix(0,0,0), numTotalFields, &ipiv(0), &rhs_and_soln_vec(0,0), 
                        numTotalFields, &info);
            
            // compute interpolant; the scalar entries are last
            scalarBasis->getValues(value_of_s_basis_at_interp_points,
                                  interp_points_ref,
                                  OPERATOR_VALUE);
            for (int pt=0;pt<numInterpPoints;pt++) {
              interpolant(0,pt)=0.0;
              for (int i=0;i<numScalarFields;i++) {
                interpolant(0,pt) += rhs_and_soln_vec(0,numVectorFields+i)
                  * value_of_s_basis_at_interp_points(i,pt);
              }
            }
            
            interp_points_ref.resize(1,numInterpPoints,cellDim);
            // get exact solution for comparison
            FieldContainer<double> exact_solution(1,numInterpPoints);
            u_exact( exact_solution , interp_points_ref , x_order, y_order, z_order);
            interp_points_ref.resize(numInterpPoints,cellDim);

            RealSpaceTools<double>::add(interpolant,exact_solution);
            
            double nrm= RealSpaceTools<double>::vectorNorm(&interpolant(0,0),interpolant.dimension(1), NORM_TWO);

            *outStream << "\nNorm-2 error between scalar components of exact solution of order ("
                       << x_order << ", " << y_order << ", " << z_order
                       << ") and finite element interpolant of order " << basis_order << ": "
                       << nrm << "\n";

            if (nrm > zero) {
              *outStream << "\n\nPatch test failed for solution polynomial order ("
                         << x_order << ", " << y_order << ", " << z_order << ") and basis order (scalar, vector)  ("
                         << basis_order << ", " << basis_order+1 << ")\n\n";
              errorFlag++;
            }
            
          }
        }
      }
    }
    
  }
  
  catch (std::logic_error err) {
    *outStream << err.what() << "\n\n";
    errorFlag = -1000;
  };
  
  if (errorFlag != 0)
    std::cout << "End Result: TEST FAILED\n";
  else
    std::cout << "End Result: TEST PASSED\n";
  
  // reset format state of std::cout
  std::cout.copyfmt(oldFormatState);
  Kokkos::finalize();
  return errorFlag;
}

//.........这里部分代码省略.........
    if (nIter == accelerationStartIteration) {
      // Initialize
      nStore = 0;
      rMat.shape(0,0);
      oldPrecF->update(1.0, *precF, -1.0);
      qrAdd(*oldPrecF);
      xMat[0]->update(1.0, solnPtr->getX(), -1.0, oldSolnPtr->getX(), 0.0);
    }
    else if (nIter > accelerationStartIteration) {
      if (nStore == storeParam) {
        Teuchos::RCP<NOX::Abstract::Vector> tempPtr = xMat[0];
        for (int ii = 0; ii<nStore-1; ii++)
          xMat[ii] = xMat[ii+1];
        xMat[nStore-1] = tempPtr;
        qrDelete();
      }
      oldPrecF->update(1.0, *precF, -1.0);
      qrAdd(*oldPrecF);
      xMat[nStore-1]->update(1.0, solnPtr->getX(), -1.0, oldSolnPtr->getX(), 0.0);
    }
  }

  // Reorthogonalize 
  if ( (nStore > 1) && (orthoFrequency > 0) )
    if (nIter % orthoFrequency == 0)
      reorthogonalize();

  // Copy current soln to the old soln.
  *oldSolnPtr = *solnPtr;
  *oldPrecF = *precF;

  // Adjust for condition number
  if (nStore > 0) {
    Teuchos::LAPACK<int,double> lapack;
    char normType = '1';
    double invCondNum = 0.0;
    int info = 0;
    if ( WORK.size() < static_cast<std::size_t>(4*nStore) )
      WORK.resize(4*nStore);
    if (IWORK.size() < static_cast<std::size_t>(nStore))
      IWORK.resize(nStore);
    lapack.GECON(normType,nStore,rMat.values(),nStore,rMat.normOne(),&invCondNum,&WORK[0],&IWORK[0],&info);
    if (utilsPtr->isPrintType(Utils::Details))
      utilsPtr->out() << "    R condition number estimate ("<< nStore << ") = " << 1.0/invCondNum << std::endl;

    if (adjustForConditionNumber) {
      while ( (1.0/invCondNum > dropTolerance) && (nStore > 1)  ) {
        Teuchos::RCP<NOX::Abstract::Vector> tempPtr = xMat[0];
        for (int ii = 0; ii<nStore-1; ii++)
          xMat[ii] = xMat[ii+1];
        xMat[nStore-1] = tempPtr;
        qrDelete();
        lapack.GECON(normType,nStore,rMat.values(),nStore,rMat.normOne(),&invCondNum,&WORK[0],&IWORK[0],&info);
        if (utilsPtr->isPrintType(Utils::Details))
          utilsPtr->out() << "    Adjusted R condition number estimate ("<< nStore << ") = " << 1.0/invCondNum << std::endl;
      }
    }
  }

  // Solve the least-squares problem.
  Teuchos::SerialDenseMatrix<int,double> gamma(nStore,1), RHS(nStore,1), Rgamma(nStore,1);
  for (int ii = 0; ii<nStore; ii++)
    RHS(ii,0) = precF->innerProduct( *(qMat[ii]) );

  //Back-solve for gamma
  for (int ii = nStore-1; ii>=0; ii--) {