10 #ifndef _TEUCHOS_SERIALDENSEMATRIX_HPP_
11 #define _TEUCHOS_SERIALDENSEMATRIX_HPP_
21 #include "Teuchos_Assert.hpp"
36 template<
typename OrdinalType,
typename ScalarType>
198 ScalarType&
operator () (OrdinalType rowIndex, OrdinalType colIndex);
208 const ScalarType&
operator () (OrdinalType rowIndex, OrdinalType colIndex)
const;
224 const ScalarType*
operator [] (OrdinalType colIndex)
const;
228 ScalarType*
values()
const {
return values_; }
258 int scale (
const ScalarType alpha );
323 OrdinalType
numRows()
const {
return(numRows_); }
326 OrdinalType
numCols()
const {
return(numCols_); }
329 OrdinalType
stride()
const {
return(stride_); }
332 bool empty()
const {
return(numRows_ == 0 || numCols_ == 0); }
350 virtual std::ostream&
print(std::ostream& os)
const;
355 void copyMat(ScalarType* inputMatrix, OrdinalType strideInput,
356 OrdinalType
numRows, OrdinalType
numCols, ScalarType* outputMatrix,
357 OrdinalType strideOutput, OrdinalType startRow, OrdinalType startCol,
360 void checkIndex( OrdinalType rowIndex, OrdinalType colIndex = 0 )
const;
363 allocateValues(
const OrdinalType numRows,
364 const OrdinalType numCols)
366 const size_t size = size_t(numRows) * size_t(numCols);
367 return new ScalarType[size];
370 OrdinalType numRows_ = 0;
371 OrdinalType numCols_ = 0;
372 OrdinalType stride_ = 0;
373 bool valuesCopied_ =
false;
374 ScalarType* values_ =
nullptr;
381 template<
typename OrdinalType,
typename ScalarType>
383 OrdinalType numRows_in, OrdinalType numCols_in,
bool zeroOut
385 : numRows_(numRows_in),
386 numCols_(numCols_in),
389 values_(allocateValues(numRows_in, numCols_in))
396 template<
typename OrdinalType,
typename ScalarType>
398 DataAccess CV, ScalarType* values_in, OrdinalType stride_in, OrdinalType numRows_in,
399 OrdinalType numCols_in
401 : numRows_(numRows_in),
402 numCols_(numCols_in),
404 valuesCopied_(false),
410 values_ = allocateValues(stride_, numCols_);
411 copyMat(values_in, stride_in, numRows_, numCols_, values_, stride_, 0, 0);
412 valuesCopied_ =
true;
416 template<
typename OrdinalType,
typename ScalarType>
418 : valuesCopied_(true)
422 numRows_ = Source.numRows_;
423 numCols_ = Source.numCols_;
425 if (!Source.valuesCopied_)
427 stride_ = Source.stride_;
428 values_ = Source.values_;
429 valuesCopied_ =
false;
434 if(stride_ > 0 && numCols_ > 0) {
435 values_ = allocateValues(stride_, numCols_);
436 copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0);
439 numRows_ = 0; numCols_ = 0; stride_ = 0;
440 valuesCopied_ =
false;
446 numRows_ = Source.numCols_;
447 numCols_ = Source.numRows_;
449 values_ = allocateValues(stride_, numCols_);
450 for (OrdinalType j=0; j<numCols_; j++) {
451 for (OrdinalType i=0; i<numRows_; i++) {
458 numRows_ = Source.numCols_;
459 numCols_ = Source.numRows_;
461 values_ = allocateValues(stride_, numCols_);
462 for (OrdinalType j=0; j<numCols_; j++) {
463 for (OrdinalType i=0; i<numRows_; i++) {
464 values_[j*stride_ + i] = Source.values_[i*Source.stride_ + j];
471 template<
typename OrdinalType,
typename ScalarType>
475 : numRows_(Source.numRows_), numCols_(Source.numCols_), stride_(Source.stride_),
476 valuesCopied_(false), values_(Source.values_)
481 values_ = allocateValues(stride_, numCols_);
482 copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0);
483 valuesCopied_ =
true;
488 template<
typename OrdinalType,
typename ScalarType>
491 OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType startRow,
494 :
CompObject(), numRows_(numRows_in), numCols_(numCols_in), stride_(Source.stride_),
495 valuesCopied_(false), values_(Source.values_)
499 stride_ = numRows_in;
500 values_ = allocateValues(stride_, numCols_in);
501 copyMat(Source.values_, Source.stride_, numRows_in, numCols_in, values_, stride_, startRow, startCol);
502 valuesCopied_ =
true;
506 values_ = values_ + (stride_ * startCol) + startRow;
510 template<
typename OrdinalType,
typename ScalarType>
520 template<
typename OrdinalType,
typename ScalarType>
522 OrdinalType numRows_in, OrdinalType numCols_in
526 numRows_ = numRows_in;
527 numCols_ = numCols_in;
529 values_ = allocateValues(stride_, numCols_);
531 valuesCopied_ =
true;
535 template<
typename OrdinalType,
typename ScalarType>
537 OrdinalType numRows_in, OrdinalType numCols_in
541 numRows_ = numRows_in;
542 numCols_ = numCols_in;
544 values_ = allocateValues(stride_, numCols_);
545 valuesCopied_ =
true;
549 template<
typename OrdinalType,
typename ScalarType>
551 OrdinalType numRows_in, OrdinalType numCols_in
555 ScalarType* values_tmp = allocateValues(numRows_in, numCols_in);
557 for(OrdinalType k = 0; k < numRows_in * numCols_in; k++)
559 values_tmp[k] = zero;
561 OrdinalType numRows_tmp = TEUCHOS_MIN(numRows_, numRows_in);
562 OrdinalType numCols_tmp = TEUCHOS_MIN(numCols_, numCols_in);
565 copyMat(values_, stride_, numRows_tmp, numCols_tmp, values_tmp,
569 numRows_ = numRows_in;
570 numCols_ = numCols_in;
572 values_ = values_tmp;
573 valuesCopied_ =
true;
581 template<
typename OrdinalType,
typename ScalarType>
585 for(OrdinalType j = 0; j < numCols_; j++)
587 for(OrdinalType i = 0; i < numRows_; i++)
589 values_[i + j*stride_] = value_in;
595 template<
typename OrdinalType,
typename ScalarType>
void
613 std::swap(values_ , B.values_);
614 std::swap(numRows_, B.numRows_);
615 std::swap(numCols_, B.numCols_);
616 std::swap(stride_, B.stride_);
617 std::swap(valuesCopied_, B.valuesCopied_);
620 template<
typename OrdinalType,
typename ScalarType>
624 for(OrdinalType j = 0; j < numCols_; j++)
626 for(OrdinalType i = 0; i < numRows_; i++)
634 template<
typename OrdinalType,
typename ScalarType>
642 if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
646 if (!Source.valuesCopied_) {
651 numRows_ = Source.numRows_;
652 numCols_ = Source.numCols_;
653 stride_ = Source.stride_;
654 values_ = Source.values_;
659 numRows_ = Source.numRows_;
660 numCols_ = Source.numCols_;
661 stride_ = Source.numRows_;
662 if(stride_ > 0 && numCols_ > 0) {
663 values_ = allocateValues(stride_, numCols_);
664 valuesCopied_ =
true;
672 if((Source.numRows_ <= stride_) && (Source.numCols_ == numCols_)) {
673 numRows_ = Source.numRows_;
674 numCols_ = Source.numCols_;
678 numRows_ = Source.numRows_;
679 numCols_ = Source.numCols_;
680 stride_ = Source.numRows_;
681 if(stride_ > 0 && numCols_ > 0) {
682 values_ = allocateValues(stride_, numCols_);
683 valuesCopied_ =
true;
687 copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0);
692 template<
typename OrdinalType,
typename ScalarType>
696 if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_))
698 TEUCHOS_CHK_REF(*
this);
704 template<
typename OrdinalType,
typename ScalarType>
708 if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_))
710 TEUCHOS_CHK_REF(*
this);
716 template<
typename OrdinalType,
typename ScalarType>
720 if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
724 if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_))
726 TEUCHOS_CHK_REF(*
this);
728 copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0);
736 template<
typename OrdinalType,
typename ScalarType>
739 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
740 checkIndex( rowIndex, colIndex );
742 return(values_[colIndex * stride_ + rowIndex]);
745 template<
typename OrdinalType,
typename ScalarType>
748 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
749 checkIndex( rowIndex, colIndex );
751 return(values_[colIndex * stride_ + rowIndex]);
754 template<
typename OrdinalType,
typename ScalarType>
757 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
758 checkIndex( 0, colIndex );
760 return(values_ + colIndex * stride_);
763 template<
typename OrdinalType,
typename ScalarType>
766 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
767 checkIndex( 0, colIndex );
769 return(values_ + colIndex * stride_);
776 template<
typename OrdinalType,
typename ScalarType>
783 for(j = 0; j < numCols_; j++)
786 ptr = values_ + j * stride_;
787 for(i = 0; i < numRows_; i++)
798 if (flopCounter_!=0) updateFlops(numRows_ * numCols_);
802 template<
typename OrdinalType,
typename ScalarType>
808 for (i = 0; i < numRows_; i++) {
810 for (j=0; j< numCols_; j++) {
813 anorm = TEUCHOS_MAX( anorm, sum );
816 if (flopCounter_!=0) updateFlops(numRows_ * numCols_);
820 template<
typename OrdinalType,
typename ScalarType>
825 for (j = 0; j < numCols_; j++) {
826 for (i = 0; i < numRows_; i++) {
832 if (flopCounter_!=0) updateFlops(numRows_ * numCols_);
840 template<
typename OrdinalType,
typename ScalarType>
844 if((numRows_ != Operand.numRows_) || (numCols_ != Operand.numCols_))
851 for(i = 0; i < numRows_; i++)
853 for(j = 0; j < numCols_; j++)
855 if((*
this)(i, j) != Operand(i, j))
865 template<
typename OrdinalType,
typename ScalarType>
868 return !((*this) == Operand);
875 template<
typename OrdinalType,
typename ScalarType>
878 this->scale( alpha );
882 template<
typename OrdinalType,
typename ScalarType>
888 for (j=0; j<numCols_; j++) {
889 ptr = values_ + j*stride_;
890 for (i=0; i<numRows_; i++) { *ptr = alpha * (*ptr); ptr++; }
893 if (flopCounter_!=0) updateFlops( numRows_*numCols_ );
897 template<
typename OrdinalType,
typename ScalarType>
904 if ((numRows_ != A.numRows_) || (numCols_ != A.numCols_))
908 for (j=0; j<numCols_; j++) {
909 ptr = values_ + j*stride_;
910 for (i=0; i<numRows_; i++) { *ptr = A(i,j) * (*ptr); ptr++; }
913 if (flopCounter_!=0) updateFlops( numRows_*numCols_ );
917 template<
typename OrdinalType,
typename ScalarType>
921 OrdinalType A_nrows = (ETranspChar[transa]!=
'N') ? A.
numCols() : A.
numRows();
922 OrdinalType A_ncols = (ETranspChar[transa]!=
'N') ? A.
numRows() : A.
numCols();
923 OrdinalType B_nrows = (ETranspChar[transb]!=
'N') ? B.
numCols() : B.
numRows();
924 OrdinalType B_ncols = (ETranspChar[transb]!=
'N') ? B.
numRows() : B.
numCols();
925 if ((numRows_ != A_nrows) || (A_ncols != B_nrows) || (numCols_ != B_ncols))
930 this->GEMM(transa, transb, numRows_, numCols_, A_ncols, alpha, A.
values(), A.
stride(), B.
values(), B.
stride(), beta, values_, stride_);
933 if (flopCounter_!=0) {
934 double nflops = 2 * numRows_;
942 template<
typename OrdinalType,
typename ScalarType>
949 if (ESideChar[sideA]==
'L') {
950 if ((numRows_ != A_nrows) || (A_ncols != B_nrows) || (numCols_ != B_ncols)) {
954 if ((numRows_ != B_nrows) || (B_ncols != A_nrows) || (numCols_ != A_ncols)) {
961 this->SYMM(sideA, uplo, numRows_, numCols_, alpha, A.
values(), A.
stride(), B.
values(), B.
stride(), beta, values_, stride_);
964 if (flopCounter_!=0) {
965 double nflops = 2 * numRows_;
973 template<
typename OrdinalType,
typename ScalarType>
978 os <<
"Values_copied : yes" << std::endl;
980 os <<
"Values_copied : no" << std::endl;
981 os <<
"Rows : " << numRows_ << std::endl;
982 os <<
"Columns : " << numCols_ << std::endl;
983 os <<
"LDA : " << stride_ << std::endl;
984 if(numRows_ == 0 || numCols_ == 0) {
985 os <<
"(matrix is empty, no values to display)" << std::endl;
987 for(OrdinalType i = 0; i < numRows_; i++) {
988 for(OrdinalType j = 0; j < numCols_; j++){
989 os << (*this)(i,j) <<
" ";
1001 template<
typename OrdinalType,
typename ScalarType>
1004 "SerialDenseMatrix<T>::checkIndex: "
1005 "Row index " << rowIndex <<
" out of range [0, "<< numRows_ <<
")");
1007 "SerialDenseMatrix<T>::checkIndex: "
1008 "Col index " << colIndex <<
" out of range [0, "<< numCols_ <<
")");
1011 template<
typename OrdinalType,
typename ScalarType>
1012 void SerialDenseMatrix<OrdinalType, ScalarType>::deleteArrays(
void)
1018 valuesCopied_ =
false;
1022 template<
typename OrdinalType,
typename ScalarType>
1023 void SerialDenseMatrix<OrdinalType, ScalarType>::copyMat(
1024 ScalarType* inputMatrix, OrdinalType strideInput, OrdinalType numRows_in,
1025 OrdinalType numCols_in, ScalarType* outputMatrix, OrdinalType strideOutput,
1026 OrdinalType startRow, OrdinalType startCol, ScalarType alpha
1030 ScalarType* ptr1 = 0;
1031 ScalarType* ptr2 = 0;
1032 for(j = 0; j < numCols_in; j++) {
1033 ptr1 = outputMatrix + (j * strideOutput);
1034 ptr2 = inputMatrix + (j + startCol) * strideInput + startRow;
1036 for(i = 0; i < numRows_in; i++)
1038 *ptr1++ += alpha*(*ptr2++);
1041 for(i = 0; i < numRows_in; i++)
1050 template<
typename OrdinalType,
typename ScalarType>
1060 template<
typename OrdinalType,
typename ScalarType>
1062 operator<<(std::ostream &out,
1065 printer.obj.print(out);
1070 template<
typename OrdinalType,
typename ScalarType>
1071 SerialDenseMatrixPrinter<OrdinalType,ScalarType>
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.
Templated interface class to BLAS routines.
SerialDenseMatrix< OrdinalType, ScalarType > & operator-=(const SerialDenseMatrix< OrdinalType, ScalarType > &Source)
Subtract another matrix from this matrix.
SerialBandDenseMatrixPrinter< OrdinalType, ScalarType > printMat(const SerialBandDenseMatrix< OrdinalType, ScalarType > &obj)
Return SerialBandDenseMatrix ostream manipulator Use as:
SerialDenseMatrix< OrdinalType, ScalarType > & operator*=(const ScalarType alpha)
Scale this matrix by alpha; *this = alpha**this.
void swap(SerialDenseMatrix< OrdinalType, ScalarType > &B)
Swap values between this matrix and incoming matrix.
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.
virtual std::ostream & print(std::ostream &os) const
Print method. Defines the behavior of the std::ostream << operator.
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.
SerialDenseMatrix()=default
Default Constructor.
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.
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.
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.
Ostream manipulator for SerialDenseMatrix.
int reshape(OrdinalType numRows, OrdinalType numCols)
Reshaping method for changing the size of a SerialDenseMatrix, keeping the entries.
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...