Anasazi Version of the Day
Tsqr_DistTsqr.hpp
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                 Anasazi: Block Eigensolvers Package
00005 //                 Copyright (2010) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #ifndef __TSQR_Tsqr_DistTsqr_hpp
00030 #define __TSQR_Tsqr_DistTsqr_hpp
00031 
00032 #include <Tsqr_DistTsqrHelper.hpp>
00033 #include <Tsqr_DistTsqrRB.hpp>
00034 
00035 #include <utility> // std::pair
00036 
00039 
00040 namespace TSQR {
00041 
00055   template< class LocalOrdinal, class Scalar >
00056   class DistTsqr {
00057   public:
00058     typedef Scalar scalar_type;
00059     typedef LocalOrdinal ordinal_type;
00060     typedef MatView< ordinal_type, scalar_type > matview_type;
00061     typedef std::vector< std::vector< scalar_type > > VecVec;
00062     typedef std::pair< VecVec, VecVec > FactorOutput;
00063     typedef int rank_type;
00064 
00069     DistTsqr (const Teuchos::RCP< MessengerBase< scalar_type > >& messenger) :
00070       messenger_ (messenger),
00071       reduceBroadcastImpl_ (messenger)
00072     {}
00073 
00076     rank_type rank() const { return messenger_->rank(); }
00077 
00080     rank_type size() const { return messenger_->size(); }
00081 
00086     ~DistTsqr () {}
00087 
00090     bool QR_produces_R_factor_with_nonnegative_diagonal () const {
00091       typedef Combine< ordinal_type, scalar_type > combine_type;
00092       return combine_type::QR_produces_R_factor_with_nonnegative_diagonal() &&
00093   reduceBroadcastImpl_.QR_produces_R_factor_with_nonnegative_diagonal();
00094     }
00095 
00115     void
00116     factorExplicit (matview_type R_mine, matview_type Q_mine)
00117     {
00118       reduceBroadcastImpl_.factorExplicit (R_mine, Q_mine);
00119     }
00120 
00124     void 
00125     getFactorExplicitTimings (std::vector< TimeStats >& stats) const
00126     {
00127       reduceBroadcastImpl_.getStats (stats);
00128     }
00129 
00133     void
00134     getFactorExplicitTimingLabels (std::vector< std::string >& labels) const
00135     {
00136       reduceBroadcastImpl_.getStatsLabels (labels);
00137     }
00138 
00162     FactorOutput
00163     factor (matview_type R_mine)
00164     {
00165       VecVec Q_factors, tau_arrays;
00166       DistTsqrHelper< ordinal_type, scalar_type > helper;
00167       const ordinal_type ncols = R_mine.ncols();
00168 
00169       std::vector< scalar_type > R_local (ncols*ncols);
00170       copy_matrix (ncols, ncols, &R_local[0], ncols, R_mine.get(), R_mine.lda());
00171 
00172       const int P = messenger_->size();
00173       const int my_rank = messenger_->rank();
00174       const int first_tag = 0;
00175       std::vector< scalar_type > work (ncols);
00176       helper.factor_helper (ncols, R_local, my_rank, 0, P-1, first_tag, 
00177           messenger_.get(), Q_factors, tau_arrays, work);
00178       copy_matrix (ncols, ncols, R_mine.get(), R_mine.lda(), &R_local[0], ncols);
00179       return std::make_pair (Q_factors, tau_arrays);
00180     }
00181 
00182     void
00183     apply (const ApplyType& apply_type,
00184      const ordinal_type ncols_C,
00185      const ordinal_type ncols_Q,
00186      scalar_type C_mine[],
00187      const ordinal_type ldc_mine,
00188      const FactorOutput& factor_output)
00189     {
00190       const bool transposed = apply_type.transposed();
00191 
00192       if (transposed)
00193   throw std::logic_error("DistTsqr: Applying Q^T or Q^H "
00194              "not yet implemented");
00195 
00196       const int P = messenger_->size();
00197       const int my_rank = messenger_->rank();
00198       const int first_tag = 0;
00199       std::vector< scalar_type > C_other (ncols_C * ncols_C);
00200       std::vector< scalar_type > work (ncols_C);
00201   
00202       const VecVec& Q_factors = factor_output.first;
00203       const VecVec& tau_arrays = factor_output.second;
00204 
00205       // assert (Q_factors.size() == tau_arrays.size());
00206       const int cur_pos = Q_factors.size() - 1;
00207       DistTsqrHelper< ordinal_type, scalar_type > helper;
00208       helper.apply_helper (apply_type, ncols_C, ncols_Q, C_mine, ldc_mine, 
00209          &C_other[0], my_rank, 0, P-1, first_tag, 
00210          messenger_.get(), Q_factors, tau_arrays, cur_pos, 
00211          work);
00212     }
00213 
00214     void
00215     explicit_Q (const ordinal_type ncols_Q,
00216     scalar_type Q_mine[],
00217     const ordinal_type ldq_mine,
00218     const FactorOutput& factor_output)
00219     {
00220       const int my_rank = messenger_->rank ();
00221       fill_matrix (ncols_Q, ncols_Q, Q_mine, ldq_mine, Scalar(0));
00222       if (my_rank == 0)
00223   {
00224     for (ordinal_type j = 0; j < ncols_Q; ++j)
00225       Q_mine[j + j*ldq_mine] = Scalar (1);
00226   }
00227       apply (ApplyType::NoTranspose, ncols_Q, ncols_Q, 
00228        Q_mine, ldq_mine, factor_output);
00229     }
00230 
00231   private:
00232     Teuchos::RCP< MessengerBase< scalar_type > > messenger_;
00233     DistTsqrRB< ordinal_type, scalar_type > reduceBroadcastImpl_;
00234   };
00235 
00236 } // namespace TSQR
00237 
00238 #endif // __TSQR_Tsqr_DistTsqr_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends