Xpetra Version of the Day
Loading...
Searching...
No Matches
Xpetra_IteratorOps.cpp
Go to the documentation of this file.
1// @HEADER
2//
3// ***********************************************************************
4//
5// Xpetra: A linear algebra interface package
6// Copyright 2012 Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact
39// Jonathan Hu (jhu@sandia.gov)
40// Andrey Prokopenko (aprokop@sandia.gov)
41// Ray Tuminaro (rstumin@sandia.gov)
42// Tobias Wiesner (tawiesn@sandia.gov)
43//
44// ***********************************************************************
45//
46// @HEADER
48
49namespace Xpetra {
50
51#if defined(HAVE_XPETRA_EPETRA) && !defined(XPETRA_EPETRA_NO_32BIT_GLOBAL_INDICES)
52 template<>
58 bool call_FillComplete_on_result,
59 bool doOptimizeStorage,
60 const std::string & label,
62 typedef double SC;
63 typedef int LO;
64 typedef int GO;
65 typedef EpetraNode NO;
66
68 "XpetraExt::MatrixMatrix::Jacobi: row map of C is not same as row map of A")
70 "XpetraExt::MatrixMatrix::Jacobi: row map of C is not same as row map of B");
73
74 bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
75
76 if (C.getRowMap()->lib() == Xpetra::UseEpetra) {
77#ifndef HAVE_XPETRA_EPETRAEXT
78 throw(Xpetra::Exceptions::RuntimeError("Xpetra::IteratorOps::Jacobi requires EpetraExt to be compiled."));
79#else
83 // FIXME
84 XPETRA_DYNAMIC_CAST(const EpetraVectorT<GO XPETRA_COMMA NO>, Dinv, epD, "Xpetra::IteratorOps::Jacobi() only accepts Xpetra::EpetraVector as input argument.");
85
86 int i = EpetraExt::MatrixMatrix::Jacobi(omega, *epD.getEpetra_Vector(), epA, epB, epC, haveMultiplyDoFillComplete);
87 if (haveMultiplyDoFillComplete) {
88 // Due to Epetra wrapper intricacies, we need to explicitly call
89 // fillComplete on Xpetra matrix here. Specifically, EpetraCrsMatrix
90 // only keeps an internal variable to check whether we are in resumed
91 // state or not, but never touches the underlying Epetra object. As
92 // such, we need to explicitly update the state of Xpetra matrix to
93 // that of Epetra one afterwords
94 C.fillComplete();
95 }
96
97 if (i != 0) {
98 std::ostringstream buf;
99 buf << i;
100 std::string msg = "EpetraExt::MatrixMatrix::Jacobi return value of " + buf.str();
101 throw(Exceptions::RuntimeError(msg));
102 }
103#endif
104 } else if (C.getRowMap()->lib() == Xpetra::UseTpetra) {
105#ifdef HAVE_XPETRA_TPETRA
106# if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
107 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
108 throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra GO=<double,int,int> enabled."));
109# else
110 const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
111 const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpB = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(B);
112 Tpetra::CrsMatrix<SC,LO,GO,NO> & tpC = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraCrs(C);
113 const RCP<Tpetra::Vector<SC,LO,GO,NO> > & tpD = toTpetra(Dinv);
114 Tpetra::MatrixMatrix::Jacobi(omega, *tpD, tpA, tpB, tpC, haveMultiplyDoFillComplete, label, params);
115# endif
116#else
117 throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra."));
118#endif
119 }
120
121 if (call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
123 ppp->set("Optimize Storage", doOptimizeStorage);
124 C.fillComplete(B.getDomainMap(), B.getRangeMap(), ppp);
125 }
126
127 // transfer striding information
128 Teuchos::RCP<Xpetra::Matrix<SC,LO,GO,NO> > rcpA = Teuchos::rcp_const_cast<Xpetra::Matrix<SC,LO,GO,NO> >(Teuchos::rcpFromRef(A));
129 Teuchos::RCP<Xpetra::Matrix<SC,LO,GO,NO> > rcpB = Teuchos::rcp_const_cast<Xpetra::Matrix<SC,LO,GO,NO> >(Teuchos::rcpFromRef(B));
130 C.CreateView("stridedMaps", rcpA, false, rcpB, false); // TODO use references instead of RCPs
131 }
132#endif
133
134#if defined(HAVE_XPETRA_EPETRA) && !defined(XPETRA_EPETRA_NO_64BIT_GLOBAL_INDICES)
135 template<>
141 bool call_FillComplete_on_result,
142 bool doOptimizeStorage,
143 const std::string & label,
145 typedef double SC;
146 typedef int LO;
147 typedef long long GO;
148 typedef EpetraNode NO;
149
151 "XpetraExt::MatrixMatrix::Jacobi: row map of C is not same as row map of A")
153 "XpetraExt::MatrixMatrix::Jacobi: row map of C is not same as row map of B");
156
157 bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
158
159 if (C.getRowMap()->lib() == Xpetra::UseEpetra) {
160#ifndef HAVE_XPETRA_EPETRAEXT
161 throw(Xpetra::Exceptions::RuntimeError("Xpetra::IteratorOps::Jacobi requires EpetraExt to be compiled."));
162#else
166 // FIXME
167 XPETRA_DYNAMIC_CAST(const EpetraVectorT<GO XPETRA_COMMA NO>, Dinv, epD, "Xpetra::IteratorOps::Jacobi() only accepts Xpetra::EpetraVector as input argument.");
168
169 int i = EpetraExt::MatrixMatrix::Jacobi(omega, *epD.getEpetra_Vector(), epA, epB, epC, haveMultiplyDoFillComplete);
170 if (haveMultiplyDoFillComplete) {
171 // Due to Epetra wrapper intricacies, we need to explicitly call
172 // fillComplete on Xpetra matrix here. Specifically, EpetraCrsMatrix
173 // only keeps an internal variable to check whether we are in resumed
174 // state or not, but never touches the underlying Epetra object. As
175 // such, we need to explicitly update the state of Xpetra matrix to
176 // that of Epetra one afterwords
177 C.fillComplete();
178 }
179
180 if (i != 0) {
181 std::ostringstream buf;
182 buf << i;
183 std::string msg = "EpetraExt::MatrixMatrix::Jacobi return value of " + buf.str();
184 throw(Exceptions::RuntimeError(msg));
185 }
186#endif
187 } else if (C.getRowMap()->lib() == Xpetra::UseTpetra) {
188#ifdef HAVE_XPETRA_TPETRA
189# if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))) || \
190 (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))))
191 throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra GO=<double,int,long long> enabled."));
192# else
193 const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpA = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(A);
194 const Tpetra::CrsMatrix<SC,LO,GO,NO> & tpB = Xpetra::Helpers<SC,LO,GO,NO>::Op2TpetraCrs(B);
195 Tpetra::CrsMatrix<SC,LO,GO,NO> & tpC = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraCrs(C);
196 const RCP<Tpetra::Vector<SC,LO,GO,NO> > & tpD = toTpetra(Dinv);
197 Tpetra::MatrixMatrix::Jacobi(omega, *tpD, tpA, tpB, tpC, haveMultiplyDoFillComplete, label, params);
198# endif
199#else
200 throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra."));
201#endif
202 }
203
204 if (call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
206 ppp->set("Optimize Storage", doOptimizeStorage);
207 C.fillComplete(B.getDomainMap(), B.getRangeMap(), ppp);
208 }
209
210 // transfer striding information
211 Teuchos::RCP<Xpetra::Matrix<SC,LO,GO,NO> > rcpA = Teuchos::rcp_const_cast<Xpetra::Matrix<SC,LO,GO,NO> >(Teuchos::rcpFromRef(A));
212 Teuchos::RCP<Xpetra::Matrix<SC,LO,GO,NO> > rcpB = Teuchos::rcp_const_cast<Xpetra::Matrix<SC,LO,GO,NO> >(Teuchos::rcpFromRef(B));
213 C.CreateView("stridedMaps", rcpA, false, rcpB, false); // TODO use references instead of RCPs
214 }
215#endif
216
217} // namespace Xpetra
#define XPETRA_DYNAMIC_CAST(type, obj, newObj, exceptionMsg)
LocalOrdinal LO
GlobalOrdinal GO
Exception throws to report errors in the internal logical of the program.
static RCP< Tpetra::CrsMatrix< SC, LO, GO, NO > > Op2NonConstTpetraCrs(RCP< Matrix > Op)
static RCP< Epetra_CrsMatrix > Op2NonConstEpetraCrs(RCP< Matrix > Op)
static RCP< const Tpetra::CrsMatrix< SC, LO, GO, NO > > Op2TpetraCrs(RCP< Matrix > Op)
Xpetra-specific matrix class.
void CreateView(viewLabel_t viewLabel, const RCP< const Map > &rowMap, const RCP< const Map > &colMap)
virtual bool isFillComplete() const =0
Returns true if fillComplete() has been called and the matrix is in compute mode.
virtual void fillComplete(const RCP< const Map > &domainMap, const RCP< const Map > &rangeMap, const RCP< ParameterList > &params=null)=0
Signal that data entry is complete, specifying domain and range maps.
virtual const RCP< const Map > & getRowMap() const
Returns the Map that describes the row distribution in this matrix.
virtual Teuchos::RCP< const Map > getRangeMap() const =0
The Map associated with the range of this operator, which must be compatible with Y....
virtual Teuchos::RCP< const Map > getDomainMap() const =0
The Map associated with the domain of this operator, which must be compatible with X....
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Xpetra namespace
void Jacobi< double, int, long long, EpetraNode >(double omega, const Xpetra::Vector< double, int, long long, EpetraNode > &Dinv, const Xpetra::Matrix< double, int, long long, EpetraNode > &A, const Xpetra::Matrix< double, int, long long, EpetraNode > &B, Xpetra::Matrix< double, int, long long, EpetraNode > &C, bool call_FillComplete_on_result, bool doOptimizeStorage, const std::string &label, const Teuchos::RCP< Teuchos::ParameterList > &params)
void Jacobi< double, int, int, EpetraNode >(double omega, const Xpetra::Vector< double, int, int, EpetraNode > &Dinv, const Xpetra::Matrix< double, int, int, EpetraNode > &A, const Xpetra::Matrix< double, int, int, EpetraNode > &B, Xpetra::Matrix< double, int, int, EpetraNode > &C, bool call_FillComplete_on_result, bool doOptimizeStorage, const std::string &label, const Teuchos::RCP< Teuchos::ParameterList > &params)
RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > toTpetra(const RCP< const CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > &graph)