Teuchos - Trilinos Tools Package Version of the Day
Loading...
Searching...
No Matches
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
52#include "Teuchos_Assert.hpp"
57
58#include "Teuchos_CommHelpers.hpp"
59#include "Teuchos_DefaultComm.hpp"
60#include "Teuchos_DefaultSerialComm.hpp"
61
62namespace Teuchos {
63
75template<typename OrdinalType, typename ScalarType>
76void 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.
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
161template<typename OrdinalType, typename ScalarType>
164{
166}
167
177template<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
194template <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
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
237template<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
273template<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_ */
Teuchos header file which uses auto-configuration information to include necessary C++ headers.
Teuchos::DataAccess Mode enumerable type.
Defines basic traits for the scalar field type.
Templated serial dense matrix class.
Templated serial dense matrix class.
Templated serial dense vector class.
Templated serial, dense, symmetric matrix class.
Templated BLAS wrapper.
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'*x+beta*y where A is a ge...
static Teuchos::RCP< const Comm< OrdinalType > > getComm()
Return the default global communicator.
Smart reference counting pointer class for automatic garbage collection.
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...
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...
Concrete serial communicator subclass.
This class creates and provides basic support for dense rectangular matrix of templated type.
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero())
Set all values in the matrix to a constant value.
void randomSyncedMatrix(Teuchos::SerialDenseMatrix< OrdinalType, ScalarType > &A)
A templated, non-member, helper function for generating a random SerialDenseMatrix that is synchroniz...
SerialDenseMatrix< OrdinalType, ScalarType > & assign(const SerialDenseMatrix< OrdinalType, ScalarType > &Source)
Copies values from one matrix to another.
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...
OrdinalType stride() const
Returns the stride between the columns of this matrix in memory.
int random()
Set all values in the matrix to be random numbers.
OrdinalType numRows() const
Returns the row dimension of this matrix.
ScalarType * values() const
Data array access method.
int shapeUninitialized(OrdinalType numRows, OrdinalType numCols)
Same as shape() except leaves uninitialized.
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.
OrdinalType numCols() const
Returns the column dimension of this matrix.
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 ...
This class creates and provides basic support for dense vectors of templated type as a specialization...
OrdinalType length() const
Returns the length of this vector.
This class creates and provides basic support for symmetric, positive-definite dense matrices of temp...
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero(), bool fullMatrix=false)
Set all values in the matrix to a constant value.
bool upper() const
Returns true if upper triangular part of this matrix has and will be used.
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...
OrdinalType numRows() const
Returns the row dimension of this matrix.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Macro for throwing an exception with breakpointing to ease debugging.
The Teuchos namespace contains all of the classes, structs and enums used by Teuchos,...
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Deprecated.
This structure defines some basic traits for a scalar field type.
static T one()
Returns representation of one for this scalar type.
static T zero()
Returns representation of zero for this scalar type.