Teuchos - Trilinos Tools Package  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 #include "Teuchos_CommHelpers.hpp"
59 #include "Teuchos_DefaultComm.hpp"
60 #include "Teuchos_DefaultSerialComm.hpp"
61 
62 namespace Teuchos {
63 
75 template<typename OrdinalType, typename ScalarType>
76 void symMatTripleProduct( ETransp transw, const ScalarType alpha, const SerialSymDenseMatrix<OrdinalType, ScalarType>& A,
78 {
79  // Local variables.
80  // Note: dimemensions of W are obtained so we can compute W^T*A*W for either cases.
81  OrdinalType A_nrowcols = A.numRows(); // A is a symmetric matrix and is assumed square.
82  OrdinalType B_nrowcols = (ETranspChar[transw]!='N') ? W.numCols() : W.numRows();
83  OrdinalType W_nrows = (ETranspChar[transw]!='N') ? W.numRows() : W.numCols();
84  OrdinalType W_ncols = (ETranspChar[transw]!='N') ? W.numCols() : W.numRows();
85 
86  bool isBUpper = B.upper();
87 
88  // Check for consistent dimensions.
89  TEUCHOS_TEST_FOR_EXCEPTION( B_nrowcols != B.numRows(), std::out_of_range,
90  "Teuchos::symMatTripleProduct<>() : "
91  "Num Rows/Cols B (" << B.numRows() << ") inconsistent with W ("<< B_nrowcols << ")");
92  TEUCHOS_TEST_FOR_EXCEPTION( A_nrowcols != W_nrows, std::out_of_range,
93  "Teuchos::symMatTripleProduct<>() : "
94  "Num Rows/Cols A (" << A_nrowcols << ") inconsistent with W ("<< W_nrows << ")");
95 
96  // Scale by zero, initialized B to zeros and return.
97  if ( alpha == ScalarTraits<ScalarType>::zero() )
98  {
99  B.putScalar();
100  return;
101  }
102 
103  // Workspace.
105 
106  // BLAS class.
108  ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
109  ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
110 
111  // Separate two cases because BLAS only supports symmetric matrix-matrix multply w/o transposes.
112  if (ETranspChar[transw]!='N') {
113  // Size AW to compute A*W
114  AW.shapeUninitialized(A_nrowcols,W_ncols);
115 
116  // A*W
118 
119  // B = W^T*A*W
120  if (isBUpper) {
121  for (OrdinalType j=0; j<B_nrowcols; ++j)
122  blas.GEMV( transw, W_nrows, j+1, one, W.values(), W.stride(), AW[j], 1, zero, &B(0,j), 1 );
123  }
124  else {
125  for (OrdinalType j=0; j<B_nrowcols; ++j)
126  blas.GEMV( transw, W_nrows, B_nrowcols-j, one, W[j], W.stride(), AW[j], 1, zero, &B(j,j), 1 );
127  }
128  }
129  else {
130  // Size AW to compute W*A
131  AW.shapeUninitialized(W_ncols, A_nrowcols);
132 
133  // W*A
135 
136  // B = W*A*W^T
137  if (isBUpper) {
138  for (OrdinalType j=0; j<B_nrowcols; ++j)
139  for (OrdinalType i=0; i<=j; ++i)
140  blas.GEMV( transw, 1, A_nrowcols, one, &AW(i,0), AW.stride(), &W(j,0), W.stride(), zero, &B(i,j), 1 );
141  }
142  else {
143  for (OrdinalType j=0; j<B_nrowcols; ++j)
144  for (OrdinalType i=j; i<B_nrowcols; ++i)
145  blas.GEMV( transw, 1, A_nrowcols, one, &AW(i,0), AW.stride(), &W(j,0), W.stride(), zero, &B(i,j), 1 );
146  }
147  }
148 
149  return;
150 }
151 
161 template<typename OrdinalType, typename ScalarType>
164 {
165  return SerialDenseVector<OrdinalType, ScalarType>(CV, A[col], A.numRows());
166 }
167 
177 template<typename OrdinalType, typename ScalarType>
179  const OrdinalType col,
181 {
182  if (v.length() != A.numRows()) return false;
183 
184  std::copy(v.values(),v.values()+v.length(),A[col]);
185 
186  return true;
187 }
188 
194 template <typename OrdinalType, typename ScalarType>
196 {
198 
199 #ifdef HAVE_MPI
200  int mpiStarted = 0;
201  MPI_Initialized(&mpiStarted);
202  if (mpiStarted)
204  else
206 #else
208 #endif
209 
210  const OrdinalType procRank = rank(*comm);
211 
212  // Construct a separate serial dense matrix and synchronize it to get around
213  // input matrices that are subviews of a larger matrix.
215  if (procRank == 0)
216  newMatrix.random();
217  else
218  newMatrix.putScalar( Teuchos::ScalarTraits<ScalarType>::zero() );
219 
220  broadcast(*comm, 0, A.numRows()*A.numCols(), newMatrix.values());
221 
222  // Assign the synchronized matrix to the input.
223  A.assign( newMatrix );
224 }
225 
226 
237 template<typename OrdinalType, typename ScalarType>
240  const OrdinalType kl, const OrdinalType ku,
241  const bool factorFormat)
242 {
243  OrdinalType m = A->numRows();
244  OrdinalType n = A->numCols();
245 
246  // Check that the new matrix is consistent.
247  TEUCHOS_TEST_FOR_EXCEPTION(A->values()==0, std::invalid_argument,
248  "SerialBandDenseSolver<T>::generalToBanded: A is an empty SerialDenseMatrix<T>!");
249  TEUCHOS_TEST_FOR_EXCEPTION(kl<0 || kl>m, std::invalid_argument,
250  "SerialBandDenseSolver<T>::generalToBanded: The lower bandwidth kl is invalid!");
251  TEUCHOS_TEST_FOR_EXCEPTION(ku<0 || ku>n, std::invalid_argument,
252  "SerialBandDenseSolver<T>::generalToBanded: The upper bandwidth ku is invalid!");
253 
254  OrdinalType extraBands = (factorFormat ? kl : 0);
256  rcp( new SerialBandDenseMatrix<OrdinalType,ScalarType>(m,n,kl,extraBands+ku,true));
257 
258  for (OrdinalType j = 0; j < n; j++) {
259  for (OrdinalType i=TEUCHOS_MAX(0,j-ku); i<=TEUCHOS_MIN(m-1,j+kl); i++) {
260  (*AB)(i,j) = (*A)(i,j);
261  }
262  }
263  return(AB);
264 }
265 
273 template<typename OrdinalType, typename ScalarType>
276 {
277 
278  OrdinalType m = AB->numRows();
279  OrdinalType n = AB->numCols();
280  OrdinalType kl = AB->lowerBandwidth();
281  OrdinalType ku = AB->upperBandwidth();
282 
283  // Check that the new matrix is consistent.
284  TEUCHOS_TEST_FOR_EXCEPTION(AB->values()==0, std::invalid_argument,
285  "SerialBandDenseSolver<T>::bandedToGeneral: AB is an empty SerialBandDenseMatrix<T>!");
286 
288  for (OrdinalType j = 0; j < n; j++) {
289  for (OrdinalType i=TEUCHOS_MAX(0,j-ku); i<=TEUCHOS_MIN(m-1,j+kl); i++) {
290  (*A)(i,j) = (*AB)(i,j);
291  }
292  }
293  return(A);
294 }
295 
296 } // namespace Teuchos
297 
298 #endif /* _TEUCHOS_SERIALDENSEHELPERS_HPP_ */
ScalarType * values() const
Data array access method.
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...
bool upper() const
Returns true if upper triangular part of this matrix has and will be used.
Templated serial dense matrix class.
void randomSyncedMatrix(Teuchos::SerialDenseMatrix< OrdinalType, ScalarType > &A)
A templated, non-member, helper function for generating a random SerialDenseMatrix that is synchroniz...
int shapeUninitialized(OrdinalType numRows, OrdinalType numCols)
Same as shape() except leaves uninitialized.
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 &lt;- alpha*A*x+beta*y or y &lt;- alpha*A&#39;*x+beta*y where A is a gene...
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...
Templated BLAS wrapper.
This class creates and provides basic support for symmetric, positive-definite dense matrices of temp...
OrdinalType numRows() const
Returns the row dimension of this matrix.
Concrete serial communicator subclass.
static Teuchos::RCP< const Comm< OrdinalType > > getComm()
Return the default global communicator.
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.
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.
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...
OrdinalType length() const
Returns the length of this vector.
Templated serial, dense, symmetric matrix class.
OrdinalType numCols() const
Returns the column dimension of this matrix.
int random()
Set all values in the matrix to be random numbers.
Templated serial dense vector class.
Defines basic traits for the scalar field type.
static T zero()
Returns representation of zero for this scalar type.
Smart reference counting pointer class for automatic garbage collection.
SerialDenseMatrix< OrdinalType, ScalarType > & assign(const SerialDenseMatrix< OrdinalType, ScalarType > &Source)
Copies values from one matrix to another.
static T one()
Returns representation of one for this scalar type.
Teuchos::DataAccess Mode enumerable type.
OrdinalType stride() const
Returns the stride between the columns of this matrix in memory.
OrdinalType numRows() const
Returns the row dimension of this matrix.
This class creates and provides basic support for dense rectangular matrix of templated type...