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 
56 namespace 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
Xpetra::CrsMatrixWrap
Concrete implementation of Xpetra::Matrix.
Definition: Xpetra_CrsMatrixWrap_fwd.hpp:51