40#ifndef _TEUCHOS_SERIALDENSEMATRIX_HPP_
41#define _TEUCHOS_SERIALDENSEMATRIX_HPP_
66template<
typename OrdinalType,
typename ScalarType>
228 ScalarType&
operator () (OrdinalType rowIndex, OrdinalType colIndex);
238 const ScalarType&
operator () (OrdinalType rowIndex, OrdinalType colIndex)
const;
254 const ScalarType*
operator [] (OrdinalType colIndex)
const;
288 int scale (
const ScalarType alpha );
381 virtual std::ostream&
print(std::ostream& os)
const;
385 void copyMat(ScalarType* inputMatrix, OrdinalType strideInput,
386 OrdinalType
numRows, OrdinalType
numCols, ScalarType* outputMatrix,
387 OrdinalType strideOutput, OrdinalType startRow, OrdinalType startCol,
390 void checkIndex( OrdinalType rowIndex, OrdinalType colIndex = 0 )
const;
397#pragma GCC diagnostic push
398#pragma GCC diagnostic ignored "-Wvla"
399 return new ScalarType[size];
400#pragma GCC diagnostic pop
414template<
typename OrdinalType,
typename ScalarType>
416 OrdinalType numRows_in, OrdinalType numCols_in,
bool zeroOut
418 : numRows_(numRows_in),
419 numCols_(numCols_in),
422 values_(allocateValues(numRows_in, numCols_in))
429template<
typename OrdinalType,
typename ScalarType>
431 DataAccess CV, ScalarType* values_in, OrdinalType stride_in, OrdinalType numRows_in,
432 OrdinalType numCols_in
434 : numRows_(numRows_in),
435 numCols_(numCols_in),
437 valuesCopied_(false),
449template<
typename OrdinalType,
typename ScalarType>
451 : valuesCopied_(true)
483 for (OrdinalType j=0; j<
numCols_; j++) {
484 for (OrdinalType i=0; i<
numRows_; i++) {
495 for (OrdinalType j=0; j<
numCols_; j++) {
496 for (OrdinalType i=0; i<
numRows_; i++) {
504template<
typename OrdinalType,
typename ScalarType>
508 : numRows_(Source.numRows_), numCols_(Source.numCols_), stride_(Source.stride_),
509 valuesCopied_(false), values_(Source.values_)
521template<
typename OrdinalType,
typename ScalarType>
524 OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType startRow,
527 :
CompObject(), numRows_(numRows_in), numCols_(numCols_in), stride_(Source.stride_),
528 valuesCopied_(false), values_(Source.values_)
543template<
typename OrdinalType,
typename ScalarType>
553template<
typename OrdinalType,
typename ScalarType>
555 OrdinalType numRows_in, OrdinalType numCols_in
559 numRows_ = numRows_in;
560 numCols_ = numCols_in;
562 values_ = allocateValues(stride_, numCols_);
564 valuesCopied_ =
true;
568template<
typename OrdinalType,
typename ScalarType>
570 OrdinalType numRows_in, OrdinalType numCols_in
574 numRows_ = numRows_in;
575 numCols_ = numCols_in;
577 values_ = allocateValues(stride_, numCols_);
578 valuesCopied_ =
true;
582template<
typename OrdinalType,
typename ScalarType>
584 OrdinalType numRows_in, OrdinalType numCols_in
588 ScalarType* values_tmp = allocateValues(numRows_in, numCols_in);
590 for(OrdinalType k = 0; k < numRows_in * numCols_in; k++)
592 values_tmp[k] = zero;
594 OrdinalType numRows_tmp =
TEUCHOS_MIN(numRows_, numRows_in);
595 OrdinalType numCols_tmp =
TEUCHOS_MIN(numCols_, numCols_in);
598 copyMat(values_, stride_, numRows_tmp, numCols_tmp, values_tmp,
602 numRows_ = numRows_in;
603 numCols_ = numCols_in;
605 values_ = values_tmp;
606 valuesCopied_ =
true;
614template<
typename OrdinalType,
typename ScalarType>
618 for(OrdinalType j = 0; j < numCols_; j++)
620 for(OrdinalType i = 0; i < numRows_; i++)
622 values_[i + j*stride_] = value_in;
628template<
typename OrdinalType,
typename ScalarType>
void
646 std::swap(values_ ,
B.values_);
647 std::swap(numRows_,
B.numRows_);
648 std::swap(numCols_,
B.numCols_);
649 std::swap(stride_,
B.stride_);
650 std::swap(valuesCopied_,
B.valuesCopied_);
653template<
typename OrdinalType,
typename ScalarType>
657 for(OrdinalType j = 0; j < numCols_; j++)
659 for(OrdinalType i = 0; i < numRows_; i++)
667template<
typename OrdinalType,
typename ScalarType>
695 if(stride_ > 0 && numCols_ > 0) {
696 values_ = allocateValues(stride_, numCols_);
697 valuesCopied_ =
true;
714 if(stride_ > 0 && numCols_ > 0) {
715 values_ = allocateValues(stride_, numCols_);
716 valuesCopied_ =
true;
720 copyMat(Source.
values_, Source.
stride_, numRows_, numCols_, values_, stride_, 0, 0);
725template<
typename OrdinalType,
typename ScalarType>
737template<
typename OrdinalType,
typename ScalarType>
749template<
typename OrdinalType,
typename ScalarType>
761 copyMat(Source.
values_, Source.
stride_, numRows_, numCols_, values_, stride_, 0, 0);
769template<
typename OrdinalType,
typename ScalarType>
772#ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
773 checkIndex( rowIndex, colIndex );
775 return(values_[colIndex * stride_ + rowIndex]);
778template<
typename OrdinalType,
typename ScalarType>
781#ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
782 checkIndex( rowIndex, colIndex );
784 return(values_[colIndex * stride_ + rowIndex]);
787template<
typename OrdinalType,
typename ScalarType>
790#ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
791 checkIndex( 0, colIndex );
793 return(values_ + colIndex * stride_);
796template<
typename OrdinalType,
typename ScalarType>
799#ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
800 checkIndex( 0, colIndex );
802 return(values_ + colIndex * stride_);
809template<
typename OrdinalType,
typename ScalarType>
816 for(j = 0; j < numCols_; j++)
819 ptr = values_ + j * stride_;
820 for(i = 0; i < numRows_; i++)
831 if (flopCounter_!=0) updateFlops(numRows_ * numCols_);
835template<
typename OrdinalType,
typename ScalarType>
841 for (i = 0; i < numRows_; i++) {
843 for (j=0; j< numCols_; j++) {
849 if (flopCounter_!=0) updateFlops(numRows_ * numCols_);
853template<
typename OrdinalType,
typename ScalarType>
858 for (j = 0; j < numCols_; j++) {
859 for (i = 0; i < numRows_; i++) {
865 if (flopCounter_!=0) updateFlops(numRows_ * numCols_);
873template<
typename OrdinalType,
typename ScalarType>
884 for(i = 0; i < numRows_; i++)
886 for(j = 0; j < numCols_; j++)
888 if((*
this)(i, j) != Operand(i, j))
898template<
typename OrdinalType,
typename ScalarType>
901 return !((*this) == Operand);
908template<
typename OrdinalType,
typename ScalarType>
911 this->scale( alpha );
915template<
typename OrdinalType,
typename ScalarType>
921 for (j=0; j<numCols_; j++) {
922 ptr = values_ + j*stride_;
923 for (i=0; i<numRows_; i++) { *ptr = alpha * (*ptr); ptr++; }
926 if (flopCounter_!=0) updateFlops( numRows_*numCols_ );
930template<
typename OrdinalType,
typename ScalarType>
937 if ((numRows_ !=
A.numRows_) || (numCols_ !=
A.numCols_))
941 for (j=0; j<numCols_; j++) {
942 ptr = values_ + j*stride_;
943 for (i=0; i<numRows_; i++) { *ptr =
A(i,j) * (*ptr); ptr++; }
946 if (flopCounter_!=0) updateFlops( numRows_*numCols_ );
950template<
typename OrdinalType,
typename ScalarType>
954 OrdinalType A_nrows = (
ETranspChar[transa]!=
'N') ?
A.numCols() :
A.numRows();
955 OrdinalType A_ncols = (
ETranspChar[transa]!=
'N') ?
A.numRows() :
A.numCols();
956 OrdinalType B_nrows = (
ETranspChar[transb]!=
'N') ?
B.numCols() :
B.numRows();
957 OrdinalType B_ncols = (
ETranspChar[transb]!=
'N') ?
B.numRows() :
B.numCols();
958 if ((numRows_ != A_nrows) || (A_ncols != B_nrows) || (numCols_ != B_ncols))
963 this->GEMM(transa, transb, numRows_, numCols_, A_ncols, alpha,
A.values(),
A.stride(),
B.values(),
B.stride(), beta, values_, stride_);
966 if (flopCounter_!=0) {
967 double nflops = 2 * numRows_;
975template<
typename OrdinalType,
typename ScalarType>
979 OrdinalType A_nrows =
A.numRows(), A_ncols =
A.numCols();
980 OrdinalType B_nrows =
B.numRows(), B_ncols =
B.numCols();
983 if ((numRows_ != A_nrows) || (A_ncols != B_nrows) || (numCols_ != B_ncols)) {
987 if ((numRows_ != B_nrows) || (B_ncols != A_nrows) || (numCols_ != A_ncols)) {
994 this->SYMM(sideA, uplo, numRows_, numCols_, alpha,
A.values(),
A.stride(),
B.values(),
B.stride(), beta, values_, stride_);
997 if (flopCounter_!=0) {
998 double nflops = 2 * numRows_;
1001 updateFlops(nflops);
1006template<
typename OrdinalType,
typename ScalarType>
1011 os <<
"Values_copied : yes" << std::endl;
1013 os <<
"Values_copied : no" << std::endl;
1014 os <<
"Rows : " << numRows_ << std::endl;
1015 os <<
"Columns : " << numCols_ << std::endl;
1016 os <<
"LDA : " << stride_ << std::endl;
1017 if(numRows_ == 0 || numCols_ == 0) {
1018 os <<
"(matrix is empty, no values to display)" << std::endl;
1020 for(OrdinalType i = 0; i < numRows_; i++) {
1021 for(OrdinalType j = 0; j < numCols_; j++){
1022 os << (*this)(i,j) <<
" ";
1034template<
typename OrdinalType,
typename ScalarType>
1037 "SerialDenseMatrix<T>::checkIndex: "
1038 "Row index " << rowIndex <<
" out of range [0, "<< numRows_ <<
")");
1040 "SerialDenseMatrix<T>::checkIndex: "
1041 "Col index " << colIndex <<
" out of range [0, "<< numCols_ <<
")");
1044template<
typename OrdinalType,
typename ScalarType>
1051 valuesCopied_ =
false;
1055template<
typename OrdinalType,
typename ScalarType>
1057 ScalarType* inputMatrix, OrdinalType strideInput, OrdinalType numRows_in,
1058 OrdinalType numCols_in, ScalarType* outputMatrix, OrdinalType strideOutput,
1059 OrdinalType startRow, OrdinalType startCol, ScalarType alpha
1063 ScalarType* ptr1 = 0;
1064 ScalarType* ptr2 = 0;
1065 for(j = 0; j < numCols_in; j++) {
1066 ptr1 = outputMatrix + (j * strideOutput);
1067 ptr2 = inputMatrix + (j + startCol) * strideInput + startRow;
1069 for(i = 0; i < numRows_in; i++)
1071 *ptr1++ += alpha*(*ptr2++);
1074 for(i = 0; i < numRows_in; i++)
1083template<
typename OrdinalType,
typename ScalarType>
1093template<
typename OrdinalType,
typename ScalarType>
1098 printer.
obj.print(out);
1103template<
typename OrdinalType,
typename ScalarType>
1104SerialDenseMatrixPrinter<OrdinalType,ScalarType>
Templated interface class to BLAS routines.
Object for storing data and providing functionality that is common to all computational classes.
Teuchos header file which uses auto-configuration information to include necessary C++ headers.
#define TEUCHOS_CHK_REF(a)
#define TEUCHOS_MIN(x, y)
#define TEUCHOS_MAX(x, y)
#define TEUCHOS_CHK_ERR(a)
Teuchos::DataAccess Mode enumerable type.
Defines basic traits for the scalar field type.
Templated serial, dense, symmetric matrix class.
Functionality and data that is common to all computational classes.
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.
ScalarType & operator()(OrdinalType rowIndex, OrdinalType colIndex)
Element access method (non-const).
void swap(SerialDenseMatrix< OrdinalType, ScalarType > &B)
Swap values between this matrix and incoming matrix.
ScalarTraits< ScalarType >::magnitudeType normOne() const
Returns the 1-norm of the matrix.
SerialDenseMatrix< OrdinalType, ScalarType > & assign(const SerialDenseMatrix< OrdinalType, ScalarType > &Source)
Copies values from one matrix to another.
int reshape(OrdinalType numRows, OrdinalType numCols)
Reshaping method for changing the size of a SerialDenseMatrix, keeping the entries.
bool operator!=(const SerialDenseMatrix< OrdinalType, ScalarType > &Operand) const
Inequality of two matrices.
int scale(const ScalarType alpha)
Scale this matrix by alpha; *this = alpha**this.
SerialDenseMatrix< OrdinalType, ScalarType > & operator=(const SerialDenseMatrix< OrdinalType, ScalarType > &Source)
Copies values from one matrix to another.
bool empty() const
Returns whether this matrix is empty.
virtual ~SerialDenseMatrix()
Destructor.
static ScalarType * allocateValues(const OrdinalType numRows, const OrdinalType numCols)
SerialDenseMatrix< OrdinalType, ScalarType > & operator-=(const SerialDenseMatrix< OrdinalType, ScalarType > &Source)
Subtract another matrix from this matrix.
ScalarType * operator[](OrdinalType colIndex)
Column access method (non-const).
SerialDenseMatrix< OrdinalType, ScalarType > & operator+=(const SerialDenseMatrix< OrdinalType, ScalarType > &Source)
Add another matrix to this matrix.
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.
ScalarType scalarType
Typedef for scalar type.
ScalarTraits< ScalarType >::magnitudeType normInf() const
Returns the Infinity-norm of the matrix.
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.
void copyMat(ScalarType *inputMatrix, OrdinalType strideInput, OrdinalType numRows, OrdinalType numCols, ScalarType *outputMatrix, OrdinalType strideOutput, OrdinalType startRow, OrdinalType startCol, ScalarType alpha=ScalarTraits< ScalarType >::zero())
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.
bool operator==(const SerialDenseMatrix< OrdinalType, ScalarType > &Operand) const
Equality of two matrices.
virtual std::ostream & print(std::ostream &os) const
Print method. Defines the behavior of the std::ostream << operator.
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
Returns the Frobenius-norm of the matrix.
OrdinalType numCols() const
Returns the column dimension of this matrix.
SerialDenseMatrix()=default
Default Constructor.
SerialDenseMatrix< OrdinalType, ScalarType > & operator*=(const ScalarType alpha)
Scale this matrix by alpha; *this = alpha**this.
int shape(OrdinalType numRows, OrdinalType numCols)
Shape method for changing the size of a SerialDenseMatrix, initializing entries to zero.
OrdinalType ordinalType
Typedef for ordinal type.
void checkIndex(OrdinalType rowIndex, OrdinalType colIndex=0) const
This class creates and provides basic support for symmetric, positive-definite dense matrices of temp...
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Macro for throwing an exception with breakpointing to ease debugging.
TEUCHOSNUMERICS_LIB_DLL_EXPORT const char ETranspChar[]
std::ostream & operator<<(std::ostream &os, BigUInt< n > a)
SerialBandDenseMatrixPrinter< OrdinalType, ScalarType > printMat(const SerialBandDenseMatrix< OrdinalType, ScalarType > &obj)
Return SerialBandDenseMatrix ostream manipulator Use as:
TEUCHOSNUMERICS_LIB_DLL_EXPORT const char ESideChar[]
This structure defines some basic traits for a scalar field type.
static magnitudeType magnitude(T a)
Returns the magnitudeType of the scalar type a.
T magnitudeType
Mandatory typedef for result of magnitude.
static T zero()
Returns representation of zero for this scalar type.
static T conjugate(T a)
Returns the conjugate of the scalar type a.
static T random()
Returns a random number (between -one() and +one()) of this scalar type.
Ostream manipulator for SerialDenseMatrix.
const SerialDenseMatrix< OrdinalType, ScalarType > & obj
SerialDenseMatrixPrinter(const SerialDenseMatrix< OrdinalType, ScalarType > &obj_in)