Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_InterlacedOperator.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
46#include "Epetra_Map.h"
47
48#include "EpetraExt_BlockUtility.h"
49
52 const Teuchos::RCP<const EpetraExt::MultiComm>& sg_comm_,
53 const Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> >& sg_basis_,
54 const Teuchos::RCP<const Stokhos::EpetraSparse3Tensor>& epetraCijk_,
55 const Teuchos::RCP<const Epetra_CrsGraph>& determ_graph,
56 const Teuchos::RCP<Teuchos::ParameterList>& params) :
57 EpetraExt::BlockCrsMatrix(*(epetraCijk_->getStochasticGraph()),
58 *determ_graph,
59 *sg_comm_),
60 sg_comm(sg_comm_),
61 sg_basis(sg_basis_),
62 epetraCijk(epetraCijk_),
63 Cijk(epetraCijk->getParallelCijk()),
64 block_ops(),
65 scale_op(true),
66 include_mean(true),
67 only_use_linear(false),
68 determOffset_(EpetraExt::BlockUtility::CalculateOffset(determ_graph->RowMap()))
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 block_ops = ops;
86
87 // Zero out matrix
88 this->PutScalar(0.0);
89
90 // Compute loop bounds
91 Cijk_type::k_iterator k_begin = Cijk->k_begin();
92 Cijk_type::k_iterator k_end = Cijk->k_end();
93 if (!include_mean && index(k_begin) == 0)
94 ++k_begin;
95 if (only_use_linear) {
96 int dim = sg_basis->dimension();
97 k_end = Cijk->find_k(dim+1);
98 }
99
100 // Assemble matrix
101 const Teuchos::Array<double>& norms = sg_basis->norm_squared();
102 for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
103 int k = index(k_it);
104 Teuchos::RCP<Epetra_RowMatrix> block =
105 Teuchos::rcp_dynamic_cast<Epetra_RowMatrix>(block_ops->getCoeffPtr(k),
106 true);
107 for (Cijk_type::kj_iterator j_it = Cijk->j_begin(k_it);
108 j_it != Cijk->j_end(k_it); ++j_it) {
109 int j = epetraCijk->GCID(index(j_it));
110 for (Cijk_type::kji_iterator i_it = Cijk->i_begin(j_it);
111 i_it != Cijk->i_end(j_it); ++i_it) {
112 int i = epetraCijk->GRID(index(i_it));
113 double c = value(i_it);
114 if (scale_op)
115 c /= norms[i];
116 this->SumIntoGlobalBlock_Deterministic(c, *block, i, j);
117 }
118 }
119 }
120}
121
122Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly >
125{
126 return block_ops;
127}
128
129Teuchos::RCP<const Stokhos::EpetraOperatorOrthogPoly >
131getSGPolynomial() const
132{
133 return block_ops;
134}
135
137SumIntoGlobalBlock_Deterministic(double alpha,const Epetra_RowMatrix & determBlock,int Row,int Col)
138{
139 const Epetra_Map & determMap = determBlock.RowMatrixRowMap();
140 const Epetra_Map & determColMap = determBlock.RowMatrixColMap();
141
142 // This routine copies entries of a BaseMatrix into big BlockCrsMatrix
143 // It performs the following operation on the global IDs row-by-row
144 // this->val[i+rowOffset][j+ColOffset] = BaseMatrix.val[i][j]
145
146 int MaxIndices = determBlock.MaxNumEntries();
147 std::vector<int> Indices(MaxIndices);
148 std::vector<double> Values(MaxIndices);
149 int NumIndices;
150 int ierr=0;
151
152 for (int i=0; i<determMap.NumMyElements(); i++) {
153 determBlock.ExtractMyRowCopy( i, MaxIndices, NumIndices, &Values[0], &Indices[0] );
154
155 // Convert to BlockMatrix Global numbering scheme
156 for( int l = 0; l < NumIndices; ++l ) {
157 Indices[l] = Col + COffset_*determColMap.GID(Indices[l]);
158 Values[l] *= alpha;
159 }
160
161 int BaseRow = determMap.GID(i);
162 ierr = this->SumIntoGlobalValues(ROffset_*BaseRow + Row, NumIndices, &Values[0], &Indices[0]);
163 if (ierr != 0) std::cout << "WARNING InterlacedOperator::SumIntoBlock_Deterministic SumIntoGlobalValues err = " << ierr <<
164 "\n\t Row " << ROffset_*BaseRow + Row << ", Col start " << Indices[0] << std::endl;
165
166 }
167}
int GID(int LID) const
int NumMyElements() const
virtual const Epetra_Map & RowMatrixColMap() const=0
virtual const Epetra_Map & RowMatrixRowMap() const=0
virtual int MaxNumEntries() const=0
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const=0
CRS matrix of dense blocks.
bool scale_op
Flag indicating whether operator be scaled with <\psi_i^2>
void SumIntoGlobalBlock_Deterministic(double alpha, const Epetra_RowMatrix &determBlock, int Row, int Col)
Sum into global matrix.
virtual void setupOperator(const Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > &poly)
Setup operator.
virtual Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > getSGPolynomial()
Get SG polynomial.
InterlacedOperator(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< Teuchos::ParameterList > &params)
Constructor.
bool only_use_linear
Flag indicating whether to only use linear terms.
bool include_mean
Flag indicating whether to include mean term.
Abstract base class for multivariate orthogonal polynomials.
Bi-directional iterator for traversing a sparse array.