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 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #ifndef _TEUCHOS_SERIALDENSEHELPERS_HPP_
00030 #define _TEUCHOS_SERIALDENSEHELPERS_HPP_
00031 
00036 #include "Teuchos_ScalarTraits.hpp"
00037 #include "Teuchos_DataAccess.hpp"
00038 #include "Teuchos_ConfigDefs.hpp"
00039 #include "Teuchos_TestForException.hpp"
00040 #include "Teuchos_SerialDenseMatrix.hpp"  
00041 #include "Teuchos_SerialSymDenseMatrix.hpp" 
00042 #include "Teuchos_SerialDenseVector.hpp"  
00043 
00044 namespace Teuchos {
00045 
00057 template<typename OrdinalType, typename ScalarType>
00058 void symMatTripleProduct( ETransp transw, const ScalarType alpha, const SerialSymDenseMatrix<OrdinalType, ScalarType>& A, 
00059         const SerialDenseMatrix<OrdinalType, ScalarType>& W, SerialSymDenseMatrix<OrdinalType, ScalarType>& B )
00060 {
00061   // Local variables.
00062   // Note: dimemensions of W are obtained so we can compute W^T*A*W for either cases.
00063   OrdinalType A_nrowcols = A.numRows();  // A is a symmetric matrix and is assumed square. 
00064   OrdinalType B_nrowcols = (ETranspChar[transw]!='N') ? W.numCols() : W.numRows();
00065   OrdinalType W_nrows = (ETranspChar[transw]!='N') ? W.numRows() : W.numCols();
00066   OrdinalType W_ncols = (ETranspChar[transw]!='N') ? W.numCols() : W.numRows();
00067   
00068   bool isBUpper = B.upper();
00069 
00070   // Check for consistent dimensions.
00071   TEST_FOR_EXCEPTION( B_nrowcols != B.numRows(), std::out_of_range,
00072     "Teuchos::symMatTripleProduct<>() : "
00073     "Num Rows/Cols B (" << B.numRows() << ") inconsistent with W ("<< B_nrowcols << ")");
00074   TEST_FOR_EXCEPTION( A_nrowcols != W_nrows, std::out_of_range,
00075     "Teuchos::symMatTripleProduct<>() : "
00076     "Num Rows/Cols A (" << A_nrowcols << ") inconsistent with W ("<< W_nrows << ")");
00077 
00078   // Scale by zero, initialized B to zeros and return.
00079   if ( alpha == ScalarTraits<ScalarType>::zero() )
00080   {
00081     B.putScalar();
00082     return;
00083   }
00084 
00085   // Workspace.
00086   SerialDenseMatrix<OrdinalType, ScalarType> AW;
00087 
00088   // BLAS class.
00089   BLAS<OrdinalType, ScalarType> blas;
00090   ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00091   ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00092 
00093   // Separate two cases because BLAS only supports symmetric matrix-matrix multply w/o transposes.
00094   if (ETranspChar[transw]!='N') {
00095     // Size AW to compute A*W
00096     AW.shapeUninitialized(A_nrowcols,W_ncols);
00097   
00098     // A*W
00099     AW.multiply( Teuchos::LEFT_SIDE, alpha, A, W, ScalarTraits<ScalarType>::zero() );
00100     
00101     // B = W^T*A*W
00102     if (isBUpper) {
00103       for (int j=0; j<B_nrowcols; ++j)
00104   blas.GEMV( transw, W_nrows, j+1, one, W.values(), W.stride(), AW[j], 1, zero, &B(0,j), 1 );
00105     }
00106     else {
00107       for (int j=0; j<B_nrowcols; ++j)
00108   blas.GEMV( transw, W_nrows, B_nrowcols-j, one, W[j], W.stride(), AW[j], 1, zero, &B(j,j), 1 );
00109     }
00110   }
00111   else {
00112     // Size AW to compute W*A
00113     AW.shapeUninitialized(W_ncols, A_nrowcols);
00114 
00115     // W*A
00116     AW.multiply( Teuchos::RIGHT_SIDE, alpha, A, W, ScalarTraits<ScalarType>::zero() );
00117 
00118     // B = W*A*W^T
00119     if (isBUpper) {
00120       for (int j=0; j<B_nrowcols; ++j)
00121   for (int i=0; i<=j; ++i) 
00122     blas.GEMV( transw, 1, A_nrowcols, one, &AW(i,0), AW.stride(), &W(j,0), W.stride(), zero, &B(i,j), 1 );
00123     }
00124     else {
00125       for (int j=0; j<B_nrowcols; ++j)
00126   for (int i=j; i<B_nrowcols; ++i) 
00127     blas.GEMV( transw, 1, A_nrowcols, one, &AW(i,0), AW.stride(), &W(j,0), W.stride(), zero, &B(i,j), 1 );
00128     }
00129   }
00130 
00131   return;
00132 }
00133   
00143 template<typename OrdinalType, typename ScalarType>
00144 SerialDenseVector<OrdinalType,ScalarType>
00145 getCol( DataAccess CV, SerialDenseMatrix<OrdinalType, ScalarType>& A, const OrdinalType col )
00146 {
00147   return SerialDenseVector<OrdinalType, ScalarType>(CV, A[col], A.numRows());
00148 }
00149 
00159 template<typename OrdinalType, typename ScalarType>
00160 bool setCol( const SerialDenseVector<OrdinalType, ScalarType>& v,
00161              const OrdinalType col,
00162              SerialDenseMatrix<OrdinalType, ScalarType>& A )
00163 {
00164   if (v.length() != A.numRows()) return false;
00165 
00166   std::copy(v.values(),v.values()+v.length(),A[col]);
00167 
00168   return true;
00169 }
00170 
00171 } // namespace Teuchos
00172 
00173 #endif /* _TEUCHOS_SERIALDENSEHELPERS_HPP_ */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 09:57:44 2011 for Teuchos - Trilinos Tools Package by  doxygen 1.6.3