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     
00109 
00110 #ifdef HAVE_AMESOS2_EPETRA
00111 
00116     template <typename LO, typename GO, typename GS, typename Node>
00117     RCP<Tpetra::Map<LO,GO,Node> >
00118     epetra_map_to_tpetra_map(const Epetra_BlockMap& map);
00119 
00125     template <typename LO, typename GO, typename GS, typename Node>
00126     RCP<Epetra_Map>
00127     tpetra_map_to_epetra_map(const Tpetra::Map<LO,GO,Node>& map);
00128 
00134     const RCP<const Teuchos::Comm<int> > to_teuchos_comm(RCP<const Epetra_Comm> c);
00135 
00141     const RCP<const Epetra_Comm> to_epetra_comm(RCP<const Teuchos::Comm<int> > c);
00142 #endif  // HAVE_AMESOS2_EPETRA
00143 
00149     template <typename Scalar,
00150         typename GlobalOrdinal,
00151         typename GlobalSizeT>
00152     void transpose(ArrayView<Scalar> vals,
00153        ArrayView<GlobalOrdinal> indices,
00154        ArrayView<GlobalSizeT> ptr,
00155        ArrayView<Scalar> trans_vals,
00156        ArrayView<GlobalOrdinal> trans_indices,
00157        ArrayView<GlobalSizeT> trans_ptr);
00158 
00172     template <typename Scalar1, typename Scalar2>
00173     void scale(ArrayView<Scalar1> vals, size_t l,
00174          size_t ld, ArrayView<Scalar2> s);
00175 
00194     template <typename Scalar1, typename Scalar2, class BinaryOp>
00195     void scale(ArrayView<Scalar1> vals, size_t l,
00196          size_t ld, ArrayView<Scalar2> s, BinaryOp binary_op);
00197 
00198 
00200     void printLine( Teuchos::FancyOStream &out );
00201 
00202 
00204     // Matrix/MultiVector Utilities //
00206     
00207 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00208     /*
00209      * The following represents a general way of getting a CRS or CCS
00210      * representation of a matrix with implicit type conversions.  The
00211      * \c get_crs_helper and \c get_ccs_helper classes are templated
00212      * on 4 types:
00213      *
00214      * - A Matrix type (conforming to the Amesos2 MatrixAdapter interface)
00215      * - A scalar type
00216      * - A global ordinal type
00217      * - A global size type
00218      *
00219      * The last three template types correspond to the input argument
00220      * types.  For example, if the scalar type is \c double , then we
00221      * require that the \c nzvals argument is a \c
00222      * Teuchos::ArrayView<double> type.
00223      *
00224      * These helpers perform any type conversions that must be
00225      * performed to go between the Matrix's types and the input types.
00226      * If no conversions are necessary the static functions can be
00227      * effectively inlined.
00228      */
00229 
00230     template <class M, typename S, typename GO, typename GS, class Op>
00231     struct same_gs_helper
00232     {
00233       static void do_get(const Teuchos::Ptr<const M> mat,
00234        const ArrayView<typename M::scalar_t> nzvals,
00235        const ArrayView<typename M::global_ordinal_t> indices,
00236        const ArrayView<GS> pointers,
00237        GS& nnz,
00238        const Teuchos::Ptr<
00239          const Tpetra::Map<typename M::local_ordinal_t,
00240                                              typename M::global_ordinal_t,
00241                            typename M::node_t> > map,
00242        EStorage_Ordering ordering)
00243       {
00244   Op::apply(mat, nzvals, indices, pointers, nnz, map, ordering);
00245       }
00246     };
00247 
00248     template <class M, typename S, typename GO, typename GS, class Op>
00249     struct diff_gs_helper
00250     {
00251       static void do_get(const Teuchos::Ptr<const M> mat,
00252        const ArrayView<typename M::scalar_t> nzvals,
00253        const ArrayView<typename M::global_ordinal_t> indices,
00254        const ArrayView<GS> pointers,
00255        GS& nnz,
00256        const Teuchos::Ptr<
00257          const Tpetra::Map<typename M::local_ordinal_t,
00258                                              typename M::global_ordinal_t,
00259                            typename M::node_t> > map,
00260        EStorage_Ordering ordering)
00261       {
00262   typedef typename M::global_size_t mat_gs_t;
00263   typename ArrayView<GS>::size_type i, size = pointers.size();
00264   Teuchos::Array<mat_gs_t> pointers_tmp(size);
00265   mat_gs_t nnz_tmp = 0;
00266 
00267   Op::apply(mat, nzvals, indices, pointers_tmp, nnz_tmp, map, ordering);
00268 
00269   for (i = 0; i < size; ++i){
00270     pointers[i] = Teuchos::as<GS>(pointers_tmp[i]);
00271   }
00272   nnz = Teuchos::as<GS>(nnz_tmp);
00273       }
00274     };
00275 
00276     template <class M, typename S, typename GO, typename GS, class Op>
00277     struct same_go_helper
00278     {
00279       static void do_get(const Teuchos::Ptr<const M> mat,
00280        const ArrayView<typename M::scalar_t> nzvals,
00281        const ArrayView<GO> indices,
00282        const ArrayView<GS> pointers,
00283        GS& nnz,
00284        const Teuchos::Ptr<
00285          const Tpetra::Map<typename M::local_ordinal_t,
00286                                              typename M::global_ordinal_t,
00287                            typename M::node_t> > map,
00288        EStorage_Ordering ordering)
00289       {
00290   typedef typename M::global_size_t mat_gs_t;
00291   if_then_else<is_same<GS,mat_gs_t>::value,
00292     same_gs_helper<M,S,GO,GS,Op>,
00293     diff_gs_helper<M,S,GO,GS,Op> >::type::do_get(mat, nzvals, indices,
00294                    pointers, nnz, map,
00295                    ordering);
00296       }
00297     };
00298 
00299     template <class M, typename S, typename GO, typename GS, class Op>
00300     struct diff_go_helper
00301     {
00302       static void do_get(const Teuchos::Ptr<const M> mat,
00303        const ArrayView<typename M::scalar_t> nzvals,
00304        const ArrayView<GO> indices,
00305        const ArrayView<GS> pointers,
00306        GS& nnz,
00307        const Teuchos::Ptr<
00308          const Tpetra::Map<typename M::local_ordinal_t,
00309                                              typename M::global_ordinal_t,
00310                            typename M::node_t> > map,
00311        EStorage_Ordering ordering)
00312       {
00313   typedef typename M::global_ordinal_t mat_go_t;
00314   typedef typename M::global_size_t mat_gs_t;
00315   typename ArrayView<GO>::size_type i, size = indices.size();
00316   Teuchos::Array<mat_go_t> indices_tmp(size);
00317 
00318   if_then_else<is_same<GS,mat_gs_t>::value,
00319     same_gs_helper<M,S,GO,GS,Op>,
00320     diff_gs_helper<M,S,GO,GS,Op> >::type::do_get(mat, nzvals, indices_tmp,
00321                    pointers, nnz, map,
00322                    ordering);
00323 
00324   for (i = 0; i < size; ++i){
00325     indices[i] = Teuchos::as<GO>(indices_tmp[i]);
00326   }
00327       }
00328     };
00329 
00330     template <class M, typename S, typename GO, typename GS, class Op>
00331     struct same_scalar_helper
00332     {
00333       static void do_get(const Teuchos::Ptr<const M> mat,
00334        const ArrayView<S> nzvals,
00335        const ArrayView<GO> indices,
00336        const ArrayView<GS> pointers,
00337        GS& nnz,
00338        const Teuchos::Ptr<
00339          const Tpetra::Map<typename M::local_ordinal_t,
00340                                              typename M::global_ordinal_t,
00341                            typename M::node_t> > map,
00342        EStorage_Ordering ordering)
00343       {
00344   typedef typename M::global_ordinal_t mat_go_t;
00345   if_then_else<is_same<GO,mat_go_t>::value,
00346     same_go_helper<M,S,GO,GS,Op>,
00347     diff_go_helper<M,S,GO,GS,Op> >::type::do_get(mat, nzvals, indices,
00348                    pointers, nnz, map,
00349                    ordering);
00350       }
00351     };
00352 
00353     template <class M, typename S, typename GO, typename GS, class Op>
00354     struct diff_scalar_helper
00355     {
00356       static void do_get(const Teuchos::Ptr<const M> mat,
00357        const ArrayView<S> nzvals,
00358        const ArrayView<GO> indices,
00359        const ArrayView<GS> pointers,
00360        GS& nnz,
00361        const Teuchos::Ptr<
00362          const Tpetra::Map<typename M::local_ordinal_t,
00363                                              typename M::global_ordinal_t,
00364                            typename M::node_t> > map,
00365        EStorage_Ordering ordering)
00366       {
00367   typedef typename M::scalar_t mat_scalar_t;
00368   typedef typename M::global_ordinal_t mat_go_t;
00369   typename ArrayView<S>::size_type i, size = nzvals.size();
00370   Teuchos::Array<mat_scalar_t> nzvals_tmp(size);
00371 
00372   if_then_else<is_same<GO,mat_go_t>::value,
00373     same_go_helper<M,S,GO,GS,Op>,
00374     diff_go_helper<M,S,GO,GS,Op> >::type::do_get(mat, nzvals_tmp, indices,
00375                    pointers, nnz, map,
00376                    ordering);
00377 
00378   for (i = 0; i < size; ++i){
00379     nzvals[i] = Teuchos::as<S>(nzvals_tmp[i]);
00380   }
00381       }
00382     };
00383 #endif  // DOXYGEN_SHOULD_SKIP_THIS
00384 
00397     template<class Matrix, typename S, typename GO, typename GS, class Op>
00398     struct get_cxs_helper
00399     {
00400       static void do_get(const Teuchos::Ptr<const Matrix> mat,
00401        const ArrayView<S> nzvals,
00402        const ArrayView<GO> indices,
00403        const ArrayView<GS> pointers,
00404        GS& nnz, EDistribution distribution,
00405        EStorage_Ordering ordering=ARBITRARY)
00406       {
00407   typedef typename Matrix::local_ordinal_t lo_t;
00408   typedef typename Matrix::global_ordinal_t go_t;
00409   typedef typename Matrix::global_size_t gs_t;
00410   typedef typename Matrix::node_t node_t;
00411   
00412   const Teuchos::RCP<const Tpetra::Map<lo_t,go_t,node_t> > map
00413     = getDistributionMap<lo_t,go_t,gs_t,node_t>(distribution,
00414                   Op::get_dimension(mat),
00415                   mat->getComm());
00416   do_get(mat, nzvals, indices, pointers, nnz, Teuchos::ptrInArg(*map), ordering);
00417       }
00418 
00423       static void do_get(const Teuchos::Ptr<const Matrix> mat,
00424        const ArrayView<S> nzvals,
00425        const ArrayView<GO> indices,
00426        const ArrayView<GS> pointers,
00427        GS& nnz, EStorage_Ordering ordering=ARBITRARY)
00428       {
00429   const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
00430                                        typename Matrix::global_ordinal_t,
00431                                        typename Matrix::node_t> > map
00432     = Op::getMap(mat);
00433   do_get(mat, nzvals, indices, pointers, nnz, Teuchos::ptrInArg(*map), ordering);
00434       }
00435 
00440       static void do_get(const Teuchos::Ptr<const Matrix> mat,
00441        const ArrayView<S> nzvals,
00442        const ArrayView<GO> indices,
00443        const ArrayView<GS> pointers,
00444        GS& nnz,
00445        const Teuchos::Ptr<
00446          const Tpetra::Map<typename Matrix::local_ordinal_t,
00447                                              typename Matrix::global_ordinal_t,
00448                            typename Matrix::node_t> > map,
00449        EStorage_Ordering ordering=ARBITRARY)
00450       {
00451   typedef typename Matrix::scalar_t mat_scalar;
00452   if_then_else<is_same<mat_scalar,S>::value,
00453     same_scalar_helper<Matrix,S,GO,GS,Op>,
00454     diff_scalar_helper<Matrix,S,GO,GS,Op> >::type::do_get(mat,
00455                 nzvals, indices,
00456                 pointers, nnz,
00457                 map, ordering);
00458       }
00459     };
00460 
00461 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00462     /*
00463      * These two function-like classes are meant to be used as the \c
00464      * Op template parameter for the \c get_cxs_helper template class.
00465      */
00466     template<class Matrix>
00467     struct get_ccs_func
00468     {
00469       static void apply(const Teuchos::Ptr<const Matrix> mat,
00470       const ArrayView<typename Matrix::scalar_t> nzvals,
00471       const ArrayView<typename Matrix::global_ordinal_t> rowind,
00472       const ArrayView<typename Matrix::global_size_t> colptr,
00473       typename Matrix::global_size_t& nnz,
00474       const Teuchos::Ptr<
00475         const Tpetra::Map<typename Matrix::local_ordinal_t,
00476                                             typename Matrix::global_ordinal_t,
00477                           typename Matrix::node_t> > map,
00478       EStorage_Ordering ordering)
00479       {
00480   mat->getCcs(nzvals, rowind, colptr, nnz, map, ordering);
00481       }
00482 
00483       static
00484       const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
00485              typename Matrix::global_ordinal_t,
00486              typename Matrix::node_t> >
00487       getMap(const Teuchos::Ptr<const Matrix> mat)
00488       {
00489   return mat->getColMap();
00490       }
00491 
00492       static
00493       typename Matrix::global_size_t
00494       get_dimension(const Teuchos::Ptr<const Matrix> mat)
00495       {
00496   return mat->getGlobalNumCols();
00497       }
00498     };
00499 
00500     template<class Matrix>
00501     struct get_crs_func
00502     {
00503       static void apply(const Teuchos::Ptr<const Matrix> mat,
00504       const ArrayView<typename Matrix::scalar_t> nzvals,
00505       const ArrayView<typename Matrix::global_ordinal_t> colind,
00506       const ArrayView<typename Matrix::global_size_t> rowptr,
00507       typename Matrix::global_size_t& nnz,
00508       const Teuchos::Ptr<
00509         const Tpetra::Map<typename Matrix::local_ordinal_t,
00510                           typename Matrix::global_ordinal_t,
00511                           typename Matrix::node_t> > map,
00512       EStorage_Ordering ordering)
00513       {
00514   mat->getCrs(nzvals, colind, rowptr, nnz, map, ordering);
00515       }
00516 
00517       static
00518       const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
00519              typename Matrix::global_ordinal_t,
00520              typename Matrix::node_t> >
00521       getMap(const Teuchos::Ptr<const Matrix> mat)
00522       {
00523   return mat->getRowMap();
00524       }
00525 
00526       static
00527       typename Matrix::global_size_t
00528       get_dimension(const Teuchos::Ptr<const Matrix> mat)
00529       {
00530   return mat->getGlobalNumRows();
00531       }
00532     };
00533 #endif  // DOXYGEN_SHOULD_SKIP_THIS
00534 
00572     template<class Matrix, typename S, typename GO, typename GS>
00573     struct get_ccs_helper : get_cxs_helper<Matrix,S,GO,GS,get_ccs_func<Matrix> >
00574     {};
00575 
00583     template<class Matrix, typename S, typename GO, typename GS>
00584     struct get_crs_helper : get_cxs_helper<Matrix,S,GO,GS,get_crs_func<Matrix> >
00585     {};
00586 
00587     /* End Matrix/MultiVector Utilities */
00588 
00589 
00591     //           Definitions              //
00593 
00594     template <typename LO, typename GO, typename GS, typename Node>
00595     const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
00596     getDistributionMap(EDistribution distribution,
00597            GS num_global_elements,
00598            const Teuchos::RCP<const Teuchos::Comm<int> >& comm)
00599     {
00600       switch( distribution ){
00601       case DISTRIBUTED:
00602       case DISTRIBUTED_NO_OVERLAP:
00603   return Tpetra::createUniformContigMap<LO,GO>(num_global_elements, comm);
00604   break;
00605       case GLOBALLY_REPLICATED:
00606   return Tpetra::createLocalMap<LO,GO>(num_global_elements, comm);
00607   break;
00608       case ROOTED:
00609   {
00610     int rank = Teuchos::rank(*comm);
00611     size_t my_num_elems = Teuchos::OrdinalTraits<size_t>::zero();
00612     if( rank == 0 ) my_num_elems = num_global_elements;
00613     return Tpetra::createContigMap<LO,GO>(num_global_elements,
00614             my_num_elems, comm);
00615     break;
00616   }
00617       default:
00618   TEUCHOS_TEST_FOR_EXCEPTION( true,
00619           std::logic_error,
00620           "Control should never reach this point.  "
00621           "Please contact the Amesos2 developers." );
00622   break;
00623       }
00624     }
00625 
00626 #ifdef HAVE_AMESOS2_EPETRA
00627     template <typename LO, typename GO, typename GS, typename Node>
00628     Teuchos::RCP<Tpetra::Map<LO,GO,Node> >
00629     epetra_map_to_tpetra_map(const Epetra_BlockMap& map)
00630     {
00631       using Teuchos::as;
00632       using Teuchos::rcp;
00633 
00634       int num_my_elements = map.NumMyElements();
00635       Teuchos::Array<int> my_global_elements(num_my_elements);
00636       map.MyGlobalElements(my_global_elements.getRawPtr());
00637 
00638       typedef Tpetra::Map<LO,GO,Node> map_t;
00639       RCP<map_t> tmap = rcp(new map_t(Teuchos::OrdinalTraits<GS>::invalid(),
00640               my_global_elements(),
00641               as<GO>(map.IndexBase()),
00642               to_teuchos_comm(Teuchos::rcpFromRef(map.Comm()))));
00643       return tmap;
00644     }
00645 
00646     template <typename LO, typename GO, typename GS, typename Node>
00647     Teuchos::RCP<Epetra_Map>
00648     tpetra_map_to_epetra_map(const Tpetra::Map<LO,GO,Node>& map)
00649     {
00650       using Teuchos::as;
00651       
00652       Teuchos::Array<GO> elements_tmp;
00653       elements_tmp = map.getNodeElementList();
00654       int num_my_elements = elements_tmp.size();
00655       Teuchos::Array<int> my_global_elements(num_my_elements);
00656       for (int i = 0; i < num_my_elements; ++i){
00657   my_global_elements[i] = as<int>(elements_tmp[i]);
00658       }
00659 
00660       using Teuchos::rcp;
00661       RCP<Epetra_Map> emap = rcp(new Epetra_Map(-1,
00662             num_my_elements,
00663             my_global_elements.getRawPtr(),
00664             as<GO>(map.getIndexBase()),
00665             *to_epetra_comm(map.getComm())));
00666       return emap;
00667     }
00668 #endif  // HAVE_AMESOS2_EPETRA
00669 
00670     template <typename Scalar,
00671         typename GlobalOrdinal,
00672         typename GlobalSizeT>
00673     void transpose(Teuchos::ArrayView<Scalar> vals,
00674        Teuchos::ArrayView<GlobalOrdinal> indices,
00675        Teuchos::ArrayView<GlobalSizeT> ptr,
00676        Teuchos::ArrayView<Scalar> trans_vals,
00677        Teuchos::ArrayView<GlobalOrdinal> trans_indices,
00678        Teuchos::ArrayView<GlobalSizeT> trans_ptr)
00679     {
00680       /* We have a compressed-row storage format of this matrix.  We
00681        * transform this into a compressed-column format using a
00682        * distribution-counting sort algorithm, which is described by
00683        * D.E. Knuth in TAOCP Vol 3, 2nd ed pg 78.
00684        */
00685 
00686 #ifdef HAVE_AMESOS2_DEBUG
00687       typename Teuchos::ArrayView<GlobalOrdinal>::iterator ind_it, ind_begin, ind_end;
00688       ind_begin = indices.begin();
00689       ind_end = indices.end();
00690       size_t min_trans_ptr_size = *std::max_element(ind_begin, ind_end) + 1;
00691       TEUCHOS_TEST_FOR_EXCEPTION( Teuchos::as<size_t>(trans_ptr.size()) < min_trans_ptr_size,
00692         std::invalid_argument,
00693         "Transpose pointer size not large enough." );
00694       TEUCHOS_TEST_FOR_EXCEPTION( trans_vals.size() < vals.size(),
00695         std::invalid_argument,
00696         "Transpose values array not large enough." );
00697       TEUCHOS_TEST_FOR_EXCEPTION( trans_indices.size() < indices.size(),
00698         std::invalid_argument,
00699         "Transpose indices array not large enough." );
00700 #else
00701       typename Teuchos::ArrayView<GlobalOrdinal>::iterator ind_it, ind_end;
00702 #endif
00703 
00704       // Count the number of entries in each column
00705       Teuchos::Array<GlobalSizeT> count(trans_ptr.size(), 0);
00706       ind_end = indices.end();
00707       for( ind_it = indices.begin(); ind_it != ind_end; ++ind_it ){
00708   ++(count[(*ind_it) + 1]);
00709       }
00710 
00711       // Accumulate
00712       typename Teuchos::Array<GlobalSizeT>::iterator cnt_it, cnt_end;
00713       cnt_end = count.end();
00714       for( cnt_it = count.begin() + 1; cnt_it != cnt_end; ++cnt_it ){
00715   *cnt_it = *cnt_it + *(cnt_it - 1);
00716       }
00717       // This becomes the array of column pointers
00718       trans_ptr.assign(count);
00719 
00720       /* Move the nonzero values into their final place in nzval, based on the
00721        * counts found previously.
00722        *
00723        * This sequence deviates from Knuth's algorithm a bit, following more
00724        * closely the description presented in Gustavson, Fred G. "Two Fast
00725        * Algorithms for Sparse Matrices: Multiplication and Permuted
00726        * Transposition" ACM Trans. Math. Softw. volume 4, number 3, 1978, pages
00727        * 250--269, http://doi.acm.org/10.1145/355791.355796.
00728        *
00729        * The output indices end up in sorted order
00730        */
00731 
00732       GlobalSizeT size = ptr.size();
00733       for( GlobalSizeT i = 0; i < size - 1; ++i ){
00734   GlobalOrdinal u = ptr[i];
00735   GlobalOrdinal v = ptr[i + 1];
00736   for( GlobalOrdinal j = u; j < v; ++j ){
00737     GlobalOrdinal k = count[indices[j]];
00738     trans_vals[k] = vals[j];
00739     trans_indices[k] = i;
00740     ++(count[indices[j]]);
00741   }
00742       }
00743     }
00744 
00745 
00746     template <typename Scalar1, typename Scalar2>
00747     void
00748     scale(Teuchos::ArrayView<Scalar1> vals, size_t l,
00749     size_t ld, Teuchos::ArrayView<Scalar2> s)
00750     {
00751       size_t vals_size = vals.size();
00752 #ifdef HAVE_AMESOS2_DEBUG
00753       size_t s_size = s.size();
00754       TEUCHOS_TEST_FOR_EXCEPTION( s_size < l,
00755         std::invalid_argument,
00756         "Scale vector must have length at least that of the vector" );
00757 #endif
00758       size_t i, s_i;
00759       for( i = 0, s_i = 0; i < vals_size; ++i, ++s_i ){
00760   if( s_i == l ){
00761     // bring i to the next multiple of ld
00762     i += ld - s_i;
00763     s_i = 0;
00764   }
00765   vals[i] *= s[s_i];
00766       }
00767     }
00768 
00769     template <typename Scalar1, typename Scalar2, class BinaryOp>
00770     void
00771     scale(Teuchos::ArrayView<Scalar1> vals, size_t l,
00772     size_t ld, Teuchos::ArrayView<Scalar2> s,
00773     BinaryOp binary_op)
00774     {
00775       size_t vals_size = vals.size();
00776 #ifdef HAVE_AMESOS2_DEBUG
00777       size_t s_size = s.size();
00778       TEUCHOS_TEST_FOR_EXCEPTION( s_size < l,
00779         std::invalid_argument,
00780         "Scale vector must have length at least that of the vector" );
00781 #endif
00782       size_t i, s_i;
00783       for( i = 0, s_i = 0; i < vals_size; ++i, ++s_i ){
00784   if( s_i == l ){
00785     // bring i to the next multiple of ld
00786     i += ld - s_i;
00787     s_i = 0;
00788   }
00789   vals[i] = binary_op(vals[i], s[s_i]);
00790       }
00791     }
00792 
00795   } // end namespace Util
00796 
00797 } // end namespace Amesos2
00798 
00799 #endif  // #ifndef AMESOS2_UTIL_HPP