Teuchos Package Browser (Single Doxygen Collection) Version of the Day
Teuchos_SerialDenseHelpers.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                    Teuchos: Common Tools Package
00005 //                 Copyright (2004) 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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #ifndef _TEUCHOS_SERIALDENSEHELPERS_HPP_
00043 #define _TEUCHOS_SERIALDENSEHELPERS_HPP_
00044 
00049 #include "Teuchos_ScalarTraits.hpp"
00050 #include "Teuchos_DataAccess.hpp"
00051 #include "Teuchos_ConfigDefs.hpp"
00052 #include "Teuchos_Assert.hpp"
00053 #include "Teuchos_SerialDenseMatrix.hpp"
00054 #include "Teuchos_SerialSymDenseMatrix.hpp"
00055 #include "Teuchos_SerialBandDenseMatrix.hpp"
00056 #include "Teuchos_SerialDenseVector.hpp"
00057 
00058 namespace Teuchos {
00059 
00071 template<typename OrdinalType, typename ScalarType>
00072 void symMatTripleProduct( ETransp transw, const ScalarType alpha, const SerialSymDenseMatrix<OrdinalType, ScalarType>& A,
00073         const SerialDenseMatrix<OrdinalType, ScalarType>& W, SerialSymDenseMatrix<OrdinalType, ScalarType>& B )
00074 {
00075   // Local variables.
00076   // Note: dimemensions of W are obtained so we can compute W^T*A*W for either cases.
00077   OrdinalType A_nrowcols = A.numRows();  // A is a symmetric matrix and is assumed square.
00078   OrdinalType B_nrowcols = (ETranspChar[transw]!='N') ? W.numCols() : W.numRows();
00079   OrdinalType W_nrows = (ETranspChar[transw]!='N') ? W.numRows() : W.numCols();
00080   OrdinalType W_ncols = (ETranspChar[transw]!='N') ? W.numCols() : W.numRows();
00081 
00082   bool isBUpper = B.upper();
00083 
00084   // Check for consistent dimensions.
00085   TEUCHOS_TEST_FOR_EXCEPTION( B_nrowcols != B.numRows(), std::out_of_range,
00086     "Teuchos::symMatTripleProduct<>() : "
00087     "Num Rows/Cols B (" << B.numRows() << ") inconsistent with W ("<< B_nrowcols << ")");
00088   TEUCHOS_TEST_FOR_EXCEPTION( A_nrowcols != W_nrows, std::out_of_range,
00089     "Teuchos::symMatTripleProduct<>() : "
00090     "Num Rows/Cols A (" << A_nrowcols << ") inconsistent with W ("<< W_nrows << ")");
00091 
00092   // Scale by zero, initialized B to zeros and return.
00093   if ( alpha == ScalarTraits<ScalarType>::zero() )
00094   {
00095     B.putScalar();
00096     return;
00097   }
00098 
00099   // Workspace.
00100   SerialDenseMatrix<OrdinalType, ScalarType> AW;
00101 
00102   // BLAS class.
00103   BLAS<OrdinalType, ScalarType> blas;
00104   ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00105   ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00106 
00107   // Separate two cases because BLAS only supports symmetric matrix-matrix multply w/o transposes.
00108   if (ETranspChar[transw]!='N') {
00109     // Size AW to compute A*W
00110     AW.shapeUninitialized(A_nrowcols,W_ncols);
00111 
00112     // A*W
00113     AW.multiply( Teuchos::LEFT_SIDE, alpha, A, W, ScalarTraits<ScalarType>::zero() );
00114 
00115     // B = W^T*A*W
00116     if (isBUpper) {
00117       for (int j=0; j<B_nrowcols; ++j)
00118   blas.GEMV( transw, W_nrows, j+1, one, W.values(), W.stride(), AW[j], 1, zero, &B(0,j), 1 );
00119     }
00120     else {
00121       for (int j=0; j<B_nrowcols; ++j)
00122   blas.GEMV( transw, W_nrows, B_nrowcols-j, one, W[j], W.stride(), AW[j], 1, zero, &B(j,j), 1 );
00123     }
00124   }
00125   else {
00126     // Size AW to compute W*A
00127     AW.shapeUninitialized(W_ncols, A_nrowcols);
00128 
00129     // W*A
00130     AW.multiply( Teuchos::RIGHT_SIDE, alpha, A, W, ScalarTraits<ScalarType>::zero() );
00131 
00132     // B = W*A*W^T
00133     if (isBUpper) {
00134       for (int j=0; j<B_nrowcols; ++j)
00135   for (int i=0; i<=j; ++i)
00136     blas.GEMV( transw, 1, A_nrowcols, one, &AW(i,0), AW.stride(), &W(j,0), W.stride(), zero, &B(i,j), 1 );
00137     }
00138     else {
00139       for (int j=0; j<B_nrowcols; ++j)
00140   for (int i=j; i<B_nrowcols; ++i)
00141     blas.GEMV( transw, 1, A_nrowcols, one, &AW(i,0), AW.stride(), &W(j,0), W.stride(), zero, &B(i,j), 1 );
00142     }
00143   }
00144 
00145   return;
00146 }
00147 
00157 template<typename OrdinalType, typename ScalarType>
00158 SerialDenseVector<OrdinalType,ScalarType>
00159 getCol( DataAccess CV, SerialDenseMatrix<OrdinalType, ScalarType>& A, const OrdinalType col )
00160 {
00161   return SerialDenseVector<OrdinalType, ScalarType>(CV, A[col], A.numRows());
00162 }
00163 
00173 template<typename OrdinalType, typename ScalarType>
00174 bool setCol( const SerialDenseVector<OrdinalType, ScalarType>& v,
00175        const OrdinalType col,
00176        SerialDenseMatrix<OrdinalType, ScalarType>& A )
00177 {
00178   if (v.length() != A.numRows()) return false;
00179 
00180   std::copy(v.values(),v.values()+v.length(),A[col]);
00181 
00182   return true;
00183 }
00184 
00195 template<typename OrdinalType, typename ScalarType>
00196 Teuchos::RCP<SerialBandDenseMatrix<OrdinalType, ScalarType> >
00197 generalToBanded(const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& A,
00198     const OrdinalType kl, const OrdinalType ku,
00199     const bool factorFormat)
00200 {
00201   OrdinalType m = A->numRows();
00202   OrdinalType n = A->numCols();
00203 
00204   // Check that the new matrix is consistent.
00205   TEUCHOS_TEST_FOR_EXCEPTION(A->values()==0, std::invalid_argument,
00206          "SerialBandDenseSolver<T>::generalToBanded: A is an empty SerialDenseMatrix<T>!");
00207   TEUCHOS_TEST_FOR_EXCEPTION(kl<0 || kl>m, std::invalid_argument,
00208          "SerialBandDenseSolver<T>::generalToBanded: The lower bandwidth kl is invalid!");
00209   TEUCHOS_TEST_FOR_EXCEPTION(ku<0 || ku>n, std::invalid_argument,
00210          "SerialBandDenseSolver<T>::generalToBanded: The upper bandwidth ku is invalid!");
00211 
00212   OrdinalType extraBands = (factorFormat ? kl : 0);
00213   Teuchos::RCP<SerialBandDenseMatrix<OrdinalType, ScalarType> > AB =
00214     rcp( new SerialBandDenseMatrix<OrdinalType,ScalarType>(m,n,kl,extraBands+ku,true));
00215 
00216   for (OrdinalType j = 0; j < n; j++) {
00217     for (OrdinalType i=TEUCHOS_MAX(0,j-ku); i<=TEUCHOS_MIN(m-1,j+kl); i++) {
00218       (*AB)(i,j) = (*A)(i,j);
00219     }
00220   }
00221   return(AB);
00222 }
00223 
00231 template<typename OrdinalType, typename ScalarType>
00232 Teuchos::RCP<SerialDenseMatrix<OrdinalType, ScalarType> >
00233 bandedToGeneral(const RCP<SerialBandDenseMatrix<OrdinalType,ScalarType> >& AB)
00234 {
00235 
00236   OrdinalType m = AB->numRows();
00237   OrdinalType n = AB->numCols();
00238   OrdinalType kl = AB->lowerBandwidth();
00239   OrdinalType ku = AB->upperBandwidth();
00240 
00241   // Check that the new matrix is consistent.
00242   TEUCHOS_TEST_FOR_EXCEPTION(AB->values()==0, std::invalid_argument,
00243          "SerialBandDenseSolver<T>::bandedToGeneral: AB is an empty SerialBandDenseMatrix<T>!");
00244 
00245   Teuchos::RCP<SerialDenseMatrix<OrdinalType, ScalarType> > A = rcp( new SerialDenseMatrix<OrdinalType,ScalarType>(m,n) );
00246   for (OrdinalType j = 0; j < n; j++) {
00247     for (OrdinalType i=TEUCHOS_MAX(0,j-ku); i<=TEUCHOS_MIN(m-1,j+kl); i++) {
00248       (*A)(i,j) = (*AB)(i,j);
00249     }
00250   }
00251   return(A);
00252 }
00253 
00254 } // namespace Teuchos
00255 
00256 #endif /* _TEUCHOS_SERIALDENSEHELPERS_HPP_ */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines