本文整理汇总了C++中teuchos::LAPACK类 的典型用法代码示例。如果您正苦于以下问题:C++ LAPACK类的具体用法?C++ LAPACK怎么用?C++ LAPACK使用的例子?那么恭喜您, 这里精选的类代码示例或许可以为您提供帮助。
在下文中一共展示了LAPACK类 的20个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于我们的系统推荐出更棒的C++代码示例。
示例1: main
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;
}
开发者ID:00liujj, 项目名称:trilinos, 代码行数:32, 代码来源:ex3.cpp
示例2: solve_state
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());
}
}
开发者ID:Russell-Jones-OxPhys, 项目名称:Trilinos, 代码行数:28, 代码来源:example_01.cpp
示例3: solve_adjoint_sensitivity
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());
}
}
开发者ID:Russell-Jones-OxPhys, 项目名称:Trilinos, 代码行数:31, 代码来源:example_01.cpp
示例4:
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;
}
开发者ID:agrippa, 项目名称:Trilinos, 代码行数:33, 代码来源:Tacho_Chol_Upper_ExternalLapack.hpp
示例5: solnT
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();
}
开发者ID:achellies, 项目名称:maxwell, 代码行数:15, 代码来源:MxDimMatrix.hpp
示例6: levecs
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;
}
开发者ID:achellies, 项目名称:maxwell, 代码行数:45, 代码来源:MxDimMatrix.hpp
示例7: BcRow
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);
}
}
开发者ID:, 项目名称:, 代码行数:44, 代码来源:
示例8: copy
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, ©(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, ©(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;
}
开发者ID:achellies, 项目名称:maxwell, 代码行数:20, 代码来源:MxDimMatrix.hpp
示例9: updateGuess
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;
}
开发者ID:MDBrothers, 项目名称:NewtonRaphsonExamples, 代码行数:23, 代码来源:complex_step.cpp
示例10: U
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.");
}
开发者ID:abhishek4747, 项目名称:trilinos, 代码行数:23, 代码来源:Ifpack2_DenseContainer_def.hpp
示例11:
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;
}
开发者ID:KineticTheory, 项目名称:Trilinos, 代码行数:37, 代码来源:LOCA_EigenvalueSort_Strategies.C
示例12: gauss
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);
}
}
开发者ID:rainiscold, 项目名称:trilinos, 代码行数:42, 代码来源:OrthogonalPolynomials.hpp
示例13: K
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);
//.........这里部分代码省略.........
开发者ID:00liujj, 项目名称:trilinos, 代码行数:101, 代码来源:Stokhos_MonoProjPCEBasisImp.hpp
示例14: main
//.........这里部分代码省略.........
// 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;
}
开发者ID:KineticTheory, 项目名称:Trilinos, 代码行数:101, 代码来源:test_02.cpp
示例15: if
//.........这里部分代码省略.........
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--) {
开发者ID:OpenModelica, 项目名称:OMCompiler-3rdParty, 代码行数:67, 代码来源:NOX_Solver_AndersonAcceleration.C
示例16: Compute
//.........这里部分代码省略.........
std::vector< std::complex<double> > cpts;
for( int jj=0; jj<ny; jj++ ) {
for( int ii=0; ii<nx; ii++ ) {
std::complex<double> cpt(xs[ii],ys[jj]);
cpts.push_back(cpt);
}
}
cpts.push_back(zero);
#ifdef HAVE_TEUCHOS_COMPLEX
const std::complex<double> one(1.0,0.0);
// Construct overdetermined Vandermonde matrix
Teuchos::SerialDenseMatrix<int, std::complex<double> > Vmatrix(cpts.size(),PolyDegree_+1);
Vmatrix.putScalar(zero);
for (int jj = 0; jj <= PolyDegree_; ++jj) {
for (int ii = 0; ii < static_cast<int> (cpts.size ()) - 1; ++ii) {
if (jj > 0) {
Vmatrix(ii,jj) = pow(cpts[ii],jj);
}
else {
Vmatrix(ii,jj) = one;
}
}
}
Vmatrix(cpts.size()-1,0)=one;
// Right hand side: all zero except last entry
Teuchos::SerialDenseMatrix< int,std::complex<double> > RHS(cpts.size(),1);
RHS.putScalar(zero);
RHS(cpts.size()-1,0)=one;
// Solve least squares problem using LAPACK
Teuchos::LAPACK< int, std::complex<double> > lapack;
const int N = Vmatrix.numCols();
Teuchos::Array<double> singularValues(N);
Teuchos::Array<double> rwork(1);
rwork.resize (std::max (1, 5 * N));
std::complex<double> lworkScalar(1.0,0.0);
int info = 0;
lapack.GELS('N', Vmatrix.numRows(), Vmatrix.numCols(), RHS.numCols(),
Vmatrix.values(), Vmatrix.numRows(), RHS.values(), RHS.numRows(),
&lworkScalar, -1, &info);
TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
"_GELSS workspace query returned INFO = "
<< info << " != 0.");
const int lwork = static_cast<int> (real(lworkScalar));
TEUCHOS_TEST_FOR_EXCEPTION(lwork < 0, std::logic_error,
"_GELSS workspace query returned LWORK = "
<< lwork << " < 0.");
// Allocate workspace. Size > 0 means &work[0] makes sense.
Teuchos::Array< std::complex<double> > work (std::max (1, lwork));
// Solve the least-squares problem.
lapack.GELS('N', Vmatrix.numRows(), Vmatrix.numCols(), RHS.numCols(),
Vmatrix.values(), Vmatrix.numRows(), RHS.values(), RHS.numRows(),
&work[0], lwork, &info);
coeff_.resize(PolyDegree_+1);
std::complex<double> c0=RHS(0,0);
for(int ii=0; ii<=PolyDegree_; ii++) {
// test that the imaginary part is nonzero
//TEUCHOS_TEST_FOR_EXCEPTION(abs(imag(RHS(ii,0))) > 1e-8, std::logic_error,
// "imaginary part of polynomial coefficients is nonzero! coeff = "
// << RHS(ii,0));
coeff_[ii]=real(RHS(ii,0)/c0);
//std::cout<<"coeff["<<ii<<"]="<<coeff_[ii]<<std::endl;
开发者ID:brian-kelley, 项目名称:Trilinos, 代码行数:67, 代码来源:Ifpack_Polynomial.cpp
示例17: main
//.........这里部分代码省略.........
rhs_at_cub_points_cell_physical,
weighted_transformed_value_of_basis_at_cub_points_cell,
COMP_BLAS);
// compute neumann b.c. contributions and adjust rhs
sideCub->getCubature(cub_points_side, cub_weights_side);
for (unsigned i=0; i<numSides; i++) {
// compute geometric cell information
CellTools<double>::mapToReferenceSubcell(cub_points_side_refcell, cub_points_side, sideDim, (int)i, cell);
CellTools<double>::setJacobian(jacobian_side_refcell, cub_points_side_refcell, cell_nodes[pcell], cell);
CellTools<double>::setJacobianDet(jacobian_det_side_refcell, jacobian_side_refcell);
// compute weighted edge measure
FunctionSpaceTools::computeEdgeMeasure<double>(weighted_measure_side_refcell,
jacobian_side_refcell,
cub_weights_side,
i,
cell);
// tabulate values of basis functions at side cubature points, in the reference parent cell domain
basis->getValues(value_of_basis_at_cub_points_side_refcell, cub_points_side_refcell, OPERATOR_VALUE);
// transform
FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_cub_points_side_refcell,
value_of_basis_at_cub_points_side_refcell);
// multiply with weighted measure
FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points_side_refcell,
weighted_measure_side_refcell,
transformed_value_of_basis_at_cub_points_side_refcell);
// compute Neumann data
// map side cubature points in reference parent cell domain to physical space
CellTools<double>::mapToPhysicalFrame(cub_points_side_physical, cub_points_side_refcell, cell_nodes[pcell], cell);
// now compute data
neumann(neumann_data_at_cub_points_side_physical, cub_points_side_physical, jacobian_side_refcell,
cell, (int)i, x_order, y_order);
FunctionSpaceTools::integrate<double>(neumann_fields_per_side,
neumann_data_at_cub_points_side_physical,
weighted_transformed_value_of_basis_at_cub_points_side_refcell,
COMP_BLAS);
// adjust RHS
RealSpaceTools<double>::add(rhs_and_soln_vector, neumann_fields_per_side);;
}
///////////////////////////////
/////////////////////////////
// Solution of linear system:
int info = 0;
Teuchos::LAPACK<int, double> solver;
solver.GESV(numFields, 1, &fe_matrix[0], numFields, &ipiv(0), &rhs_and_soln_vector[0], numFields, &info);
/////////////////////////////
////////////////////////
// Building interpolant:
// evaluate basis at interpolation points
basis->getValues(value_of_basis_at_interp_points, interp_points_ref, OPERATOR_VALUE);
// transform values of basis functions
FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_interp_points,
value_of_basis_at_interp_points);
FunctionSpaceTools::evaluate<double>(interpolant, rhs_and_soln_vector, transformed_value_of_basis_at_interp_points);
////////////////////////
/******************* END COMPUTATION ***********************/
RealSpaceTools<double>::subtract(interpolant, exact_solution);
*outStream << "\nRelative norm-2 error between exact solution polynomial of order ("
<< x_order << ", " << y_order << ") and finite element interpolant of order " << basis_order << ": "
<< RealSpaceTools<double>::vectorNorm(&interpolant[0], interpolant.dimension(1), NORM_TWO) /
RealSpaceTools<double>::vectorNorm(&exact_solution[0], exact_solution.dimension(1), NORM_TWO) << "\n";
if (RealSpaceTools<double>::vectorNorm(&interpolant[0], interpolant.dimension(1), NORM_TWO) /
RealSpaceTools<double>::vectorNorm(&exact_solution[0], exact_solution.dimension(1), NORM_TWO) > zero) {
*outStream << "\n\nPatch test failed for solution polynomial order ("
<< x_order << ", " << y_order << ") and basis order " << basis_order << "\n\n";
errorFlag++;
}
} // end for y_order
} // end for x_order
} // end for pcell
}
// Catch unexpected errors
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);
return errorFlag;
}
开发者ID:rainiscold, 项目名称:trilinos, 代码行数:101, 代码来源:test_02.cpp
示例18: if
void RCGIter<ScalarType,MV,OP>::iterate()
{
TEUCHOS_TEST_FOR_EXCEPTION( initialized_ == false, RCGIterFailure,
"Belos::RCGIter::iterate(): RCGIter class not initialized." );
// We'll need LAPACK
Teuchos::LAPACK<int,ScalarType> lapack;
// Create convenience variables for zero and one.
ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
// Allocate memory for scalars
std::vector<int> index(1);
Teuchos::SerialDenseMatrix<int,ScalarType> pAp(1,1), rTz(1,1);
// Get the current solution std::vector.
Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
// Check that the current solution std::vector only has one column.
TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1, RCGIterFailure,
"Belos::RCGIter::iterate(): current linear system has more than one std::vector!" );
// Compute the current search dimension.
int searchDim = numBlocks_+1;
// index of iteration within current cycle
int i_ = 0;
////////////////////////////////////////////////////////////////
// iterate until the status test tells us to stop.
//
// also break if our basis is full
//
Teuchos::RCP<const MV> p_ = Teuchos::null;
Teuchos::RCP<MV> pnext_ = Teuchos::null;
while (stest_->checkStatus(this) != Passed && curDim_+1 <= searchDim) {
// Ap = A*p;
index.resize( 1 );
index[0] = i_;
p_ = MVT::CloneView( *P_, index );
lp_->applyOp( *p_, *Ap_ );
// d = p'*Ap;
MVT::MvTransMv( one, *p_, *Ap_, pAp );
(*D_)(i_,0) = pAp(0,0);
// alpha = rTz_old / pAp
(*Alpha_)(i_,0) = (*rTz_old_)(0,0) / pAp(0,0);
// Check that alpha is a positive number
TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(pAp(0,0)) <= zero, RCGIterFailure, "Belos::RCGIter::iterate(): non-positive value for p^H*A*p encountered!" );
// x = x + (alpha * p);
MVT::MvAddMv( one, *cur_soln_vec, (*Alpha_)(i_,0), *p_, *cur_soln_vec );
lp_->updateSolution();
// r = r - (alpha * Ap);
MVT::MvAddMv( one, *r_, -(*Alpha_)(i_,0), *Ap_, *r_ );
std::vector<MagnitudeType> norm(1);
MVT::MvNorm( *r_, norm );
//printf("i = %i\tnorm(r) = %e\n",i_,norm[0]);
// z = M\r
if ( lp_->getLeftPrec() != Teuchos::null ) {
lp_->applyLeftPrec( *r_, *z_ );
}
else if ( lp_->getRightPrec() != Teuchos::null ) {
lp_->applyRightPrec( *r_, *z_ );
}
else {
z_ = r_;
}
// rTz_new = r'*z;
MVT::MvTransMv( one, *r_, *z_, rTz );
// beta = rTz_new/rTz_old;
(*Beta_)(i_,0) = rTz(0,0) / (*rTz_old_)(0,0);
// rTz_old = rTz_new;
(*rTz_old_)(0,0) = rTz(0,0);
// get pointer for next p
index.resize( 1 );
index[0] = i_+1;
pnext_ = MVT::CloneViewNonConst( *P_, index );
if (existU_) {
// mu = UTAU \ (AU'*z);
Teuchos::SerialDenseMatrix<int,ScalarType> mu( Teuchos::View, *Delta_, recycleBlocks_, 1, 0, i_ );
MVT::MvTransMv( one, *AU_, *z_, mu );
char TRANS = 'N';
int info;
lapack.GETRS( TRANS, recycleBlocks_, 1, LUUTAU_->values(), LUUTAU_->stride(), &(*ipiv_)[0], mu.values(), mu.stride(), &info );
TEUCHOS_TEST_FOR_EXCEPTION(info != 0, RCGIterLAPACKFailure,
"Belos::RCGIter::solve(): LAPACK GETRS failed to compute a solution.");
// p = -(U*mu) + (beta*p) + z (in two steps)
//.........这里部分代码省略.........
开发者ID:gitter-badger, 项目名称:quinoa, 代码行数:101, 代码来源:BelosRCGIter.hpp
六六分期app的软件客服如何联系?不知道吗?加qq群【895510560】即可!标题:六六分期
阅读:17579| 2023-10-27
今天小编告诉大家如何处理win10系统火狐flash插件总是崩溃的问题,可能很多用户都不知
阅读:9427| 2022-11-06
今天小编告诉大家如何对win10系统删除桌面回收站图标进行设置,可能很多用户都不知道
阅读:8031| 2022-11-06
今天小编告诉大家如何对win10系统电脑设置节能降温的设置方法,想必大家都遇到过需要
阅读:8399| 2022-11-06
我们在使用xp系统的过程中,经常需要对xp系统无线网络安装向导设置进行设置,可能很多
阅读:8307| 2022-11-06
今天小编告诉大家如何处理win7系统玩cf老是与主机连接不稳定的问题,可能很多用户都不
阅读:9160| 2022-11-06
电脑对日常生活的重要性小编就不多说了,可是一旦碰到win7系统设置cf烟雾头的问题,很
阅读:8277| 2022-11-06
我们在日常使用电脑的时候,有的小伙伴们可能在打开应用的时候会遇见提示应用程序无法
阅读:7699| 2022-11-06
今天小编告诉大家如何对win7系统打开vcf文件进行设置,可能很多用户都不知道怎么对win
阅读:8247| 2022-11-06
今天小编告诉大家如何对win10系统s4开启USB调试模式进行设置,可能很多用户都不知道怎
阅读:7271| 2022-11-06
请发表评论