Anasazi Version of the Day
Tsqr_printGlobalMatrix.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_printGlobalMatrix_hpp
00030 #define __Tsqr_printGlobalMatrix_hpp
00031 
00032 #include "Tsqr_MessengerBase.hpp"
00033 #include "Tsqr_Util.hpp"
00034 #include "Tsqr_Matrix.hpp"
00035 
00036 #include <limits>
00037 #include <ostream>
00038 #include <stdexcept>
00039 
00042 
00043 namespace TSQR {
00044 
00048   template< class ConstMatrixViewType >
00049   void
00050   printGlobalMatrix (std::ostream& out,
00051          const ConstMatrixViewType& A_local,
00052          MessengerBase< typename ConstMatrixViewType::scalar_type >* const scalarComm,
00053          MessengerBase< typename ConstMatrixViewType::ordinal_type >* const ordinalComm)
00054     {
00055       typedef typename ConstMatrixViewType::ordinal_type LocalOrdinal;
00056       typedef typename ConstMatrixViewType::scalar_type Scalar;
00057       using std::endl;
00058 
00059       const int myRank = scalarComm->rank ();
00060       const int nprocs = scalarComm->size ();
00061       const LocalOrdinal nrowsLocal = A_local.nrows();
00062       const LocalOrdinal ncols = A_local.ncols();
00063       const Scalar quiet_NaN = std::numeric_limits< Scalar >::quiet_NaN();
00064 
00065       if (myRank == 0)
00066   {
00067     // Print the remote matrix data
00068     // out << "Processor " << my_rank << ":" << endl;
00069     print_local_matrix (out, A_local.nrows(), A_local.ncols(), 
00070             A_local.get(), A_local.lda());
00071 
00072     // Space for remote matrix data.  Other processors are allowed
00073     // to have different nrows_local values; we make space as
00074     // necessary.
00075     Matrix< LocalOrdinal, Scalar > A_remote (nrowsLocal, ncols, quiet_NaN);
00076 
00077     // Loop through all the other processors in order.
00078     // Fetch their matrix data and print it.
00079     for (int srcProc = 1; srcProc < nprocs; ++srcProc)
00080       {
00081         // Get processor proc's local matrix dimensions
00082         LocalOrdinal dims[2];
00083         ordinalComm->recv (&dims[0], 2, srcProc, 0);
00084 
00085         // Make space for the remote matrix data
00086         if (std::numeric_limits< LocalOrdinal >::is_signed)
00087     {
00088       if (dims[0] <= 0 || dims[1] <= 0)
00089         throw std::runtime_error ("Invalid dimensions of remote matrix");
00090     }
00091         else
00092     {
00093       if (dims[0] == 0 || dims[1] == 0)
00094         throw std::runtime_error ("Invalid dimensions of remote matrix");
00095     }
00096         A_remote.reshape (dims[0], dims[1]);
00097 
00098         // Receive the remote matrix data, which we assume is
00099         // stored contiguously.
00100         scalarComm->recv (A_remote.get(), dims[0]*dims[1], srcProc, 0);
00101 
00102         // Print the remote matrix data
00103         // out << "Processor " << proc << ":" << endl;
00104         print_local_matrix (out, dims[0], dims[0], A_remote.get(), A_remote.lda());
00105       }
00106   }
00107       else
00108   {
00109     // Send my local matrix dimensions to proc 0.
00110     int rootProc = 0;
00111     LocalOrdinal dims[2];
00112 
00113     dims[0] = nrowsLocal;
00114     dims[1] = ncols;
00115     ordinalComm->send (dims, 2, rootProc, 0);
00116 
00117     // Create a (contiguous) buffer and copy the data into it.
00118     Matrix< LocalOrdinal, Scalar > A_buf (nrowsLocal, ncols, quiet_NaN);
00119     A_buf.copy (A_local);
00120 
00121     // Send the actual data to proc 0.
00122     scalarComm->send (A_buf.get(), nrowsLocal*ncols, rootProc, 0);
00123   }
00124       scalarComm->barrier ();
00125     }
00126 
00127 } // namespace TSQR
00128 
00129 #endif // __Tsqr_printGlobalMatrix_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends