43 #ifndef _TEUCHOS_SERIALSYMDENSEMATRIX_HPP_ 44 #define _TEUCHOS_SERIALSYMDENSEMATRIX_HPP_ 54 #include "Teuchos_Assert.hpp" 120 template<
typename OrdinalType,
typename ScalarType>
202 int shape(OrdinalType numRowsCols);
227 int reshape(OrdinalType numRowsCols);
292 ScalarType&
operator () (OrdinalType rowIndex, OrdinalType colIndex);
302 const ScalarType&
operator () (OrdinalType rowIndex, OrdinalType colIndex)
const;
307 ScalarType*
values()
const {
return(values_); }
315 bool upper()
const {
return(upper_);};
318 char UPLO()
const {
return(UPLO_);};
368 OrdinalType
numRows()
const {
return(numRowCols_); }
371 OrdinalType
numCols()
const {
return(numRowCols_); }
374 OrdinalType
stride()
const {
return(stride_); }
377 bool empty()
const {
return(numRowCols_ == 0); }
396 virtual void print(std::ostream& os)
const;
404 void scale(
const ScalarType alpha );
407 void copyMat(
bool inputUpper, ScalarType* inputMatrix, OrdinalType inputStride,
408 OrdinalType numRowCols,
bool outputUpper, ScalarType* outputMatrix,
409 OrdinalType outputStride, OrdinalType startRowCol,
413 void copyUPLOMat(
bool inputUpper, ScalarType* inputMatrix,
414 OrdinalType inputStride, OrdinalType inputRows);
417 void checkIndex( OrdinalType rowIndex, OrdinalType colIndex = 0 )
const;
418 OrdinalType numRowCols_;
432 template<
typename OrdinalType,
typename ScalarType>
434 :
CompObject(), numRowCols_(0), stride_(0), valuesCopied_(false), values_(0), upper_(false), UPLO_(
'L')
437 template<
typename OrdinalType,
typename ScalarType>
439 :
CompObject(), numRowCols_(numRowCols_in), stride_(numRowCols_in), valuesCopied_(false), upper_(false), UPLO_(
'L')
441 values_ =
new ScalarType[stride_*numRowCols_];
442 valuesCopied_ =
true;
447 template<
typename OrdinalType,
typename ScalarType>
449 DataAccess CV,
bool upper_in, ScalarType* values_in, OrdinalType stride_in, OrdinalType numRowCols_in
451 :
CompObject(), numRowCols_(numRowCols_in), stride_(stride_in), valuesCopied_(false),
452 values_(values_in), upper_(upper_in)
461 stride_ = numRowCols_;
462 values_ =
new ScalarType[stride_*numRowCols_];
463 copyMat(upper_in, values_in, stride_in, numRowCols_, upper_, values_, stride_, 0);
464 valuesCopied_ =
true;
468 template<
typename OrdinalType,
typename ScalarType>
471 if (!Source.valuesCopied_)
473 stride_ = Source.stride_;
474 values_ = Source.values_;
475 valuesCopied_ =
false;
479 stride_ = numRowCols_;
480 const OrdinalType newsize = stride_ * numRowCols_;
482 values_ =
new ScalarType[newsize];
483 copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0);
486 numRowCols_ = 0; stride_ = 0;
487 valuesCopied_ =
false;
492 template<
typename OrdinalType,
typename ScalarType>
495 ScalarType> &Source, OrdinalType numRowCols_in, OrdinalType startRowCol )
496 :
CompObject(), numRowCols_(numRowCols_in), stride_(Source.stride_), valuesCopied_(false), upper_(Source.upper_), UPLO_(Source.UPLO_)
500 stride_ = numRowCols_in;
501 values_ =
new ScalarType[stride_ * numRowCols_in];
502 copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_in, upper_, values_, stride_, startRowCol);
503 valuesCopied_ =
true;
507 values_ = Source.values_ + (stride_ * startRowCol) + startRowCol;
511 template<
typename OrdinalType,
typename ScalarType>
521 template<
typename OrdinalType,
typename ScalarType>
525 numRowCols_ = numRowCols_in;
526 stride_ = numRowCols_;
527 values_ =
new ScalarType[stride_*numRowCols_];
529 valuesCopied_ =
true;
533 template<
typename OrdinalType,
typename ScalarType>
537 numRowCols_ = numRowCols_in;
538 stride_ = numRowCols_;
539 values_ =
new ScalarType[stride_*numRowCols_];
540 valuesCopied_ =
true;
544 template<
typename OrdinalType,
typename ScalarType>
548 ScalarType* values_tmp =
new ScalarType[numRowCols_in * numRowCols_in];
550 for(OrdinalType k = 0; k < numRowCols_in * numRowCols_in; k++)
552 values_tmp[k] = zero;
554 OrdinalType numRowCols_tmp = TEUCHOS_MIN(numRowCols_, numRowCols_in);
557 copyMat(upper_, values_, stride_, numRowCols_tmp, upper_, values_tmp, numRowCols_in, 0);
560 numRowCols_ = numRowCols_in;
561 stride_ = numRowCols_;
562 values_ = values_tmp;
563 valuesCopied_ =
true;
571 template<
typename OrdinalType,
typename ScalarType>
575 if (upper_ !=
false) {
576 copyUPLOMat(
true, values_, stride_, numRowCols_ );
582 template<
typename OrdinalType,
typename ScalarType>
586 if (upper_ ==
false) {
587 copyUPLOMat(
false, values_, stride_, numRowCols_ );
593 template<
typename OrdinalType,
typename ScalarType>
597 if (fullMatrix ==
true) {
598 for(OrdinalType j = 0; j < numRowCols_; j++)
600 for(OrdinalType i = 0; i < numRowCols_; i++)
602 values_[i + j*stride_] = value_in;
609 for(OrdinalType j = 0; j < numRowCols_; j++) {
610 for(OrdinalType i = 0; i <= j; i++) {
611 values_[i + j*stride_] = value_in;
616 for(OrdinalType j = 0; j < numRowCols_; j++) {
617 for(OrdinalType i = j; i < numRowCols_; i++) {
618 values_[i + j*stride_] = value_in;
626 template<
typename OrdinalType,
typename ScalarType>
634 for(OrdinalType j = 0; j < numRowCols_; j++) {
635 for(OrdinalType i = 0; i < j; i++) {
643 for(OrdinalType j = 0; j < numRowCols_; j++) {
644 for(OrdinalType i = j+1; i < numRowCols_; i++) {
653 for(OrdinalType i = 0; i < numRowCols_; i++) {
654 values_[i + i*stride_] = diagSum[i] + bias;
659 template<
typename OrdinalType,
typename ScalarType>
665 if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_)) {
666 upper_ = Source.upper_;
671 if (!Source.valuesCopied_) {
676 numRowCols_ = Source.numRowCols_;
677 stride_ = Source.stride_;
678 values_ = Source.values_;
679 upper_ = Source.upper_;
680 UPLO_ = Source.UPLO_;
685 numRowCols_ = Source.numRowCols_;
686 stride_ = Source.numRowCols_;
687 upper_ = Source.upper_;
688 UPLO_ = Source.UPLO_;
689 const OrdinalType newsize = stride_ * numRowCols_;
691 values_ =
new ScalarType[newsize];
692 valuesCopied_ =
true;
700 if((Source.numRowCols_ <= stride_) && (Source.numRowCols_ == numRowCols_)) {
701 numRowCols_ = Source.numRowCols_;
702 upper_ = Source.upper_;
703 UPLO_ = Source.UPLO_;
707 numRowCols_ = Source.numRowCols_;
708 stride_ = Source.numRowCols_;
709 upper_ = Source.upper_;
710 UPLO_ = Source.UPLO_;
711 const OrdinalType newsize = stride_ * numRowCols_;
713 values_ =
new ScalarType[newsize];
714 valuesCopied_ =
true;
718 copyMat(Source.upper_, Source.values_, Source.stride_, Source.numRowCols_, upper_, values_, stride_, 0);
723 template<
typename OrdinalType,
typename ScalarType>
730 template<
typename OrdinalType,
typename ScalarType>
734 if ((numRowCols_ != Source.numRowCols_))
736 TEUCHOS_CHK_REF(*
this);
742 template<
typename OrdinalType,
typename ScalarType>
746 if ((numRowCols_ != Source.numRowCols_))
748 TEUCHOS_CHK_REF(*
this);
754 template<
typename OrdinalType,
typename ScalarType>
758 if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_)) {
759 upper_ = Source.upper_;
764 if ((numRowCols_ != Source.numRowCols_))
766 TEUCHOS_CHK_REF(*
this);
768 copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0 );
776 template<
typename OrdinalType,
typename ScalarType>
779 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK 780 checkIndex( rowIndex, colIndex );
782 if ( rowIndex <= colIndex ) {
785 return(values_[colIndex * stride_ + rowIndex]);
787 return(values_[rowIndex * stride_ + colIndex]);
792 return(values_[rowIndex * stride_ + colIndex]);
794 return(values_[colIndex * stride_ + rowIndex]);
798 template<
typename OrdinalType,
typename ScalarType>
801 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK 802 checkIndex( rowIndex, colIndex );
804 if ( rowIndex <= colIndex ) {
807 return(values_[colIndex * stride_ + rowIndex]);
809 return(values_[rowIndex * stride_ + colIndex]);
814 return(values_[rowIndex * stride_ + colIndex]);
816 return(values_[colIndex * stride_ + rowIndex]);
824 template<
typename OrdinalType,
typename ScalarType>
830 template<
typename OrdinalType,
typename ScalarType>
841 for (j=0; j<numRowCols_; j++) {
843 ptr = values_ + j*stride_;
844 for (i=0; i<j; i++) {
847 ptr = values_ + j + j*stride_;
848 for (i=j; i<numRowCols_; i++) {
852 anorm = TEUCHOS_MAX( anorm, sum );
856 for (j=0; j<numRowCols_; j++) {
858 ptr = values_ + j + j*stride_;
859 for (i=j; i<numRowCols_; i++) {
863 for (i=0; i<j; i++) {
867 anorm = TEUCHOS_MAX( anorm, sum );
873 template<
typename OrdinalType,
typename ScalarType>
883 for (j = 0; j < numRowCols_; j++) {
884 for (i = 0; i < j; i++) {
891 for (j = 0; j < numRowCols_; j++) {
893 for (i = j+1; i < numRowCols_; i++) {
906 template<
typename OrdinalType,
typename ScalarType>
910 if((numRowCols_ != Operand.numRowCols_)) {
915 for(i = 0; i < numRowCols_; i++) {
916 for(j = 0; j < numRowCols_; j++) {
917 if((*
this)(i, j) != Operand(i, j)) {
926 template<
typename OrdinalType,
typename ScalarType>
929 return !((*this) == Operand);
936 template<
typename OrdinalType,
typename ScalarType>
943 for (j=0; j<numRowCols_; j++) {
944 ptr = values_ + j*stride_;
945 for (i=0; i<= j; i++) { *ptr = alpha * (*ptr); ptr++; }
949 for (j=0; j<numRowCols_; j++) {
950 ptr = values_ + j*stride_ + j;
951 for (i=j; i<numRowCols_; i++) { *ptr = alpha * (*ptr); ptr++; }
985 template<
typename OrdinalType,
typename ScalarType>
990 os <<
"Values_copied : yes" << std::endl;
992 os <<
"Values_copied : no" << std::endl;
993 os <<
"Rows / Columns : " << numRowCols_ << std::endl;
994 os <<
"LDA : " << stride_ << std::endl;
996 os <<
"Storage: Upper" << std::endl;
998 os <<
"Storage: Lower" << std::endl;
999 if(numRowCols_ == 0) {
1000 os <<
"(matrix is empty, no values to display)" << std::endl;
1002 for(OrdinalType i = 0; i < numRowCols_; i++) {
1003 for(OrdinalType j = 0; j < numRowCols_; j++){
1004 os << (*this)(i,j) <<
" ";
1015 template<
typename OrdinalType,
typename ScalarType>
1018 "SerialSymDenseMatrix<T>::checkIndex: " 1019 "Row index " << rowIndex <<
" out of range [0, "<< numRowCols_ <<
")");
1021 "SerialSymDenseMatrix<T>::checkIndex: " 1022 "Col index " << colIndex <<
" out of range [0, "<< numRowCols_ <<
")");
1025 template<
typename OrdinalType,
typename ScalarType>
1032 valuesCopied_ =
false;
1036 template<
typename OrdinalType,
typename ScalarType>
1038 bool inputUpper, ScalarType* inputMatrix,
1039 OrdinalType inputStride, OrdinalType numRowCols_in,
1040 bool outputUpper, ScalarType* outputMatrix,
1041 OrdinalType outputStride, OrdinalType startRowCol,
1046 ScalarType* ptr1 = 0;
1047 ScalarType* ptr2 = 0;
1049 for (j = 0; j < numRowCols_in; j++) {
1050 if (inputUpper ==
true) {
1052 ptr2 = inputMatrix + (j + startRowCol) * inputStride + startRowCol;
1053 if (outputUpper ==
true) {
1055 ptr1 = outputMatrix + j*outputStride;
1057 for(i = 0; i <= j; i++) {
1058 *ptr1++ += alpha*(*ptr2++);
1061 for(i = 0; i <= j; i++) {
1069 ptr1 = outputMatrix + j;
1071 for(i = 0; i <= j; i++) {
1072 *ptr1 += alpha*(*ptr2++);
1073 ptr1 += outputStride;
1076 for(i = 0; i <= j; i++) {
1078 ptr1 += outputStride;
1085 ptr2 = inputMatrix + (startRowCol+j) * inputStride + startRowCol + j;
1086 if (outputUpper ==
true) {
1089 ptr1 = outputMatrix + j*outputStride + j;
1091 for(i = j; i < numRowCols_in; i++) {
1092 *ptr1 += alpha*(*ptr2++);
1093 ptr1 += outputStride;
1096 for(i = j; i < numRowCols_in; i++) {
1098 ptr1 += outputStride;
1104 ptr1 = outputMatrix + j*outputStride + j;
1106 for(i = j; i < numRowCols_in; i++) {
1107 *ptr1++ += alpha*(*ptr2++);
1110 for(i = j; i < numRowCols_in; i++) {
1119 template<
typename OrdinalType,
typename ScalarType>
1121 bool inputUpper, ScalarType* inputMatrix,
1122 OrdinalType inputStride, OrdinalType inputRows
1126 ScalarType * ptr1 = 0;
1127 ScalarType * ptr2 = 0;
1130 for (j=1; j<inputRows; j++) {
1131 ptr1 = inputMatrix + j;
1132 ptr2 = inputMatrix + j*inputStride;
1133 for (i=0; i<j; i++) {
1140 for (i=1; i<inputRows; i++) {
1141 ptr1 = inputMatrix + i;
1142 ptr2 = inputMatrix + i*inputStride;
1143 for (j=0; j<i; j++) {
bool upper() const
Returns true if upper triangular part of this matrix has and will be used.
T magnitudeType
Mandatory typedef for result of magnitude.
SerialSymDenseMatrix()
Default constructor; defines a zero size object.
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.
ScalarType * values() const
Returns the pointer to the ScalarType data array contained in the object.
ScalarType & operator()(OrdinalType rowIndex, OrdinalType colIndex)
Element access method (non-const).
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Macro for throwing an exception with breakpointing to ease debugging.
int shape(OrdinalType numRowsCols)
Set dimensions of a Teuchos::SerialSymDenseMatrix object; init values to zero.
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...
OrdinalType numRows() const
Returns the row dimension of this matrix.
ScalarTraits< ScalarType >::magnitudeType normInf() const
Returns the Infinity-norm of the matrix.
Object for storing data and providing functionality that is common to all computational classes...
int reshape(OrdinalType numRowsCols)
Reshape a Teuchos::SerialSymDenseMatrix object.
This structure defines some basic traits for a scalar field type.
void setLower()
Specify that the lower triangle of the this matrix should be used.
void setUpper()
Specify that the upper triangle of the this matrix should be used.
Functionality and data that is common to all computational classes.
SerialSymDenseMatrix< OrdinalType, ScalarType > & operator*=(const ScalarType alpha)
Inplace scalar-matrix product A = alpha*A.
bool operator==(const SerialSymDenseMatrix< OrdinalType, ScalarType > &Operand) const
Equality of two matrices.
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero(), bool fullMatrix=false)
Set all values in the matrix to a constant value.
SerialSymDenseMatrix< OrdinalType, ScalarType > & operator+=(const SerialSymDenseMatrix< OrdinalType, ScalarType > &Source)
Add another matrix to this matrix.
ScalarTraits< ScalarType >::magnitudeType normOne() const
Returns the 1-norm of the matrix.
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...
char UPLO() const
Returns character value of UPLO used by LAPACK routines.
virtual ~SerialSymDenseMatrix()
Teuchos::SerialSymDenseMatrix destructor.
static magnitudeType magnitude(T a)
Returns the magnitudeType of the scalar type a.
SerialSymDenseMatrix< OrdinalType, ScalarType > & operator-=(const SerialSymDenseMatrix< OrdinalType, ScalarType > &Source)
Subtract another matrix from this matrix.
The Teuchos namespace contains all of the classes, structs and enums used by Teuchos, as well as a number of utility routines.
OrdinalType ordinalType
Typedef for ordinal type.
SerialSymDenseMatrix< OrdinalType, ScalarType > & operator=(const SerialSymDenseMatrix< OrdinalType, ScalarType > &Source)
Copies values from one matrix to another.
int shapeUninitialized(OrdinalType numRowsCols)
Set dimensions of a Teuchos::SerialSymDenseMatrix object; don't initialize values.
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.
OrdinalType numCols() const
Returns the column dimension of this matrix.
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
Returns the Frobenius-norm of the matrix.
SerialSymDenseMatrix< OrdinalType, ScalarType > & assign(const SerialSymDenseMatrix< OrdinalType, ScalarType > &Source)
Copies values from one matrix to another.
bool operator!=(const SerialSymDenseMatrix< OrdinalType, ScalarType > &Operand) const
Inequality of two matrices.
Teuchos::DataAccess Mode enumerable type.
bool empty() const
Returns whether this matrix is empty.
ScalarType scalarType
Typedef for scalar type.