Xpetra Version of the Day
Loading...
Searching...
No Matches
Xpetra_ThyraUtils.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
47
48#include "Xpetra_ConfigDefs.hpp"
49#ifdef HAVE_XPETRA_THYRA
50
51#include "Xpetra_BlockedCrsMatrix.hpp"
52
53//#include "Xpetra_ThyraUtils.hpp"
54
55
56namespace Xpetra {
57
58#ifdef HAVE_XPETRA_EPETRA
59
60#ifndef XPETRA_EPETRA_NO_32BIT_GLOBAL_INDICES
61 // implementation of "toThyra" for full specialization on SC=double, LO=GO=int and NO=EpetraNode
62 // We need the specialization in the cpp file due to a circle dependency in the .hpp files for BlockedCrsMatrix
63 Teuchos::RCP<Thyra::LinearOpBase<double> >
64 ThyraUtils<double, int, int, EpetraNode>::toThyra(const Teuchos::RCP<Xpetra::BlockedCrsMatrix<double, int, int, EpetraNode> >& mat) {
65
66 int nRows = mat->Rows();
67 int nCols = mat->Cols();
68
69 Teuchos::RCP<Xpetra::Matrix<double, int, int, EpetraNode> > Ablock = mat->getInnermostCrsMatrix();
70 Teuchos::RCP<Xpetra::CrsMatrixWrap<double, int, int, EpetraNode> > Ablock_wrap = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<double, int, int, EpetraNode>>(Ablock);
71 TEUCHOS_TEST_FOR_EXCEPT(Ablock_wrap.is_null() == true);
72
73 bool bTpetra = false;
74 bool bEpetra = false;
75#ifdef HAVE_XPETRA_TPETRA
76 // Note: Epetra is enabled
77#if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
78 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
79 Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpetraMat = Teuchos::rcp_dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(Ablock_wrap->getCrsMatrix());
80 if(tpetraMat!=Teuchos::null) bTpetra = true;
81#else
82 bTpetra = false;
83#endif
84#endif
85
86#ifdef HAVE_XPETRA_EPETRA
87 Teuchos::RCP<Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> > epetraMat = Teuchos::rcp_dynamic_cast<Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >(Ablock_wrap->getCrsMatrix());
88 if(epetraMat!=Teuchos::null) bEpetra = true;
89#endif
90
91 TEUCHOS_TEST_FOR_EXCEPT(bTpetra == bEpetra); // we only allow Epetra OR Tpetra
92
93 // create new Thyra blocked operator
94 Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<Scalar> > blockMat =
95 Thyra::defaultBlockedLinearOp<Scalar>();
96
97 blockMat->beginBlockFill(nRows,nCols);
98
99 for (int r=0; r<nRows; ++r) {
100 for (int c=0; c<nCols; ++c) {
101 Teuchos::RCP<Matrix> xpmat = mat->getMatrix(r,c);
102
103 if(xpmat == Teuchos::null) continue; // shortcut for empty blocks
104
105 Teuchos::RCP<Thyra::LinearOpBase<Scalar> > thBlock = Teuchos::null;
106
107 // check whether the subblock is again a blocked operator
108 Teuchos::RCP<BlockedCrsMatrix> xpblock = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(xpmat);
109 if(xpblock != Teuchos::null) {
110 if(xpblock->Rows() == 1 && xpblock->Cols() == 1) {
111 // If it is a single block operator, unwrap it
112 Teuchos::RCP<CrsMatrixWrap> xpwrap = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(xpblock->getCrsMatrix());
113 TEUCHOS_TEST_FOR_EXCEPT(xpwrap.is_null() == true);
114 thBlock = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(xpwrap->getCrsMatrix());
115 } else {
116 // recursive call for general blocked operators
117 thBlock = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(xpblock);
118 }
119 } else {
120 // check whether it is a CRSMatrix object
121 Teuchos::RCP<CrsMatrixWrap> xpwrap = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(xpmat);
122 TEUCHOS_TEST_FOR_EXCEPT(xpwrap.is_null() == true);
123 thBlock = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(xpwrap->getCrsMatrix());
124 }
125
126 blockMat->setBlock(r,c,thBlock);
127 }
128 }
129
130 blockMat->endBlockFill();
131
132 return blockMat;
133 }
134# endif //#ifndef XPETRA_EPETRA_NO_32BIT_GLOBAL_INDICES
135
136#ifndef XPETRA_EPETRA_NO_64BIT_GLOBAL_INDICES
137 // implementation of "toThyra" for full specialization on SC=double, LO=int, GO=long long and NO=EpetraNode
138 // We need the specialization in the cpp file due to a circle dependency in the .hpp files for BlockedCrsMatrix
139 Teuchos::RCP<Thyra::LinearOpBase<double> >
140 ThyraUtils<double, int, long long, EpetraNode>::toThyra(const Teuchos::RCP<Xpetra::BlockedCrsMatrix<double, int, long long, EpetraNode> >& mat) {
141
142 int nRows = mat->Rows();
143 int nCols = mat->Cols();
144
145 Teuchos::RCP<Xpetra::Matrix<double, int, long long, EpetraNode> > Ablock = mat->getInnermostCrsMatrix();
146 Teuchos::RCP<Xpetra::CrsMatrixWrap<double, int, long long, EpetraNode> > Ablock_wrap = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<double, int, long long, EpetraNode>>(Ablock);
147 TEUCHOS_TEST_FOR_EXCEPT(Ablock_wrap.is_null() == true);
148
149 bool bTpetra = false;
150 bool bEpetra = false;
151#ifdef HAVE_XPETRA_TPETRA
152 // Note: Epetra is enabled
153#if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
154 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)))
155 Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpetraMat = Teuchos::rcp_dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(Ablock_wrap->getCrsMatrix());
156 if(tpetraMat!=Teuchos::null) bTpetra = true;
157#else
158 bTpetra = false;
159#endif
160#endif
161
162#ifdef HAVE_XPETRA_EPETRA
163 Teuchos::RCP<Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> > epetraMat = Teuchos::rcp_dynamic_cast<Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >(Ablock_wrap->getCrsMatrix());
164 if(epetraMat!=Teuchos::null) bEpetra = true;
165#endif
166
167 TEUCHOS_TEST_FOR_EXCEPT(bTpetra == bEpetra); // we only allow Epetra OR Tpetra
168
169 // create new Thyra blocked operator
170 Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<Scalar> > blockMat =
171 Thyra::defaultBlockedLinearOp<Scalar>();
172
173 blockMat->beginBlockFill(nRows,nCols);
174
175 for (int r=0; r<nRows; ++r) {
176 for (int c=0; c<nCols; ++c) {
177 Teuchos::RCP<Matrix> xpmat = mat->getMatrix(r,c);
178
179 if(xpmat == Teuchos::null) continue; // shortcut for empty blocks
180
181 Teuchos::RCP<Thyra::LinearOpBase<Scalar> > thBlock = Teuchos::null;
182
183 // check whether the subblock is again a blocked operator
184 Teuchos::RCP<BlockedCrsMatrix> xpblock = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(xpmat);
185 if(xpblock != Teuchos::null) {
186 if(xpblock->Rows() == 1 && xpblock->Cols() == 1) {
187 // If it is a single block operator, unwrap it
188 Teuchos::RCP<CrsMatrixWrap> xpwrap = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(xpblock->getCrsMatrix());
189 TEUCHOS_TEST_FOR_EXCEPT(xpwrap.is_null() == true);
190 thBlock = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(xpwrap->getCrsMatrix());
191 } else {
192 // recursive call for general blocked operators
193 thBlock = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(xpblock);
194 }
195 } else {
196 // check whether it is a CRSMatrix object
197 Teuchos::RCP<CrsMatrixWrap> xpwrap = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(xpmat);
198 TEUCHOS_TEST_FOR_EXCEPT(xpwrap.is_null() == true);
199 thBlock = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(xpwrap->getCrsMatrix());
200 }
201
202 blockMat->setBlock(r,c,thBlock);
203 }
204 }
205
206 blockMat->endBlockFill();
207
208 return blockMat;
209 }
210#endif // #ifndef XPETRA_EPETRA_NO_64BIT_GLOBAL_INDICES
211
212#endif
213
214} // namespace Xpetra
215
216#endif // HAVE_XPETRA_THYRA
Xpetra_BlockedCrsMatrix.hpp
Xpetra_ConfigDefs.hpp
Teuchos::RCP::is_null
bool is_null() const
TEUCHOS_TEST_FOR_EXCEPT
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
Xpetra
Xpetra namespace
Definition Xpetra_BlockedCrsMatrix.hpp:89