Amesos2 - Direct Sparse Solver Interfaces Version of the Day
Amesos2_Util.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 //
00003 // ***********************************************************************
00004 //
00005 //           Amesos2: Templated Direct Sparse Solver Package
00006 //                  Copyright 2011 Sandia Corporation
00007 //
00008 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00009 // the U.S. Government retains certain rights in this software.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00039 //
00040 // ***********************************************************************
00041 //
00042 // @HEADER
00043 
00052 #ifndef AMESOS2_UTIL_HPP
00053 #define AMESOS2_UTIL_HPP
00054 
00055 #include "Amesos2_config.h"
00056 
00057 #include <Teuchos_RCP.hpp>
00058 #include <Teuchos_BLAS_types.hpp>
00059 #include <Teuchos_ArrayView.hpp>
00060 #include <Teuchos_FancyOStream.hpp>
00061 
00062 #ifdef HAVE_AMESOS2_EPETRA
00063 #  include <Epetra_Map.h>
00064 #endif
00065 
00066 #include <Tpetra_Map.hpp>
00067 
00068 #include "Amesos2_TypeDecl.hpp"
00069 #include "Amesos2_Meta.hpp"
00070 
00071 
00072 namespace Amesos2 {
00073 
00074   namespace Util {
00075 
00082     using Teuchos::RCP;
00083     using Teuchos::ArrayView;
00084 
00085     using Meta::is_same;
00086     using Meta::if_then_else;
00087 
00103     template <typename LO, typename GO, typename GS, typename Node>
00104     const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
00105     getDistributionMap(EDistribution distribution,
00106            GS num_global_elements,
00107            const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
00108            GO indexBase = 0);
00109     
00110 
00111 #ifdef HAVE_AMESOS2_EPETRA
00112 
00117     template <typename LO, typename GO, typename GS, typename Node>
00118     RCP<Tpetra::Map<LO,GO,Node> >
00119     epetra_map_to_tpetra_map(const Epetra_BlockMap& map);
00120 
00126     template <typename LO, typename GO, typename GS, typename Node>
00127     RCP<Epetra_Map>
00128     tpetra_map_to_epetra_map(const Tpetra::Map<LO,GO,Node>& map);
00129 
00135     const RCP<const Teuchos::Comm<int> > to_teuchos_comm(RCP<const Epetra_Comm> c);
00136 
00142     const RCP<const Epetra_Comm> to_epetra_comm(RCP<const Teuchos::Comm<int> > c);
00143 #endif  // HAVE_AMESOS2_EPETRA
00144 
00150     template <typename Scalar,
00151         typename GlobalOrdinal,
00152         typename GlobalSizeT>
00153     void transpose(ArrayView<Scalar> vals,
00154        ArrayView<GlobalOrdinal> indices,
00155        ArrayView<GlobalSizeT> ptr,
00156        ArrayView<Scalar> trans_vals,
00157        ArrayView<GlobalOrdinal> trans_indices,
00158        ArrayView<GlobalSizeT> trans_ptr);
00159 
00173     template <typename Scalar1, typename Scalar2>
00174     void scale(ArrayView<Scalar1> vals, size_t l,
00175          size_t ld, ArrayView<Scalar2> s);
00176 
00195     template <typename Scalar1, typename Scalar2, class BinaryOp>
00196     void scale(ArrayView<Scalar1> vals, size_t l,
00197          size_t ld, ArrayView<Scalar2> s, BinaryOp binary_op);
00198 
00199 
00201     void printLine( Teuchos::FancyOStream &out );
00202 
00203 
00205     // Matrix/MultiVector Utilities //
00207     
00208 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00209     /*
00210      * The following represents a general way of getting a CRS or CCS
00211      * representation of a matrix with implicit type conversions.  The
00212      * \c get_crs_helper and \c get_ccs_helper classes are templated
00213      * on 4 types:
00214      *
00215      * - A Matrix type (conforming to the Amesos2 MatrixAdapter interface)
00216      * - A scalar type
00217      * - A global ordinal type
00218      * - A global size type
00219      *
00220      * The last three template types correspond to the input argument
00221      * types.  For example, if the scalar type is \c double , then we
00222      * require that the \c nzvals argument is a \c
00223      * Teuchos::ArrayView<double> type.
00224      *
00225      * These helpers perform any type conversions that must be
00226      * performed to go between the Matrix's types and the input types.
00227      * If no conversions are necessary the static functions can be
00228      * effectively inlined.
00229      */
00230 
00231     template <class M, typename S, typename GO, typename GS, class Op>
00232     struct same_gs_helper
00233     {
00234       static void do_get(const Teuchos::Ptr<const M> mat,
00235        const ArrayView<typename M::scalar_t> nzvals,
00236        const ArrayView<typename M::global_ordinal_t> indices,
00237        const ArrayView<GS> pointers,
00238        GS& nnz,
00239        const Teuchos::Ptr<
00240          const Tpetra::Map<typename M::local_ordinal_t,
00241                                              typename M::global_ordinal_t,
00242                            typename M::node_t> > map,
00243        EStorage_Ordering ordering)
00244       {
00245   Op::apply(mat, nzvals, indices, pointers, nnz, map, ordering);
00246       }
00247     };
00248 
00249     template <class M, typename S, typename GO, typename GS, class Op>
00250     struct diff_gs_helper
00251     {
00252       static void do_get(const Teuchos::Ptr<const M> mat,
00253        const ArrayView<typename M::scalar_t> nzvals,
00254        const ArrayView<typename M::global_ordinal_t> indices,
00255        const ArrayView<GS> pointers,
00256        GS& nnz,
00257        const Teuchos::Ptr<
00258          const Tpetra::Map<typename M::local_ordinal_t,
00259                                              typename M::global_ordinal_t,
00260                            typename M::node_t> > map,
00261        EStorage_Ordering ordering)
00262       {
00263   typedef typename M::global_size_t mat_gs_t;
00264   typename ArrayView<GS>::size_type i, size = pointers.size();
00265   Teuchos::Array<mat_gs_t> pointers_tmp(size);
00266   mat_gs_t nnz_tmp = 0;
00267 
00268   Op::apply(mat, nzvals, indices, pointers_tmp, nnz_tmp, map, ordering);
00269 
00270   for (i = 0; i < size; ++i){
00271     pointers[i] = Teuchos::as<GS>(pointers_tmp[i]);
00272   }
00273   nnz = Teuchos::as<GS>(nnz_tmp);
00274       }
00275     };
00276 
00277     template <class M, typename S, typename GO, typename GS, class Op>
00278     struct same_go_helper
00279     {
00280       static void do_get(const Teuchos::Ptr<const M> mat,
00281        const ArrayView<typename M::scalar_t> nzvals,
00282        const ArrayView<GO> indices,
00283        const ArrayView<GS> pointers,
00284        GS& nnz,
00285        const Teuchos::Ptr<
00286          const Tpetra::Map<typename M::local_ordinal_t,
00287                                              typename M::global_ordinal_t,
00288                            typename M::node_t> > map,
00289        EStorage_Ordering ordering)
00290       {
00291   typedef typename M::global_size_t mat_gs_t;
00292   if_then_else<is_same<GS,mat_gs_t>::value,
00293     same_gs_helper<M,S,GO,GS,Op>,
00294     diff_gs_helper<M,S,GO,GS,Op> >::type::do_get(mat, nzvals, indices,
00295                    pointers, nnz, map,
00296                    ordering);
00297       }
00298     };
00299 
00300     template <class M, typename S, typename GO, typename GS, class Op>
00301     struct diff_go_helper
00302     {
00303       static void do_get(const Teuchos::Ptr<const M> mat,
00304        const ArrayView<typename M::scalar_t> nzvals,
00305        const ArrayView<GO> indices,
00306        const ArrayView<GS> pointers,
00307        GS& nnz,
00308        const Teuchos::Ptr<
00309          const Tpetra::Map<typename M::local_ordinal_t,
00310                                              typename M::global_ordinal_t,
00311                            typename M::node_t> > map,
00312        EStorage_Ordering ordering)
00313       {
00314   typedef typename M::global_ordinal_t mat_go_t;
00315   typedef typename M::global_size_t mat_gs_t;
00316   typename ArrayView<GO>::size_type i, size = indices.size();
00317   Teuchos::Array<mat_go_t> indices_tmp(size);
00318 
00319   if_then_else<is_same<GS,mat_gs_t>::value,
00320     same_gs_helper<M,S,GO,GS,Op>,
00321     diff_gs_helper<M,S,GO,GS,Op> >::type::do_get(mat, nzvals, indices_tmp,
00322                    pointers, nnz, map,
00323                    ordering);
00324 
00325   for (i = 0; i < size; ++i){
00326     indices[i] = Teuchos::as<GO>(indices_tmp[i]);
00327   }
00328       }
00329     };
00330 
00331     template <class M, typename S, typename GO, typename GS, class Op>
00332     struct same_scalar_helper
00333     {
00334       static void do_get(const Teuchos::Ptr<const M> mat,
00335        const ArrayView<S> nzvals,
00336        const ArrayView<GO> indices,
00337        const ArrayView<GS> pointers,
00338        GS& nnz,
00339        const Teuchos::Ptr<
00340          const Tpetra::Map<typename M::local_ordinal_t,
00341                                              typename M::global_ordinal_t,
00342                            typename M::node_t> > map,
00343        EStorage_Ordering ordering)
00344       {
00345   typedef typename M::global_ordinal_t mat_go_t;
00346   if_then_else<is_same<GO,mat_go_t>::value,
00347     same_go_helper<M,S,GO,GS,Op>,
00348     diff_go_helper<M,S,GO,GS,Op> >::type::do_get(mat, nzvals, indices,
00349                    pointers, nnz, map,
00350                    ordering);
00351       }
00352     };
00353 
00354     template <class M, typename S, typename GO, typename GS, class Op>
00355     struct diff_scalar_helper
00356     {
00357       static void do_get(const Teuchos::Ptr<const M> mat,
00358        const ArrayView<S> nzvals,
00359        const ArrayView<GO> indices,
00360        const ArrayView<GS> pointers,
00361        GS& nnz,
00362        const Teuchos::Ptr<
00363          const Tpetra::Map<typename M::local_ordinal_t,
00364                                              typename M::global_ordinal_t,
00365                            typename M::node_t> > map,
00366        EStorage_Ordering ordering)
00367       {
00368   typedef typename M::scalar_t mat_scalar_t;
00369   typedef typename M::global_ordinal_t mat_go_t;
00370   typename ArrayView<S>::size_type i, size = nzvals.size();
00371   Teuchos::Array<mat_scalar_t> nzvals_tmp(size);
00372 
00373   if_then_else<is_same<GO,mat_go_t>::value,
00374     same_go_helper<M,S,GO,GS,Op>,
00375     diff_go_helper<M,S,GO,GS,Op> >::type::do_get(mat, nzvals_tmp, indices,
00376                    pointers, nnz, map,
00377                    ordering);
00378 
00379   for (i = 0; i < size; ++i){
00380     nzvals[i] = Teuchos::as<S>(nzvals_tmp[i]);
00381   }
00382       }
00383     };
00384 #endif  // DOXYGEN_SHOULD_SKIP_THIS
00385 
00398     template<class Matrix, typename S, typename GO, typename GS, class Op>
00399     struct get_cxs_helper
00400     {
00401       static void do_get(const Teuchos::Ptr<const Matrix> mat,
00402        const ArrayView<S> nzvals,
00403        const ArrayView<GO> indices,
00404        const ArrayView<GS> pointers,
00405        GS& nnz, EDistribution distribution,
00406        EStorage_Ordering ordering=ARBITRARY,
00407        GO indexBase = 0)
00408       {
00409   typedef typename Matrix::local_ordinal_t lo_t;
00410   typedef typename Matrix::global_ordinal_t go_t;
00411   typedef typename Matrix::global_size_t gs_t;
00412   typedef typename Matrix::node_t node_t;
00413   
00414   const Teuchos::RCP<const Tpetra::Map<lo_t,go_t,node_t> > map
00415     = getDistributionMap<lo_t,go_t,gs_t,node_t>(distribution,
00416                   Op::get_dimension(mat),
00417                   mat->getComm(), indexBase);
00418   do_get(mat, nzvals, indices, pointers, nnz, Teuchos::ptrInArg(*map), ordering);
00419       }
00420 
00425       static void do_get(const Teuchos::Ptr<const Matrix> mat,
00426        const ArrayView<S> nzvals,
00427        const ArrayView<GO> indices,
00428        const ArrayView<GS> pointers,
00429        GS& nnz, EStorage_Ordering ordering=ARBITRARY)
00430       {
00431   const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
00432                                        typename Matrix::global_ordinal_t,
00433                                        typename Matrix::node_t> > map
00434     = Op::getMap(mat);
00435   do_get(mat, nzvals, indices, pointers, nnz, Teuchos::ptrInArg(*map), ordering);
00436       }
00437 
00442       static void do_get(const Teuchos::Ptr<const Matrix> mat,
00443        const ArrayView<S> nzvals,
00444        const ArrayView<GO> indices,
00445        const ArrayView<GS> pointers,
00446        GS& nnz,
00447        const Teuchos::Ptr<
00448          const Tpetra::Map<typename Matrix::local_ordinal_t,
00449                                              typename Matrix::global_ordinal_t,
00450                            typename Matrix::node_t> > map,
00451        EStorage_Ordering ordering=ARBITRARY)
00452       {
00453   typedef typename Matrix::scalar_t mat_scalar;
00454   if_then_else<is_same<mat_scalar,S>::value,
00455     same_scalar_helper<Matrix,S,GO,GS,Op>,
00456     diff_scalar_helper<Matrix,S,GO,GS,Op> >::type::do_get(mat,
00457                 nzvals, indices,
00458                 pointers, nnz,
00459                 map, ordering);
00460       }
00461     };
00462 
00463 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00464     /*
00465      * These two function-like classes are meant to be used as the \c
00466      * Op template parameter for the \c get_cxs_helper template class.
00467      */
00468     template<class Matrix>
00469     struct get_ccs_func
00470     {
00471       static void apply(const Teuchos::Ptr<const Matrix> mat,
00472       const ArrayView<typename Matrix::scalar_t> nzvals,
00473       const ArrayView<typename Matrix::global_ordinal_t> rowind,
00474       const ArrayView<typename Matrix::global_size_t> colptr,
00475       typename Matrix::global_size_t& nnz,
00476       const Teuchos::Ptr<
00477         const Tpetra::Map<typename Matrix::local_ordinal_t,
00478                                             typename Matrix::global_ordinal_t,
00479                           typename Matrix::node_t> > map,
00480       EStorage_Ordering ordering)
00481       {
00482   mat->getCcs(nzvals, rowind, colptr, nnz, map, ordering);
00483       }
00484 
00485       static
00486       const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
00487              typename Matrix::global_ordinal_t,
00488              typename Matrix::node_t> >
00489       getMap(const Teuchos::Ptr<const Matrix> mat)
00490       {
00491   return mat->getColMap();
00492       }
00493 
00494       static
00495       typename Matrix::global_size_t
00496       get_dimension(const Teuchos::Ptr<const Matrix> mat)
00497       {
00498   return mat->getGlobalNumCols();
00499       }
00500     };
00501 
00502     template<class Matrix>
00503     struct get_crs_func
00504     {
00505       static void apply(const Teuchos::Ptr<const Matrix> mat,
00506       const ArrayView<typename Matrix::scalar_t> nzvals,
00507       const ArrayView<typename Matrix::global_ordinal_t> colind,
00508       const ArrayView<typename Matrix::global_size_t> rowptr,
00509       typename Matrix::global_size_t& nnz,
00510       const Teuchos::Ptr<
00511         const Tpetra::Map<typename Matrix::local_ordinal_t,
00512                           typename Matrix::global_ordinal_t,
00513                           typename Matrix::node_t> > map,
00514       EStorage_Ordering ordering)
00515       {
00516   mat->getCrs(nzvals, colind, rowptr, nnz, map, ordering);
00517       }
00518 
00519       static
00520       const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
00521              typename Matrix::global_ordinal_t,
00522              typename Matrix::node_t> >
00523       getMap(const Teuchos::Ptr<const Matrix> mat)
00524       {
00525   return mat->getRowMap();
00526       }
00527 
00528       static
00529       typename Matrix::global_size_t
00530       get_dimension(const Teuchos::Ptr<const Matrix> mat)
00531       {
00532   return mat->getGlobalNumRows();
00533       }
00534     };
00535 #endif  // DOXYGEN_SHOULD_SKIP_THIS
00536 
00574     template<class Matrix, typename S, typename GO, typename GS>
00575     struct get_ccs_helper : get_cxs_helper<Matrix,S,GO,GS,get_ccs_func<Matrix> >
00576     {};
00577 
00585     template<class Matrix, typename S, typename GO, typename GS>
00586     struct get_crs_helper : get_cxs_helper<Matrix,S,GO,GS,get_crs_func<Matrix> >
00587     {};
00588 
00589     /* End Matrix/MultiVector Utilities */
00590 
00591 
00593     //           Definitions              //
00595 
00596     template <typename LO, typename GO, typename GS, typename Node>
00597     const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
00598     getDistributionMap(EDistribution distribution,
00599            GS num_global_elements,
00600            const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
00601            GO indexBase)
00602     {
00603         // TODO: Need to add indexBase to cases other than ROOTED
00604         //  We do not support these maps in any solver now.
00605       switch( distribution ){
00606       case DISTRIBUTED:
00607       case DISTRIBUTED_NO_OVERLAP:
00608   return Tpetra::createUniformContigMapWithNode<LO,GO, Node>(num_global_elements, comm);
00609   break;
00610       case GLOBALLY_REPLICATED:
00611   return Tpetra::createLocalMapWithNode<LO,GO, Node>(num_global_elements, comm);
00612   break;
00613       case ROOTED:
00614   {
00615     int rank = Teuchos::rank(*comm);
00616     size_t my_num_elems = Teuchos::OrdinalTraits<size_t>::zero();
00617     if( rank == 0 ) my_num_elems = num_global_elements;
00618     return rcp(new Tpetra::Map<LO,GO, Node>(num_global_elements,
00619             my_num_elems, indexBase, comm));
00620     break;
00621   }
00622       default:
00623   TEUCHOS_TEST_FOR_EXCEPTION( true,
00624           std::logic_error,
00625           "Control should never reach this point.  "
00626           "Please contact the Amesos2 developers." );
00627   break;
00628       }
00629     }
00630 
00631 #ifdef HAVE_AMESOS2_EPETRA
00632     template <typename LO, typename GO, typename GS, typename Node>
00633     Teuchos::RCP<Tpetra::Map<LO,GO,Node> >
00634     epetra_map_to_tpetra_map(const Epetra_BlockMap& map)
00635     {
00636       using Teuchos::as;
00637       using Teuchos::rcp;
00638 
00639       int num_my_elements = map.NumMyElements();
00640       Teuchos::Array<int> my_global_elements(num_my_elements);
00641       map.MyGlobalElements(my_global_elements.getRawPtr());
00642 
00643       typedef Tpetra::Map<LO,GO,Node> map_t;
00644       RCP<map_t> tmap = rcp(new map_t(Teuchos::OrdinalTraits<GS>::invalid(),
00645               my_global_elements(),
00646               as<GO>(map.IndexBase()),
00647               to_teuchos_comm(Teuchos::rcpFromRef(map.Comm()))));
00648       return tmap;
00649     }
00650 
00651     template <typename LO, typename GO, typename GS, typename Node>
00652     Teuchos::RCP<Epetra_Map>
00653     tpetra_map_to_epetra_map(const Tpetra::Map<LO,GO,Node>& map)
00654     {
00655       using Teuchos::as;
00656       
00657       Teuchos::Array<GO> elements_tmp;
00658       elements_tmp = map.getNodeElementList();
00659       int num_my_elements = elements_tmp.size();
00660       Teuchos::Array<int> my_global_elements(num_my_elements);
00661       for (int i = 0; i < num_my_elements; ++i){
00662   my_global_elements[i] = as<int>(elements_tmp[i]);
00663       }
00664 
00665       using Teuchos::rcp;
00666       RCP<Epetra_Map> emap = rcp(new Epetra_Map(-1,
00667             num_my_elements,
00668             my_global_elements.getRawPtr(),
00669             as<GO>(map.getIndexBase()),
00670             *to_epetra_comm(map.getComm())));
00671       return emap;
00672     }
00673 #endif  // HAVE_AMESOS2_EPETRA
00674 
00675     template <typename Scalar,
00676         typename GlobalOrdinal,
00677         typename GlobalSizeT>
00678     void transpose(Teuchos::ArrayView<Scalar> vals,
00679        Teuchos::ArrayView<GlobalOrdinal> indices,
00680        Teuchos::ArrayView<GlobalSizeT> ptr,
00681        Teuchos::ArrayView<Scalar> trans_vals,
00682        Teuchos::ArrayView<GlobalOrdinal> trans_indices,
00683        Teuchos::ArrayView<GlobalSizeT> trans_ptr)
00684     {
00685       /* We have a compressed-row storage format of this matrix.  We
00686        * transform this into a compressed-column format using a
00687        * distribution-counting sort algorithm, which is described by
00688        * D.E. Knuth in TAOCP Vol 3, 2nd ed pg 78.
00689        */
00690 
00691 #ifdef HAVE_AMESOS2_DEBUG
00692       typename Teuchos::ArrayView<GlobalOrdinal>::iterator ind_it, ind_begin, ind_end;
00693       ind_begin = indices.begin();
00694       ind_end = indices.end();
00695       size_t min_trans_ptr_size = *std::max_element(ind_begin, ind_end) + 1;
00696       TEUCHOS_TEST_FOR_EXCEPTION( Teuchos::as<size_t>(trans_ptr.size()) < min_trans_ptr_size,
00697         std::invalid_argument,
00698         "Transpose pointer size not large enough." );
00699       TEUCHOS_TEST_FOR_EXCEPTION( trans_vals.size() < vals.size(),
00700         std::invalid_argument,
00701         "Transpose values array not large enough." );
00702       TEUCHOS_TEST_FOR_EXCEPTION( trans_indices.size() < indices.size(),
00703         std::invalid_argument,
00704         "Transpose indices array not large enough." );
00705 #else
00706       typename Teuchos::ArrayView<GlobalOrdinal>::iterator ind_it, ind_end;
00707 #endif
00708 
00709       // Count the number of entries in each column
00710       Teuchos::Array<GlobalSizeT> count(trans_ptr.size(), 0);
00711       ind_end = indices.end();
00712       for( ind_it = indices.begin(); ind_it != ind_end; ++ind_it ){
00713   ++(count[(*ind_it) + 1]);
00714       }
00715 
00716       // Accumulate
00717       typename Teuchos::Array<GlobalSizeT>::iterator cnt_it, cnt_end;
00718       cnt_end = count.end();
00719       for( cnt_it = count.begin() + 1; cnt_it != cnt_end; ++cnt_it ){
00720   *cnt_it = *cnt_it + *(cnt_it - 1);
00721       }
00722       // This becomes the array of column pointers
00723       trans_ptr.assign(count);
00724 
00725       /* Move the nonzero values into their final place in nzval, based on the
00726        * counts found previously.
00727        *
00728        * This sequence deviates from Knuth's algorithm a bit, following more
00729        * closely the description presented in Gustavson, Fred G. "Two Fast
00730        * Algorithms for Sparse Matrices: Multiplication and Permuted
00731        * Transposition" ACM Trans. Math. Softw. volume 4, number 3, 1978, pages
00732        * 250--269, http://doi.acm.org/10.1145/355791.355796.
00733        *
00734        * The output indices end up in sorted order
00735        */
00736 
00737       GlobalSizeT size = ptr.size();
00738       for( GlobalSizeT i = 0; i < size - 1; ++i ){
00739   GlobalOrdinal u = ptr[i];
00740   GlobalOrdinal v = ptr[i + 1];
00741   for( GlobalOrdinal j = u; j < v; ++j ){
00742     GlobalOrdinal k = count[indices[j]];
00743     trans_vals[k] = vals[j];
00744     trans_indices[k] = i;
00745     ++(count[indices[j]]);
00746   }
00747       }
00748     }
00749 
00750 
00751     template <typename Scalar1, typename Scalar2>
00752     void
00753     scale(Teuchos::ArrayView<Scalar1> vals, size_t l,
00754     size_t ld, Teuchos::ArrayView<Scalar2> s)
00755     {
00756       size_t vals_size = vals.size();
00757 #ifdef HAVE_AMESOS2_DEBUG
00758       size_t s_size = s.size();
00759       TEUCHOS_TEST_FOR_EXCEPTION( s_size < l,
00760         std::invalid_argument,
00761         "Scale vector must have length at least that of the vector" );
00762 #endif
00763       size_t i, s_i;
00764       for( i = 0, s_i = 0; i < vals_size; ++i, ++s_i ){
00765   if( s_i == l ){
00766     // bring i to the next multiple of ld
00767     i += ld - s_i;
00768     s_i = 0;
00769   }
00770   vals[i] *= s[s_i];
00771       }
00772     }
00773 
00774     template <typename Scalar1, typename Scalar2, class BinaryOp>
00775     void
00776     scale(Teuchos::ArrayView<Scalar1> vals, size_t l,
00777     size_t ld, Teuchos::ArrayView<Scalar2> s,
00778     BinaryOp binary_op)
00779     {
00780       size_t vals_size = vals.size();
00781 #ifdef HAVE_AMESOS2_DEBUG
00782       size_t s_size = s.size();
00783       TEUCHOS_TEST_FOR_EXCEPTION( s_size < l,
00784         std::invalid_argument,
00785         "Scale vector must have length at least that of the vector" );
00786 #endif
00787       size_t i, s_i;
00788       for( i = 0, s_i = 0; i < vals_size; ++i, ++s_i ){
00789   if( s_i == l ){
00790     // bring i to the next multiple of ld
00791     i += ld - s_i;
00792     s_i = 0;
00793   }
00794   vals[i] = binary_op(vals[i], s[s_i]);
00795       }
00796     }
00797 
00800   } // end namespace Util
00801 
00802 } // end namespace Amesos2
00803 
00804 #endif  // #ifndef AMESOS2_UTIL_HPP