42#ifndef _TEUCHOS_SERIALBANDDENSEMATRIX_HPP_
43#define _TEUCHOS_SERIALBANDDENSEMATRIX_HPP_
54#include "Teuchos_Assert.hpp"
132template<
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); }
426 virtual std::ostream& print(std::ostream& os)
const;
431 void copyMat(ScalarType* inputMatrix, OrdinalType strideInput,
432 OrdinalType
numRows, OrdinalType
numCols, ScalarType* outputMatrix,
435 void checkIndex( OrdinalType rowIndex, OrdinalType colIndex = 0 )
const;
436 OrdinalType numRows_;
437 OrdinalType numCols_;
450template<
typename OrdinalType,
typename ScalarType>
453 BLAS<OrdinalType,ScalarType>(),
459 valuesCopied_ (false),
463template<
typename OrdinalType,
typename ScalarType>
466 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_];
486template<
typename OrdinalType,
typename ScalarType>
489 ScalarType* values_in,
490 OrdinalType stride_in,
491 OrdinalType numRows_in,
492 OrdinalType numCols_in,
496 BLAS<OrdinalType,ScalarType>(),
497 numRows_ (numRows_in),
498 numCols_ (numCols_in),
502 valuesCopied_ (false),
507 values_ =
new ScalarType[stride_*numCols_];
508 copyMat (values_in, stride_in, numRows_, numCols_, values_, stride_, 0);
509 valuesCopied_ =
true;
513template<
typename OrdinalType,
typename ScalarType>
517 BLAS<OrdinalType,ScalarType>(),
523 valuesCopied_ (true),
527 numRows_ = Source.numRows_;
528 numCols_ = Source.numCols_;
532 values_ =
new ScalarType[stride_*numCols_];
533 copyMat (Source.values_, Source.stride_, numRows_, numCols_,
534 values_, stride_, 0);
537 numRows_ = Source.numCols_;
538 numCols_ = Source.numRows_;
542 values_ =
new ScalarType[stride_*numCols_];
543 for (OrdinalType j = 0; j < numCols_; ++j) {
544 for (OrdinalType i = TEUCHOS_MAX(0,j-ku_);
545 i <= TEUCHOS_MIN(numRows_-1,j+kl_); ++i) {
546 values_[j*stride_ + (ku_+i-j)] =
552 numRows_ = Source.numCols_;
553 numCols_ = Source.numRows_;
557 values_ =
new ScalarType[stride_*numCols_];
558 for (OrdinalType j=0; j<numCols_; j++) {
559 for (OrdinalType i = TEUCHOS_MAX(0,j-ku_);
560 i <= TEUCHOS_MIN(numRows_-1,j+kl_); ++i) {
561 values_[j*stride_ + (ku_+i-j)] = Source.values_[i*Source.stride_ + (Source.ku_+j-i)];
567template<
typename OrdinalType,
typename ScalarType>
571 OrdinalType numRows_in,
572 OrdinalType numCols_in,
573 OrdinalType startCol)
575 BLAS<OrdinalType,ScalarType>(),
576 numRows_ (numRows_in),
577 numCols_ (numCols_in),
578 stride_ (Source.stride_),
581 valuesCopied_ (false),
582 values_ (Source.values_)
585 values_ =
new ScalarType[stride_ * numCols_in];
586 copyMat (Source.values_, Source.stride_, numRows_in, numCols_in,
587 values_, stride_, startCol);
588 valuesCopied_ =
true;
590 values_ = values_ + (stride_ * startCol);
594template<
typename OrdinalType,
typename ScalarType>
604template<
typename OrdinalType,
typename ScalarType>
606 OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType kl_in, OrdinalType ku_in
610 numRows_ = numRows_in;
611 numCols_ = numCols_in;
615 values_ =
new ScalarType[stride_*numCols_];
617 valuesCopied_ =
true;
622template<
typename OrdinalType,
typename ScalarType>
624 OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType kl_in, OrdinalType ku_in
628 numRows_ = numRows_in;
629 numCols_ = numCols_in;
633 values_ =
new ScalarType[stride_*numCols_];
634 valuesCopied_ =
true;
639template<
typename OrdinalType,
typename ScalarType>
641 OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType kl_in, OrdinalType ku_in
646 ScalarType* values_tmp =
new ScalarType[(kl_in+ku_in+1) * numCols_in];
648 for(OrdinalType k = 0; k < (kl_in+ku_in+1) * numCols_in; k++) {
649 values_tmp[k] = zero;
651 OrdinalType numRows_tmp = TEUCHOS_MIN(numRows_, numRows_in);
652 OrdinalType numCols_tmp = TEUCHOS_MIN(numCols_, numCols_in);
654 copyMat(values_, stride_, numRows_tmp, numCols_tmp, values_tmp,
658 numRows_ = numRows_in;
659 numCols_ = numCols_in;
663 values_ = values_tmp;
664 valuesCopied_ =
true;
673template<
typename OrdinalType,
typename ScalarType>
678 for(OrdinalType j = 0; j < numCols_; j++) {
679 for (OrdinalType i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
680 values_[(ku_+i-j) + j*stride_] = value_in;
687template<
typename OrdinalType,
typename ScalarType>
692 for(OrdinalType j = 0; j < numCols_; j++) {
693 for (OrdinalType i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
701template<
typename OrdinalType,
typename ScalarType>
710 if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
714 if (!Source.valuesCopied_) {
719 numRows_ = Source.numRows_;
720 numCols_ = Source.numCols_;
723 stride_ = Source.stride_;
724 values_ = Source.values_;
728 numRows_ = Source.numRows_;
729 numCols_ = Source.numCols_;
733 const OrdinalType newsize = stride_ * numCols_;
735 values_ =
new ScalarType[newsize];
736 valuesCopied_ =
true;
742 if((Source.numRows_ <= stride_) && (Source.numCols_ == numCols_)) {
743 numRows_ = Source.numRows_;
744 numCols_ = Source.numCols_;
750 numRows_ = Source.numRows_;
751 numCols_ = Source.numCols_;
755 const OrdinalType newsize = stride_ * numCols_;
757 values_ =
new ScalarType[newsize];
758 valuesCopied_ =
true;
762 copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0);
768template<
typename OrdinalType,
typename ScalarType>
773 if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_) || (kl_ != Source.kl_) || (ku_ != Source.ku_)) {
774 TEUCHOS_CHK_REF(*
this);
781template<
typename OrdinalType,
typename ScalarType>
786 if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_) || (kl_ != Source.kl_) || (ku_ != Source.ku_)) {
787 TEUCHOS_CHK_REF(*
this);
794template<
typename OrdinalType,
typename ScalarType>
799 if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
803 if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_) || (kl_ != Source.kl_) || (ku_ != Source.ku_)) {
804 TEUCHOS_CHK_REF(*
this);
806 copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0);
815template<
typename OrdinalType,
typename ScalarType>
818#ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
819 checkIndex( rowIndex, colIndex );
821 return(values_[colIndex * stride_ + ku_+rowIndex-colIndex]);
824template<
typename OrdinalType,
typename ScalarType>
827#ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
828 checkIndex( rowIndex, colIndex );
830 return(values_[colIndex * stride_ + ku_+rowIndex-colIndex]);
833template<
typename OrdinalType,
typename ScalarType>
836#ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
837 checkIndex( 0, colIndex );
839 return(values_ + colIndex * stride_);
842template<
typename OrdinalType,
typename ScalarType>
845#ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
846 checkIndex( 0, colIndex );
848 return(values_ + colIndex * stride_);
855template<
typename OrdinalType,
typename ScalarType>
863 for(j = 0; j < numCols_; j++) {
865 ptr = values_ + j * stride_ + TEUCHOS_MAX(0, ku_-j);
866 for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
874 updateFlops((kl_+ku_+1) * numCols_);
879template<
typename OrdinalType,
typename ScalarType>
885 for (i = 0; i < numRows_; i++) {
887 for (j=TEUCHOS_MAX(0,i-kl_); j<=TEUCHOS_MIN(numCols_-1,i+ku_); j++) {
890 anorm = TEUCHOS_MAX( anorm, sum );
892 updateFlops((kl_+ku_+1) * numCols_);
897template<
typename OrdinalType,
typename ScalarType>
903 for (j = 0; j < numCols_; j++) {
904 for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
909 updateFlops((kl_+ku_+1) * numCols_);
918template<
typename OrdinalType,
typename ScalarType>
923 if((numRows_ != Operand.numRows_) || (numCols_ != Operand.numCols_) || (kl_ != Operand.kl_) || (ku_ != Operand.ku_)) {
927 for(j = 0; j < numCols_; j++) {
928 for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
929 if((*
this)(i, j) != Operand(i, j)) {
939template<
typename OrdinalType,
typename ScalarType>
942 return !((*this) == Operand);
949template<
typename OrdinalType,
typename ScalarType>
952 this->scale( alpha );
956template<
typename OrdinalType,
typename ScalarType>
963 for (j=0; j<numCols_; j++) {
964 ptr = values_ + j*stride_ + TEUCHOS_MAX(0, ku_-j);
965 for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
966 *ptr = alpha * (*ptr); ptr++;
969 updateFlops( (kl_+ku_+1)*numCols_ );
974template<
typename OrdinalType,
typename ScalarType>
982 if ((numRows_ != A.numRows_) || (numCols_ != A.numCols_) || (kl_ != A.kl_) || (ku_ != A.ku_)) {
985 for (j=0; j<numCols_; j++) {
986 ptr = values_ + j*stride_ + TEUCHOS_MAX(0, ku_-j);
987 for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
988 *ptr = A(i,j) * (*ptr); ptr++;
991 updateFlops( (kl_+ku_+1)*numCols_ );
996template<
typename OrdinalType,
typename ScalarType>
1001 os <<
"Values_copied : yes" << std::endl;
1003 os <<
"Values_copied : no" << std::endl;
1004 os <<
"Rows : " << numRows_ << std::endl;
1005 os <<
"Columns : " << numCols_ << std::endl;
1006 os <<
"Lower Bandwidth : " << kl_ << std::endl;
1007 os <<
"Upper Bandwidth : " << ku_ << std::endl;
1008 os <<
"LDA : " << stride_ << std::endl;
1009 if(numRows_ == 0 || numCols_ == 0) {
1010 os <<
"(matrix is empty, no values to display)" << std::endl;
1013 for(OrdinalType i = 0; i < numRows_; i++) {
1014 for (OrdinalType j=TEUCHOS_MAX(0,i-kl_); j<=TEUCHOS_MIN(numCols_-1,i+ku_); j++) {
1015 os << (*this)(i,j) <<
" ";
1027template<
typename OrdinalType,
typename ScalarType>
1031 rowIndex < TEUCHOS_MAX(0,colIndex-ku_) || rowIndex > TEUCHOS_MIN(numRows_-1,colIndex+kl_),
1033 "SerialBandDenseMatrix<T>::checkIndex: "
1034 "Row index " << rowIndex <<
" out of range [0, "<< numRows_ <<
")");
1036 "SerialBandDenseMatrix<T>::checkIndex: "
1037 "Col index " << colIndex <<
" out of range [0, "<< numCols_ <<
")");
1041template<
typename OrdinalType,
typename ScalarType>
1042void SerialBandDenseMatrix<OrdinalType, ScalarType>::deleteArrays(
void)
1044 if (valuesCopied_) {
1047 valuesCopied_ =
false;
1051template<
typename OrdinalType,
typename ScalarType>
1052void SerialBandDenseMatrix<OrdinalType, ScalarType>::copyMat(
1053 ScalarType* inputMatrix, OrdinalType strideInput, OrdinalType numRows_in,
1054 OrdinalType numCols_in, ScalarType* outputMatrix, OrdinalType strideOutput, OrdinalType startCol, ScalarType alpha
1058 ScalarType* ptr1 = 0;
1059 ScalarType* ptr2 = 0;
1061 for(j = 0; j < numCols_in; j++) {
1062 ptr1 = outputMatrix + (j * strideOutput) + TEUCHOS_MAX(0, ku_-j);
1063 ptr2 = inputMatrix + (j + startCol) * strideInput + TEUCHOS_MAX(0, ku_-j);
1065 for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_in-1,j+kl_); i++) {
1066 *ptr1++ += alpha*(*ptr2++);
1069 for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_in-1,j+kl_); i++) {
1077template<
typename OrdinalType,
typename ScalarType>
1087template<
typename OrdinalType,
typename ScalarType>
1089operator<<(std::ostream &out,
1092 printer.obj.print(out);
1097template<
typename OrdinalType,
typename ScalarType>
1098SerialBandDenseMatrixPrinter<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.
Reference-counted pointer class and non-member templated function implementations.
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 banded dense matrices of templated type.
SerialBandDenseMatrix< OrdinalType, ScalarType > & operator+=(const SerialBandDenseMatrix< OrdinalType, ScalarType > &Source)
Add another matrix to this matrix.
virtual ~SerialBandDenseMatrix()
Destructor.
int shapeUninitialized(OrdinalType numRows, OrdinalType numCols, OrdinalType kl, OrdinalType ku)
Same as shape() except leaves uninitialized.
SerialBandDenseMatrix()
Default Constructor.
int random()
Set all values in the matrix to be random numbers.
OrdinalType numCols() const
Returns the column dimension of this matrix.
virtual std::ostream & print(std::ostream &os) const
Print method. Defines the behavior of the std::ostream << operator.
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
Returns the Frobenius-norm of the matrix.
ScalarType scalarType
Typedef for scalar type.
SerialBandDenseMatrix< OrdinalType, ScalarType > & operator=(const SerialBandDenseMatrix< OrdinalType, ScalarType > &Source)
Copies values from one matrix to another.
bool empty() const
Returns whether this matrix is empty.
OrdinalType ordinalType
Typedef for ordinal type.
ScalarType * values() const
Data array access method.
OrdinalType numRows() const
Returns the row dimension of this matrix.
int scale(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.
ScalarType * operator[](OrdinalType colIndex)
Column access method (non-const).
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero())
Set all values in the matrix to a constant value.
SerialBandDenseMatrix< OrdinalType, ScalarType > & operator*=(const ScalarType alpha)
Scale this matrix by alpha; *this = alpha**this.
OrdinalType lowerBandwidth() const
Returns the lower bandwidth of this matrix.
bool operator==(const SerialBandDenseMatrix< OrdinalType, ScalarType > &Operand) const
Equality of two matrices.
int reshape(OrdinalType numRows, OrdinalType numCols, OrdinalType kl, OrdinalType ku)
Reshaping method for changing the size of a SerialBandDenseMatrix, keeping the entries.
OrdinalType stride() const
Returns the stride between the columns of this matrix in memory.
SerialBandDenseMatrix< OrdinalType, ScalarType > & assign(const SerialBandDenseMatrix< OrdinalType, ScalarType > &Source)
Copies values from one matrix to another.
bool operator!=(const SerialBandDenseMatrix< OrdinalType, ScalarType > &Operand) const
Inequality of two matrices.
ScalarType & operator()(OrdinalType rowIndex, OrdinalType colIndex)
Element access method (non-const).
int shape(OrdinalType numRows, OrdinalType numCols, OrdinalType kl, OrdinalType ku)
Shape method for changing the size of a SerialBandDenseMatrix, initializing entries to zero.
OrdinalType upperBandwidth() const
Returns the upper bandwidth of this matrix.
ScalarTraits< ScalarType >::magnitudeType normInf() const
Returns the Infinity-norm of the matrix.
ScalarTraits< ScalarType >::magnitudeType normOne() const
Returns the 1-norm of the matrix.
#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.
Ostream manipulator for SerialBandDenseMatrix.