43#ifndef _TEUCHOS_SERIALTRIDIMATRIX_HPP_
44#define _TEUCHOS_SERIALTRIDIMATRIX_HPP_
54#include "Teuchos_Assert.hpp"
66template<
typename OrdinalType,
typename ScalarType>
141 int shape(OrdinalType numRows);
209 ScalarType&
operator () (OrdinalType rowIndex, OrdinalType colIndex);
219 const ScalarType&
operator () (OrdinalType rowIndex, OrdinalType colIndex)
const;
243 ScalarType*
values()
const {
return(values_); }
245 ScalarType* D()
const {
return D_;}
246 ScalarType* DL()
const {
return DL_;}
247 ScalarType* DU()
const {
return DU_;}
248 ScalarType* DU2()
const {
return DU2_;}
259 SerialTriDiMatrix<OrdinalType, ScalarType>&
operator+= (
const SerialTriDiMatrix<OrdinalType, ScalarType>& Source);
265 SerialTriDiMatrix<OrdinalType, ScalarType>&
operator-= (
const SerialTriDiMatrix<OrdinalType, ScalarType>& Source);
271 SerialTriDiMatrix<OrdinalType, ScalarType>&
operator*= (
const ScalarType alpha);
278 int scale (
const ScalarType alpha );
287 int scale (
const SerialTriDiMatrix<OrdinalType, ScalarType>& A );
329 bool operator== (
const SerialTriDiMatrix<OrdinalType, ScalarType> &Operand)
const;
335 bool operator!= (
const SerialTriDiMatrix<OrdinalType, ScalarType> &Operand)
const;
352 bool empty()
const {
return(numRowsCols_ == 0); }
371 virtual void print(std::ostream& os)
const;
376 OrdinalType startCol,
379 void checkIndex( OrdinalType rowIndex, OrdinalType colIndex = 0 )
const;
380 OrdinalType numRowsCols_;
395template<
typename OrdinalType,
typename ScalarType>
400 valuesCopied_(false),
408template<
typename OrdinalType,
typename ScalarType>
410 :
CompObject(), numRowsCols_(numRowsCols_in) {
412 OrdinalType numvals = (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
413 values_ =
new ScalarType [numvals];
415 D_ = DL_ + (numRowsCols_-1);
416 DU_ = D_ + numRowsCols_;
417 DU2_ = DU_ + (numRowsCols_-1);
419 valuesCopied_ =
true;
424template<
typename OrdinalType,
typename ScalarType>
427 valuesCopied_(false), values_(values_in)
429 const OrdinalType numvals = (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
431 values_ =
new ScalarType[numvals];
432 valuesCopied_ =
true;
437 valuesCopied_ =
false;
440 D_ = DL_ + (numRowsCols_-1);
441 DU_ = D_ + numRowsCols_;
442 DU2_ = DU_ + (numRowsCols_-1);
444 for(OrdinalType i = 0 ; i < numRowsCols_ ; ++i )
445 values_[i] = values_in[i];
449template<
typename OrdinalType,
typename ScalarType>
453 numRowsCols_ = Source.numRowsCols_;
455 const OrdinalType numvals = (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
456 values_ =
new ScalarType[numvals];
458 D_ = DL_+ (numRowsCols_-1);
459 DU_ = D_ + numRowsCols_;
460 DU2_ = DU_ + (numRowsCols_-1);
462 copyMat(Source, 0, 0);
466 numRowsCols_ = Source.numRowsCols_;
467 const OrdinalType numvals = (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
468 values_ =
new ScalarType[numvals];
470 D_ = DL_+(numRowsCols_-1);
471 DU_ = D_ + numRowsCols_;
472 DU2_ = DU_ + (numRowsCols_-1);
474 OrdinalType min = numRowsCols_;
475 if(min > Source.numRowsCols_) min = Source.numRowsCols_;
477 for(OrdinalType i = 0 ; i< min ; ++i) {
490 numRowsCols_ = Source.numRowsCols_;
491 const OrdinalType numvals = (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
492 values_ =
new ScalarType[numvals];
493 OrdinalType min = numRowsCols_;
494 if(min > Source.numRowsCols_) min = Source.numRowsCols_;
495 for(OrdinalType i = 0 ; i< min ; ++i) {
496 D_[i] = Source.D_[i];
498 DL_[i] = Source.DL_[i];
499 DU_[i] = Source.DU_[i];
502 DU2_[i] = Source.DU2_[i];
508template<
typename OrdinalType,
typename ScalarType>
511 OrdinalType numRowsCols_in, OrdinalType startRow )
513 valuesCopied_(false), values_(Source.values_) {
516 const OrdinalType numvals = (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
517 values_ =
new ScalarType[numvals];
518 copyMat(Source, startRow);
519 valuesCopied_ =
true;
528template<
typename OrdinalType,
typename ScalarType>
538template<
typename OrdinalType,
typename ScalarType>
542 numRowsCols_ = numRowsCols_in;
543 const OrdinalType numvals = ( numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
544 values_ =
new ScalarType[numvals];
547 valuesCopied_ =
true;
551template<
typename OrdinalType,
typename ScalarType>
555 numRowsCols_ = numRowsCols_in;
556 const OrdinalType numvals = ( numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
557 values_ =
new ScalarType[numvals];
558 valuesCopied_ =
true;
562template<
typename OrdinalType,
typename ScalarType>
565 if(numRowsCols_in <1 ) {
570 const OrdinalType numvals = ( numRowsCols_in == 1) ? 1 : 4*(numRowsCols_in - 1);
571 ScalarType* values_tmp =
new ScalarType[numvals];
573 for(OrdinalType i= 0; i<numvals ; ++i)
574 values_tmp[i] = zero;
576 OrdinalType min = TEUCHOS_MIN(numRowsCols_, numRowsCols_in);
577 ScalarType* dl = values_tmp;
578 ScalarType* d = values_tmp + (numRowsCols_in-1);
579 ScalarType* du = d+(numRowsCols_in);
580 ScalarType* du2 = du+(numRowsCols_in - 1);
583 for(OrdinalType i = 0 ; i< min ; ++i) {
592 numRowsCols_ = numRowsCols_in;
594 values_ = values_tmp;
600 valuesCopied_ =
true;
608template<
typename OrdinalType,
typename ScalarType>
611 const OrdinalType numvals = (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
613 for(OrdinalType i = 0; i<numvals ; ++i)
615 values_[i] = value_in;
634template<
typename OrdinalType,
typename ScalarType>
640 if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
644 if (!Source.valuesCopied_) {
649 numRowsCols_ = Source.numRowsCols_;
650 values_ = Source.values_;
655 numRowsCols_ = Source.numRowsCols_;
656 const OrdinalType numvals = ( Source.numRowsCols_ == 1) ? 1 : 4*(Source.numRowsCols_ - 1);
658 values_ =
new ScalarType[numvals];
659 valuesCopied_ =
true;
669 numRowsCols_ = Source.numRowsCols_;
670 const OrdinalType numvals = ( Source.numRowsCols_ == 1) ? 1 : 4*(Source.numRowsCols_ - 1);
672 values_ =
new ScalarType[numvals];
673 valuesCopied_ =
true;
678 D_ = DL_ + (numRowsCols_-1);
679 DU_ = D_ + numRowsCols_;
680 DU2_ = DU_ + (numRowsCols_-1);
682 OrdinalType min = TEUCHOS_MIN(numRowsCols_, Source.numRowsCols_);
683 for(OrdinalType i = 0 ; i < min ; ++i ) {
684 D_[i] = Source.D()[i];
687 DL_[i] = Source.DL()[i];
688 DU_[i] = Source.DU()[i];
691 DU2_[i] = Source.DU2()[i];
702template<
typename OrdinalType,
typename ScalarType>
706 if ((numRowsCols_ != Source.numRowsCols_) )
708 TEUCHOS_CHK_REF(*
this);
714template<
typename OrdinalType,
typename ScalarType>
718 if ((numRowsCols_ != Source.numRowsCols_) )
720 TEUCHOS_CHK_REF(*
this);
726template<
typename OrdinalType,
typename ScalarType>
731 if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
735 if ((numRowsCols_ != Source.numRowsCols_) )
737 TEUCHOS_CHK_REF(*
this);
739 copyMat(Source, 0, 0);
747template<
typename OrdinalType,
typename ScalarType>
750 OrdinalType diff = colIndex - rowIndex;
753 checkIndex( rowIndex, colIndex );
757 return DL_[colIndex];
761 return DU_[rowIndex];
763 return DU2_[rowIndex];
766 "SerialTriDiMatrix<T>::operator (row,col) "
767 "Index (" << rowIndex <<
","<<colIndex<<
") out of range ");
771template<
typename OrdinalType,
typename ScalarType>
774 OrdinalType diff = colIndex - rowIndex;
776 checkIndex( rowIndex, colIndex );
780 return DL_[colIndex];
784 return DU_[rowIndex];
786 return DU2_[rowIndex];
789 "SerialTriDiMatrix<T>::operator (row,col) "
790 "Index (" << rowIndex <<
","<<colIndex<<
") out of range ");
798template<
typename OrdinalType,
typename ScalarType>
807 for(j = 0; j < numRowsCols_; j++)
819 updateFlops(numRowsCols_ * numRowsCols_);
823template<
typename OrdinalType,
typename ScalarType>
829 for (i = 0; i < numRowsCols_; i++) {
831 for (j=i-1; j<= i+1; j++) {
834 anorm = TEUCHOS_MAX( anorm, sum );
836 updateFlops(numRowsCols_ * numRowsCols_);
840template<
typename OrdinalType,
typename ScalarType>
846 for (j = 0; j < numRowsCols_; j++) {
847 for (i = j-1; i <= j+1; i++) {
852 updateFlops( (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1) );
860template<
typename OrdinalType,
typename ScalarType>
863 bool allmatch =
true;
865 if((numRowsCols_ != Operand.numRowsCols_) )
871 OrdinalType numvals = (numRowsCols_ == 1)? 1 : 4*(numRowsCols_ -1 );
873 for(OrdinalType i = 0; i< numvals; ++i)
874 allmatch &= (Operand.values_[i] == values_[i]);
879template<
typename OrdinalType,
typename ScalarType>
881 return !((*this) == Operand);
888template<
typename OrdinalType,
typename ScalarType>
891 this->scale( alpha );
895template<
typename OrdinalType,
typename ScalarType>
899 OrdinalType numvals = (numRowsCols_ == 1)? 1 : 4*(numRowsCols_ -1 );
900 for (i=0; i < numvals ; ++i ) {
906template<
typename OrdinalType,
typename ScalarType>
912 if ((numRowsCols_ != A.numRowsCols_) )
916 OrdinalType numvals = (numRowsCols_ == 1)? 1 : 4*(numRowsCols_ -1 );
917 for (j=0; j<numvals; j++) {
918 values_[j] = A.values_ * values_[j];
920 updateFlops( numvals );
924template<
typename OrdinalType,
typename ScalarType>
929 os <<
"A_Copied: yes" << std::endl;
931 os <<
"A_Copied: no" << std::endl;
932 os <<
"Rows and Columns: " << numRowsCols_ << std::endl;
933 if(numRowsCols_ == 0) {
934 os <<
"(matrix is empty, no values to display)" << std::endl;
939 os <<
"DL: "<<std::endl;
940 for(
int i=0;i<numRowsCols_-1;++i)
943 os <<
"D: "<<std::endl;
944 for(
int i=0;i<numRowsCols_;++i)
947 os <<
"DU: "<<std::endl;
948 for(
int i=0;i<numRowsCols_-1;++i)
951 os <<
"DU2: "<<std::endl;
952 for(
int i=0;i<numRowsCols_-2;++i)
957 os <<
" square format:"<<std::endl;
958 for(
int i=0 ; i < numRowsCols_ ; ++i ) {
959 for(
int j=0;j<numRowsCols_;++j) {
960 if ( j >= i-1 && j <= i+1) {
961 os << (*this)(i,j)<<
" ";
975template<
typename OrdinalType,
typename ScalarType>
978 OrdinalType diff = colIndex - rowIndex;
980 "SerialTriDiMatrix<T>::checkIndex: "
981 "Row index " << rowIndex <<
" out of range [0, "<< numRowsCols_ <<
"]");
984 "SerialTriDiMatrix<T>::checkIndex: "
985 "Col index " << colIndex <<
" out of range [0, "<< numRowsCols_ <<
"]");
987 "SerialTriDiMatrix<T>::checkIndex: "
988 "index difference " << diff <<
" out of range [-1, 2]");
991template<
typename OrdinalType,
typename ScalarType>
992void SerialTriDiMatrix<OrdinalType, ScalarType>::deleteArrays(
void)
998 valuesCopied_ =
false;
1002template<
typename OrdinalType,
typename ScalarType>
1003void SerialTriDiMatrix<OrdinalType, ScalarType>::copyMat(SerialTriDiMatrix<OrdinalType, ScalarType> inputMatrix,
1004 OrdinalType startRowCol,
1008 OrdinalType max = inputMatrix.numRowsCols_;
1009 if(max > numRowsCols_ ) max = numRowsCols_;
1010 if(startRowCol > max )
return;
1012 for(i = startRowCol ; i < max ; ++i) {
1016 D()[i] += inputMatrix.D()[i];
1017 if(i<(max-1) && (i-1) >= startRowCol) {
1018 DL()[i] += inputMatrix.DL()[i];
1019 DU()[i] += inputMatrix.DU()[i];
1021 if(i<(max-2) && (i-2) >= startRowCol) {
1022 DU2()[i] += inputMatrix.DU2()[i];
1027 D()[i] = inputMatrix.D()[i];
1028 if(i<(max-1) && (i-1) >= startRowCol) {
1029 DL()[i] = inputMatrix.DL()[i];
1030 DU()[i] = inputMatrix.DU()[i];
1032 if(i<(max-2) && (i-2) >= startRowCol) {
1033 DU2()[i] = inputMatrix.DU2()[i];
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 TriDi matrix of templated type.
int reshape(OrdinalType numRowsCols)
Reshaping method for changing the size of a SerialTriDiMatrix, keeping the entries.
virtual void print(std::ostream &os) const
Print method. Defines the behavior of the std::ostream << operator inherited from the Object class.
ScalarType scalarType
Typedef for scalar type.
bool operator==(const SerialTriDiMatrix< OrdinalType, ScalarType > &Operand) const
Equality of two matrices.
SerialTriDiMatrix< OrdinalType, ScalarType > & assign(const SerialTriDiMatrix< OrdinalType, ScalarType > &Source)
Copies values from one matrix to another.
ScalarTraits< ScalarType >::magnitudeType normOne() const
Returns the 1-norm of the matrix.
int shape(OrdinalType numRows)
Shape method for changing the size of a SerialTriDiMatrix, initializing entries to zero.
int scale(const ScalarType alpha)
Scale this matrix by alpha; *this = alpha**this.
SerialTriDiMatrix()
Default Constructor.
SerialTriDiMatrix< OrdinalType, ScalarType > & operator+=(const SerialTriDiMatrix< OrdinalType, ScalarType > &Source)
Add another matrix to this matrix.
SerialTriDiMatrix< OrdinalType, ScalarType > & operator-=(const SerialTriDiMatrix< OrdinalType, ScalarType > &Source)
Subtract another matrix from this matrix.
bool empty() const
Returns the column dimension of this matrix.
OrdinalType numRowsCols() const
Returns the row dimension of this matrix.
ScalarTraits< ScalarType >::magnitudeType normInf() const
Returns the Infinity-norm of the matrix.
int shapeUninitialized(OrdinalType numRows)
Same as shape() except leaves uninitialized.
virtual ~SerialTriDiMatrix()
Destructor.
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero())
Set all values in the matrix to a constant value.
SerialTriDiMatrix< OrdinalType, ScalarType > & operator*=(const ScalarType alpha)
Scale this matrix by alpha; *this = alpha**this.
SerialTriDiMatrix< OrdinalType, ScalarType > & operator=(const SerialTriDiMatrix< OrdinalType, ScalarType > &Source)
Copies values from one matrix to another.
ScalarType & operator()(OrdinalType rowIndex, OrdinalType colIndex)
Element access method (non-const).
OrdinalType ordinalType
Typedef for ordinal type.
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
Returns the Frobenius-norm of the matrix.
bool operator!=(const SerialTriDiMatrix< OrdinalType, ScalarType > &Operand) const
Inequality of two matrices.
ScalarType * values() const
Column access method (non-const).
#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,...
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 conjugate(T a)
Returns the conjugate of the scalar type a.