RTOpPack_LapackWrappers.hpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 // RTOp: Interfaces and Support Software for Vector Reduction Transformation
00005 //       Operations
00006 //                Copyright (2006) Sandia Corporation
00007 // 
00008 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00009 // license for use of this work by or on behalf of the U.S. Government.
00010 // 
00011 // This library is free software; you can redistribute it and/or modify
00012 // it under the terms of the GNU Lesser General Public License as
00013 // published by the Free Software Foundation; either version 2.1 of the
00014 // License, or (at your option) any later version.
00015 //  
00016 // This library is distributed in the hope that it will be useful, but
00017 // WITHOUT ANY WARRANTY; without even the implied warranty of
00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00019 // Lesser General Public License for more details.
00020 //  
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License along with this library; if not, write to the Free Software
00023 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00024 // USA
00025 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00026 // 
00027 // ***********************************************************************
00028 // @HEADER
00029 
00030 #ifndef RTOPPACK_LAPACK_WRAPPERS_HPP
00031 #define RTOPPACK_LAPACK_WRAPPERS_HPP
00032 
00033 
00034 #include "RTOpPack_Types.hpp"
00035 #include "Teuchos_LAPACK.hpp"
00036 #include "Teuchos_as.hpp"
00037 
00038 
00039 namespace RTOpPack {
00040 
00041 
00043 enum ETransp {
00044   NOTRANS, TRANS, CONJTRANS
00045 };
00046 const int NUM_ETRANS_ARGS = 3;
00047 
00049 extern const Teuchos::Tuple<char,NUM_ETRANS_ARGS> transpMap;
00050 
00051 
00062 template<class Scalar>
00063 void getrf(
00064   const SubMultiVectorView<Scalar> &A,
00065   const ArrayView<int> &ipiv,
00066   const Ptr<int> &rank
00067   );
00068 
00069 
00071 template<class Scalar>
00072 void getrs(
00073   const ConstSubMultiVectorView<Scalar> &A,
00074   const ArrayView<const int> &ipiv,
00075   const ETransp transp,
00076   const Ptr<const SubMultiVectorView<Scalar> > &BX
00077   );
00078 
00079 
00080 } // namespace RTOpPack
00081 
00082 
00083 //
00084 // Implementations
00085 //
00086 
00087 
00088 template<class Scalar>
00089 void RTOpPack::getrf(
00090   const SubMultiVectorView<Scalar> &A,
00091   const ArrayView<int> &ipiv,
00092   const Ptr<int> &rank
00093   )
00094 {
00095   using Teuchos::as;
00096   const int maxRank = TEUCHOS_MIN( A.subDim(), A.numSubCols() );
00097 #ifdef TEUCHOS_DEBUG
00098   TEST_FOR_EXCEPT( A.subDim() == 0  );
00099   TEST_FOR_EXCEPT( A.numSubCols() == 0  );
00100   TEST_FOR_EXCEPT( is_null(A.values()) );
00101   TEUCHOS_ASSERT_EQUALITY( as<int>(ipiv.size()), maxRank );
00102 #endif
00103 
00104   Teuchos::LAPACK<int, Scalar> lapack;
00105   int info = -1;
00106   lapack.GETRF( A.subDim(), A.numSubCols(), A.values().get(), A.leadingDim(),
00107     &ipiv[0], &info );
00108   *rank = maxRank;
00109   TEST_FOR_EXCEPTION(
00110     info < 0, std::invalid_argument
00111     ,"getrf(...): Error, Invalid argument "
00112     << -info << " sent to LAPACK function xGETRF(...)" );
00113   if (info > 0)
00114     *rank = info - 1;
00115 }
00116 
00117 
00118 template<class Scalar>
00119 void RTOpPack::getrs(
00120   const ConstSubMultiVectorView<Scalar> &A,
00121   const ArrayView<const int> &ipiv,
00122   const ETransp transp,
00123   const Ptr<const SubMultiVectorView<Scalar> > &BX
00124   )
00125 {
00126   using Teuchos::as;
00127 #ifdef TEUCHOS_DEBUG
00128   TEUCHOS_ASSERT( !is_null(BX) );
00129   TEUCHOS_ASSERT_EQUALITY( A.subDim(), BX->subDim() );
00130   TEUCHOS_ASSERT_EQUALITY( A.subDim(), A.numSubCols() );
00131   TEST_FOR_EXCEPT( A.subDim() == 0  );
00132   TEST_FOR_EXCEPT( A.numSubCols() == 0  );
00133   TEST_FOR_EXCEPT( is_null(A.values()) );
00134   TEUCHOS_ASSERT_EQUALITY( A.subDim(), ipiv.size() );
00135 #endif
00136   Teuchos::LAPACK<int, Scalar> lapack;
00137   int info = -1;
00138   lapack.GETRS(
00139     transpMap[transp],
00140     A.subDim(), BX->numSubCols(), A.values().get(), A.leadingDim(),
00141     &ipiv[0], BX->values().get(), BX->leadingDim(), &info
00142     );
00143   TEST_FOR_EXCEPTION(
00144     info < 0, std::invalid_argument
00145     ,"getrs(...): Error, Invalid argument "
00146     << -info << " sent to LAPACK function xGETRS(...)" );
00147   // If we get here B is the solution to the linear system.
00148 }
00149 
00150 
00151 #endif // RTOPPACK_LAPACK_WRAPPERS_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 09:59:06 2011 for RTOp Package Browser (Single Doxygen Collection) by  doxygen 1.6.3