Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_FullyAssembledOperator.cpp
Go to the documentation of this file.
1// $Id$
2// $Source$
3// @HEADER
4// ***********************************************************************
5//
6// Stokhos Package
7// Copyright (2009) Sandia Corporation
8//
9// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
10// license for use of this work by or on behalf of the U.S. Government.
11//
12// Redistribution and use in source and binary forms, with or without
13// modification, are permitted provided that the following conditions are
14// met:
15//
16// 1. Redistributions of source code must retain the above copyright
17// notice, this list of conditions and the following disclaimer.
18//
19// 2. Redistributions in binary form must reproduce the above copyright
20// notice, this list of conditions and the following disclaimer in the
21// documentation and/or other materials provided with the distribution.
22//
23// 3. Neither the name of the Corporation nor the names of the
24// contributors may be used to endorse or promote products derived from
25// this software without specific prior written permission.
26//
27// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38//
39// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
40//
41// ***********************************************************************
42// @HEADER
43
45#include "Teuchos_TimeMonitor.hpp"
46
49 const Teuchos::RCP<const EpetraExt::MultiComm>& sg_comm_,
50 const Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> >& sg_basis_,
51 const Teuchos::RCP<const Stokhos::EpetraSparse3Tensor>& epetraCijk_,
52 const Teuchos::RCP<const Epetra_CrsGraph>& base_graph,
53 const Teuchos::RCP<const Epetra_Map>& domain_sg_map_,
54 const Teuchos::RCP<const Epetra_Map>& range_sg_map_,
55 const Teuchos::RCP<Teuchos::ParameterList>& params) :
56 EpetraExt::BlockCrsMatrix(*base_graph,
57 *(epetraCijk_->getStochasticGraph()),
58 *sg_comm_),
59 sg_comm(sg_comm_),
60 sg_basis(sg_basis_),
61 epetraCijk(epetraCijk_),
62 domain_sg_map(domain_sg_map_),
63 range_sg_map(range_sg_map_),
64 Cijk(epetraCijk->getParallelCijk()),
65 block_ops(),
66 scale_op(true),
67 include_mean(true),
68 only_use_linear(false)
69{
70 scale_op = params->get("Scale Operator by Inverse Basis Norms", true);
71 include_mean = params->get("Include Mean", true);
72 only_use_linear = params->get("Only Use Linear Terms", false);
73}
74
77{
78}
79
80void
83 const Teuchos::RCP<Stokhos::EpetraOperatorOrthogPoly >& ops)
84{
85#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
86 TEUCHOS_FUNC_TIME_MONITOR("Stokhos: SG Fully Assembled Operator Assembly");
87#endif
88
89 block_ops = ops;
90
91 // Zero out matrix
92 this->PutScalar(0.0);
93
94 // Compute loop bounds
95 Cijk_type::k_iterator k_begin = Cijk->k_begin();
96 Cijk_type::k_iterator k_end = Cijk->k_end();
97 if (!include_mean && index(k_begin) == 0)
98 ++k_begin;
99 if (only_use_linear) {
100 int dim = sg_basis->dimension();
101 k_end = Cijk->find_k(dim+1);
102 }
103
104 // Assemble matrix
105 const Teuchos::Array<double>& norms = sg_basis->norm_squared();
106 for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
107 int k = index(k_it);
108 Teuchos::RCP<Epetra_RowMatrix> block =
109 Teuchos::rcp_dynamic_cast<Epetra_RowMatrix>(block_ops->getCoeffPtr(k),
110 true);
111 for (Cijk_type::kj_iterator j_it = Cijk->j_begin(k_it);
112 j_it != Cijk->j_end(k_it); ++j_it) {
113 int j = epetraCijk->GCID(index(j_it));
114 for (Cijk_type::kji_iterator i_it = Cijk->i_begin(j_it);
115 i_it != Cijk->i_end(j_it); ++i_it) {
116 int i = epetraCijk->GRID(index(i_it));
117 double c = value(i_it);
118 if (scale_op)
119 c /= norms[i];
120 this->SumIntoGlobalBlock(c, *block, i, j);
121 }
122 }
123 }
124
125 this->FillComplete(*domain_sg_map, *range_sg_map);
126}
127
128Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly >
131{
132 return block_ops;
133}
134
135Teuchos::RCP<const Stokhos::EpetraOperatorOrthogPoly >
137getSGPolynomial() const
138{
139 return block_ops;
140}
141
142int
144Apply(const Epetra_MultiVector& Input, Epetra_MultiVector& Result) const
145{
146#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
147 TEUCHOS_FUNC_TIME_MONITOR("Stokhos: SG Operator Apply");
148#endif
149 return Epetra_CrsMatrix::Apply(Input, Result);
150}
int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
CRS matrix of dense blocks.
bool scale_op
Flag indicating whether operator be scaled with <\psi_i^2>
virtual int Apply(const Epetra_MultiVector &Input, Epetra_MultiVector &Result) const
Wrap Apply() to add a timer.
bool include_mean
Flag indicating whether to include mean term.
virtual Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > getSGPolynomial()
Get SG polynomial.
bool only_use_linear
Flag indicating whether to only use linear terms.
virtual void setupOperator(const Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > &poly)
Setup operator.
FullyAssembledOperator(const Teuchos::RCP< const EpetraExt::MultiComm > &sg_comm, const Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > &sg_basis, const Teuchos::RCP< const Stokhos::EpetraSparse3Tensor > &epetraCijk, const Teuchos::RCP< const Epetra_CrsGraph > &base_graph, const Teuchos::RCP< const Epetra_Map > &domain_sg_map, const Teuchos::RCP< const Epetra_Map > &range_sg_map, const Teuchos::RCP< Teuchos::ParameterList > &params)
Constructor.
Abstract base class for multivariate orthogonal polynomials.
Bi-directional iterator for traversing a sparse array.