42 #ifndef _TEUCHOS_SERIALBANDDENSEMATRIX_HPP_ 43 #define _TEUCHOS_SERIALBANDDENSEMATRIX_HPP_ 54 #include "Teuchos_Assert.hpp" 132 template<
typename OrdinalType,
typename ScalarType>
293 ScalarType&
operator () (OrdinalType rowIndex, OrdinalType colIndex);
303 const ScalarType&
operator () (OrdinalType rowIndex, OrdinalType colIndex)
const;
323 const ScalarType*
operator [] (OrdinalType colIndex)
const;
327 ScalarType*
values()
const {
return(values_); }
357 int scale (
const ScalarType alpha );
391 OrdinalType
numRows()
const {
return(numRows_); }
394 OrdinalType
numCols()
const {
return(numCols_); }
403 OrdinalType
stride()
const {
return(stride_); }
406 bool empty()
const {
return(numRows_ == 0 || numCols_ == 0); }
424 virtual void print(std::ostream& os)
const;
429 void copyMat(ScalarType* inputMatrix, OrdinalType strideInput,
430 OrdinalType
numRows, OrdinalType
numCols, ScalarType* outputMatrix,
433 void checkIndex( OrdinalType rowIndex, OrdinalType colIndex = 0 )
const;
434 OrdinalType numRows_;
435 OrdinalType numCols_;
448 template<
typename OrdinalType,
typename ScalarType>
456 valuesCopied_ (false),
460 template<
typename OrdinalType,
typename ScalarType>
463 OrdinalType numCols_in,
468 numRows_ (numRows_in),
469 numCols_ (numCols_in),
470 stride_ (kl_in+ku_in+1),
473 valuesCopied_ (true),
476 values_ =
new ScalarType[stride_ * numCols_];
482 template<
typename OrdinalType,
typename ScalarType>
485 ScalarType* values_in,
486 OrdinalType stride_in,
487 OrdinalType numRows_in,
488 OrdinalType numCols_in,
492 numRows_ (numRows_in),
493 numCols_ (numCols_in),
497 valuesCopied_ (false),
502 values_ =
new ScalarType[stride_*numCols_];
503 copyMat (values_in, stride_in, numRows_, numCols_, values_, stride_, 0);
504 valuesCopied_ =
true;
508 template<
typename OrdinalType,
typename ScalarType>
517 valuesCopied_ (true),
521 numRows_ = Source.numRows_;
522 numCols_ = Source.numCols_;
526 values_ =
new ScalarType[stride_*numCols_];
527 copyMat (Source.values_, Source.stride_, numRows_, numCols_,
528 values_, stride_, 0);
531 numRows_ = Source.numCols_;
532 numCols_ = Source.numRows_;
536 values_ =
new ScalarType[stride_*numCols_];
537 for (OrdinalType j = 0; j < numCols_; ++j) {
538 for (OrdinalType i = TEUCHOS_MAX(0,j-ku_);
539 i <= TEUCHOS_MIN(numRows_-1,j+kl_); ++i) {
540 values_[j*stride_ + (ku_+i-j)] =
546 numRows_ = Source.numCols_;
547 numCols_ = Source.numRows_;
551 values_ =
new ScalarType[stride_*numCols_];
552 for (OrdinalType j=0; j<numCols_; j++) {
553 for (OrdinalType i = TEUCHOS_MAX(0,j-ku_);
554 i <= TEUCHOS_MIN(numRows_-1,j+kl_); ++i) {
555 values_[j*stride_ + (ku_+i-j)] = Source.values_[i*Source.stride_ + (Source.ku_+j-i)];
561 template<
typename OrdinalType,
typename ScalarType>
565 OrdinalType numRows_in,
566 OrdinalType numCols_in,
567 OrdinalType startCol)
569 numRows_ (numRows_in),
570 numCols_ (numCols_in),
571 stride_ (Source.stride_),
574 valuesCopied_ (false),
575 values_ (Source.values_)
578 values_ =
new ScalarType[stride_ * numCols_in];
579 copyMat (Source.values_, Source.stride_, numRows_in, numCols_in,
580 values_, stride_, startCol);
581 valuesCopied_ =
true;
583 values_ = values_ + (stride_ * startCol);
587 template<
typename OrdinalType,
typename ScalarType>
597 template<
typename OrdinalType,
typename ScalarType>
599 OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType kl_in, OrdinalType ku_in
603 numRows_ = numRows_in;
604 numCols_ = numCols_in;
608 values_ =
new ScalarType[stride_*numCols_];
610 valuesCopied_ =
true;
615 template<
typename OrdinalType,
typename ScalarType>
617 OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType kl_in, OrdinalType ku_in
621 numRows_ = numRows_in;
622 numCols_ = numCols_in;
626 values_ =
new ScalarType[stride_*numCols_];
627 valuesCopied_ =
true;
632 template<
typename OrdinalType,
typename ScalarType>
634 OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType kl_in, OrdinalType ku_in
639 ScalarType* values_tmp =
new ScalarType[(kl_in+ku_in+1) * numCols_in];
641 for(OrdinalType k = 0; k < (kl_in+ku_in+1) * numCols_in; k++) {
642 values_tmp[k] = zero;
644 OrdinalType numRows_tmp = TEUCHOS_MIN(numRows_, numRows_in);
645 OrdinalType numCols_tmp = TEUCHOS_MIN(numCols_, numCols_in);
647 copyMat(values_, stride_, numRows_tmp, numCols_tmp, values_tmp,
651 numRows_ = numRows_in;
652 numCols_ = numCols_in;
656 values_ = values_tmp;
657 valuesCopied_ =
true;
666 template<
typename OrdinalType,
typename ScalarType>
671 for(OrdinalType j = 0; j < numCols_; j++) {
672 for (OrdinalType i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
673 values_[(ku_+i-j) + j*stride_] = value_in;
680 template<
typename OrdinalType,
typename ScalarType>
685 for(OrdinalType j = 0; j < numCols_; j++) {
686 for (OrdinalType i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
694 template<
typename OrdinalType,
typename ScalarType>
703 if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
707 if (!Source.valuesCopied_) {
712 numRows_ = Source.numRows_;
713 numCols_ = Source.numCols_;
716 stride_ = Source.stride_;
717 values_ = Source.values_;
721 numRows_ = Source.numRows_;
722 numCols_ = Source.numCols_;
726 const OrdinalType newsize = stride_ * numCols_;
728 values_ =
new ScalarType[newsize];
729 valuesCopied_ =
true;
735 if((Source.numRows_ <= stride_) && (Source.numCols_ == numCols_)) {
736 numRows_ = Source.numRows_;
737 numCols_ = Source.numCols_;
743 numRows_ = Source.numRows_;
744 numCols_ = Source.numCols_;
748 const OrdinalType newsize = stride_ * numCols_;
750 values_ =
new ScalarType[newsize];
751 valuesCopied_ =
true;
755 copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0);
761 template<
typename OrdinalType,
typename ScalarType>
766 if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_) || (kl_ != Source.kl_) || (ku_ != Source.ku_)) {
767 TEUCHOS_CHK_REF(*
this);
774 template<
typename OrdinalType,
typename ScalarType>
779 if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_) || (kl_ != Source.kl_) || (ku_ != Source.ku_)) {
780 TEUCHOS_CHK_REF(*
this);
787 template<
typename OrdinalType,
typename ScalarType>
792 if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
796 if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_) || (kl_ != Source.kl_) || (ku_ != Source.ku_)) {
797 TEUCHOS_CHK_REF(*
this);
799 copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0);
808 template<
typename OrdinalType,
typename ScalarType>
811 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK 812 checkIndex( rowIndex, colIndex );
814 return(values_[colIndex * stride_ + ku_+rowIndex-colIndex]);
817 template<
typename OrdinalType,
typename ScalarType>
820 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK 821 checkIndex( rowIndex, colIndex );
823 return(values_[colIndex * stride_ + ku_+rowIndex-colIndex]);
826 template<
typename OrdinalType,
typename ScalarType>
829 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK 830 checkIndex( 0, colIndex );
832 return(values_ + colIndex * stride_);
835 template<
typename OrdinalType,
typename ScalarType>
838 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK 839 checkIndex( 0, colIndex );
841 return(values_ + colIndex * stride_);
848 template<
typename OrdinalType,
typename ScalarType>
856 for(j = 0; j < numCols_; j++) {
858 ptr = values_ + j * stride_ + TEUCHOS_MAX(0, ku_-j);
859 for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
872 template<
typename OrdinalType,
typename ScalarType>
878 for (i = 0; i < numRows_; i++) {
880 for (j=TEUCHOS_MAX(0,i-kl_); j<=TEUCHOS_MIN(numCols_-1,i+ku_); j++) {
883 anorm = TEUCHOS_MAX( anorm, sum );
890 template<
typename OrdinalType,
typename ScalarType>
896 for (j = 0; j < numCols_; j++) {
897 for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
911 template<
typename OrdinalType,
typename ScalarType>
916 if((numRows_ != Operand.numRows_) || (numCols_ != Operand.numCols_) || (kl_ != Operand.kl_) || (ku_ != Operand.ku_)) {
920 for(j = 0; j < numCols_; j++) {
921 for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
922 if((*
this)(i, j) != Operand(i, j)) {
932 template<
typename OrdinalType,
typename ScalarType>
935 return !((*this) == Operand);
942 template<
typename OrdinalType,
typename ScalarType>
945 this->
scale( alpha );
949 template<
typename OrdinalType,
typename ScalarType>
956 for (j=0; j<numCols_; j++) {
957 ptr = values_ + j*stride_ + TEUCHOS_MAX(0, ku_-j);
958 for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
959 *ptr = alpha * (*ptr); ptr++;
967 template<
typename OrdinalType,
typename ScalarType>
975 if ((numRows_ != A.numRows_) || (numCols_ != A.numCols_) || (kl_ != A.kl_) || (ku_ != A.ku_)) {
978 for (j=0; j<numCols_; j++) {
979 ptr = values_ + j*stride_ + TEUCHOS_MAX(0, ku_-j);
980 for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
981 *ptr = A(i,j) * (*ptr); ptr++;
989 template<
typename OrdinalType,
typename ScalarType>
994 os <<
"Values_copied : yes" << std::endl;
996 os <<
"Values_copied : no" << std::endl;
997 os <<
"Rows : " << numRows_ << std::endl;
998 os <<
"Columns : " << numCols_ << std::endl;
999 os <<
"Lower Bandwidth : " << kl_ << std::endl;
1000 os <<
"Upper Bandwidth : " << ku_ << std::endl;
1001 os <<
"LDA : " << stride_ << std::endl;
1002 if(numRows_ == 0 || numCols_ == 0) {
1003 os <<
"(matrix is empty, no values to display)" << std::endl;
1006 for(OrdinalType i = 0; i < numRows_; i++) {
1007 for (OrdinalType j=TEUCHOS_MAX(0,i-kl_); j<=TEUCHOS_MIN(numCols_-1,i+ku_); j++) {
1008 os << (*this)(i,j) <<
" ";
1019 template<
typename OrdinalType,
typename ScalarType>
1023 rowIndex < TEUCHOS_MAX(0,colIndex-ku_) || rowIndex > TEUCHOS_MIN(numRows_-1,colIndex+kl_),
1025 "SerialBandDenseMatrix<T>::checkIndex: " 1026 "Row index " << rowIndex <<
" out of range [0, "<< numRows_ <<
")");
1028 "SerialBandDenseMatrix<T>::checkIndex: " 1029 "Col index " << colIndex <<
" out of range [0, "<< numCols_ <<
")");
1033 template<
typename OrdinalType,
typename ScalarType>
1036 if (valuesCopied_) {
1039 valuesCopied_ =
false;
1043 template<
typename OrdinalType,
typename ScalarType>
1045 ScalarType* inputMatrix, OrdinalType strideInput, OrdinalType numRows_in,
1046 OrdinalType numCols_in, ScalarType* outputMatrix, OrdinalType strideOutput, OrdinalType startCol, ScalarType alpha
1050 ScalarType* ptr1 = 0;
1051 ScalarType* ptr2 = 0;
1053 for(j = 0; j < numCols_in; j++) {
1054 ptr1 = outputMatrix + (j * strideOutput) + TEUCHOS_MAX(0, ku_-j);
1055 ptr2 = inputMatrix + (j + startCol) * strideInput + TEUCHOS_MAX(0, ku_-j);
1057 for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_in-1,j+kl_); i++) {
1058 *ptr1++ += alpha*(*ptr2++);
1061 for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_in-1,j+kl_); i++) {
bool empty() const
Returns whether this matrix is empty.
OrdinalType numRows() const
Returns the row dimension of this matrix.
SerialBandDenseMatrix< OrdinalType, ScalarType > & operator+=(const SerialBandDenseMatrix< OrdinalType, ScalarType > &Source)
Add another matrix to this matrix.
ScalarType & operator()(OrdinalType rowIndex, OrdinalType colIndex)
Element access method (non-const).
int scale(const ScalarType alpha)
Scale this matrix by alpha; *this = alpha**this.
int shape(OrdinalType numRows, OrdinalType numCols, OrdinalType kl, OrdinalType ku)
Shape method for changing the size of a SerialBandDenseMatrix, initializing entries to zero...
ScalarTraits< ScalarType >::magnitudeType normInf() const
Returns the Infinity-norm of the matrix.
Templated interface class to BLAS routines.
virtual ~SerialBandDenseMatrix()
Destructor.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Macro for throwing an exception with breakpointing to ease debugging.
int reshape(OrdinalType numRows, OrdinalType numCols, OrdinalType kl, OrdinalType ku)
Reshaping method for changing the size of a SerialBandDenseMatrix, keeping the entries.
Teuchos header file which uses auto-configuration information to include necessary C++ headers...
ScalarTraits< ScalarType >::magnitudeType normOne() const
Returns the 1-norm of the matrix.
SerialBandDenseMatrix< OrdinalType, ScalarType > & assign(const SerialBandDenseMatrix< OrdinalType, ScalarType > &Source)
Copies values from one matrix to another.
Object for storing data and providing functionality that is common to all computational classes...
This structure defines some basic traits for a scalar field type.
OrdinalType stride() const
Returns the stride between the columns of this matrix in memory.
virtual void print(std::ostream &os) const
Print method. Defines the behavior of the std::ostream << operator inherited from the Object class...
ScalarType * values() const
Data array access method.
Functionality and data that is common to all computational classes.
bool operator!=(const SerialBandDenseMatrix< OrdinalType, ScalarType > &Operand) const
Inequality of two matrices.
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
Returns the Frobenius-norm of the matrix.
This class creates and provides basic support for banded dense matrices of templated type...
SerialBandDenseMatrix< OrdinalType, ScalarType > & operator*=(const ScalarType alpha)
Scale this matrix by alpha; *this = alpha**this.
SerialBandDenseMatrix< OrdinalType, ScalarType > & operator-=(const SerialBandDenseMatrix< OrdinalType, ScalarType > &Source)
Subtract another matrix from this matrix.
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero())
Set all values in the matrix to a constant value.
int random()
Set all values in the matrix to be random numbers.
SerialBandDenseMatrix()
Default Constructor.
bool operator==(const SerialBandDenseMatrix< OrdinalType, ScalarType > &Operand) const
Equality of two matrices.
static magnitudeType magnitude(T a)
Returns the magnitudeType of the scalar type a.
OrdinalType ordinalType
Typedef for ordinal type.
SerialBandDenseMatrix< OrdinalType, ScalarType > & operator=(const SerialBandDenseMatrix< OrdinalType, ScalarType > &Source)
Copies values from one matrix to another.
OrdinalType upperBandwidth() const
Returns the upper bandwidth of this matrix.
The Teuchos namespace contains all of the classes, structs and enums used by Teuchos, as well as a number of utility routines.
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.
OrdinalType lowerBandwidth() const
Returns the lower bandwidth of this matrix.
ScalarType * operator[](OrdinalType colIndex)
Column access method (non-const).
ScalarType scalarType
Typedef for scalar type.
OrdinalType numCols() const
Returns the column dimension of this matrix.
Reference-counted pointer class and non-member templated function implementations.
int shapeUninitialized(OrdinalType numRows, OrdinalType numCols, OrdinalType kl, OrdinalType ku)
Same as shape() except leaves uninitialized.
Teuchos::DataAccess Mode enumerable type.