Teuchos_SerialTriDiMatrix.hpp

Go to the documentation of this file.

// @HEADER
// ***********************************************************************
//
//                    Teuchos: Common Tools Package
//                 Copyright (2004) Sandia Corporation
//
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
//
// ***********************************************************************
// @HEADER


#ifndef _TEUCHOS_SERIALTRIDIMATRIX_HPP_
#define _TEUCHOS_SERIALTRIDIMATRIX_HPP_

#include "Teuchos_CompObject.hpp"
#include "Teuchos_BLAS.hpp"
#include "Teuchos_ScalarTraits.hpp"
#include "Teuchos_DataAccess.hpp"
#include "Teuchos_ConfigDefs.hpp"
#include "Teuchos_Assert.hpp"

namespace Teuchos {

template<typename OrdinalType, typename ScalarType>
class SerialTriDiMatrix : public CompObject, public Object, public BLAS<OrdinalType, ScalarType > {
public:

  typedef OrdinalType ordinalType;
  typedef ScalarType scalarType;




  SerialTriDiMatrix();


  SerialTriDiMatrix(OrdinalType numRows, OrdinalType numCols, bool zeroOut = true);


  SerialTriDiMatrix(DataAccess CV, ScalarType* values, OrdinalType numRowsCols);


  SerialTriDiMatrix(const SerialTriDiMatrix<OrdinalType, ScalarType> &Source, ETransp trans = Teuchos::NO_TRANS);


  SerialTriDiMatrix(DataAccess CV, const SerialTriDiMatrix<OrdinalType, ScalarType> &Source, OrdinalType numRowsCols, OrdinalType startRowCols=0);

  virtual ~SerialTriDiMatrix();



  int shape(OrdinalType numRows);

  int shapeUninitialized(OrdinalType numRows);


  int reshape(OrdinalType numRowsCols);





  SerialTriDiMatrix<OrdinalType, ScalarType>& operator= (const SerialTriDiMatrix<OrdinalType, ScalarType>& Source);


  SerialTriDiMatrix<OrdinalType, ScalarType>& assign (const SerialTriDiMatrix<OrdinalType, ScalarType>& Source);


  SerialTriDiMatrix<OrdinalType, ScalarType>& operator= (const ScalarType value) { putScalar(value); return(*this); }


  int putScalar( const ScalarType value = Teuchos::ScalarTraits<ScalarType>::zero() );

  //  int random();





  ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex);


  const ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex) const;


  //  ScalarType* operator [] (OrdinalType colIndex);


  //  const ScalarType* operator [] (OrdinalType colIndex) const;


  ScalarType* values() const { return(values_); }

  ScalarType* D() const { return D_;}
  ScalarType* DL() const { return DL_;}
  ScalarType* DU() const { return DU_;}
  ScalarType* DU2() const { return DU2_;}





  SerialTriDiMatrix<OrdinalType, ScalarType>& operator+= (const SerialTriDiMatrix<OrdinalType, ScalarType>& Source);


  SerialTriDiMatrix<OrdinalType, ScalarType>& operator-= (const SerialTriDiMatrix<OrdinalType, ScalarType>& Source);


  SerialTriDiMatrix<OrdinalType, ScalarType>& operator*= (const ScalarType alpha);


  int scale ( const ScalarType alpha );


  int scale ( const SerialTriDiMatrix<OrdinalType, ScalarType>& A );


  //int multiply (ETransp transa, ETransp transb, ScalarType alpha, const SerialTriDiMatrix<OrdinalType, ScalarType> &A, const SerialTriDiMatrix<OrdinalType, ScalarType> &B, ScalarType beta);


  //int multiply (ESide sideA, ScalarType alpha, const SerialSymTriDiMatrix<OrdinalType, ScalarType> &A, const SerialTriDiMatrix<OrdinalType, ScalarType> &B, ScalarType beta);





  bool operator== (const SerialTriDiMatrix<OrdinalType, ScalarType> &Operand) const;


  bool operator!= (const SerialTriDiMatrix<OrdinalType, ScalarType> &Operand) const;




  OrdinalType numRowsCols() const { return(numRowsCols_); }

  //  OrdinalType numCols() const { return(numRowsCols_); }

  //  OrdinalType stride() const { return(stride_); }

  bool empty() const { return(numRowsCols_ == 0); }



  typename ScalarTraits<ScalarType>::magnitudeType normOne() const;

  typename ScalarTraits<ScalarType>::magnitudeType normInf() const;

  typename ScalarTraits<ScalarType>::magnitudeType normFrobenius() const;


  virtual void print(std::ostream& os) const;

protected:
  void copyMat(SerialTriDiMatrix<OrdinalType, ScalarType> matrix,
               OrdinalType startCol,
               ScalarType alpha = ScalarTraits<ScalarType>::zero() );
  void deleteArrays();
  void checkIndex( OrdinalType rowIndex, OrdinalType colIndex = 0 ) const;
  OrdinalType numRowsCols_;

  bool valuesCopied_;
  ScalarType* values_;
  ScalarType* DL_;
  ScalarType* D_;
  ScalarType* DU_;
  ScalarType* DU2_;

}; // class Teuchos_SerialTriDiMatrix

//----------------------------------------------------------------------------------------------------
//  Constructors and Destructor
//----------------------------------------------------------------------------------------------------

template<typename OrdinalType, typename ScalarType>
SerialTriDiMatrix<OrdinalType, ScalarType>::SerialTriDiMatrix()
  :
    CompObject(),
    numRowsCols_(0),
    valuesCopied_(false),
    values_(0),
    DL_(NULL),
    D_(NULL),
    DU_(NULL),
    DU2_(NULL)
{}

template<typename OrdinalType, typename ScalarType>
SerialTriDiMatrix<OrdinalType, ScalarType>::SerialTriDiMatrix( OrdinalType numRowsCols_in, OrdinalType /* numCols_in */, bool zeroOut)
  : CompObject(), numRowsCols_(numRowsCols_in) {

  OrdinalType numvals = (numRowsCols_ == 1) ? 1 :  4*(numRowsCols_-1);
  values_ = new ScalarType [numvals];
  DL_  = values_;
  D_   = DL_  + (numRowsCols_-1);
  DU_  = D_   +  numRowsCols_;
  DU2_ = DU_  + (numRowsCols_-1);

  valuesCopied_ = true;
  if (zeroOut == true)
    putScalar();
}

template<typename OrdinalType, typename ScalarType>
SerialTriDiMatrix<OrdinalType, ScalarType>::SerialTriDiMatrix(DataAccess CV, ScalarType* values_in, OrdinalType numRowsCols_in  )
  : CompObject(), numRowsCols_(numRowsCols_in),
    valuesCopied_(false), values_(values_in)
{
  const OrdinalType numvals = (numRowsCols_ == 1) ? 1 :  4*(numRowsCols_-1);
  if(CV == Copy) {
    values_ = new ScalarType[numvals];
    valuesCopied_ = true;
  }
  else //CV == View
  {
    values_ = values_in;
    valuesCopied_ = false;
  }
  DL_ = values_;
  D_  = DL_    + (numRowsCols_-1);
  DU_ = D_     + numRowsCols_;
  DU2_ = DU_   + (numRowsCols_-1);
  if(CV == Copy) {
    for(OrdinalType i = 0 ; i < numRowsCols_ ; ++i )
      values_[i] = values_in[i];
  }
}

template<typename OrdinalType, typename ScalarType>
SerialTriDiMatrix<OrdinalType, ScalarType>::SerialTriDiMatrix(const SerialTriDiMatrix<OrdinalType, ScalarType> &Source, ETransp trans) : CompObject(), BLAS<OrdinalType,ScalarType>(), numRowsCols_(0), valuesCopied_(true), values_(0)
{
  if ( trans == Teuchos::NO_TRANS )  {
    numRowsCols_ = Source.numRowsCols_;

    const OrdinalType numvals = (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
    values_ = new ScalarType[numvals];
    DL_ = values_;
    D_  = DL_+ (numRowsCols_-1);
    DU_ = D_ + numRowsCols_;
    DU2_ = DU_ + (numRowsCols_-1);

    copyMat(Source, 0, 0);
  }
  else if ( trans == Teuchos::CONJ_TRANS && ScalarTraits<ScalarType>::isComplex )
    {
      numRowsCols_ = Source.numRowsCols_;
      const OrdinalType numvals = (numRowsCols_ == 1) ? 1 :  4*(numRowsCols_-1);
      values_ = new ScalarType[numvals];
      DL_ = values_;
      D_  = DL_+(numRowsCols_-1);
      DU_ = D_ + numRowsCols_;
      DU2_ = DU_ + (numRowsCols_-1);

      OrdinalType min = numRowsCols_;
      if(min > Source.numRowsCols_) min = Source.numRowsCols_;

      for(OrdinalType i = 0 ; i< min ; ++i) {
        D_[i] = Teuchos::ScalarTraits<ScalarType>::conjugate(Source.D_[i]);
        if(i < (min-1)) {
          DL_[i] = Teuchos::ScalarTraits<ScalarType>::conjugate(Source.DL_[i]);
          DU_[i] = Teuchos::ScalarTraits<ScalarType>::conjugate(Source.DU_[i]);
        }
        if(i < (min-2)) {
          DU2_[i] = Teuchos::ScalarTraits<ScalarType>::conjugate(Source.DU2_[i]);
        }
      }
    }
  else
    {
      numRowsCols_ = Source.numRowsCols_;
      const OrdinalType numvals = (numRowsCols_  == 1) ? 1 : 4*(numRowsCols_-1);
      values_ = new ScalarType[numvals];
      OrdinalType min = numRowsCols_;
      if(min > Source.numRowsCols_) min = Source.numRowsCols_;
      for(OrdinalType i = 0 ; i< min ; ++i) {
        D_[i] = Source.D_[i];
        if(i < (min-1)) {
          DL_[i] = Source.DL_[i];
          DU_[i] = Source.DU_[i];
        }
        if(i < (min-2)) {
          DU2_[i] = Source.DU2_[i];
        }
    }
  }
}

template<typename OrdinalType, typename ScalarType>
SerialTriDiMatrix<OrdinalType, ScalarType>::SerialTriDiMatrix(
  DataAccess CV, const SerialTriDiMatrix<OrdinalType, ScalarType> &Source,
  OrdinalType numRowsCols_in, OrdinalType startRow )
  : CompObject(), numRowsCols_(numRowsCols_in),
    valuesCopied_(false), values_(Source.values_) {
  if(CV == Copy)
    {
      const OrdinalType numvals = (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
      values_ = new ScalarType[numvals];
      copyMat(Source, startRow);
      valuesCopied_ = true;
    }
  else // CV == View
    {
      // \todo ???
      //    values_ = values_ + (stride_ * startCol) + startRow;
    }
}

template<typename OrdinalType, typename ScalarType>
SerialTriDiMatrix<OrdinalType, ScalarType>::~SerialTriDiMatrix()
{
  deleteArrays();
}

//----------------------------------------------------------------------------------------------------
//  Shape methods
//----------------------------------------------------------------------------------------------------

template<typename OrdinalType, typename ScalarType>
int SerialTriDiMatrix<OrdinalType, ScalarType>::shape( OrdinalType numRowsCols_in  )
{
  deleteArrays(); // Get rid of anything that might be already allocated
  numRowsCols_ = numRowsCols_in;
  const OrdinalType numvals = ( numRowsCols_ == 1) ? 1 :  4*(numRowsCols_-1);
  values_ = new ScalarType[numvals];

  putScalar();
  valuesCopied_ = true;
  return(0);
}

template<typename OrdinalType, typename ScalarType>
int SerialTriDiMatrix<OrdinalType, ScalarType>::shapeUninitialized(  OrdinalType numRowsCols_in  )
{
  deleteArrays(); // Get rid of anything that might be already allocated
  numRowsCols_ = numRowsCols_in;
  const OrdinalType numvals = ( numRowsCols_ == 1) ? 1 :  4*(numRowsCols_-1);
  values_ = new ScalarType[numvals];
  valuesCopied_ = true;
  return(0);
}

template<typename OrdinalType, typename ScalarType>
int SerialTriDiMatrix<OrdinalType, ScalarType>::reshape(  OrdinalType numRowsCols_in  )
{
  if(numRowsCols_in <1 ) {
    deleteArrays();
    return 0;
  }
  // Allocate space for new matrix
  const OrdinalType numvals = ( numRowsCols_in == 1) ? 1 :  4*(numRowsCols_in - 1);
  ScalarType* values_tmp = new ScalarType[numvals];
  ScalarType zero = ScalarTraits<ScalarType>::zero();
  for(OrdinalType i= 0; i<numvals ; ++i)
    values_tmp[i] = zero;

  OrdinalType min = TEUCHOS_MIN(numRowsCols_, numRowsCols_in);
  ScalarType* dl = values_tmp;
  ScalarType* d = values_tmp + (numRowsCols_in-1);
  ScalarType* du = d+(numRowsCols_in);
  ScalarType* du2 = du+(numRowsCols_in - 1);

  if(values_ != 0) {
   for(OrdinalType i = 0 ; i< min ; ++i) {
      dl[i] = DL_[i];
      d[i]  = D_[i];
      du[i] = DU_[i];
      du2[i] = DU2_[i];
    }
  }

  deleteArrays(); // Get rid of anything that might be already allocated
  numRowsCols_ = numRowsCols_in;

  values_ = values_tmp; // Set pointer to new A
  DL_ = dl;
  D_  = d;
  DU_ = du;
  DU2_ = du2;

  valuesCopied_ = true;
  return(0);
}

//----------------------------------------------------------------------------------------------------
//   Set methods
//----------------------------------------------------------------------------------------------------

template<typename OrdinalType, typename ScalarType>
int SerialTriDiMatrix<OrdinalType, ScalarType>::putScalar( const ScalarType value_in ) {
  // Set each value of the TriDi matrix to "value".
  const OrdinalType numvals = (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);

  for(OrdinalType i = 0; i<numvals ; ++i)
    {
      values_[i] = value_in;
    }
  return 0;
}

// template<typename OrdinalType, typename ScalarType>
// int SerialTriDiMatrix<OrdinalType, ScalarType>::random()
// {
//   // Set each value of the TriDi matrix to a random value.
//   for(OrdinalType j = 0; j < numCols_; j++)
//   {
//     for(OrdinalType i = 0; i < numRowsCols_; i++)
//        {
//          values_[i + j*stride_] = ScalarTraits<ScalarType>::random();
//        }
//   }
//   return 0;
// }

template<typename OrdinalType, typename ScalarType>
SerialTriDiMatrix<OrdinalType,ScalarType>&
SerialTriDiMatrix<OrdinalType, ScalarType>::operator=(const SerialTriDiMatrix<OrdinalType,ScalarType>& Source )
{
  if(this == &Source)
    return(*this); // Special case of source same as target
  if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
    return(*this); // Special case of both are views to same data.

  // If the source is a view then we will return a view, else we will return a copy.
  if (!Source.valuesCopied_) {
    if(valuesCopied_)       {
      // Clean up stored data if this was previously a copy.
      deleteArrays();
    }
    numRowsCols_ = Source.numRowsCols_;
    values_ = Source.values_;
  }
  else {
    // If we were a view, we will now be a copy.
    if(!valuesCopied_) {
      numRowsCols_ = Source.numRowsCols_;
      const OrdinalType numvals = ( Source.numRowsCols_ == 1) ? 1 :  4*(Source.numRowsCols_ - 1);
      if(numvals > 0) {
        values_ = new ScalarType[numvals];
        valuesCopied_ = true;
      }
      else {
        values_ = 0;
      }
    }
    // If we were a copy, we will stay a copy.
    else {
      // we need to allocate more space (or less space)
      deleteArrays();
      numRowsCols_ = Source.numRowsCols_;
      const OrdinalType numvals = ( Source.numRowsCols_ == 1) ? 1 :  4*(Source.numRowsCols_ - 1);
      if(numvals > 0) {
        values_ = new ScalarType[numvals];
        valuesCopied_ = true;
      }
    }
    DL_ = values_;
    if(values_ != 0) {
      D_  = DL_    + (numRowsCols_-1);
      DU_ = D_     +  numRowsCols_;
      DU2_ = DU_   + (numRowsCols_-1);

      OrdinalType min = TEUCHOS_MIN(numRowsCols_, Source.numRowsCols_);
      for(OrdinalType i = 0 ; i < min ; ++i ) {
        D_[i] = Source.D()[i];
        if(i< (min-1 ) )
          {
            DL_[i] = Source.DL()[i];
            DU_[i] = Source.DU()[i];
          }
        if(i< (min-2))
          DU2_[i] = Source.DU2()[i];
      }

    }
    else {
      D_ = DU_ = DU2_ = 0;
    }
  }
  return(*this);
}

template<typename OrdinalType, typename ScalarType>
SerialTriDiMatrix<OrdinalType, ScalarType>& SerialTriDiMatrix<OrdinalType, ScalarType>::operator+= (const SerialTriDiMatrix<OrdinalType,ScalarType>& Source )
{
  // Check for compatible dimensions
  if ((numRowsCols_ != Source.numRowsCols_) )
    {
      TEUCHOS_CHK_REF(*this); // Return *this without altering it.
    }
  copyMat(Source, 0, ScalarTraits<ScalarType>::one());
  return(*this);
}

template<typename OrdinalType, typename ScalarType>
SerialTriDiMatrix<OrdinalType, ScalarType>& SerialTriDiMatrix<OrdinalType, ScalarType>::operator-= (const SerialTriDiMatrix<OrdinalType,ScalarType>& Source )
{
  // Check for compatible dimensions
  if ((numRowsCols_ != Source.numRowsCols_) )
    {
      TEUCHOS_CHK_REF(*this); // Return *this without altering it.
    }
  copyMat(Source, 0, -ScalarTraits<ScalarType>::one());
  return(*this);
}

template<typename OrdinalType,typename ScalarType>
SerialTriDiMatrix<OrdinalType,ScalarType> & SerialTriDiMatrix<OrdinalType,ScalarType>::assign(const SerialTriDiMatrix<OrdinalType,ScalarType> & Source)
{
  if(this == &Source)
    return(*this); // Special case of source same as target
  if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
    return(*this); // Special case of both are views to same data.

  // Check for compatible dimensions
  if ((numRowsCols_ != Source.numRowsCols_) )
    {
      TEUCHOS_CHK_REF(*this); // Return *this without altering it.
    }
  copyMat(Source,  0, 0);
  return(*this);
}

//----------------------------------------------------------------------------------------------------
//   Accessor methods
//----------------------------------------------------------------------------------------------------

template<typename OrdinalType,typename ScalarType>
inline const ScalarType& SerialTriDiMatrix<OrdinalType,ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex) const
{
  OrdinalType diff = colIndex - rowIndex;

  //#ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
  checkIndex( rowIndex, colIndex );
  //#endif
  switch (diff) {
  case -1:
    return DL_[colIndex];
  case 0:
    return D_[colIndex];
  case 1:
    return DU_[rowIndex];
  case 2:
    return DU2_[rowIndex];
  default:
    TEUCHOS_TEST_FOR_EXCEPTION(true, std::out_of_range,
                               "SerialTriDiMatrix<T>::operator (row,col) "
                               "Index (" << rowIndex <<","<<colIndex<<") out of range ");
  }
}

template<typename OrdinalType,typename ScalarType>
inline ScalarType& Teuchos::SerialTriDiMatrix<OrdinalType,ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex)
{
  OrdinalType diff = colIndex - rowIndex;
  //#ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
  checkIndex( rowIndex, colIndex );
  //#endif
  switch (diff) {
  case -1:
    return DL_[colIndex];
  case 0:
    return D_[colIndex];
  case 1:
    return DU_[rowIndex];
  case 2:
    return DU2_[rowIndex];
  default:
    TEUCHOS_TEST_FOR_EXCEPTION(true, std::out_of_range,
                               "SerialTriDiMatrix<T>::operator (row,col) "
                               "Index (" << rowIndex <<","<<colIndex<<") out of range ");
  }
}

//----------------------------------------------------------------------------------------------------
//   Norm methods
//----------------------------------------------------------------------------------------------------

template<typename OrdinalType,typename ScalarType>
typename ScalarTraits<ScalarType>::magnitudeType SerialTriDiMatrix<OrdinalType,ScalarType>::normOne() const
{
  OrdinalType i, j;
  typename ScalarTraits<ScalarType>::magnitudeType anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
  typename ScalarTraits<ScalarType>::magnitudeType absSum = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());

  // Fix this for Tri DI

  for(j = 0; j < numRowsCols_; j++)
    {
      ScalarType sum = 0;
      if(j-1>=0) sum +=  ScalarTraits<ScalarType>::magnitude((*this)(j-1,j));
      sum+= ScalarTraits<ScalarType>::magnitude((*this)(j,j));
      if(j+1<numRowsCols_) sum+= ScalarTraits<ScalarType>::magnitude((*this)(j+1,j));
      absSum = ScalarTraits<ScalarType>::magnitude(sum);
      if(absSum > anorm)
        {
          anorm = absSum;
        }
    }
  updateFlops(numRowsCols_ * numRowsCols_);
  return(anorm);
}

template<typename OrdinalType, typename ScalarType>
typename ScalarTraits<ScalarType>::magnitudeType SerialTriDiMatrix<OrdinalType, ScalarType>::normInf() const
{
  OrdinalType i,j;
  typename ScalarTraits<ScalarType>::magnitudeType sum, anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());

  for (i = 0; i < numRowsCols_; i++) {
    sum = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
    for (j=i-1; j<= i+1; j++) {
      if(j >= 0 && j < numRowsCols_) sum += ScalarTraits<ScalarType>::magnitude((*this)(i,j));
    }
    anorm = TEUCHOS_MAX( anorm, sum );
  }
  updateFlops(numRowsCols_ * numRowsCols_);
  return(anorm);
}

template<typename OrdinalType, typename ScalarType>
typename ScalarTraits<ScalarType>::magnitudeType SerialTriDiMatrix<OrdinalType, ScalarType>::normFrobenius() const
{
  // \todo fix this
  OrdinalType i, j;
  typename ScalarTraits<ScalarType>::magnitudeType anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
  for (j = 0; j < numRowsCols_; j++) {
    for (i = j-1; i <= j+1; i++) {
      if(i >= 0 && i < numRowsCols_ ) anorm += ScalarTraits<ScalarType>::magnitude((*this)(i,j));
    }
  }
  anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::squareroot(anorm));
  updateFlops( (numRowsCols_ == 1) ? 1 :  4*(numRowsCols_-1) );
  return(anorm);
}

//----------------------------------------------------------------------------------------------------
//   Comparison methods
//----------------------------------------------------------------------------------------------------

template<typename OrdinalType, typename ScalarType>
bool SerialTriDiMatrix<OrdinalType, ScalarType>::operator== (const SerialTriDiMatrix<OrdinalType, ScalarType> &Operand) const
{
  bool allmatch = true;
  // bool result = 1; // unused
  if((numRowsCols_ != Operand.numRowsCols_) )
    {
      // result = 0; // unused
    }
  else
    {
      OrdinalType numvals = (numRowsCols_ == 1)? 1 : 4*(numRowsCols_ -1 );

      for(OrdinalType i = 0; i< numvals; ++i)
        allmatch &= (Operand.values_[i] == values_[i]);
    }
  return allmatch;
}

template<typename OrdinalType, typename ScalarType>
bool SerialTriDiMatrix<OrdinalType, ScalarType>::operator!= (const SerialTriDiMatrix<OrdinalType, ScalarType> &Operand) const {
  return !((*this) == Operand);
}

//----------------------------------------------------------------------------------------------------
//   Multiplication method
//----------------------------------------------------------------------------------------------------

template<typename OrdinalType, typename ScalarType>
SerialTriDiMatrix<OrdinalType, ScalarType>& SerialTriDiMatrix<OrdinalType, ScalarType>::operator*= (const ScalarType alpha )
{
  this->scale( alpha );
  return(*this);
}

template<typename OrdinalType, typename ScalarType>
int SerialTriDiMatrix<OrdinalType, ScalarType>::scale( const ScalarType alpha )
{
  OrdinalType i;
  OrdinalType numvals = (numRowsCols_ == 1)? 1 : 4*(numRowsCols_ -1 );
  for (i=0; i < numvals ; ++i ) {
    values_[i] *= alpha;
  }
  return(0);
}

template<typename OrdinalType, typename ScalarType>
int SerialTriDiMatrix<OrdinalType, ScalarType>::scale( const SerialTriDiMatrix<OrdinalType,ScalarType>& A )
{
  OrdinalType j;

  // Check for compatible dimensions
  if ((numRowsCols_ != A.numRowsCols_) )
    {
      TEUCHOS_CHK_ERR(-1); // Return error
    }
  OrdinalType numvals = (numRowsCols_ == 1)? 1 : 4*(numRowsCols_ -1 );
  for (j=0; j<numvals; j++) {
    values_[j] = A.values_ * values_[j];
  }
  updateFlops( numvals );
  return(0);
}

template<typename OrdinalType, typename ScalarType>
void SerialTriDiMatrix<OrdinalType, ScalarType>::print(std::ostream& os) const
{
  os << std::endl;
  if(valuesCopied_)
    os << "A_Copied: yes" << std::endl;
  else
    os << "A_Copied: no" << std::endl;
  os << "Rows and Columns: " << numRowsCols_ << std::endl;
  if(numRowsCols_ == 0) {
    os << "(matrix is empty, no values to display)" << std::endl;
    return;
  }
  else
    {
      os << "DL: "<<std::endl;
      for(int i=0;i<numRowsCols_-1;++i)
        os << DL_[i]<<" ";
      os << std::endl;
      os << "D: "<<std::endl;
      for(int i=0;i<numRowsCols_;++i)
        os << D_[i]<<" ";
      os << std::endl;
      os << "DU: "<<std::endl;
      for(int i=0;i<numRowsCols_-1;++i)
        os << DU_[i]<<" ";
      os << std::endl;
      os << "DU2: "<<std::endl;
      for(int i=0;i<numRowsCols_-2;++i)
        os << DU2_[i]<<" ";
      os << std::endl;
    }

  os <<" square format:"<<std::endl;
  for(int i=0 ; i < numRowsCols_ ; ++i )  {
    for(int j=0;j<numRowsCols_;++j)  {
      if ( j >= i-1  && j <= i+1) {
        os << (*this)(i,j)<<" ";
      }
      else
        os <<". ";
    }
    os << std::endl;
  }

}

//----------------------------------------------------------------------------------------------------
//   Protected methods
//----------------------------------------------------------------------------------------------------

template<typename OrdinalType, typename ScalarType>
inline void SerialTriDiMatrix<OrdinalType, ScalarType>::checkIndex( OrdinalType rowIndex, OrdinalType colIndex ) const
{
  OrdinalType diff = colIndex - rowIndex;
  TEUCHOS_TEST_FOR_EXCEPTION(rowIndex < 0 || rowIndex >= numRowsCols_, std::out_of_range,
                             "SerialTriDiMatrix<T>::checkIndex: "
                             "Row index " << rowIndex << " out of range [0, "<< numRowsCols_ << "]");
  TEUCHOS_TEST_FOR_EXCEPTION(colIndex < 0 || colIndex >= numRowsCols_,
                             std::out_of_range,
                             "SerialTriDiMatrix<T>::checkIndex: "
                             "Col index " << colIndex << " out of range [0, "<< numRowsCols_ << "]");
  TEUCHOS_TEST_FOR_EXCEPTION(diff > 2 || diff < -1 , std::out_of_range,
                             "SerialTriDiMatrix<T>::checkIndex: "
                             "index difference " << diff << " out of range [-1, 2]");
}

template<typename OrdinalType, typename ScalarType>
void SerialTriDiMatrix<OrdinalType, ScalarType>::deleteArrays(void)
{
  if (valuesCopied_)
    {
      delete [] values_;
      values_ = 0;
      valuesCopied_ = false;
    }
}

template<typename OrdinalType, typename ScalarType>
void SerialTriDiMatrix<OrdinalType, ScalarType>::copyMat(SerialTriDiMatrix<OrdinalType, ScalarType> inputMatrix,
                                                         OrdinalType startRowCol,
                                                         ScalarType alpha )
{
  OrdinalType i;
  OrdinalType max = inputMatrix.numRowsCols_;
  if(max > numRowsCols_ ) max = numRowsCols_;
  if(startRowCol > max ) return; //

  for(i = startRowCol ; i < max ; ++i) {

    if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
      // diagonal
      D()[i] += inputMatrix.D()[i];
      if(i<(max-1) && (i-1) >= startRowCol) {
        DL()[i] += inputMatrix.DL()[i];
        DU()[i] += inputMatrix.DU()[i];
      }
      if(i<(max-2) && (i-2) >= startRowCol) {
        DU2()[i] += inputMatrix.DU2()[i];
      }
    }
    else {
      // diagonal
      D()[i] = inputMatrix.D()[i];
      if(i<(max-1) && (i-1) >= startRowCol) {
        DL()[i] = inputMatrix.DL()[i];
        DU()[i] = inputMatrix.DU()[i];
      }
      if(i<(max-2) && (i-2) >= startRowCol) {
        DU2()[i] = inputMatrix.DU2()[i];
      }
    }
  }
}

}


#endif /* _TEUCHOS_SERIALTRIDIMATRIX_HPP_ */