Amesos2 - Direct Sparse Solver Interfaces Version of the Day
Amesos2_Lapack_def.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 
00044 
00053 #ifndef AMESOS2_LAPACK_DEF_HPP
00054 #define AMESOS2_LAPACK_DEF_HPP
00055 
00056 #include <Teuchos_RCP.hpp>
00057 
00058 #include "Amesos2_SolverCore_def.hpp"
00059 #include "Amesos2_Lapack_decl.hpp"
00060 
00061 namespace Amesos2 {
00062 
00063 
00064   template <class Matrix, class Vector>
00065   Lapack<Matrix,Vector>::Lapack(Teuchos::RCP<const Matrix> A,
00066         Teuchos::RCP<Vector>       X,
00067         Teuchos::RCP<const Vector> B)
00068     : SolverCore<Amesos2::Lapack,Matrix,Vector>(A, X, B) // instantiate superclass
00069     , nzvals_()
00070     , rowind_()
00071     , colptr_()
00072   {
00073     // Set default parameters
00074     Teuchos::RCP<Teuchos::ParameterList> default_params
00075       = Teuchos::parameterList( *(this->getValidParameters()) );
00076     this->setParameters(default_params);
00077   }
00078 
00079 
00080   template <class Matrix, class Vector>
00081   Lapack<Matrix,Vector>::~Lapack( )
00082   {
00083     /*
00084      * Free any memory allocated by the Lapack library functions (i.e. none)
00085      */
00086   }
00087 
00088 
00089   template<class Matrix, class Vector>
00090   int
00091   Lapack<Matrix,Vector>::preOrdering_impl()
00092   {
00093     return(0);
00094   }
00095 
00096 
00097   template <class Matrix, class Vector>
00098   int
00099   Lapack<Matrix,Vector>::symbolicFactorization_impl()
00100   {
00101     return(0);
00102   }
00103 
00104 
00105   template <class Matrix, class Vector>
00106   int
00107   Lapack<Matrix,Vector>::numericFactorization_impl()
00108   {
00109     int factor_ierr = 0;
00110 
00111     if( this->root_ ){
00112       // Set here so that solver_ can refresh it's internal state
00113       solver_.setMatrix( Teuchos::rcpFromRef(lu_) );
00114 
00115       {
00116 #ifdef HAVE_AMESOS2_TIMERS
00117   Teuchos::TimeMonitor numFactTimer( this->timers_.numFactTime_ );
00118 #endif
00119   factor_ierr = solver_.factor();
00120       }
00121     }
00122 
00123     Teuchos::broadcast(*(this->getComm()), 0, &factor_ierr);
00124     TEUCHOS_TEST_FOR_EXCEPTION( factor_ierr != 0,
00125         std::runtime_error,
00126         "Lapack factor routine returned error code "
00127         << factor_ierr );
00128     return( 0 );
00129   }
00130 
00131 
00132   template <class Matrix, class Vector>
00133   int
00134   Lapack<Matrix,Vector>::solve_impl(const Teuchos::Ptr<MultiVecAdapter<Vector> >       X,
00135             const Teuchos::Ptr<const MultiVecAdapter<Vector> > B) const
00136   {
00137     using Teuchos::as;
00138 
00139     // Convert X and B to SerialDenseMatrix's
00140     const global_size_type ld_rhs = this->root_ ? X->getGlobalLength() : 0;
00141     const size_t nrhs = X->getGlobalNumVectors();
00142 
00143     const size_t val_store_size = as<size_t>(ld_rhs * nrhs);
00144     if( this->root_ ){
00145       rhsvals_.resize(val_store_size);
00146     }
00147 
00148     {                             // Get values from RHS B
00149 #ifdef HAVE_AMESOS2_TIMERS
00150       Teuchos::TimeMonitor mvConvTimer( this->timers_.vecConvTime_ );
00151       Teuchos::TimeMonitor redistTimer( this->timers_.vecRedistTime_ );
00152 #endif
00153       typedef Util::get_1d_copy_helper<MultiVecAdapter<Vector>,
00154                scalar_type> copy_helper;
00155       copy_helper::do_get(B, rhsvals_(), as<size_t>(ld_rhs), ROOTED);
00156     }
00157 
00158     int solve_ierr         = 0;
00159     int unequilibrate_ierr = 0;
00160 
00161     if( this->root_ ){
00162 #ifdef HAVE_AMESOS2_TIMERS
00163       Teuchos::TimeMonitor solveTimer( this->timers_.solveTime_ );
00164 #endif
00165 
00166       using Teuchos::rcpFromRef;
00167       typedef Teuchos::SerialDenseMatrix<int,scalar_type> DenseMat;
00168 
00169       DenseMat rhs_dense_mat(Teuchos::View, rhsvals_.getRawPtr(),
00170           as<int>(ld_rhs), as<int>(ld_rhs), as<int>(nrhs));
00171 
00172       solver_.setVectors( rcpFromRef(rhs_dense_mat),
00173         rcpFromRef(rhs_dense_mat) );
00174 
00175       solve_ierr = solver_.solve();
00176 
00177       // Solution is found in rhsvals_
00178     }
00179 
00180     // Consolidate and check error codes
00181     Teuchos::broadcast(*(this->getComm()), 0, &solve_ierr);
00182     TEUCHOS_TEST_FOR_EXCEPTION( solve_ierr != 0,
00183         std::runtime_error,
00184         "Lapack solver solve method returned with error code "
00185         << solve_ierr );
00186 
00187     /* Update X's global values */
00188     {
00189 #ifdef HAVE_AMESOS2_TIMERS
00190       Teuchos::TimeMonitor redistTimer( this->timers_.vecRedistTime_ );
00191 #endif
00192 
00193       Util::put_1d_data_helper<
00194       MultiVecAdapter<Vector>,scalar_type>::do_put(X, rhsvals_(),
00195                as<size_t>(ld_rhs),
00196                ROOTED);
00197     }
00198 
00199     return( 0 );
00200   }
00201 
00202 
00203   template <class Matrix, class Vector>
00204   bool
00205   Lapack<Matrix,Vector>::matrixShapeOK_impl() const
00206   {
00207     // Factorization of rectangular matrices is supported, but not
00208     // their solution.  For solution we can have square matrices.
00209 
00210     return( this->globalNumCols_ == this->globalNumRows_ );
00211   }
00212 
00213 
00214   template <class Matrix, class Vector>
00215   void
00216   Lapack<Matrix,Vector>::setParameters_impl(const Teuchos::RCP<Teuchos::ParameterList> & parameterList )
00217   {
00218     solver_.solveWithTranspose( parameterList->get<bool>("Transpose",
00219                this->control_.useTranspose_) );
00220 
00221     solver_.factorWithEquilibration( parameterList->get<bool>("Equilibrate", true) );
00222 
00223     // solver_.solveToRefinedSolution( parameterList->get<bool>("Refine", false) );
00224   }
00225 
00226   template <class Matrix, class Vector>
00227   Teuchos::RCP<const Teuchos::ParameterList>
00228   Lapack<Matrix,Vector>::getValidParameters_impl() const
00229   {
00230     using Teuchos::ParameterList;
00231 
00232     static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
00233 
00234     if( is_null(valid_params) ){
00235       Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
00236 
00237       pl->set("Equilibrate", true, "Whether to equilibrate the input matrix");
00238 
00239       // TODO: Refinement does not seem to be working with the SerialDenseSolver.  Not sure why.
00240       // pl->set("Refine", false, "Whether to apply iterative refinement");
00241 
00242       valid_params = pl;
00243     }
00244 
00245     return valid_params;
00246   }
00247 
00248   template <class Matrix, class Vector>
00249   bool
00250   Lapack<Matrix,Vector>::loadA_impl(EPhase current_phase)
00251   {
00252 #ifdef HAVE_AMESOS2_TIMERS
00253     Teuchos::TimeMonitor convTimer(this->timers_.mtxConvTime_);
00254 #endif
00255 
00256     // We only load the matrix when numericFactorization is called
00257     if( current_phase < NUMFACT ) return( false );
00258 
00259     if( this->root_ ){
00260       nzvals_.resize(this->globalNumNonZeros_);
00261       rowind_.resize(this->globalNumNonZeros_);
00262       colptr_.resize(this->globalNumCols_ + 1);
00263     }
00264 
00265     // global_size_type nnz_ret = 0;
00266     int nnz_ret = 0;
00267     {
00268 #ifdef HAVE_AMESOS2_TIMERS
00269       Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
00270 #endif
00271 
00272       // typedef Util::get_ccs_helper<MatrixAdapter<Matrix>,
00273       //  scalar_type, global_ordinal_type, global_size_type> ccs_helper;
00274       typedef Util::get_ccs_helper<MatrixAdapter<Matrix>,
00275   scalar_type, int, int> ccs_helper;
00276       ccs_helper::do_get(this->matrixA_.ptr(),
00277        nzvals_(), rowind_(), colptr_(),
00278        nnz_ret, ROOTED, SORTED_INDICES);
00279   }
00280 
00281     if( this->root_ ){
00282       // entries are initialized to zero in here:
00283       lu_.shape(this->globalNumRows_, this->globalNumCols_);
00284 
00285       // Put entries of ccs representation into the dense matrix
00286       global_size_type end_col = this->globalNumCols_;
00287       for( global_size_type col = 0; col < end_col; ++col ){
00288   global_ordinal_type ptr = colptr_[col];
00289   global_ordinal_type end_ptr = colptr_[col+1];
00290   for( ; ptr < end_ptr; ++ptr ){
00291     lu_(rowind_[ptr], col) = nzvals_[ptr];
00292   }
00293       }
00294 
00295       // lu_.print(std::cout);
00296     }
00297 
00298   return( true );
00299 }
00300 
00301 
00302   template<class Matrix, class Vector>
00303   const char* Lapack<Matrix,Vector>::name = "LAPACK";
00304 
00305 
00306 } // end namespace Amesos2
00307 
00308 #endif  // AMESOS2_LAPACK_DEF_HPP