41#ifndef _TEUCHOS_SERIALSYMDENSEMATRIX_HPP_
42#define _TEUCHOS_SERIALSYMDENSEMATRIX_HPP_
52#include "Teuchos_Assert.hpp"
120template<
typename OrdinalType,
typename ScalarType>
202 int shape(OrdinalType numRowsCols);
227 int reshape(OrdinalType numRowsCols);
298 ScalarType&
operator () (OrdinalType rowIndex, OrdinalType colIndex);
308 const ScalarType&
operator () (OrdinalType rowIndex, OrdinalType colIndex)
const;
313 ScalarType*
values()
const {
return(values_); }
321 bool upper()
const {
return(upper_);};
324 char UPLO()
const {
return(UPLO_);};
374 OrdinalType
numRows()
const {
return(numRowCols_); }
377 OrdinalType
numCols()
const {
return(numRowCols_); }
380 OrdinalType
stride()
const {
return(stride_); }
383 bool empty()
const {
return(numRowCols_ == 0); }
403 virtual std::ostream& print(std::ostream& os)
const;
410 void scale(
const ScalarType alpha );
413 void copyMat(
bool inputUpper, ScalarType* inputMatrix, OrdinalType inputStride,
414 OrdinalType numRowCols,
bool outputUpper, ScalarType* outputMatrix,
415 OrdinalType outputStride, OrdinalType startRowCol,
419 void copyUPLOMat(
bool inputUpper, ScalarType* inputMatrix,
420 OrdinalType inputStride, OrdinalType inputRows);
423 void checkIndex( OrdinalType rowIndex, OrdinalType colIndex = 0 )
const;
426 allocateValues(
const OrdinalType
numRows,
430#pragma GCC diagnostic push
431#pragma GCC diagnostic ignored "-Wvla"
432 return new ScalarType[size];
433#pragma GCC diagnostic pop
436 OrdinalType numRowCols_ = 0;
437 OrdinalType stride_ = 0;
438 bool valuesCopied_ =
false;
439 ScalarType* values_ =
nullptr;
448template<
typename OrdinalType,
typename ScalarType>
450 : numRowCols_(numRowCols_in), stride_(numRowCols_in), valuesCopied_(false), upper_(false), UPLO_(
'L')
452 values_ = allocateValues(stride_, numRowCols_);
453 valuesCopied_ =
true;
460template<
typename OrdinalType,
typename ScalarType>
462 DataAccess CV,
bool upper_in, ScalarType* values_in, OrdinalType stride_in, OrdinalType numRowCols_in
464 : numRowCols_(numRowCols_in), stride_(stride_in), valuesCopied_(false),
465 values_(values_in), upper_(upper_in)
474 stride_ = numRowCols_;
475 values_ = allocateValues(stride_, numRowCols_);
476 copyMat(upper_in, values_in, stride_in, numRowCols_, upper_, values_, stride_, 0);
477 valuesCopied_ =
true;
481template<
typename OrdinalType,
typename ScalarType>
483 : numRowCols_(Source.numRowCols_), stride_(0), valuesCopied_(true),
484 values_(0), upper_(Source.upper_), UPLO_(Source.UPLO_)
486 if (!Source.valuesCopied_)
488 stride_ = Source.stride_;
489 values_ = Source.values_;
490 valuesCopied_ =
false;
494 stride_ = numRowCols_;
495 if(stride_ > 0 && numRowCols_ > 0) {
496 values_ = allocateValues(stride_, numRowCols_);
497 copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0);
502 valuesCopied_ =
false;
507template<
typename OrdinalType,
typename ScalarType>
510 ScalarType> &Source, OrdinalType numRowCols_in, OrdinalType startRowCol )
512 numRowCols_(numRowCols_in), stride_(Source.stride_), valuesCopied_(false), upper_(Source.upper_), UPLO_(Source.UPLO_)
516 stride_ = numRowCols_in;
517 values_ = allocateValues(stride_, numRowCols_in);
518 copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_in, upper_, values_, stride_, startRowCol);
519 valuesCopied_ =
true;
523 values_ = Source.values_ + (stride_ * startRowCol) + startRowCol;
527template<
typename OrdinalType,
typename ScalarType>
537template<
typename OrdinalType,
typename ScalarType>
541 numRowCols_ = numRowCols_in;
542 stride_ = numRowCols_;
543 values_ = allocateValues(stride_, numRowCols_);
545 valuesCopied_ =
true;
549template<
typename OrdinalType,
typename ScalarType>
553 numRowCols_ = numRowCols_in;
554 stride_ = numRowCols_;
555 values_ = allocateValues(stride_, numRowCols_);
556 valuesCopied_ =
true;
560template<
typename OrdinalType,
typename ScalarType>
564 ScalarType* values_tmp = allocateValues(numRowCols_in, numRowCols_in);
566 for(OrdinalType k = 0; k < numRowCols_in * numRowCols_in; k++)
568 values_tmp[k] = zero;
570 OrdinalType numRowCols_tmp = TEUCHOS_MIN(numRowCols_, numRowCols_in);
573 copyMat(upper_, values_, stride_, numRowCols_tmp, upper_, values_tmp, numRowCols_in, 0);
576 numRowCols_ = numRowCols_in;
577 stride_ = numRowCols_;
578 values_ = values_tmp;
579 valuesCopied_ =
true;
587template<
typename OrdinalType,
typename ScalarType>
591 if (upper_ !=
false) {
592 copyUPLOMat(
true, values_, stride_, numRowCols_ );
598template<
typename OrdinalType,
typename ScalarType>
602 if (upper_ ==
false) {
603 copyUPLOMat(
false, values_, stride_, numRowCols_ );
609template<
typename OrdinalType,
typename ScalarType>
613 if (fullMatrix ==
true) {
614 for(OrdinalType j = 0; j < numRowCols_; j++)
616 for(OrdinalType i = 0; i < numRowCols_; i++)
618 values_[i + j*stride_] = value_in;
625 for(OrdinalType j = 0; j < numRowCols_; j++) {
626 for(OrdinalType i = 0; i <= j; i++) {
627 values_[i + j*stride_] = value_in;
632 for(OrdinalType j = 0; j < numRowCols_; j++) {
633 for(OrdinalType i = j; i < numRowCols_; i++) {
634 values_[i + j*stride_] = value_in;
642template<
typename OrdinalType,
typename ScalarType>
void
646 std::swap(values_ , B.values_);
647 std::swap(numRowCols_, B.numRowCols_);
648 std::swap(stride_, B.stride_);
649 std::swap(valuesCopied_, B.valuesCopied_);
650 std::swap(upper_, B.upper_);
651 std::swap(UPLO_, B.UPLO_);
654template<
typename OrdinalType,
typename ScalarType>
662 for(OrdinalType j = 0; j < numRowCols_; j++) {
663 for(OrdinalType i = 0; i < j; i++) {
671 for(OrdinalType j = 0; j < numRowCols_; j++) {
672 for(OrdinalType i = j+1; i < numRowCols_; i++) {
681 for(OrdinalType i = 0; i < numRowCols_; i++) {
682 values_[i + i*stride_] = diagSum[i] + bias;
687template<
typename OrdinalType,
typename ScalarType>
693 if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_)) {
694 upper_ = Source.upper_;
699 if (!Source.valuesCopied_) {
704 numRowCols_ = Source.numRowCols_;
705 stride_ = Source.stride_;
706 values_ = Source.values_;
707 upper_ = Source.upper_;
708 UPLO_ = Source.UPLO_;
713 numRowCols_ = Source.numRowCols_;
714 stride_ = Source.numRowCols_;
715 upper_ = Source.upper_;
716 UPLO_ = Source.UPLO_;
717 if(stride_ > 0 && numRowCols_ > 0) {
718 values_ = allocateValues(stride_, numRowCols_);
719 valuesCopied_ =
true;
727 if((Source.numRowCols_ <= stride_) && (Source.numRowCols_ == numRowCols_)) {
728 numRowCols_ = Source.numRowCols_;
729 upper_ = Source.upper_;
730 UPLO_ = Source.UPLO_;
734 numRowCols_ = Source.numRowCols_;
735 stride_ = Source.numRowCols_;
736 upper_ = Source.upper_;
737 UPLO_ = Source.UPLO_;
738 if(stride_ > 0 && numRowCols_ > 0) {
739 values_ = allocateValues(stride_, numRowCols_);
740 valuesCopied_ =
true;
744 copyMat(Source.upper_, Source.values_, Source.stride_, Source.numRowCols_, upper_, values_, stride_, 0);
749template<
typename OrdinalType,
typename ScalarType>
756template<
typename OrdinalType,
typename ScalarType>
760 if ((numRowCols_ != Source.numRowCols_))
762 TEUCHOS_CHK_REF(*
this);
768template<
typename OrdinalType,
typename ScalarType>
772 if ((numRowCols_ != Source.numRowCols_))
774 TEUCHOS_CHK_REF(*
this);
780template<
typename OrdinalType,
typename ScalarType>
784 if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_)) {
785 upper_ = Source.upper_;
790 if ((numRowCols_ != Source.numRowCols_))
792 TEUCHOS_CHK_REF(*
this);
794 copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0 );
802template<
typename OrdinalType,
typename ScalarType>
805#ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
806 checkIndex( rowIndex, colIndex );
808 if ( rowIndex <= colIndex ) {
811 return(values_[colIndex * stride_ + rowIndex]);
813 return(values_[rowIndex * stride_ + colIndex]);
818 return(values_[rowIndex * stride_ + colIndex]);
820 return(values_[colIndex * stride_ + rowIndex]);
824template<
typename OrdinalType,
typename ScalarType>
827#ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
828 checkIndex( rowIndex, colIndex );
830 if ( rowIndex <= colIndex ) {
833 return(values_[colIndex * stride_ + rowIndex]);
835 return(values_[rowIndex * stride_ + colIndex]);
840 return(values_[rowIndex * stride_ + colIndex]);
842 return(values_[colIndex * stride_ + rowIndex]);
850template<
typename OrdinalType,
typename ScalarType>
856template<
typename OrdinalType,
typename ScalarType>
867 for (j=0; j<numRowCols_; j++) {
869 ptr = values_ + j*stride_;
870 for (i=0; i<j; i++) {
873 ptr = values_ + j + j*stride_;
874 for (i=j; i<numRowCols_; i++) {
878 anorm = TEUCHOS_MAX( anorm, sum );
882 for (j=0; j<numRowCols_; j++) {
884 ptr = values_ + j + j*stride_;
885 for (i=j; i<numRowCols_; i++) {
889 for (i=0; i<j; i++) {
893 anorm = TEUCHOS_MAX( anorm, sum );
899template<
typename OrdinalType,
typename ScalarType>
909 for (j = 0; j < numRowCols_; j++) {
910 for (i = 0; i < j; i++) {
917 for (j = 0; j < numRowCols_; j++) {
919 for (i = j+1; i < numRowCols_; i++) {
932template<
typename OrdinalType,
typename ScalarType>
936 if((numRowCols_ != Operand.numRowCols_)) {
941 for(i = 0; i < numRowCols_; i++) {
942 for(j = 0; j < numRowCols_; j++) {
943 if((*
this)(i, j) != Operand(i, j)) {
952template<
typename OrdinalType,
typename ScalarType>
955 return !((*this) == Operand);
962template<
typename OrdinalType,
typename ScalarType>
969 for (j=0; j<numRowCols_; j++) {
970 ptr = values_ + j*stride_;
971 for (i=0; i<= j; i++) { *ptr = alpha * (*ptr); ptr++; }
975 for (j=0; j<numRowCols_; j++) {
976 ptr = values_ + j*stride_ + j;
977 for (i=j; i<numRowCols_; i++) { *ptr = alpha * (*ptr); ptr++; }
1011template<
typename OrdinalType,
typename ScalarType>
1016 os <<
"Values_copied : yes" << std::endl;
1018 os <<
"Values_copied : no" << std::endl;
1019 os <<
"Rows / Columns : " << numRowCols_ << std::endl;
1020 os <<
"LDA : " << stride_ << std::endl;
1022 os <<
"Storage: Upper" << std::endl;
1024 os <<
"Storage: Lower" << std::endl;
1025 if(numRowCols_ == 0) {
1026 os <<
"(matrix is empty, no values to display)" << std::endl;
1028 for(OrdinalType i = 0; i < numRowCols_; i++) {
1029 for(OrdinalType j = 0; j < numRowCols_; j++){
1030 os << (*this)(i,j) <<
" ";
1042template<
typename OrdinalType,
typename ScalarType>
1045 "SerialSymDenseMatrix<T>::checkIndex: "
1046 "Row index " << rowIndex <<
" out of range [0, "<< numRowCols_ <<
")");
1048 "SerialSymDenseMatrix<T>::checkIndex: "
1049 "Col index " << colIndex <<
" out of range [0, "<< numRowCols_ <<
")");
1052template<
typename OrdinalType,
typename ScalarType>
1053void SerialSymDenseMatrix<OrdinalType, ScalarType>::deleteArrays(
void)
1059 valuesCopied_ =
false;
1063template<
typename OrdinalType,
typename ScalarType>
1064void SerialSymDenseMatrix<OrdinalType, ScalarType>::copyMat(
1065 bool inputUpper, ScalarType* inputMatrix,
1066 OrdinalType inputStride, OrdinalType numRowCols_in,
1067 bool outputUpper, ScalarType* outputMatrix,
1068 OrdinalType outputStride, OrdinalType startRowCol,
1073 ScalarType* ptr1 = 0;
1074 ScalarType* ptr2 = 0;
1076 for (j = 0; j < numRowCols_in; j++) {
1077 if (inputUpper ==
true) {
1079 ptr2 = inputMatrix + (j + startRowCol) * inputStride + startRowCol;
1080 if (outputUpper ==
true) {
1082 ptr1 = outputMatrix + j*outputStride;
1084 for(i = 0; i <= j; i++) {
1085 *ptr1++ += alpha*(*ptr2++);
1088 for(i = 0; i <= j; i++) {
1096 ptr1 = outputMatrix + j;
1098 for(i = 0; i <= j; i++) {
1099 *ptr1 += alpha*(*ptr2++);
1100 ptr1 += outputStride;
1103 for(i = 0; i <= j; i++) {
1105 ptr1 += outputStride;
1112 ptr2 = inputMatrix + (startRowCol+j) * inputStride + startRowCol + j;
1113 if (outputUpper ==
true) {
1116 ptr1 = outputMatrix + j*outputStride + j;
1118 for(i = j; i < numRowCols_in; i++) {
1119 *ptr1 += alpha*(*ptr2++);
1120 ptr1 += outputStride;
1123 for(i = j; i < numRowCols_in; i++) {
1125 ptr1 += outputStride;
1131 ptr1 = outputMatrix + j*outputStride + j;
1133 for(i = j; i < numRowCols_in; i++) {
1134 *ptr1++ += alpha*(*ptr2++);
1137 for(i = j; i < numRowCols_in; i++) {
1146template<
typename OrdinalType,
typename ScalarType>
1147void SerialSymDenseMatrix<OrdinalType, ScalarType>::copyUPLOMat(
1148 bool inputUpper, ScalarType* inputMatrix,
1149 OrdinalType inputStride, OrdinalType inputRows
1153 ScalarType * ptr1 = 0;
1154 ScalarType * ptr2 = 0;
1157 for (j=1; j<inputRows; j++) {
1158 ptr1 = inputMatrix + j;
1159 ptr2 = inputMatrix + j*inputStride;
1160 for (i=0; i<j; i++) {
1167 for (i=1; i<inputRows; i++) {
1168 ptr1 = inputMatrix + i;
1169 ptr2 = inputMatrix + i*inputStride;
1170 for (j=0; j<i; j++) {
1179template<
typename OrdinalType,
typename ScalarType>
1189template<
typename OrdinalType,
typename ScalarType>
1191operator<<(std::ostream &out,
1194 printer.obj.print(out);
1199template<
typename OrdinalType,
typename ScalarType>
1200SerialSymDenseMatrixPrinter<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.
Teuchos::DataAccess Mode enumerable type.
Defines basic traits for the scalar field type.
Functionality and data that is common to all computational classes.
This class creates and provides basic support for symmetric, positive-definite dense matrices of temp...
SerialSymDenseMatrix()=default
Default constructor; defines a zero size object.
bool operator==(const SerialSymDenseMatrix< OrdinalType, ScalarType > &Operand) const
Equality of two matrices.
int shapeUninitialized(OrdinalType numRowsCols)
Set dimensions of a Teuchos::SerialSymDenseMatrix object; don't initialize values.
int shape(OrdinalType numRowsCols)
Set dimensions of a Teuchos::SerialSymDenseMatrix object; init values to zero.
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
Returns the Frobenius-norm of the matrix.
ScalarType & operator()(OrdinalType rowIndex, OrdinalType colIndex)
Element access method (non-const).
OrdinalType stride() const
Returns the stride between the columns of this matrix in memory.
SerialSymDenseMatrix< OrdinalType, ScalarType > & operator-=(const SerialSymDenseMatrix< OrdinalType, ScalarType > &Source)
Subtract another matrix from this matrix.
virtual std::ostream & print(std::ostream &os) const
Print method. Defines the behavior of the std::ostream << operator.
SerialSymDenseMatrix< OrdinalType, ScalarType > & operator+=(const SerialSymDenseMatrix< OrdinalType, ScalarType > &Source)
Add another matrix to this matrix.
bool operator!=(const SerialSymDenseMatrix< OrdinalType, ScalarType > &Operand) const
Inequality of two matrices.
void swap(SerialSymDenseMatrix< OrdinalType, ScalarType > &B)
Swap values between this matrix and incoming matrix.
OrdinalType numCols() const
Returns the column dimension of this matrix.
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero(), bool fullMatrix=false)
Set all values in the matrix to a constant value.
virtual ~SerialSymDenseMatrix()
Teuchos::SerialSymDenseMatrix destructor.
char UPLO() const
Returns character value of UPLO used by LAPACK routines.
SerialSymDenseMatrix< OrdinalType, ScalarType > & operator*=(const ScalarType alpha)
Inplace scalar-matrix product A = alpha*A.
int random(const ScalarType bias=0.1 *Teuchos::ScalarTraits< ScalarType >::one())
Set all values in the active area (upper/lower triangle) of this matrix to be random numbers.
ScalarType scalarType
Typedef for scalar type.
bool upper() const
Returns true if upper triangular part of this matrix has and will be used.
ScalarTraits< ScalarType >::magnitudeType normInf() const
Returns the Infinity-norm of the matrix.
ScalarType * values() const
Returns the pointer to the ScalarType data array contained in the object.
SerialSymDenseMatrix< OrdinalType, ScalarType > & operator=(const SerialSymDenseMatrix< OrdinalType, ScalarType > &Source)
Copies values from one matrix to another.
bool empty() const
Returns whether this matrix is empty.
SerialSymDenseMatrix< OrdinalType, ScalarType > & assign(const SerialSymDenseMatrix< OrdinalType, ScalarType > &Source)
Copies values from one matrix to another.
void setLower()
Specify that the lower triangle of the this matrix should be used.
int reshape(OrdinalType numRowsCols)
Reshape a Teuchos::SerialSymDenseMatrix object.
OrdinalType numRows() const
Returns the row dimension of this matrix.
void setUpper()
Specify that the upper triangle of the this matrix should be used.
ScalarTraits< ScalarType >::magnitudeType normOne() const
Returns the 1-norm of the matrix.
OrdinalType ordinalType
Typedef for ordinal type.
#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,...
SerialBandDenseMatrixPrinter< OrdinalType, ScalarType > printMat(const SerialBandDenseMatrix< OrdinalType, ScalarType > &obj)
Return SerialBandDenseMatrix ostream manipulator Use as:
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 random()
Returns a random number (between -one() and +one()) of this scalar type.
Ostream manipulator for SerialSymDenseMatrix.