Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_DiagEpetraOp.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
44#include "Epetra_config.h"
45#include "EpetraExt_BlockMultiVector.h"
46#include "EpetraExt_MatrixMatrix.h"
48
50 const Teuchos::RCP<const Epetra_Map>& domain_base_map_,
51 const Teuchos::RCP<const Epetra_Map>& range_base_map_,
52 const Teuchos::RCP<const Epetra_Map>& domain_sg_map_,
53 const Teuchos::RCP<const Epetra_Map>& range_sg_map_,
54 const Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> >& sg_basis_,
55 const Teuchos::RCP<const Stokhos::Sparse3Tensor<int,double> >& Cijk_,
56 const Teuchos::RCP<Stokhos::EpetraOperatorOrthogPoly >& ops_)
57 : label("Stokhos Diagonal Operator"),
58 domain_base_map(domain_base_map_),
59 range_base_map(range_base_map_),
60 domain_sg_map(domain_sg_map_),
61 range_sg_map(range_sg_map_),
62 sg_basis(sg_basis_),
63 Cijk(Cijk_),
64 block_ops(ops_),
65 useTranspose(false),
66 expansion_size(sg_basis->size()),
67 num_blocks(block_ops->size()),
68 input_block(expansion_size),
69 result_block(expansion_size),
70 tmp(),
71 tmp_trans()
72{
73}
74
76{
77}
78
79void
81 const Teuchos::RCP<Stokhos::EpetraOperatorOrthogPoly >& ops)
82{
83 block_ops = ops;
84}
85
86Teuchos::RCP<const Stokhos::EpetraOperatorOrthogPoly >
88{
89 return block_ops;
90}
91
92Teuchos::RCP<Stokhos::EpetraOperatorOrthogPoly >
94{
95 return block_ops;
96}
97
98int
100{
101 useTranspose = UseTranspose;
102 for (int i=0; i<num_blocks; i++)
103 (*block_ops)[i].SetUseTranspose(useTranspose);
104
105 return 0;
106}
107
108int
109Stokhos::DiagEpetraOp::Apply(std::vector< Teuchos::RCP< const Epetra_CrsMatrix> >& sg_J_all, std::vector< Teuchos::RCP< Epetra_CrsMatrix> >& sg_Kkk_all) const
110{
111
112// std::vector< Teuchos::RCP< Epetra_CrsMatrix> > sg_Kkk_all ;
113// std::vector< Teuchos::RCP< const Epetra_CrsMatrix> > sg_J_all;
114
115 for (int i=0;i<num_blocks+1;i++) {
116 Teuchos::RCP<Epetra_CrsMatrix> sg_J_poly_Crs =
117 Teuchos::rcp_dynamic_cast< Epetra_CrsMatrix>((*block_ops).getCoeffPtr(i),true);
118 sg_J_all.push_back(sg_J_poly_Crs);
119 }
120
121/* Teuchos::RCP<Epetra_CrsMatrix> sg_J_poly_Crs =
122 Teuchos::rcp_dynamic_cast< Epetra_CrsMatrix>((*sg_J_poly).getCoeffPtr(0),true);
123
124 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::is_null(sg_J_poly_Crs), std::runtime_error,
125 "Dynamic cast of sg_J_poly failed!");
126*/
127
128 for(int k=0;k<expansion_size;k++) {
129 Teuchos::RCP<Epetra_CrsMatrix> Kkk =
130 Teuchos::rcp(new Epetra_CrsMatrix(*(sg_J_all[0])));
131 sg_Kkk_all.push_back(Kkk);
132 sg_Kkk_all[k]->PutScalar(0.0);
133 }
134
135 //Compute diagonal blocks of SG matrix
136 for (int k=0; k<expansion_size; k++) {
137 int nj = Cijk->num_j(k);
138 const Teuchos::Array<int>& j_indices = Cijk->Jindices(k);
139 for (int jj=0; jj<nj; jj++) {
140 int j = j_indices[jj];
141 if (j==k) {
142 const Teuchos::Array<double>& cijk_values = Cijk->values(k,jj);
143 const Teuchos::Array<int>& i_indices = Cijk->Iindices(k,jj);
144 int ni = i_indices.size();
145 for (int ii=0; ii<ni; ii++) {
146 int i = i_indices[ii];
147 if (i<num_blocks+1) {
148 double cikk = cijk_values[ii]; // C(i,j,k)
149 EpetraExt::MatrixMatrix::Add((*sg_J_all[i]), false, cikk, *(sg_Kkk_all[k]), 1.0);
150 }
151 }
152 }
153 }
154 }
155
156 /* // Apply block SG operator via
157 // w_i =
158 // \sum_{j=0}^P \sum_{k=0}^L J_k v_j < \psi_i \psi_j \psi_k > / <\psi_i^2>
159 // for i=0,...,P where P = expansion_size, L = num_blocks, w_j is the jth
160 // input block, w_i is the ith result block, and J_k is the kth block operator
161 const Teuchos::Array<double>& norms = sg_basis->norm_squared();
162 for (int k=0; k<num_blocks; k++) {
163 int nj = Cijk->num_j(k);
164 const Teuchos::Array<int>& j_indices = Cijk->Jindices(k);
165 Teuchos::Array<double*> j_ptr(nj*m);
166 Teuchos::Array<int> mj_indices(nj*m);
167 for (int l=0; l<nj; l++) {
168 for (int mm=0; mm<m; mm++) {
169 j_ptr[l*m+mm] = input_block[j_indices[l]]->Values()+mm*N;
170 mj_indices[l*m+mm] = j_indices[l]*m+mm;
171 }
172 }
173 Epetra_MultiVector input_tmp(View, *input_base_map, &j_ptr[0], nj*m);
174 Epetra_MultiVector result_tmp(View, *tmp_result, &mj_indices[0], nj*m);
175 (*block_ops)[k].Apply(input_tmp, result_tmp);
176 for (int l=0; l<nj; l++) {
177 const Teuchos::Array<int>& i_indices = Cijk->Iindices(k,l);
178 const Teuchos::Array<double>& c_values = Cijk->values(k,l);
179 for (int i=0; i<i_indices.size(); i++) {
180 int ii = i_indices[i];
181 for (int mm=0; mm<m; mm++)
182 (*result_block[ii])(mm)->Update(c_values[i]/norms[ii],
183 *result_tmp(l*m+mm), 1.0);
184 }
185 }
186 }
187
188 // Destroy blocks
189 for (int i=0; i<expansion_size; i++) {
190 input_block[i] = Teuchos::null;
191 result_block[i] = Teuchos::null;
192 }
193
194 if (made_copy)
195 delete input;
196 */
197 return 0;
198}
199
200int
202 Epetra_MultiVector& Result) const
203{
204 throw "DiagEpetraOp::ApplyInverse not defined!";
205 return -1;
206}
207
208double
210{
211 return 1.0;
212}
213
214
215const char*
217{
218 return const_cast<char*>(label.c_str());
219}
220
221bool
223{
224 return useTranspose;
225}
226
227bool
229{
230 return false;
231}
232
233const Epetra_Comm &
235{
236 return domain_base_map->Comm();
237}
238const Epetra_Map&
240{
241 return *domain_sg_map;
242}
243
244const Epetra_Map&
246{
247 return *range_sg_map;
248}
virtual void reset(const Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > &ops)
Reset operator blocks.
virtual bool UseTranspose() const
Returns the current UseTranspose setting.
virtual Teuchos::RCP< const Stokhos::EpetraOperatorOrthogPoly > getOperatorBlocks() const
Get operator blocks.
DiagEpetraOp(const Teuchos::RCP< const Epetra_Map > &domain_base_map_, const Teuchos::RCP< const Epetra_Map > &range_base_map_, const Teuchos::RCP< const Epetra_Map > &domain_sg_map_, const Teuchos::RCP< const Epetra_Map > &range_sg_map_, const Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > &sg_basis, const Teuchos::RCP< const Stokhos::Sparse3Tensor< int, double > > &Cijk, const Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > &ops)
Constructor.
virtual int SetUseTranspose(bool UseTranspose)
Set to true if the transpose of the operator is requested.
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this matrix operator.
virtual const char * Label() const
Returns a character string describing the operator.
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of the inverse of the operator applied to a Epetra_MultiVector Input in Result as ...
virtual ~DiagEpetraOp()
Destructor.
virtual double NormInf() const
Returns an approximate infinity norm of the operator matrix.
virtual int Apply(std::vector< Teuchos::RCP< const Epetra_CrsMatrix > > &sg_J_all, std::vector< Teuchos::RCP< Epetra_CrsMatrix > > &sg_Kkk_all) const
Returns Diagonal blocks of SG matrix when PC coefficients of the SG matrix are given.
virtual const Epetra_Comm & Comm() const
Returns a reference to the Epetra_Comm communicator associated with this operator.
virtual bool HasNormInf() const
Returns true if the this object can provide an approximate Inf-norm, false otherwise.
virtual const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this matrix operator.
Abstract base class for multivariate orthogonal polynomials.
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.