43 #ifndef _TEUCHOS_SERIALTRIDIMATRIX_HPP_ 44 #define _TEUCHOS_SERIALTRIDIMATRIX_HPP_ 54 #include "Teuchos_Assert.hpp" 66 template<
typename OrdinalType,
typename ScalarType>
142 int shape(OrdinalType numRows);
210 ScalarType&
operator () (OrdinalType rowIndex, OrdinalType colIndex);
220 const ScalarType&
operator () (OrdinalType rowIndex, OrdinalType colIndex)
const;
244 ScalarType*
values()
const {
return(values_); }
246 ScalarType* D()
const {
return D_;}
247 ScalarType* DL()
const {
return DL_;}
248 ScalarType* DU()
const {
return DU_;}
249 ScalarType* DU2()
const {
return DU2_;}
279 int scale (
const ScalarType alpha );
353 bool empty()
const {
return(numRowsCols_ == 0); }
371 virtual void print(std::ostream& os)
const;
377 OrdinalType startCol,
380 void checkIndex( OrdinalType rowIndex, OrdinalType colIndex = 0 )
const;
381 OrdinalType numRowsCols_;
396 template<
typename OrdinalType,
typename ScalarType>
401 valuesCopied_(false),
409 template<
typename OrdinalType,
typename ScalarType>
411 :
CompObject(), numRowsCols_(numRowsCols_in) {
413 OrdinalType numvals = (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
414 values_ =
new ScalarType [numvals];
416 D_ = DL_ + (numRowsCols_-1);
417 DU_ = D_ + numRowsCols_;
418 DU2_ = DU_ + (numRowsCols_-1);
420 valuesCopied_ =
true;
425 template<
typename OrdinalType,
typename ScalarType>
428 valuesCopied_(false), values_(values_in)
431 const OrdinalType numvals = (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
432 values_ =
new ScalarType[numvals];
434 D_ = DL_ + (numRowsCols_-1);
435 DU_ = D_ + numRowsCols_;
436 DU2_ = DU_ + (numRowsCols_-1);
437 valuesCopied_ =
true;
439 for(OrdinalType i = 0 ; i < numRowsCols_ ; ++i )
440 values_[i] = values_in[i];
444 template<
typename OrdinalType,
typename ScalarType>
448 numRowsCols_ = Source.numRowsCols_;
450 const OrdinalType numvals = (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
451 values_ =
new ScalarType[numvals];
453 D_ = DL_+ (numRowsCols_-1);
454 DU_ = D_ + numRowsCols_;
455 DU2_ = DU_ + (numRowsCols_-1);
457 copyMat(Source, 0, 0);
461 numRowsCols_ = Source.numRowsCols_;
462 const OrdinalType numvals = (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
463 values_ =
new ScalarType[numvals];
465 D_ = DL_+(numRowsCols_-1);
466 DU_ = D_ + numRowsCols_;
467 DU2_ = DU_ + (numRowsCols_-1);
469 OrdinalType min = numRowsCols_;
470 if(min > Source.numRowsCols_) min = Source.numRowsCols_;
472 for(OrdinalType i = 0 ; i< min ; ++i) {
485 numRowsCols_ = Source.numCols_;
486 const OrdinalType numvals = (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
487 values_ =
new ScalarType[numvals];
488 OrdinalType min = numRowsCols_;
489 if(min > Source.numRowsCols_) min = Source.numRowsCols_;
490 for(OrdinalType i = 0 ; i< min ; ++i) {
497 DU2_ = Source.DU2_[i];
503 template<
typename OrdinalType,
typename ScalarType>
506 OrdinalType numRowsCols_in, OrdinalType startRow )
508 valuesCopied_(false), values_(Source.values_) {
511 const OrdinalType numvals = (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
512 values_ =
new ScalarType[numvals];
513 copyMat(Source, startRow);
514 valuesCopied_ =
true;
523 template<
typename OrdinalType,
typename ScalarType>
533 template<
typename OrdinalType,
typename ScalarType>
537 numRowsCols_ = numRowsCols_in;
538 const OrdinalType numvals = ( numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
539 values_ =
new ScalarType[numvals];
542 valuesCopied_ =
true;
546 template<
typename OrdinalType,
typename ScalarType>
550 numRowsCols_ = numRowsCols_in;
551 const OrdinalType numvals = ( numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
552 values_ =
new ScalarType[numvals];
553 valuesCopied_ =
true;
557 template<
typename OrdinalType,
typename ScalarType>
560 if(numRowsCols_in <1 ) {
565 const OrdinalType numvals = ( numRowsCols_in == 1) ? 1 : 4*(numRowsCols_in - 1);
566 ScalarType* values_tmp =
new ScalarType[numvals];
568 for(OrdinalType i= 0; i<numvals ; ++i)
569 values_tmp[i] = zero;
571 OrdinalType min = TEUCHOS_MIN(numRowsCols_, numRowsCols_in);
572 ScalarType* dl = values_tmp;
573 ScalarType* d = values_tmp + (numRowsCols_in-1);
574 ScalarType* du = d+(numRowsCols_in);
575 ScalarType* du2 = du+(numRowsCols_in - 1);
578 for(OrdinalType i = 0 ; i< min ; ++i) {
587 numRowsCols_ = numRowsCols_in;
589 values_ = values_tmp;
595 valuesCopied_ =
true;
603 template<
typename OrdinalType,
typename ScalarType>
606 const OrdinalType numvals = (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
608 for(OrdinalType i = 0; i<numvals ; ++i)
610 values_[i] = value_in;
629 template<
typename OrdinalType,
typename ScalarType>
635 if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
639 if (!Source.valuesCopied_) {
644 numRowsCols_ = Source.numRowsCols_;
645 values_ = Source.values_;
650 numRowsCols_ = Source.numRowsCols_;
651 const OrdinalType numvals = ( Source.numRowsCols_ == 1) ? 1 : 4*(Source.numRowsCols_ - 1);
653 values_ =
new ScalarType[numvals];
654 valuesCopied_ =
true;
664 numRowsCols_ = Source.numRowsCols_;
665 const OrdinalType numvals = ( Source.numRowsCols_ == 1) ? 1 : 4*(Source.numRowsCols_ - 1);
667 values_ =
new ScalarType[numvals];
668 valuesCopied_ =
true;
673 D_ = DL_ + (numRowsCols_-1);
674 DU_ = D_ + numRowsCols_;
675 DU2_ = DU_ + (numRowsCols_-1);
677 OrdinalType min = TEUCHOS_MIN(numRowsCols_, Source.numRowsCols_);
678 for(OrdinalType i = 0 ; i < min ; ++i ) {
679 D_[i] = Source.D()[i];
682 DL_[i] = Source.DL()[i];
683 DU_[i] = Source.DU()[i];
686 DU2_[i] = Source.DU2()[i];
697 template<
typename OrdinalType,
typename ScalarType>
701 if ((numRowsCols_ != Source.numRowsCols_) )
703 TEUCHOS_CHK_REF(*
this);
709 template<
typename OrdinalType,
typename ScalarType>
713 if ((numRowsCols_ != Source.numRowsCols_) )
715 TEUCHOS_CHK_REF(*
this);
721 template<
typename OrdinalType,
typename ScalarType>
726 if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
730 if ((numRowsCols_ != Source.numRowsCols_) )
732 TEUCHOS_CHK_REF(*
this);
734 copyMat(Source, 0, 0);
742 template<
typename OrdinalType,
typename ScalarType>
745 OrdinalType diff = colIndex - rowIndex;
748 checkIndex( rowIndex, colIndex );
752 return DL_[colIndex];
758 return DU_[rowIndex];
761 return DU2_[rowIndex];
765 "SerialTriDiMatrix<T>::operator (row,col) " 766 "Index (" << rowIndex <<
","<<colIndex<<
") out of range ");
769 "SerialTriDiMatrix<T>::operator (row,col) " 770 "Index (" << rowIndex <<
","<<colIndex<<
") out of range ");
774 template<
typename OrdinalType,
typename ScalarType>
777 OrdinalType diff = colIndex - rowIndex;
779 checkIndex( rowIndex, colIndex );
783 return DL_[colIndex];
789 return DU_[rowIndex];
792 return DU2_[rowIndex];
796 "SerialTriDiMatrix<T>::operator (row,col) " 797 "Index (" << rowIndex <<
","<<colIndex<<
") out of range ");
800 "SerialTriDiMatrix<T>::operator (row,col) " 801 "Index (" << rowIndex <<
","<<colIndex<<
") out of range ");
810 template<
typename OrdinalType,
typename ScalarType>
820 for(j = 0; j < numRowsCols_; j++)
836 template<
typename OrdinalType,
typename ScalarType>
842 for (i = 0; i < numRowsCols_; i++) {
844 for (j=i-1; j<= i+1; j++) {
847 anorm = TEUCHOS_MAX( anorm, sum );
853 template<
typename OrdinalType,
typename ScalarType>
859 for (j = 0; j < numRowsCols_; j++) {
860 for (i = j-1; i <= j+1; i++) {
865 updateFlops( (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1) );
873 template<
typename OrdinalType,
typename ScalarType>
876 bool allmatch =
true;
878 if((numRowsCols_ != Operand.numRowsCols_) )
884 OrdinalType numvals = (numRowsCols_ == 1)? 1 : 4*(numRowsCols_ -1 );
886 for(OrdinalType i = 0; i< numvals; ++i)
887 allmatch &= (Operand.values_[i] == values_[i]);
892 template<
typename OrdinalType,
typename ScalarType>
894 return !((*this) == Operand);
901 template<
typename OrdinalType,
typename ScalarType>
904 this->
scale( alpha );
908 template<
typename OrdinalType,
typename ScalarType>
912 OrdinalType numvals = (numRowsCols_ == 1)? 1 : 4*(numRowsCols_ -1 );
913 for (i=0; i < numvals ; ++i ) {
919 template<
typename OrdinalType,
typename ScalarType>
925 if ((numRowsCols_ != A.numRowsCols_) )
929 OrdinalType numvals = (numRowsCols_ == 1)? 1 : 4*(numRowsCols_ -1 );
930 for (j=0; j<numvals; j++) {
931 values_[j] = A.values_ * values_[j];
937 template<
typename OrdinalType,
typename ScalarType>
942 os <<
"A_Copied: yes" << std::endl;
944 os <<
"A_Copied: no" << std::endl;
945 os <<
"Rows and Columns: " << numRowsCols_ << std::endl;
946 if(numRowsCols_ == 0) {
947 os <<
"(matrix is empty, no values to display)" << std::endl;
952 os <<
"DL: "<<std::endl;
953 for(
int i=0;i<numRowsCols_-1;++i)
956 os <<
"D: "<<std::endl;
957 for(
int i=0;i<numRowsCols_;++i)
960 os <<
"DU: "<<std::endl;
961 for(
int i=0;i<numRowsCols_-1;++i)
964 os <<
"DU2: "<<std::endl;
965 for(
int i=0;i<numRowsCols_-2;++i)
970 os <<
" square format:"<<std::endl;
971 for(
int i=0 ; i < numRowsCols_ ; ++i ) {
972 for(
int j=0;j<numRowsCols_;++j) {
973 if ( j >= i-1 && j <= i+1) {
974 os << (*this)(i,j)<<
" ";
988 template<
typename OrdinalType,
typename ScalarType>
991 OrdinalType diff = colIndex - rowIndex;
993 "SerialTriDiMatrix<T>::checkIndex: " 994 "Row index " << rowIndex <<
" out of range [0, "<< numRowsCols_ <<
"]");
997 "SerialTriDiMatrix<T>::checkIndex: " 998 "Col index " << colIndex <<
" out of range [0, "<< numRowsCols_ <<
"]");
1000 "SerialTriDiMatrix<T>::checkIndex: " 1001 "index difference " << diff <<
" out of range [-1, 2]");
1004 template<
typename OrdinalType,
typename ScalarType>
1011 valuesCopied_ =
false;
1015 template<
typename OrdinalType,
typename ScalarType>
1017 OrdinalType startRowCol,
1021 OrdinalType max = inputMatrix.numRowsCols_;
1022 if(max > numRowsCols_ ) max = numRowsCols_;
1023 if(startRowCol > max )
return;
1025 for(i = startRowCol ; i < max ; ++i) {
1029 D()[i] += inputMatrix.D()[i];
1030 if(i<(max-1) && (i-1) >= startRowCol) {
1031 DL()[i] += inputMatrix.DL()[i];
1032 DU()[i] += inputMatrix.DU()[i];
1034 if(i<(max-2) && (i-2) >= startRowCol) {
1035 DU2()[i] += inputMatrix.DU2()[i];
1040 D()[i] = inputMatrix.D()[i];
1041 if(i<(max-1) && (i-1) >= startRowCol) {
1042 DL()[i] = inputMatrix.DL()[i];
1043 DU()[i] = inputMatrix.DU()[i];
1045 if(i<(max-2) && (i-2) >= startRowCol) {
1046 DU2()[i] = inputMatrix.DU2()[i];
ScalarType & operator()(OrdinalType rowIndex, OrdinalType colIndex)
Element access method (non-const).
ScalarTraits< ScalarType >::magnitudeType normInf() const
Returns the Infinity-norm of the matrix.
int shape(OrdinalType numRows)
Shape method for changing the size of a SerialTriDiMatrix, initializing entries to zero...
OrdinalType numRowsCols() const
Returns the row dimension of this matrix.
ScalarType scalarType
Typedef for scalar type.
Templated interface class to BLAS routines.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Macro for throwing an exception with breakpointing to ease debugging.
virtual void print(std::ostream &os) const
Print method. Defines the behavior of the std::ostream << operator inherited from the Object class...
virtual ~SerialTriDiMatrix()
Destructor.
Teuchos header file which uses auto-configuration information to include necessary C++ headers...
SerialTriDiMatrix< OrdinalType, ScalarType > & operator=(const SerialTriDiMatrix< OrdinalType, ScalarType > &Source)
Copies values from one matrix to another.
SerialTriDiMatrix< OrdinalType, ScalarType > & assign(const SerialTriDiMatrix< OrdinalType, ScalarType > &Source)
Copies values from one matrix to another.
Object for storing data and providing functionality that is common to all computational classes...
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
Returns the Frobenius-norm of the matrix.
This structure defines some basic traits for a scalar field type.
This class creates and provides basic support for TriDi matrix of templated type. ...
ScalarType * values() const
Column access method (non-const).
Functionality and data that is common to all computational classes.
static T conjugate(T a)
Returns the conjugate of the scalar type a.
ScalarTraits< ScalarType >::magnitudeType normOne() const
Returns the 1-norm of the matrix.
int scale(const ScalarType alpha)
Scale this matrix by alpha; *this = alpha**this.
bool operator==(const SerialTriDiMatrix< OrdinalType, ScalarType > &Operand) const
Equality of two matrices.
SerialTriDiMatrix< OrdinalType, ScalarType > & operator*=(const ScalarType alpha)
Scale this matrix by alpha; *this = alpha**this.
bool empty() const
Returns the column dimension of this matrix.
static magnitudeType magnitude(T a)
Returns the magnitudeType of the scalar type a.
SerialTriDiMatrix()
Default Constructor.
SerialTriDiMatrix< OrdinalType, ScalarType > & operator+=(const SerialTriDiMatrix< OrdinalType, ScalarType > &Source)
Add another matrix to this matrix.
OrdinalType ordinalType
Typedef for ordinal type.
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.
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero())
Set all values in the matrix to a constant value.
int reshape(OrdinalType numRowsCols)
Reshaping method for changing the size of a SerialTriDiMatrix, keeping the entries.
bool operator!=(const SerialTriDiMatrix< OrdinalType, ScalarType > &Operand) const
Inequality of two matrices.
SerialTriDiMatrix< OrdinalType, ScalarType > & operator-=(const SerialTriDiMatrix< OrdinalType, ScalarType > &Source)
Subtract another matrix from this matrix.
int shapeUninitialized(OrdinalType numRows)
Same as shape() except leaves uninitialized.
Teuchos::DataAccess Mode enumerable type.