42 #ifndef _TEUCHOS_SERIALDENSEMATRIX_HPP_ 43 #define _TEUCHOS_SERIALDENSEMATRIX_HPP_ 53 #include "Teuchos_Assert.hpp" 66 template<
typename OrdinalType,
typename ScalarType>
222 ScalarType&
operator () (OrdinalType rowIndex, OrdinalType colIndex);
232 const ScalarType&
operator () (OrdinalType rowIndex, OrdinalType colIndex)
const;
252 const ScalarType*
operator [] (OrdinalType colIndex)
const;
256 ScalarType*
values()
const {
return(values_); }
286 int scale (
const ScalarType alpha );
351 OrdinalType
numRows()
const {
return(numRows_); }
354 OrdinalType
numCols()
const {
return(numCols_); }
357 OrdinalType
stride()
const {
return(stride_); }
360 bool empty()
const {
return(numRows_ == 0 || numCols_ == 0); }
378 virtual void print(std::ostream& os)
const;
383 void copyMat(ScalarType* inputMatrix, OrdinalType strideInput,
384 OrdinalType
numRows, OrdinalType
numCols, ScalarType* outputMatrix,
385 OrdinalType strideOutput, OrdinalType startRow, OrdinalType startCol,
388 void checkIndex( OrdinalType rowIndex, OrdinalType colIndex = 0 )
const;
389 OrdinalType numRows_;
390 OrdinalType numCols_;
401 template<
typename OrdinalType,
typename ScalarType>
403 :
CompObject(), numRows_(0), numCols_(0), stride_(0), valuesCopied_(false), values_(0)
406 template<
typename OrdinalType,
typename ScalarType>
408 OrdinalType numRows_in, OrdinalType numCols_in,
bool zeroOut
410 :
CompObject(), numRows_(numRows_in), numCols_(numCols_in), stride_(numRows_in)
412 values_ =
new ScalarType[stride_*numCols_];
413 valuesCopied_ =
true;
418 template<
typename OrdinalType,
typename ScalarType>
420 DataAccess CV, ScalarType* values_in, OrdinalType stride_in, OrdinalType numRows_in,
421 OrdinalType numCols_in
423 :
CompObject(), numRows_(numRows_in), numCols_(numCols_in), stride_(stride_in),
424 valuesCopied_(false), values_(values_in)
429 values_ =
new ScalarType[stride_*numCols_];
430 copyMat(values_in, stride_in, numRows_, numCols_, values_, stride_, 0, 0);
431 valuesCopied_ =
true;
435 template<
typename OrdinalType,
typename ScalarType>
436 SerialDenseMatrix<OrdinalType, ScalarType>::SerialDenseMatrix(
const SerialDenseMatrix<OrdinalType, ScalarType> &Source,
ETransp trans) :
CompObject(),
BLAS<OrdinalType,ScalarType>(), numRows_(0), numCols_(0), stride_(0), valuesCopied_(true), values_(0)
440 numRows_ = Source.numRows_;
441 numCols_ = Source.numCols_;
443 if (!Source.valuesCopied_)
445 stride_ = Source.stride_;
446 values_ = Source.values_;
447 valuesCopied_ =
false;
452 const OrdinalType newsize = stride_ * numCols_;
454 values_ =
new ScalarType[newsize];
455 copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0);
458 numRows_ = 0; numCols_ = 0; stride_ = 0;
459 valuesCopied_ =
false;
465 numRows_ = Source.numCols_;
466 numCols_ = Source.numRows_;
468 values_ =
new ScalarType[stride_*numCols_];
469 for (OrdinalType j=0; j<numCols_; j++) {
470 for (OrdinalType i=0; i<numRows_; i++) {
477 numRows_ = Source.numCols_;
478 numCols_ = Source.numRows_;
480 values_ =
new ScalarType[stride_*numCols_];
481 for (OrdinalType j=0; j<numCols_; j++) {
482 for (OrdinalType i=0; i<numRows_; i++) {
483 values_[j*stride_ + i] = Source.values_[i*Source.stride_ + j];
490 template<
typename OrdinalType,
typename ScalarType>
494 :
CompObject(), numRows_(Source.numRows_), numCols_(Source.numCols_), stride_(Source.stride_),
495 valuesCopied_(false), values_(Source.values_)
500 values_ =
new ScalarType[stride_ * numCols_];
501 copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0);
502 valuesCopied_ =
true;
507 template<
typename OrdinalType,
typename ScalarType>
510 OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType startRow,
513 :
CompObject(), numRows_(numRows_in), numCols_(numCols_in), stride_(Source.stride_),
514 valuesCopied_(false), values_(Source.values_)
518 stride_ = numRows_in;
519 values_ =
new ScalarType[stride_ * numCols_in];
520 copyMat(Source.values_, Source.stride_, numRows_in, numCols_in, values_, stride_, startRow, startCol);
521 valuesCopied_ =
true;
525 values_ = values_ + (stride_ * startCol) + startRow;
529 template<
typename OrdinalType,
typename ScalarType>
539 template<
typename OrdinalType,
typename ScalarType>
541 OrdinalType numRows_in, OrdinalType numCols_in
545 numRows_ = numRows_in;
546 numCols_ = numCols_in;
548 values_ =
new ScalarType[stride_*numCols_];
550 valuesCopied_ =
true;
554 template<
typename OrdinalType,
typename ScalarType>
556 OrdinalType numRows_in, OrdinalType numCols_in
560 numRows_ = numRows_in;
561 numCols_ = numCols_in;
563 values_ =
new ScalarType[stride_*numCols_];
564 valuesCopied_ =
true;
568 template<
typename OrdinalType,
typename ScalarType>
570 OrdinalType numRows_in, OrdinalType numCols_in
574 ScalarType* values_tmp =
new ScalarType[numRows_in * numCols_in];
576 for(OrdinalType k = 0; k < numRows_in * numCols_in; k++)
578 values_tmp[k] = zero;
580 OrdinalType numRows_tmp = TEUCHOS_MIN(numRows_, numRows_in);
581 OrdinalType numCols_tmp = TEUCHOS_MIN(numCols_, numCols_in);
584 copyMat(values_, stride_, numRows_tmp, numCols_tmp, values_tmp,
588 numRows_ = numRows_in;
589 numCols_ = numCols_in;
591 values_ = values_tmp;
592 valuesCopied_ =
true;
600 template<
typename OrdinalType,
typename ScalarType>
604 for(OrdinalType j = 0; j < numCols_; j++)
606 for(OrdinalType i = 0; i < numRows_; i++)
608 values_[i + j*stride_] = value_in;
614 template<
typename OrdinalType,
typename ScalarType>
618 for(OrdinalType j = 0; j < numCols_; j++)
620 for(OrdinalType i = 0; i < numRows_; i++)
628 template<
typename OrdinalType,
typename ScalarType>
636 if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
640 if (!Source.valuesCopied_) {
645 numRows_ = Source.numRows_;
646 numCols_ = Source.numCols_;
647 stride_ = Source.stride_;
648 values_ = Source.values_;
653 numRows_ = Source.numRows_;
654 numCols_ = Source.numCols_;
655 stride_ = Source.numRows_;
656 const OrdinalType newsize = stride_ * numCols_;
658 values_ =
new ScalarType[newsize];
659 valuesCopied_ =
true;
667 if((Source.numRows_ <= stride_) && (Source.numCols_ == numCols_)) {
668 numRows_ = Source.numRows_;
669 numCols_ = Source.numCols_;
673 numRows_ = Source.numRows_;
674 numCols_ = Source.numCols_;
675 stride_ = Source.numRows_;
676 const OrdinalType newsize = stride_ * numCols_;
678 values_ =
new ScalarType[newsize];
679 valuesCopied_ =
true;
683 copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0);
688 template<
typename OrdinalType,
typename ScalarType>
692 if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_))
694 TEUCHOS_CHK_REF(*
this);
700 template<
typename OrdinalType,
typename ScalarType>
704 if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_))
706 TEUCHOS_CHK_REF(*
this);
712 template<
typename OrdinalType,
typename ScalarType>
716 if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
720 if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_))
722 TEUCHOS_CHK_REF(*
this);
724 copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0);
732 template<
typename OrdinalType,
typename ScalarType>
735 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK 736 checkIndex( rowIndex, colIndex );
738 return(values_[colIndex * stride_ + rowIndex]);
741 template<
typename OrdinalType,
typename ScalarType>
744 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK 745 checkIndex( rowIndex, colIndex );
747 return(values_[colIndex * stride_ + rowIndex]);
750 template<
typename OrdinalType,
typename ScalarType>
753 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK 754 checkIndex( 0, colIndex );
756 return(values_ + colIndex * stride_);
759 template<
typename OrdinalType,
typename ScalarType>
762 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK 763 checkIndex( 0, colIndex );
765 return(values_ + colIndex * stride_);
772 template<
typename OrdinalType,
typename ScalarType>
779 for(j = 0; j < numCols_; j++)
782 ptr = values_ + j * stride_;
783 for(i = 0; i < numRows_; i++)
797 template<
typename OrdinalType,
typename ScalarType>
803 for (i = 0; i < numRows_; i++) {
805 for (j=0; j< numCols_; j++) {
808 anorm = TEUCHOS_MAX( anorm, sum );
814 template<
typename OrdinalType,
typename ScalarType>
819 for (j = 0; j < numCols_; j++) {
820 for (i = 0; i < numRows_; i++) {
833 template<
typename OrdinalType,
typename ScalarType>
837 if((numRows_ != Operand.numRows_) || (numCols_ != Operand.numCols_))
844 for(i = 0; i < numRows_; i++)
846 for(j = 0; j < numCols_; j++)
848 if((*
this)(i, j) != Operand(i, j))
858 template<
typename OrdinalType,
typename ScalarType>
861 return !((*this) == Operand);
868 template<
typename OrdinalType,
typename ScalarType>
871 this->
scale( alpha );
875 template<
typename OrdinalType,
typename ScalarType>
881 for (j=0; j<numCols_; j++) {
882 ptr = values_ + j*stride_;
883 for (i=0; i<numRows_; i++) { *ptr = alpha * (*ptr); ptr++; }
889 template<
typename OrdinalType,
typename ScalarType>
896 if ((numRows_ != A.numRows_) || (numCols_ != A.numCols_))
900 for (j=0; j<numCols_; j++) {
901 ptr = values_ + j*stride_;
902 for (i=0; i<numRows_; i++) { *ptr = A(i,j) * (*ptr); ptr++; }
908 template<
typename OrdinalType,
typename ScalarType>
912 OrdinalType A_nrows = (ETranspChar[transa]!=
'N') ? A.
numCols() : A.
numRows();
913 OrdinalType A_ncols = (ETranspChar[transa]!=
'N') ? A.
numRows() : A.
numCols();
914 OrdinalType B_nrows = (ETranspChar[transb]!=
'N') ? B.
numCols() : B.
numRows();
915 OrdinalType B_ncols = (ETranspChar[transb]!=
'N') ? B.
numRows() : B.
numCols();
916 if ((numRows_ != A_nrows) || (A_ncols != B_nrows) || (numCols_ != B_ncols))
921 this->
GEMM(transa, transb, numRows_, numCols_, A_ncols, alpha, A.
values(), A.
stride(), B.
values(), B.
stride(), beta, values_, stride_);
922 double nflops = 2 * numRows_;
929 template<
typename OrdinalType,
typename ScalarType>
936 if (ESideChar[sideA]==
'L') {
937 if ((numRows_ != A_nrows) || (A_ncols != B_nrows) || (numCols_ != B_ncols)) {
941 if ((numRows_ != B_nrows) || (B_ncols != A_nrows) || (numCols_ != A_ncols)) {
948 this->
SYMM(sideA, uplo, numRows_, numCols_, alpha, A.
values(), A.
stride(), B.
values(), B.
stride(), beta, values_, stride_);
949 double nflops = 2 * numRows_;
956 template<
typename OrdinalType,
typename ScalarType>
961 os <<
"Values_copied : yes" << std::endl;
963 os <<
"Values_copied : no" << std::endl;
964 os <<
"Rows : " << numRows_ << std::endl;
965 os <<
"Columns : " << numCols_ << std::endl;
966 os <<
"LDA : " << stride_ << std::endl;
967 if(numRows_ == 0 || numCols_ == 0) {
968 os <<
"(matrix is empty, no values to display)" << std::endl;
970 for(OrdinalType i = 0; i < numRows_; i++) {
971 for(OrdinalType j = 0; j < numCols_; j++){
972 os << (*this)(i,j) <<
" ";
983 template<
typename OrdinalType,
typename ScalarType>
986 "SerialDenseMatrix<T>::checkIndex: " 987 "Row index " << rowIndex <<
" out of range [0, "<< numRows_ <<
")");
989 "SerialDenseMatrix<T>::checkIndex: " 990 "Col index " << colIndex <<
" out of range [0, "<< numCols_ <<
")");
993 template<
typename OrdinalType,
typename ScalarType>
1000 valuesCopied_ =
false;
1004 template<
typename OrdinalType,
typename ScalarType>
1006 ScalarType* inputMatrix, OrdinalType strideInput, OrdinalType numRows_in,
1007 OrdinalType numCols_in, ScalarType* outputMatrix, OrdinalType strideOutput,
1008 OrdinalType startRow, OrdinalType startCol, ScalarType alpha
1012 ScalarType* ptr1 = 0;
1013 ScalarType* ptr2 = 0;
1014 for(j = 0; j < numCols_in; j++) {
1015 ptr1 = outputMatrix + (j * strideOutput);
1016 ptr2 = inputMatrix + (j + startCol) * strideInput + startRow;
1018 for(i = 0; i < numRows_in; i++)
1020 *ptr1++ += alpha*(*ptr2++);
1023 for(i = 0; i < numRows_in; i++)
ScalarTraits< ScalarType >::magnitudeType normOne() const
Returns the 1-norm of the matrix.
ScalarType * values() const
Data array access method.
bool upper() const
Returns true if upper triangular part of this matrix has and will be used.
bool operator==(const SerialDenseMatrix< OrdinalType, ScalarType > &Operand) const
Equality of two matrices.
int shapeUninitialized(OrdinalType numRows, OrdinalType numCols)
Same as shape() except leaves uninitialized.
virtual void print(std::ostream &os) const
Print method. Defines the behavior of the std::ostream << operator inherited from the Object class...
Templated interface class to BLAS routines.
SerialDenseMatrix< OrdinalType, ScalarType > & operator-=(const SerialDenseMatrix< OrdinalType, ScalarType > &Source)
Subtract another matrix from this matrix.
SerialDenseMatrix< OrdinalType, ScalarType > & operator*=(const ScalarType alpha)
Scale this matrix by alpha; *this = alpha**this.
ScalarType * values() const
Returns the pointer to the ScalarType data array contained in the object.
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.
Teuchos header file which uses auto-configuration information to include necessary C++ headers...
This class creates and provides basic support for symmetric, positive-definite dense matrices of temp...
int scale(const ScalarType alpha)
Scale this matrix by alpha; *this = alpha**this.
OrdinalType numRows() const
Returns the row dimension of this matrix.
Object for storing data and providing functionality that is common to all computational classes...
virtual ~SerialDenseMatrix()
Destructor.
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
Returns the Frobenius-norm of the matrix.
This structure defines some basic traits for a scalar field type.
SerialDenseMatrix< OrdinalType, ScalarType > & operator+=(const SerialDenseMatrix< OrdinalType, ScalarType > &Source)
Add another matrix to this matrix.
void SYMM(ESide side, EUplo uplo, const OrdinalType m, const OrdinalType n, const alpha_type alpha, const A_type *A, const OrdinalType lda, const B_type *B, const OrdinalType ldb, const beta_type beta, ScalarType *C, const OrdinalType ldc) const
Performs the matrix-matrix operation: C <- alpha*A*B+beta*C or C <- alpha*B*A+beta*C where A is an m ...
SerialDenseMatrix< OrdinalType, ScalarType > & operator=(const SerialDenseMatrix< OrdinalType, ScalarType > &Source)
Copies values from one matrix to another.
Functionality and data that is common to all computational classes.
static T conjugate(T a)
Returns the conjugate of the scalar type a.
bool empty() const
Returns whether this matrix is empty.
ScalarType scalarType
Typedef for scalar type.
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero())
Set all values in the matrix to a constant value.
void GEMM(ETransp transa, ETransp transb, const OrdinalType m, const OrdinalType n, const OrdinalType k, const alpha_type alpha, const A_type *A, const OrdinalType lda, const B_type *B, const OrdinalType ldb, const beta_type beta, ScalarType *C, const OrdinalType ldc) const
General matrix-matrix multiply.
SerialDenseMatrix()
Default Constructor.
bool operator!=(const SerialDenseMatrix< OrdinalType, ScalarType > &Operand) const
Inequality of two matrices.
ScalarType & operator()(OrdinalType rowIndex, OrdinalType colIndex)
Element access method (non-const).
Templated serial, dense, symmetric matrix class.
OrdinalType ordinalType
Typedef for ordinal type.
ScalarType * operator[](OrdinalType colIndex)
Column access method (non-const).
static magnitudeType magnitude(T a)
Returns the magnitudeType of the scalar type a.
OrdinalType numCols() const
Returns the column dimension of this matrix.
ScalarTraits< ScalarType >::magnitudeType normInf() const
Returns the Infinity-norm of the matrix.
int random()
Set all values in the matrix to be random numbers.
The Teuchos namespace contains all of the classes, structs and enums used by Teuchos, as well as a number of utility routines.
int reshape(OrdinalType numRows, OrdinalType numCols)
Reshaping method for changing the size of a SerialDenseMatrix, keeping the entries.
void updateFlops(int addflops) const
Increment Flop count for this object.
Defines basic traits for the scalar field type.
static T zero()
Returns representation of zero for this scalar type.
static T random()
Returns a random number (between -one() and +one()) of this scalar type.
OrdinalType stride() const
Returns the stride between the columns of this matrix in memory.
int shape(OrdinalType numRows, OrdinalType numCols)
Shape method for changing the size of a SerialDenseMatrix, initializing entries to zero...
OrdinalType numCols() const
Returns the column dimension of this matrix.
SerialDenseMatrix< OrdinalType, ScalarType > & assign(const SerialDenseMatrix< OrdinalType, ScalarType > &Source)
Copies values from one matrix to another.
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...