Anasazi Version of the Day
Tsqr_TwoLevelDistTsqr.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_TwoLevelDistTsqr_hpp
00030 #define __TSQR_Tsqr_TwoLevelDistTsqr_hpp
00031 
00032 #include <Tsqr_DistTsqr.hpp>
00033 #include <Tsqr_MpiCommFactory.hpp>
00034 #include <sstream>
00035 
00038 
00039 namespace TSQR {
00040 
00047   template< class LocalOrdinal, 
00048       class Scalar, 
00049       class DistTsqrType = DistTsqr< LocalOrdinal, Scalar > >
00050   class TwoLevelDistTsqr {
00051   public:
00052     typedef Scalar scalar_type;
00053     typedef LocalOrdinal ordinal_type;
00054     typedef DistTsqrType dist_tsqr_type;
00055     typedef typename dist_tsqr_type::rank_type rank_type;
00056     typedef typename Teuchos::RCP< dist_tsqr_type > dist_tsqr_ptr;
00057     typedef typename dist_tsqr_type::FactorOutput DistTsqrFactorOutput;
00058     typedef std::pair< DistTsqrFactorOutput, DistTsqrFactorOutput > FactorOutput;
00059 
00062     TwoLevelDistTsqr () :
00063       worldMess_ (TSQR::MPI::makeMpiCommWorld()),
00064       nodeDistTsqr_ (TSQR::MPI::makeMpiCommNode()),
00065       networkDistTsqr_ (TSQR::MPI::makeMpiCommNetwork())
00066     {}
00067 
00071     ~TwoLevelDistTsqr () {}
00072 
00075     bool QR_produces_R_factor_with_nonnegative_diagonal () const {
00076       return nodeDistTsqr_->QR_produces_R_factor_with_nonnegative_diagonal() &&
00077   networkDistTsqr_->QR_produces_R_factor_with_nonnegative_diagonal();
00078     }
00079 
00097     FactorOutput
00098     factor (MatView< LocalOrdinal, Scalar > R_mine)
00099     {
00100       DistTsqrFactorOutput nodeOutput = nodeDistTsqr_->factor (R_mine);
00101       DistTsqrFactorOutput networkOutput = networkDistTsqr_->factor (R_mine);
00102       return std::make_pair (nodeOutput, networkOutput);
00103     }
00104 
00105     void
00106     apply (const ApplyType& applyType,
00107      const LocalOrdinal ncols_C,
00108      const LocalOrdinal ncols_Q,
00109      Scalar C_mine[],
00110      const LocalOrdinal ldc_mine,
00111      const FactorOutput& factorOutput)
00112     {
00113       if (applyType.transposed())
00114   {
00115     nodeDistTsqr_->apply (applyType, ncols_C, ncols_Q, 
00116         C_mine, ldc_mine, factorOutput.first);
00117     networkDistTsqr_->apply (applyType, ncols_C, ncols_Q, 
00118            C_mine, ldc_mine, factorOutput.second);
00119   }
00120       else 
00121   {
00122     networkDistTsqr_->apply (applyType, ncols_C, ncols_Q, 
00123            C_mine, ldc_mine, factorOutput.second);
00124     nodeDistTsqr_->apply (applyType, ncols_C, ncols_Q, 
00125         C_mine, ldc_mine, factorOutput.first);
00126   }
00127     }
00128 
00129     void
00130     explicit_Q (const LocalOrdinal ncols_Q,
00131     Scalar Q_mine[],
00132     const LocalOrdinal ldq_mine,
00133     const FactorOutput& factorOutput)
00134     {
00135       typedef MatView< LocalOrdinal, Scalar > matview_type;
00136       matview_type Q_view (ncols_Q, ncols_Q, Q_mine, ldq_mine, Scalar(0));
00137       Q_view.fill (Scalar(0));
00138   
00139       const rank_type myRank = worldMess_->rank();
00140       if (myRank == 0)
00141   {
00142     if (networkMess_->rank() != 0)
00143       {
00144         std::ostringstream os;
00145         os << "My rank with respect to MPI_COMM_WORLD is 0, but my rank "
00146     "with respect to MPI_COMM_NETWORK is nonzero (= " 
00147      << networkMess_->rank() << ").  We could deal with this by "
00148     "swapping data between those two ranks, but we haven\'t "
00149     "implemented that fix yet.";
00150         throw std::logic_error (os.str());
00151       }
00152     for (LocalOrdinal j = 0; j < ncols_Q; ++j)
00153       Q_view(j, j) = Scalar (1);
00154   }
00155       apply (ApplyType::NoTranspose, ncols_Q, ncols_Q, 
00156        Q_mine, ldq_mine, factorOutput);
00157     }
00158 
00159   private:
00160     Teuchos::RCP< MessengerBase< Scalar > > worldMess_;
00161     dist_tsqr_ptr nodeDistTsqr_, networkDistTsqr_;
00162   };
00163 
00164 } // namespace TSQR
00165 
00166 #endif // __TSQR_Tsqr_TwoLevelDistTsqr_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends