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