Anasazi  Version of the Day
Tsqr_TwoLevelDistTsqr.hpp
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Anasazi: Block Eigensolvers Package
5 // Copyright (2010) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // This library is free software; you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as
12 // published by the Free Software Foundation; either version 2.1 of the
13 // License, or (at your option) any later version.
14 //
15 // This library is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23 // USA
24 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
25 //
26 // ***********************************************************************
27 // @HEADER
28 
29 #ifndef __TSQR_Tsqr_TwoLevelDistTsqr_hpp
30 #define __TSQR_Tsqr_TwoLevelDistTsqr_hpp
31 
32 #include <Tsqr_DistTsqr.hpp>
33 #include <Tsqr_MpiCommFactory.hpp>
34 #include <sstream>
35 
38 
39 namespace TSQR {
40 
47  template< class LocalOrdinal,
48  class Scalar,
49  class DistTsqrType = DistTsqr< LocalOrdinal, Scalar > >
51  public:
52  typedef Scalar scalar_type;
53  typedef LocalOrdinal ordinal_type;
54  typedef DistTsqrType dist_tsqr_type;
55  typedef typename dist_tsqr_type::rank_type rank_type;
56  typedef typename Teuchos::RCP< dist_tsqr_type > dist_tsqr_ptr;
57  typedef typename dist_tsqr_type::FactorOutput DistTsqrFactorOutput;
58  typedef std::pair< DistTsqrFactorOutput, DistTsqrFactorOutput > FactorOutput;
59 
63  worldMess_ (TSQR::MPI::makeMpiCommWorld()),
64  nodeDistTsqr_ (TSQR::MPI::makeMpiCommNode()),
65  networkDistTsqr_ (TSQR::MPI::makeMpiCommNetwork())
66  {}
67 
72 
76  return nodeDistTsqr_->QR_produces_R_factor_with_nonnegative_diagonal() &&
77  networkDistTsqr_->QR_produces_R_factor_with_nonnegative_diagonal();
78  }
79 
97  FactorOutput
98  factor (MatView< LocalOrdinal, Scalar > R_mine)
99  {
100  DistTsqrFactorOutput nodeOutput = nodeDistTsqr_->factor (R_mine);
101  DistTsqrFactorOutput networkOutput = networkDistTsqr_->factor (R_mine);
102  return std::make_pair (nodeOutput, networkOutput);
103  }
104 
105  void
106  apply (const ApplyType& applyType,
107  const LocalOrdinal ncols_C,
108  const LocalOrdinal ncols_Q,
109  Scalar C_mine[],
110  const LocalOrdinal ldc_mine,
111  const FactorOutput& factorOutput)
112  {
113  if (applyType.transposed())
114  {
115  nodeDistTsqr_->apply (applyType, ncols_C, ncols_Q,
116  C_mine, ldc_mine, factorOutput.first);
117  networkDistTsqr_->apply (applyType, ncols_C, ncols_Q,
118  C_mine, ldc_mine, factorOutput.second);
119  }
120  else
121  {
122  networkDistTsqr_->apply (applyType, ncols_C, ncols_Q,
123  C_mine, ldc_mine, factorOutput.second);
124  nodeDistTsqr_->apply (applyType, ncols_C, ncols_Q,
125  C_mine, ldc_mine, factorOutput.first);
126  }
127  }
128 
129  void
130  explicit_Q (const LocalOrdinal ncols_Q,
131  Scalar Q_mine[],
132  const LocalOrdinal ldq_mine,
133  const FactorOutput& factorOutput)
134  {
135  typedef MatView< LocalOrdinal, Scalar > matview_type;
136  matview_type Q_view (ncols_Q, ncols_Q, Q_mine, ldq_mine, Scalar(0));
137  Q_view.fill (Scalar(0));
138 
139  const rank_type myRank = worldMess_->rank();
140  if (myRank == 0)
141  {
142  if (networkMess_->rank() != 0)
143  {
144  std::ostringstream os;
145  os << "My rank with respect to MPI_COMM_WORLD is 0, but my rank "
146  "with respect to MPI_COMM_NETWORK is nonzero (= "
147  << networkMess_->rank() << "). We could deal with this by "
148  "swapping data between those two ranks, but we haven\'t "
149  "implemented that fix yet.";
150  throw std::logic_error (os.str());
151  }
152  for (LocalOrdinal j = 0; j < ncols_Q; ++j)
153  Q_view(j, j) = Scalar (1);
154  }
155  apply (ApplyType::NoTranspose, ncols_Q, ncols_Q,
156  Q_mine, ldq_mine, factorOutput);
157  }
158 
159  private:
160  Teuchos::RCP< MessengerBase< Scalar > > worldMess_;
161  dist_tsqr_ptr nodeDistTsqr_, networkDistTsqr_;
162  };
163 
164 } // namespace TSQR
165 
166 #endif // __TSQR_Tsqr_TwoLevelDistTsqr_hpp
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.
Interprocess part of TSQR.