Teuchos - Trilinos Tools Package  Version of the Day
Teuchos_SerialDenseHelpers.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Teuchos: Common Tools Package
5 // Copyright (2004) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #ifndef _TEUCHOS_SERIALDENSEHELPERS_HPP_
43 #define _TEUCHOS_SERIALDENSEHELPERS_HPP_
44 
49 #include "Teuchos_ScalarTraits.hpp"
50 #include "Teuchos_DataAccess.hpp"
51 #include "Teuchos_ConfigDefs.hpp"
52 #include "Teuchos_Assert.hpp"
57 
58 namespace Teuchos {
59 
71 template<typename OrdinalType, typename ScalarType>
72 void symMatTripleProduct( ETransp transw, const ScalarType alpha, const SerialSymDenseMatrix<OrdinalType, ScalarType>& A,
74 {
75  // Local variables.
76  // Note: dimemensions of W are obtained so we can compute W^T*A*W for either cases.
77  OrdinalType A_nrowcols = A.numRows(); // A is a symmetric matrix and is assumed square.
78  OrdinalType B_nrowcols = (ETranspChar[transw]!='N') ? W.numCols() : W.numRows();
79  OrdinalType W_nrows = (ETranspChar[transw]!='N') ? W.numRows() : W.numCols();
80  OrdinalType W_ncols = (ETranspChar[transw]!='N') ? W.numCols() : W.numRows();
81 
82  bool isBUpper = B.upper();
83 
84  // Check for consistent dimensions.
85  TEUCHOS_TEST_FOR_EXCEPTION( B_nrowcols != B.numRows(), std::out_of_range,
86  "Teuchos::symMatTripleProduct<>() : "
87  "Num Rows/Cols B (" << B.numRows() << ") inconsistent with W ("<< B_nrowcols << ")");
88  TEUCHOS_TEST_FOR_EXCEPTION( A_nrowcols != W_nrows, std::out_of_range,
89  "Teuchos::symMatTripleProduct<>() : "
90  "Num Rows/Cols A (" << A_nrowcols << ") inconsistent with W ("<< W_nrows << ")");
91 
92  // Scale by zero, initialized B to zeros and return.
93  if ( alpha == ScalarTraits<ScalarType>::zero() )
94  {
95  B.putScalar();
96  return;
97  }
98 
99  // Workspace.
101 
102  // BLAS class.
104  ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
105  ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
106 
107  // Separate two cases because BLAS only supports symmetric matrix-matrix multply w/o transposes.
108  if (ETranspChar[transw]!='N') {
109  // Size AW to compute A*W
110  AW.shapeUninitialized(A_nrowcols,W_ncols);
111 
112  // A*W
114 
115  // B = W^T*A*W
116  if (isBUpper) {
117  for (int j=0; j<B_nrowcols; ++j)
118  blas.GEMV( transw, W_nrows, j+1, one, W.values(), W.stride(), AW[j], 1, zero, &B(0,j), 1 );
119  }
120  else {
121  for (int j=0; j<B_nrowcols; ++j)
122  blas.GEMV( transw, W_nrows, B_nrowcols-j, one, W[j], W.stride(), AW[j], 1, zero, &B(j,j), 1 );
123  }
124  }
125  else {
126  // Size AW to compute W*A
127  AW.shapeUninitialized(W_ncols, A_nrowcols);
128 
129  // W*A
131 
132  // B = W*A*W^T
133  if (isBUpper) {
134  for (int j=0; j<B_nrowcols; ++j)
135  for (int i=0; i<=j; ++i)
136  blas.GEMV( transw, 1, A_nrowcols, one, &AW(i,0), AW.stride(), &W(j,0), W.stride(), zero, &B(i,j), 1 );
137  }
138  else {
139  for (int j=0; j<B_nrowcols; ++j)
140  for (int i=j; i<B_nrowcols; ++i)
141  blas.GEMV( transw, 1, A_nrowcols, one, &AW(i,0), AW.stride(), &W(j,0), W.stride(), zero, &B(i,j), 1 );
142  }
143  }
144 
145  return;
146 }
147 
157 template<typename OrdinalType, typename ScalarType>
160 {
161  return SerialDenseVector<OrdinalType, ScalarType>(CV, A[col], A.numRows());
162 }
163 
173 template<typename OrdinalType, typename ScalarType>
175  const OrdinalType col,
177 {
178  if (v.length() != A.numRows()) return false;
179 
180  std::copy(v.values(),v.values()+v.length(),A[col]);
181 
182  return true;
183 }
184 
195 template<typename OrdinalType, typename ScalarType>
198  const OrdinalType kl, const OrdinalType ku,
199  const bool factorFormat)
200 {
201  OrdinalType m = A->numRows();
202  OrdinalType n = A->numCols();
203 
204  // Check that the new matrix is consistent.
205  TEUCHOS_TEST_FOR_EXCEPTION(A->values()==0, std::invalid_argument,
206  "SerialBandDenseSolver<T>::generalToBanded: A is an empty SerialDenseMatrix<T>!");
207  TEUCHOS_TEST_FOR_EXCEPTION(kl<0 || kl>m, std::invalid_argument,
208  "SerialBandDenseSolver<T>::generalToBanded: The lower bandwidth kl is invalid!");
209  TEUCHOS_TEST_FOR_EXCEPTION(ku<0 || ku>n, std::invalid_argument,
210  "SerialBandDenseSolver<T>::generalToBanded: The upper bandwidth ku is invalid!");
211 
212  OrdinalType extraBands = (factorFormat ? kl : 0);
214  rcp( new SerialBandDenseMatrix<OrdinalType,ScalarType>(m,n,kl,extraBands+ku,true));
215 
216  for (OrdinalType j = 0; j < n; j++) {
217  for (OrdinalType i=TEUCHOS_MAX(0,j-ku); i<=TEUCHOS_MIN(m-1,j+kl); i++) {
218  (*AB)(i,j) = (*A)(i,j);
219  }
220  }
221  return(AB);
222 }
223 
231 template<typename OrdinalType, typename ScalarType>
234 {
235 
236  OrdinalType m = AB->numRows();
237  OrdinalType n = AB->numCols();
238  OrdinalType kl = AB->lowerBandwidth();
239  OrdinalType ku = AB->upperBandwidth();
240 
241  // Check that the new matrix is consistent.
242  TEUCHOS_TEST_FOR_EXCEPTION(AB->values()==0, std::invalid_argument,
243  "SerialBandDenseSolver<T>::bandedToGeneral: AB is an empty SerialBandDenseMatrix<T>!");
244 
246  for (OrdinalType j = 0; j < n; j++) {
247  for (OrdinalType i=TEUCHOS_MAX(0,j-ku); i<=TEUCHOS_MIN(m-1,j+kl); i++) {
248  (*A)(i,j) = (*AB)(i,j);
249  }
250  }
251  return(A);
252 }
253 
254 } // namespace Teuchos
255 
256 #endif /* _TEUCHOS_SERIALDENSEHELPERS_HPP_ */
bool setCol(const SerialDenseVector< OrdinalType, ScalarType > &v, const OrdinalType col, SerialDenseMatrix< OrdinalType, ScalarType > &A)
A templated, non-member, helper function for setting a SerialDenseMatrix column using a SerialDenseVe...
Templated serial dense matrix class.
ScalarType * values() const
Data array access method.
int shapeUninitialized(OrdinalType numRows, OrdinalType numCols)
Same as shape() except leaves uninitialized.
void symMatTripleProduct(ETransp transw, const ScalarType alpha, const SerialSymDenseMatrix< OrdinalType, ScalarType > &A, const SerialDenseMatrix< OrdinalType, ScalarType > &W, SerialSymDenseMatrix< OrdinalType, ScalarType > &B)
A templated, non-member, helper function for computing the matrix triple-product: B = alpha*W^T*A*W o...
int multiply(ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix< OrdinalType, ScalarType > &A, const SerialDenseMatrix< OrdinalType, ScalarType > &B, ScalarType beta)
Multiply A * B and add them to this; this = beta * this + alpha*A*B.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Macro for throwing an exception with breakpointing to ease debugging.
SerialDenseVector< OrdinalType, ScalarType > getCol(DataAccess CV, SerialDenseMatrix< OrdinalType, ScalarType > &A, const OrdinalType col)
A templated, non-member, helper function for viewing or copying a column of a SerialDenseMatrix as a ...
Teuchos header file which uses auto-configuration information to include necessary C++ headers...
void GEMV(ETransp trans, const OrdinalType m, const OrdinalType n, const alpha_type alpha, const A_type *A, const OrdinalType lda, const x_type *x, const OrdinalType incx, const beta_type beta, ScalarType *y, const OrdinalType incy) const
Performs the matrix-vector operation: y <- alpha*A*x+beta*y or y <- alpha*A&#39;*x+beta*y where A is a ge...
Templated BLAS wrapper.
This class creates and provides basic support for symmetric, positive-definite dense matrices of temp...
This class creates and provides basic support for dense vectors of templated type as a specialization...
This structure defines some basic traits for a scalar field type.
Templated serial dense matrix class.
OrdinalType numRows() const
Returns the row dimension of this matrix.
Teuchos::RCP< SerialDenseMatrix< OrdinalType, ScalarType > > bandedToGeneral(const RCP< SerialBandDenseMatrix< OrdinalType, ScalarType > > &AB)
A templated, non-member, helper function for converting a SerialBandDenseMatrix to a SerialDenseMatri...
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Deprecated.
bool upper() const
Returns true if upper triangular part of this matrix has and will be used.
OrdinalType numRows() const
Returns the row dimension of this matrix.
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero(), bool fullMatrix=false)
Set all values in the matrix to a constant value.
This class creates and provides basic support for banded dense matrices of templated type...
Teuchos::RCP< SerialBandDenseMatrix< OrdinalType, ScalarType > > generalToBanded(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &A, const OrdinalType kl, const OrdinalType ku, const bool factorFormat)
A templated, non-member, helper function for converting a SerialDenseMatrix to a SerialBandDenseMatri...
Templated serial, dense, symmetric matrix class.
The Teuchos namespace contains all of the classes, structs and enums used by Teuchos, as well as a number of utility routines.
Templated serial dense vector class.
Defines basic traits for the scalar field type.
static T zero()
Returns representation of zero for this scalar type.
OrdinalType stride() const
Returns the stride between the columns of this matrix in memory.
Smart reference counting pointer class for automatic garbage collection.
OrdinalType numCols() const
Returns the column dimension of this matrix.
OrdinalType length() const
Returns the length of this vector.
static T one()
Returns representation of one for this scalar type.
Teuchos::DataAccess Mode enumerable type.
This class creates and provides basic support for dense rectangular matrix of templated type...