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_TestForException.hpp"
00053 #include "Teuchos_SerialDenseMatrix.hpp"  
00054 #include "Teuchos_SerialSymDenseMatrix.hpp" 
00055 #include "Teuchos_SerialDenseVector.hpp"  
00056 
00057 namespace Teuchos {
00058 
00070 template<typename OrdinalType, typename ScalarType>
00071 void symMatTripleProduct( ETransp transw, const ScalarType alpha, const SerialSymDenseMatrix<OrdinalType, ScalarType>& A, 
00072         const SerialDenseMatrix<OrdinalType, ScalarType>& W, SerialSymDenseMatrix<OrdinalType, ScalarType>& B )
00073 {
00074   // Local variables.
00075   // Note: dimemensions of W are obtained so we can compute W^T*A*W for either cases.
00076   OrdinalType A_nrowcols = A.numRows();  // A is a symmetric matrix and is assumed square. 
00077   OrdinalType B_nrowcols = (ETranspChar[transw]!='N') ? W.numCols() : W.numRows();
00078   OrdinalType W_nrows = (ETranspChar[transw]!='N') ? W.numRows() : W.numCols();
00079   OrdinalType W_ncols = (ETranspChar[transw]!='N') ? W.numCols() : W.numRows();
00080   
00081   bool isBUpper = B.upper();
00082 
00083   // Check for consistent dimensions.
00084   TEST_FOR_EXCEPTION( B_nrowcols != B.numRows(), std::out_of_range,
00085     "Teuchos::symMatTripleProduct<>() : "
00086     "Num Rows/Cols B (" << B.numRows() << ") inconsistent with W ("<< B_nrowcols << ")");
00087   TEST_FOR_EXCEPTION( A_nrowcols != W_nrows, std::out_of_range,
00088     "Teuchos::symMatTripleProduct<>() : "
00089     "Num Rows/Cols A (" << A_nrowcols << ") inconsistent with W ("<< W_nrows << ")");
00090 
00091   // Scale by zero, initialized B to zeros and return.
00092   if ( alpha == ScalarTraits<ScalarType>::zero() )
00093   {
00094     B.putScalar();
00095     return;
00096   }
00097 
00098   // Workspace.
00099   SerialDenseMatrix<OrdinalType, ScalarType> AW;
00100 
00101   // BLAS class.
00102   BLAS<OrdinalType, ScalarType> blas;
00103   ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00104   ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00105 
00106   // Separate two cases because BLAS only supports symmetric matrix-matrix multply w/o transposes.
00107   if (ETranspChar[transw]!='N') {
00108     // Size AW to compute A*W
00109     AW.shapeUninitialized(A_nrowcols,W_ncols);
00110   
00111     // A*W
00112     AW.multiply( Teuchos::LEFT_SIDE, alpha, A, W, ScalarTraits<ScalarType>::zero() );
00113     
00114     // B = W^T*A*W
00115     if (isBUpper) {
00116       for (int j=0; j<B_nrowcols; ++j)
00117   blas.GEMV( transw, W_nrows, j+1, one, W.values(), W.stride(), AW[j], 1, zero, &B(0,j), 1 );
00118     }
00119     else {
00120       for (int j=0; j<B_nrowcols; ++j)
00121   blas.GEMV( transw, W_nrows, B_nrowcols-j, one, W[j], W.stride(), AW[j], 1, zero, &B(j,j), 1 );
00122     }
00123   }
00124   else {
00125     // Size AW to compute W*A
00126     AW.shapeUninitialized(W_ncols, A_nrowcols);
00127 
00128     // W*A
00129     AW.multiply( Teuchos::RIGHT_SIDE, alpha, A, W, ScalarTraits<ScalarType>::zero() );
00130 
00131     // B = W*A*W^T
00132     if (isBUpper) {
00133       for (int j=0; j<B_nrowcols; ++j)
00134   for (int i=0; i<=j; ++i) 
00135     blas.GEMV( transw, 1, A_nrowcols, one, &AW(i,0), AW.stride(), &W(j,0), W.stride(), zero, &B(i,j), 1 );
00136     }
00137     else {
00138       for (int j=0; j<B_nrowcols; ++j)
00139   for (int i=j; i<B_nrowcols; ++i) 
00140     blas.GEMV( transw, 1, A_nrowcols, one, &AW(i,0), AW.stride(), &W(j,0), W.stride(), zero, &B(i,j), 1 );
00141     }
00142   }
00143 
00144   return;
00145 }
00146   
00156 template<typename OrdinalType, typename ScalarType>
00157 SerialDenseVector<OrdinalType,ScalarType>
00158 getCol( DataAccess CV, SerialDenseMatrix<OrdinalType, ScalarType>& A, const OrdinalType col )
00159 {
00160   return SerialDenseVector<OrdinalType, ScalarType>(CV, A[col], A.numRows());
00161 }
00162 
00172 template<typename OrdinalType, typename ScalarType>
00173 bool setCol( const SerialDenseVector<OrdinalType, ScalarType>& v,
00174              const OrdinalType col,
00175              SerialDenseMatrix<OrdinalType, ScalarType>& A )
00176 {
00177   if (v.length() != A.numRows()) return false;
00178 
00179   std::copy(v.values(),v.values()+v.length(),A[col]);
00180 
00181   return true;
00182 }
00183 
00184 } // namespace Teuchos
00185 
00186 #endif /* _TEUCHOS_SERIALDENSEHELPERS_HPP_ */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines