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>
452 BLAS<OrdinalType,ScalarType>(),
458 valuesCopied_ (false),
462 template<
typename OrdinalType,
typename ScalarType>
465 OrdinalType numCols_in,
471 BLAS<OrdinalType,ScalarType>(),
472 numRows_ (numRows_in),
473 numCols_ (numCols_in),
474 stride_ (kl_in+ku_in+1),
477 valuesCopied_ (true),
480 values_ =
new ScalarType[stride_ * numCols_];
486 template<
typename OrdinalType,
typename ScalarType>
489 ScalarType* values_in,
490 OrdinalType stride_in,
491 OrdinalType numRows_in,
492 OrdinalType numCols_in,
497 BLAS<OrdinalType,ScalarType>(),
498 numRows_ (numRows_in),
499 numCols_ (numCols_in),
503 valuesCopied_ (false),
508 values_ =
new ScalarType[stride_*numCols_];
509 copyMat (values_in, stride_in, numRows_, numCols_, values_, stride_, 0);
510 valuesCopied_ =
true;
514 template<
typename OrdinalType,
typename ScalarType>
519 BLAS<OrdinalType,ScalarType>(),
525 valuesCopied_ (true),
529 numRows_ = Source.numRows_;
530 numCols_ = Source.numCols_;
534 values_ =
new ScalarType[stride_*numCols_];
535 copyMat (Source.values_, Source.stride_, numRows_, numCols_,
536 values_, stride_, 0);
539 numRows_ = Source.numCols_;
540 numCols_ = Source.numRows_;
544 values_ =
new ScalarType[stride_*numCols_];
545 for (OrdinalType j = 0; j < numCols_; ++j) {
546 for (OrdinalType i = TEUCHOS_MAX(0,j-ku_);
547 i <= TEUCHOS_MIN(numRows_-1,j+kl_); ++i) {
548 values_[j*stride_ + (ku_+i-j)] =
554 numRows_ = Source.numCols_;
555 numCols_ = Source.numRows_;
559 values_ =
new ScalarType[stride_*numCols_];
560 for (OrdinalType j=0; j<numCols_; j++) {
561 for (OrdinalType i = TEUCHOS_MAX(0,j-ku_);
562 i <= TEUCHOS_MIN(numRows_-1,j+kl_); ++i) {
563 values_[j*stride_ + (ku_+i-j)] = Source.values_[i*Source.stride_ + (Source.ku_+j-i)];
569 template<
typename OrdinalType,
typename ScalarType>
573 OrdinalType numRows_in,
574 OrdinalType numCols_in,
575 OrdinalType startCol)
578 BLAS<OrdinalType,ScalarType>(),
579 numRows_ (numRows_in),
580 numCols_ (numCols_in),
581 stride_ (Source.stride_),
584 valuesCopied_ (false),
585 values_ (Source.values_)
588 values_ =
new ScalarType[stride_ * numCols_in];
589 copyMat (Source.values_, Source.stride_, numRows_in, numCols_in,
590 values_, stride_, startCol);
591 valuesCopied_ =
true;
593 values_ = values_ + (stride_ * startCol);
597 template<
typename OrdinalType,
typename ScalarType>
607 template<
typename OrdinalType,
typename ScalarType>
609 OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType kl_in, OrdinalType ku_in
613 numRows_ = numRows_in;
614 numCols_ = numCols_in;
618 values_ =
new ScalarType[stride_*numCols_];
620 valuesCopied_ =
true;
625 template<
typename OrdinalType,
typename ScalarType>
627 OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType kl_in, OrdinalType ku_in
631 numRows_ = numRows_in;
632 numCols_ = numCols_in;
636 values_ =
new ScalarType[stride_*numCols_];
637 valuesCopied_ =
true;
642 template<
typename OrdinalType,
typename ScalarType>
644 OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType kl_in, OrdinalType ku_in
649 ScalarType* values_tmp =
new ScalarType[(kl_in+ku_in+1) * numCols_in];
651 for(OrdinalType k = 0; k < (kl_in+ku_in+1) * numCols_in; k++) {
652 values_tmp[k] = zero;
654 OrdinalType numRows_tmp = TEUCHOS_MIN(numRows_, numRows_in);
655 OrdinalType numCols_tmp = TEUCHOS_MIN(numCols_, numCols_in);
657 copyMat(values_, stride_, numRows_tmp, numCols_tmp, values_tmp,
661 numRows_ = numRows_in;
662 numCols_ = numCols_in;
666 values_ = values_tmp;
667 valuesCopied_ =
true;
676 template<
typename OrdinalType,
typename ScalarType>
681 for(OrdinalType j = 0; j < numCols_; j++) {
682 for (OrdinalType i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
683 values_[(ku_+i-j) + j*stride_] = value_in;
690 template<
typename OrdinalType,
typename ScalarType>
695 for(OrdinalType j = 0; j < numCols_; j++) {
696 for (OrdinalType i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
704 template<
typename OrdinalType,
typename ScalarType>
713 if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
717 if (!Source.valuesCopied_) {
722 numRows_ = Source.numRows_;
723 numCols_ = Source.numCols_;
726 stride_ = Source.stride_;
727 values_ = Source.values_;
731 numRows_ = Source.numRows_;
732 numCols_ = Source.numCols_;
736 const OrdinalType newsize = stride_ * numCols_;
738 values_ =
new ScalarType[newsize];
739 valuesCopied_ =
true;
745 if((Source.numRows_ <= stride_) && (Source.numCols_ == numCols_)) {
746 numRows_ = Source.numRows_;
747 numCols_ = Source.numCols_;
753 numRows_ = Source.numRows_;
754 numCols_ = Source.numCols_;
758 const OrdinalType newsize = stride_ * numCols_;
760 values_ =
new ScalarType[newsize];
761 valuesCopied_ =
true;
765 copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0);
771 template<
typename OrdinalType,
typename ScalarType>
776 if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_) || (kl_ != Source.kl_) || (ku_ != Source.ku_)) {
777 TEUCHOS_CHK_REF(*
this);
784 template<
typename OrdinalType,
typename ScalarType>
789 if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_) || (kl_ != Source.kl_) || (ku_ != Source.ku_)) {
790 TEUCHOS_CHK_REF(*
this);
797 template<
typename OrdinalType,
typename ScalarType>
802 if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
806 if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_) || (kl_ != Source.kl_) || (ku_ != Source.ku_)) {
807 TEUCHOS_CHK_REF(*
this);
809 copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0);
818 template<
typename OrdinalType,
typename ScalarType>
821 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK 822 checkIndex( rowIndex, colIndex );
824 return(values_[colIndex * stride_ + ku_+rowIndex-colIndex]);
827 template<
typename OrdinalType,
typename ScalarType>
830 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK 831 checkIndex( rowIndex, colIndex );
833 return(values_[colIndex * stride_ + ku_+rowIndex-colIndex]);
836 template<
typename OrdinalType,
typename ScalarType>
839 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK 840 checkIndex( 0, colIndex );
842 return(values_ + colIndex * stride_);
845 template<
typename OrdinalType,
typename ScalarType>
848 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK 849 checkIndex( 0, colIndex );
851 return(values_ + colIndex * stride_);
858 template<
typename OrdinalType,
typename ScalarType>
866 for(j = 0; j < numCols_; j++) {
868 ptr = values_ + j * stride_ + TEUCHOS_MAX(0, ku_-j);
869 for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
877 updateFlops((kl_+ku_+1) * numCols_);
882 template<
typename OrdinalType,
typename ScalarType>
888 for (i = 0; i < numRows_; i++) {
890 for (j=TEUCHOS_MAX(0,i-kl_); j<=TEUCHOS_MIN(numCols_-1,i+ku_); j++) {
893 anorm = TEUCHOS_MAX( anorm, sum );
895 updateFlops((kl_+ku_+1) * numCols_);
900 template<
typename OrdinalType,
typename ScalarType>
906 for (j = 0; j < numCols_; j++) {
907 for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
912 updateFlops((kl_+ku_+1) * numCols_);
921 template<
typename OrdinalType,
typename ScalarType>
926 if((numRows_ != Operand.numRows_) || (numCols_ != Operand.numCols_) || (kl_ != Operand.kl_) || (ku_ != Operand.ku_)) {
930 for(j = 0; j < numCols_; j++) {
931 for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
932 if((*
this)(i, j) != Operand(i, j)) {
942 template<
typename OrdinalType,
typename ScalarType>
945 return !((*this) == Operand);
952 template<
typename OrdinalType,
typename ScalarType>
955 this->scale( alpha );
959 template<
typename OrdinalType,
typename ScalarType>
966 for (j=0; j<numCols_; j++) {
967 ptr = values_ + j*stride_ + TEUCHOS_MAX(0, ku_-j);
968 for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
969 *ptr = alpha * (*ptr); ptr++;
972 updateFlops( (kl_+ku_+1)*numCols_ );
977 template<
typename OrdinalType,
typename ScalarType>
985 if ((numRows_ != A.numRows_) || (numCols_ != A.numCols_) || (kl_ != A.kl_) || (ku_ != A.ku_)) {
988 for (j=0; j<numCols_; j++) {
989 ptr = values_ + j*stride_ + TEUCHOS_MAX(0, ku_-j);
990 for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
991 *ptr = A(i,j) * (*ptr); ptr++;
994 updateFlops( (kl_+ku_+1)*numCols_ );
999 template<
typename OrdinalType,
typename ScalarType>
1004 os <<
"Values_copied : yes" << std::endl;
1006 os <<
"Values_copied : no" << std::endl;
1007 os <<
"Rows : " << numRows_ << std::endl;
1008 os <<
"Columns : " << numCols_ << std::endl;
1009 os <<
"Lower Bandwidth : " << kl_ << std::endl;
1010 os <<
"Upper Bandwidth : " << ku_ << std::endl;
1011 os <<
"LDA : " << stride_ << std::endl;
1012 if(numRows_ == 0 || numCols_ == 0) {
1013 os <<
"(matrix is empty, no values to display)" << std::endl;
1016 for(OrdinalType i = 0; i < numRows_; i++) {
1017 for (OrdinalType j=TEUCHOS_MAX(0,i-kl_); j<=TEUCHOS_MIN(numCols_-1,i+ku_); j++) {
1018 os << (*this)(i,j) <<
" ";
1029 template<
typename OrdinalType,
typename ScalarType>
1033 rowIndex < TEUCHOS_MAX(0,colIndex-ku_) || rowIndex > TEUCHOS_MIN(numRows_-1,colIndex+kl_),
1035 "SerialBandDenseMatrix<T>::checkIndex: " 1036 "Row index " << rowIndex <<
" out of range [0, "<< numRows_ <<
")");
1038 "SerialBandDenseMatrix<T>::checkIndex: " 1039 "Col index " << colIndex <<
" out of range [0, "<< numCols_ <<
")");
1043 template<
typename OrdinalType,
typename ScalarType>
1044 void SerialBandDenseMatrix<OrdinalType, ScalarType>::deleteArrays(
void)
1046 if (valuesCopied_) {
1049 valuesCopied_ =
false;
1053 template<
typename OrdinalType,
typename ScalarType>
1054 void SerialBandDenseMatrix<OrdinalType, ScalarType>::copyMat(
1055 ScalarType* inputMatrix, OrdinalType strideInput, OrdinalType numRows_in,
1056 OrdinalType numCols_in, ScalarType* outputMatrix, OrdinalType strideOutput, OrdinalType startCol, ScalarType alpha
1060 ScalarType* ptr1 = 0;
1061 ScalarType* ptr2 = 0;
1063 for(j = 0; j < numCols_in; j++) {
1064 ptr1 = outputMatrix + (j * strideOutput) + TEUCHOS_MAX(0, ku_-j);
1065 ptr2 = inputMatrix + (j + startCol) * strideInput + TEUCHOS_MAX(0, ku_-j);
1067 for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_in-1,j+kl_); i++) {
1068 *ptr1++ += alpha*(*ptr2++);
1071 for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_in-1,j+kl_); i++) {
ScalarTraits< ScalarType >::magnitudeType normInf() const
Returns the Infinity-norm of the matrix.
bool operator==(const SerialBandDenseMatrix< OrdinalType, ScalarType > &Operand) const
Equality of two matrices.
SerialBandDenseMatrix< OrdinalType, ScalarType > & operator+=(const SerialBandDenseMatrix< OrdinalType, ScalarType > &Source)
Add another matrix to this matrix.
bool operator!=(const SerialBandDenseMatrix< OrdinalType, ScalarType > &Operand) const
Inequality of two matrices.
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...
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.
OrdinalType upperBandwidth() const
Returns the upper bandwidth of this matrix.
Teuchos header file which uses auto-configuration information to include necessary C++ headers...
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.
Functionality and data that is common to all computational classes.
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
Returns the Frobenius-norm of the matrix.
ScalarTraits< ScalarType >::magnitudeType normOne() const
Returns the 1-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.
ScalarType * values() const
Data array access method.
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 empty() const
Returns whether this matrix is empty.
OrdinalType numRows() const
Returns the row dimension of this matrix.
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.
The Teuchos namespace contains all of the classes, structs and enums used by Teuchos, as well as a number of utility routines.
Defines basic traits for the scalar field type.
static T zero()
Returns representation of zero for this scalar type.
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.
OrdinalType stride() const
Returns the stride between the columns of this matrix in memory.
Reference-counted pointer class and non-member templated function implementations.
OrdinalType lowerBandwidth() const
Returns the lower bandwidth of this matrix.
virtual void print(std::ostream &os) const
Print method. Defines the behavior of the std::ostream << operator inherited from the Object class...
int shapeUninitialized(OrdinalType numRows, OrdinalType numCols, OrdinalType kl, OrdinalType ku)
Same as shape() except leaves uninitialized.
Teuchos::DataAccess Mode enumerable type.