Anasazi Version of the Day
Loading...
Searching...
No Matches
Tsqr_TwoLevelDistTsqr.hpp
1// @HEADER
2// ***********************************************************************
3//
4// Anasazi: Block Eigensolvers Package
5// Copyright 2004 Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40// @HEADER
41
42#ifndef __TSQR_Tsqr_TwoLevelDistTsqr_hpp
43#define __TSQR_Tsqr_TwoLevelDistTsqr_hpp
44
45#include <Tsqr_DistTsqr.hpp>
46#include <Tsqr_MpiCommFactory.hpp>
47#include <sstream>
48
51
52namespace TSQR {
53
60 template< class LocalOrdinal,
61 class Scalar,
62 class DistTsqrType = DistTsqr< LocalOrdinal, Scalar > >
64 public:
65 typedef Scalar scalar_type;
66 typedef LocalOrdinal ordinal_type;
67 typedef DistTsqrType dist_tsqr_type;
68 typedef typename dist_tsqr_type::rank_type rank_type;
69 typedef typename Teuchos::RCP< dist_tsqr_type > dist_tsqr_ptr;
70 typedef typename dist_tsqr_type::FactorOutput DistTsqrFactorOutput;
71 typedef std::pair< DistTsqrFactorOutput, DistTsqrFactorOutput > FactorOutput;
72
76 worldMess_ (TSQR::MPI::makeMpiCommWorld()),
77 nodeDistTsqr_ (TSQR::MPI::makeMpiCommNode()),
78 networkDistTsqr_ (TSQR::MPI::makeMpiCommNetwork())
79 {}
80
85
89 return nodeDistTsqr_->QR_produces_R_factor_with_nonnegative_diagonal() &&
90 networkDistTsqr_->QR_produces_R_factor_with_nonnegative_diagonal();
91 }
92
110 FactorOutput
111 factor (MatView< LocalOrdinal, Scalar > R_mine)
112 {
113 DistTsqrFactorOutput nodeOutput = nodeDistTsqr_->factor (R_mine);
114 DistTsqrFactorOutput networkOutput = networkDistTsqr_->factor (R_mine);
115 return std::make_pair (nodeOutput, networkOutput);
116 }
117
118 void
119 apply (const ApplyType& applyType,
120 const LocalOrdinal ncols_C,
121 const LocalOrdinal ncols_Q,
122 Scalar C_mine[],
123 const LocalOrdinal ldc_mine,
124 const FactorOutput& factorOutput)
125 {
126 if (applyType.transposed())
127 {
128 nodeDistTsqr_->apply (applyType, ncols_C, ncols_Q,
129 C_mine, ldc_mine, factorOutput.first);
130 networkDistTsqr_->apply (applyType, ncols_C, ncols_Q,
131 C_mine, ldc_mine, factorOutput.second);
132 }
133 else
134 {
135 networkDistTsqr_->apply (applyType, ncols_C, ncols_Q,
136 C_mine, ldc_mine, factorOutput.second);
137 nodeDistTsqr_->apply (applyType, ncols_C, ncols_Q,
138 C_mine, ldc_mine, factorOutput.first);
139 }
140 }
141
142 void
143 explicit_Q (const LocalOrdinal ncols_Q,
144 Scalar Q_mine[],
145 const LocalOrdinal ldq_mine,
146 const FactorOutput& factorOutput)
147 {
148 typedef MatView< LocalOrdinal, Scalar > matview_type;
149 matview_type Q_view (ncols_Q, ncols_Q, Q_mine, ldq_mine, Scalar(0));
150 Q_view.fill (Scalar(0));
151
152 const rank_type myRank = worldMess_->rank();
153 if (myRank == 0)
154 {
155 if (networkMess_->rank() != 0)
156 {
157 std::ostringstream os;
158 os << "My rank with respect to MPI_COMM_WORLD is 0, but my rank "
159 "with respect to MPI_COMM_NETWORK is nonzero (= "
160 << networkMess_->rank() << "). We could deal with this by "
161 "swapping data between those two ranks, but we haven\'t "
162 "implemented that fix yet.";
163 throw std::logic_error (os.str());
164 }
165 for (LocalOrdinal j = 0; j < ncols_Q; ++j)
166 Q_view(j, j) = Scalar (1);
167 }
168 apply (ApplyType::NoTranspose, ncols_Q, ncols_Q,
169 Q_mine, ldq_mine, factorOutput);
170 }
171
172 private:
173 Teuchos::RCP< MessengerBase< Scalar > > worldMess_;
174 dist_tsqr_ptr nodeDistTsqr_, networkDistTsqr_;
175 };
176
177} // namespace TSQR
178
179#endif // __TSQR_Tsqr_TwoLevelDistTsqr_hpp
Interprocess part of TSQR.
bool QR_produces_R_factor_with_nonnegative_diagonal() const
FactorOutput factor(MatView< LocalOrdinal, Scalar > R_mine)
Compute QR factorization of R factors, one per MPI process.