Teuchos_SerialBandDenseMatrix.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_SERIALBANDDENSEMATRIX_HPP_
#define _TEUCHOS_SERIALBANDDENSEMATRIX_HPP_

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

namespace Teuchos {

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

  typedef OrdinalType ordinalType;
  typedef ScalarType scalarType;




  SerialBandDenseMatrix();


  SerialBandDenseMatrix(OrdinalType numRows, OrdinalType numCols, OrdinalType kl, OrdinalType ku, bool zeroOut = true);


  SerialBandDenseMatrix(DataAccess CV, ScalarType* values, OrdinalType stride, OrdinalType numRows, OrdinalType numCols, OrdinalType kl, OrdinalType ku);


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


  SerialBandDenseMatrix(DataAccess CV, const SerialBandDenseMatrix<OrdinalType, ScalarType> &Source, OrdinalType numRows, OrdinalType numCols, OrdinalType startCol=0);

  virtual ~SerialBandDenseMatrix();



  int shape(OrdinalType numRows, OrdinalType numCols, OrdinalType kl, OrdinalType ku);

  int shapeUninitialized(OrdinalType numRows, OrdinalType numCols, OrdinalType kl, OrdinalType ku);


  int reshape(OrdinalType numRows, OrdinalType numCols, OrdinalType kl, OrdinalType ku);






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


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


  SerialBandDenseMatrix<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_); }





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


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


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


  int scale ( const ScalarType alpha );


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





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


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




  OrdinalType numRows() const { return(numRows_); }

  OrdinalType numCols() const { return(numCols_); }

  OrdinalType lowerBandwidth() const { return(kl_); }

  OrdinalType upperBandwidth() const { return(ku_); }

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

  bool empty() const { return(numRows_ == 0 || numCols_ == 0); }



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

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

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



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


protected:
  void copyMat(ScalarType* inputMatrix, OrdinalType strideInput,
    OrdinalType numRows, OrdinalType numCols, ScalarType* outputMatrix,
    OrdinalType strideOutput, OrdinalType startCol, ScalarType alpha = ScalarTraits<ScalarType>::zero() );
  void deleteArrays();
  void checkIndex( OrdinalType rowIndex, OrdinalType colIndex = 0 ) const;
  OrdinalType numRows_;
  OrdinalType numCols_;
  OrdinalType stride_;
  OrdinalType kl_;
  OrdinalType ku_;
  bool valuesCopied_;
  ScalarType* values_;

}; // class Teuchos_SerialBandDenseMatrix

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

template<typename OrdinalType, typename ScalarType>
SerialBandDenseMatrix<OrdinalType, ScalarType>::SerialBandDenseMatrix ()
  : CompObject (),
    BLAS<OrdinalType,ScalarType>(),
    numRows_ (0),
    numCols_ (0),
    stride_ (0),
    kl_ (0),
    ku_ (0),
    valuesCopied_ (false),
    values_ (0)
{}

template<typename OrdinalType, typename ScalarType>
SerialBandDenseMatrix<OrdinalType, ScalarType>::
SerialBandDenseMatrix (OrdinalType numRows_in,
                       OrdinalType numCols_in,
                       OrdinalType kl_in,
                       OrdinalType ku_in,
                       bool zeroOut)
  : CompObject (),
    BLAS<OrdinalType,ScalarType>(),
    numRows_ (numRows_in),
    numCols_ (numCols_in),
    stride_ (kl_in+ku_in+1),
    kl_ (kl_in),
    ku_ (ku_in),
    valuesCopied_ (true),
    values_ (NULL) // for safety, in case allocation fails below
{
  values_ = new ScalarType[stride_ * numCols_];
  if (zeroOut) {
    putScalar ();
  }
}

template<typename OrdinalType, typename ScalarType>
SerialBandDenseMatrix<OrdinalType, ScalarType>::
SerialBandDenseMatrix (DataAccess CV,
                       ScalarType* values_in,
                       OrdinalType stride_in,
                       OrdinalType numRows_in,
                       OrdinalType numCols_in,
                       OrdinalType kl_in,
                       OrdinalType ku_in)
  : CompObject (),
    BLAS<OrdinalType,ScalarType>(),
    numRows_ (numRows_in),
    numCols_ (numCols_in),
    stride_ (stride_in),
    kl_ (kl_in),
    ku_ (ku_in),
    valuesCopied_ (false),
    values_ (values_in)
{
  if (CV == Copy) {
    stride_ = kl_+ku_+1;
    values_ = new ScalarType[stride_*numCols_];
    copyMat (values_in, stride_in, numRows_, numCols_, values_, stride_, 0);
    valuesCopied_ = true;
  }
}

template<typename OrdinalType, typename ScalarType>
SerialBandDenseMatrix<OrdinalType, ScalarType>::
SerialBandDenseMatrix (const SerialBandDenseMatrix<OrdinalType, ScalarType> &Source, ETransp trans)
  : CompObject (),
    BLAS<OrdinalType,ScalarType>(),
    numRows_ (0),
    numCols_ (0),
    stride_ (0),
    kl_ (0),
    ku_ (0),
    valuesCopied_ (true),
    values_ (NULL)
{
  if (trans == NO_TRANS) {
    numRows_ = Source.numRows_;
    numCols_ = Source.numCols_;
    kl_ = Source.kl_;
    ku_ = Source.ku_;
    stride_ = kl_+ku_+1;
    values_ = new ScalarType[stride_*numCols_];
    copyMat (Source.values_, Source.stride_, numRows_, numCols_,
             values_, stride_, 0);
  }
  else if (trans == CONJ_TRANS && ScalarTraits<ScalarType>::isComplex) {
    numRows_ = Source.numCols_;
    numCols_ = Source.numRows_;
    kl_ = Source.ku_;
    ku_ = Source.kl_;
    stride_ = kl_+ku_+1;
    values_ = new ScalarType[stride_*numCols_];
    for (OrdinalType j = 0; j < numCols_; ++j) {
      for (OrdinalType i = TEUCHOS_MAX(0,j-ku_);
           i <= TEUCHOS_MIN(numRows_-1,j+kl_); ++i) {
        values_[j*stride_ + (ku_+i-j)] =
          ScalarTraits<ScalarType>::conjugate (Source.values_[i*Source.stride_ + (Source.ku_+j-i)]);
      }
    }
  }
  else {
    numRows_ = Source.numCols_;
    numCols_ = Source.numRows_;
    kl_ = Source.ku_;
    ku_ = Source.kl_;
    stride_ = kl_+ku_+1;
    values_ = new ScalarType[stride_*numCols_];
    for (OrdinalType j=0; j<numCols_; j++) {
      for (OrdinalType i = TEUCHOS_MAX(0,j-ku_);
           i <= TEUCHOS_MIN(numRows_-1,j+kl_); ++i) {
        values_[j*stride_ + (ku_+i-j)] = Source.values_[i*Source.stride_ + (Source.ku_+j-i)];
      }
    }
  }
}

template<typename OrdinalType, typename ScalarType>
SerialBandDenseMatrix<OrdinalType, ScalarType>::
SerialBandDenseMatrix (DataAccess CV,
                       const SerialBandDenseMatrix<OrdinalType, ScalarType> &Source,
                       OrdinalType numRows_in,
                       OrdinalType numCols_in,
                       OrdinalType startCol)
  : CompObject (),
    BLAS<OrdinalType,ScalarType>(),
    numRows_ (numRows_in),
    numCols_ (numCols_in),
    stride_ (Source.stride_),
    kl_ (Source.kl_),
    ku_ (Source.ku_),
    valuesCopied_ (false),
    values_ (Source.values_)
{
  if (CV == Copy) {
    values_ = new ScalarType[stride_ * numCols_in];
    copyMat (Source.values_, Source.stride_, numRows_in, numCols_in,
             values_, stride_, startCol);
    valuesCopied_ = true;
  } else { // CV = View
    values_ = values_ + (stride_ * startCol);
  }
}

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

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

template<typename OrdinalType, typename ScalarType>
int SerialBandDenseMatrix<OrdinalType, ScalarType>::shape(
  OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType kl_in, OrdinalType ku_in
  )
{
  deleteArrays(); // Get rid of anything that might be already allocated
  numRows_ = numRows_in;
  numCols_ = numCols_in;
  kl_ = kl_in;
  ku_ = ku_in;
  stride_ = kl_+ku_+1;
  values_ = new ScalarType[stride_*numCols_];
  putScalar();
  valuesCopied_ = true;

  return(0);
}

template<typename OrdinalType, typename ScalarType>
int SerialBandDenseMatrix<OrdinalType, ScalarType>::shapeUninitialized(
  OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType kl_in, OrdinalType ku_in
  )
{
  deleteArrays(); // Get rid of anything that might be already allocated
  numRows_ = numRows_in;
  numCols_ = numCols_in;
  kl_ = kl_in;
  ku_ = ku_in;
  stride_ = kl_+ku_+1;
  values_ = new ScalarType[stride_*numCols_];
  valuesCopied_ = true;

  return(0);
}

template<typename OrdinalType, typename ScalarType>
int SerialBandDenseMatrix<OrdinalType, ScalarType>::reshape(
  OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType kl_in, OrdinalType ku_in
  )
{

  // Allocate space for new matrix
  ScalarType* values_tmp = new ScalarType[(kl_in+ku_in+1) * numCols_in];
  ScalarType zero = ScalarTraits<ScalarType>::zero();
  for(OrdinalType k = 0; k < (kl_in+ku_in+1) * numCols_in; k++) {
    values_tmp[k] = zero;
  }
  OrdinalType numRows_tmp = TEUCHOS_MIN(numRows_, numRows_in);
  OrdinalType numCols_tmp = TEUCHOS_MIN(numCols_, numCols_in);
  if(values_ != 0) {
    copyMat(values_, stride_, numRows_tmp, numCols_tmp, values_tmp,
            kl_in+ku_in+1, 0); // Copy principal submatrix of A to new A
  }
  deleteArrays(); // Get rid of anything that might be already allocated
  numRows_ = numRows_in;
  numCols_ = numCols_in;
  kl_ = kl_in;
  ku_ = ku_in;
  stride_ = kl_+ku_+1;
  values_ = values_tmp; // Set pointer to new A
  valuesCopied_ = true;

  return(0);
}

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

template<typename OrdinalType, typename ScalarType>
int SerialBandDenseMatrix<OrdinalType, ScalarType>::putScalar( const ScalarType value_in )
{

  // Set each value of the dense matrix to "value".
  for(OrdinalType j = 0; j < numCols_; j++) {
    for (OrdinalType i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
      values_[(ku_+i-j) + j*stride_] = value_in;
    }
  }

  return 0;
}

template<typename OrdinalType, typename ScalarType>
int SerialBandDenseMatrix<OrdinalType, ScalarType>::random()
{

  // Set each value of the dense matrix to a random value.
  for(OrdinalType j = 0; j < numCols_; j++) {
    for (OrdinalType i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
      values_[(ku_+i-j) + j*stride_] = ScalarTraits<ScalarType>::random();
    }
  }

  return 0;
}

template<typename OrdinalType, typename ScalarType>
SerialBandDenseMatrix<OrdinalType,ScalarType>&
SerialBandDenseMatrix<OrdinalType, ScalarType>::operator=(
  const SerialBandDenseMatrix<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();
    }
    numRows_ = Source.numRows_;
    numCols_ = Source.numCols_;
    kl_ = Source.kl_;
    ku_ = Source.ku_;
    stride_ = Source.stride_;
    values_ = Source.values_;
  } else {
    // If we were a view, we will now be a copy.
    if(!valuesCopied_) {
      numRows_ = Source.numRows_;
      numCols_ = Source.numCols_;
      kl_ = Source.kl_;
      ku_ = Source.ku_;
      stride_ = kl_+ku_+1;
      const OrdinalType newsize = stride_ * numCols_;
      if(newsize > 0) {
        values_ = new ScalarType[newsize];
        valuesCopied_ = true;
      } else {
        values_ = 0;
      }
    } else {
      // If we were a copy, we will stay a copy.
      if((Source.numRows_ <= stride_) && (Source.numCols_ == numCols_)) { // we don't need to reallocate
        numRows_ = Source.numRows_;
        numCols_ = Source.numCols_;
        kl_ = Source.kl_;
        ku_ = Source.ku_;
      } else {
        // we need to allocate more space (or less space)
        deleteArrays();
        numRows_ = Source.numRows_;
        numCols_ = Source.numCols_;
        kl_ = Source.kl_;
        ku_ = Source.ku_;
        stride_ = kl_+ku_+1;
        const OrdinalType newsize = stride_ * numCols_;
        if(newsize > 0) {
          values_ = new ScalarType[newsize];
          valuesCopied_ = true;
        }
      }
    }
    copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0);
  }
  return(*this);

}

template<typename OrdinalType, typename ScalarType>
SerialBandDenseMatrix<OrdinalType, ScalarType>& SerialBandDenseMatrix<OrdinalType, ScalarType>::operator+= (const SerialBandDenseMatrix<OrdinalType,ScalarType>& Source )
{

  // Check for compatible dimensions
  if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_) || (kl_ != Source.kl_) || (ku_ != Source.ku_)) {
    TEUCHOS_CHK_REF(*this); // Return *this without altering it.
  }
  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, ScalarTraits<ScalarType>::one());
  return(*this);

}

template<typename OrdinalType, typename ScalarType>
SerialBandDenseMatrix<OrdinalType, ScalarType>& SerialBandDenseMatrix<OrdinalType, ScalarType>::operator-= (const SerialBandDenseMatrix<OrdinalType,ScalarType>& Source )
{

  // Check for compatible dimensions
  if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_) || (kl_ != Source.kl_) || (ku_ != Source.ku_)) {
    TEUCHOS_CHK_REF(*this); // Return *this without altering it.
  }
  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, -ScalarTraits<ScalarType>::one());
  return(*this);

}

template<typename OrdinalType, typename ScalarType>
SerialBandDenseMatrix<OrdinalType,ScalarType>& SerialBandDenseMatrix<OrdinalType, ScalarType>::assign (const SerialBandDenseMatrix<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 ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_) || (kl_ != Source.kl_) || (ku_ != Source.ku_)) {
    TEUCHOS_CHK_REF(*this); // Return *this without altering it.
  }
  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0);
  return(*this);

}

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

template<typename OrdinalType, typename ScalarType>
inline ScalarType& SerialBandDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex)
{
#ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
  checkIndex( rowIndex, colIndex );
#endif
  return(values_[colIndex * stride_ + ku_+rowIndex-colIndex]);
}

template<typename OrdinalType, typename ScalarType>
inline const ScalarType& SerialBandDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex) const
{
#ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
  checkIndex( rowIndex, colIndex );
#endif
  return(values_[colIndex * stride_ + ku_+rowIndex-colIndex]);
}

template<typename OrdinalType, typename ScalarType>
inline const ScalarType* SerialBandDenseMatrix<OrdinalType, ScalarType>::operator [] (OrdinalType colIndex) const
{
#ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
  checkIndex( 0, colIndex );
#endif
  return(values_ + colIndex * stride_);
}

template<typename OrdinalType, typename ScalarType>
inline ScalarType* SerialBandDenseMatrix<OrdinalType, ScalarType>::operator [] (OrdinalType colIndex)
{
#ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
  checkIndex( 0, colIndex );
#endif
  return(values_ + colIndex * stride_);
}

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

template<typename OrdinalType, typename ScalarType>
typename ScalarTraits<ScalarType>::magnitudeType SerialBandDenseMatrix<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());

  ScalarType* ptr;
  for(j = 0; j < numCols_; j++) {
    ScalarType sum = 0;
    ptr = values_ + j * stride_ + TEUCHOS_MAX(0, ku_-j);
    for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
      sum += ScalarTraits<ScalarType>::magnitude(*ptr++);
    }
    absSum = ScalarTraits<ScalarType>::magnitude(sum);
    if(absSum > anorm) {
      anorm = absSum;
    }
  }
  updateFlops((kl_+ku_+1) * numCols_);

  return(anorm);
}

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

  for (i = 0; i < numRows_; i++) {
    sum = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
    for (j=TEUCHOS_MAX(0,i-kl_); j<=TEUCHOS_MIN(numCols_-1,i+ku_); j++) {
      sum += ScalarTraits<ScalarType>::magnitude(*(values_+(ku_+i-j)+j*stride_));
    }
    anorm = TEUCHOS_MAX( anorm, sum );
  }
  updateFlops((kl_+ku_+1) * numCols_);

  return(anorm);
}

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

  for (j = 0; j < numCols_; j++) {
    for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
      anorm += ScalarTraits<ScalarType>::magnitude(values_[(ku_+i-j)+j*stride_]*values_[(ku_+i-j)+j*stride_]);
    }
  }
  anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::squareroot(anorm));
  updateFlops((kl_+ku_+1) * numCols_);

  return(anorm);
}

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

template<typename OrdinalType, typename ScalarType>
bool SerialBandDenseMatrix<OrdinalType, ScalarType>::operator== (const SerialBandDenseMatrix<OrdinalType, ScalarType> &Operand) const
{
  bool result = 1;

  if((numRows_ != Operand.numRows_) || (numCols_ != Operand.numCols_) || (kl_ != Operand.kl_) || (ku_ != Operand.ku_)) {
    result = 0;
  } else {
    OrdinalType i, j;
    for(j = 0; j < numCols_; j++) {
      for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
        if((*this)(i, j) != Operand(i, j)) {
          return 0;
        }
      }
    }
  }

  return result;
}

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

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

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

template<typename OrdinalType, typename ScalarType>
int SerialBandDenseMatrix<OrdinalType, ScalarType>::scale( const ScalarType alpha )
{

  OrdinalType i, j;
  ScalarType* ptr;

  for (j=0; j<numCols_; j++) {
    ptr = values_ + j*stride_ + TEUCHOS_MAX(0, ku_-j);
    for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
      *ptr = alpha * (*ptr); ptr++;
    }
  }
  updateFlops( (kl_+ku_+1)*numCols_ );

  return(0);
}

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

  OrdinalType i, j;
  ScalarType* ptr;

  // Check for compatible dimensions
  if ((numRows_ != A.numRows_) || (numCols_ != A.numCols_) || (kl_ != A.kl_) || (ku_ != A.ku_)) {
    TEUCHOS_CHK_ERR(-1); // Return error
  }
  for (j=0; j<numCols_; j++) {
    ptr = values_ + j*stride_ + TEUCHOS_MAX(0, ku_-j);
    for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
      *ptr = A(i,j) * (*ptr); ptr++;
    }
  }
  updateFlops( (kl_+ku_+1)*numCols_ );

  return(0);
}

template<typename OrdinalType, typename ScalarType>
std::ostream& SerialBandDenseMatrix<OrdinalType, ScalarType>::print(std::ostream& os) const
{
  os << std::endl;
  if(valuesCopied_)
    os << "Values_copied : yes" << std::endl;
  else
    os << "Values_copied : no" << std::endl;
  os << "Rows : " << numRows_ << std::endl;
  os << "Columns : " << numCols_ << std::endl;
  os << "Lower Bandwidth : " << kl_ << std::endl;
  os << "Upper Bandwidth : " << ku_ << std::endl;
  os << "LDA : " << stride_ << std::endl;
  if(numRows_ == 0 || numCols_ == 0) {
    os << "(matrix is empty, no values to display)" << std::endl;
  } else {

    for(OrdinalType i = 0; i < numRows_; i++) {
      for (OrdinalType j=TEUCHOS_MAX(0,i-kl_); j<=TEUCHOS_MIN(numCols_-1,i+ku_); j++) {
        os << (*this)(i,j) << " ";
      }
      os << std::endl;
    }
  }
  return os;
}

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

template<typename OrdinalType, typename ScalarType>
inline void SerialBandDenseMatrix<OrdinalType, ScalarType>::checkIndex( OrdinalType rowIndex, OrdinalType colIndex ) const {

  TEUCHOS_TEST_FOR_EXCEPTION(rowIndex < 0 || rowIndex >= numRows_ ||
                             rowIndex < TEUCHOS_MAX(0,colIndex-ku_) || rowIndex > TEUCHOS_MIN(numRows_-1,colIndex+kl_),
                             std::out_of_range,
                             "SerialBandDenseMatrix<T>::checkIndex: "
                             "Row index " << rowIndex << " out of range [0, "<< numRows_ << ")");
  TEUCHOS_TEST_FOR_EXCEPTION(colIndex < 0 || colIndex >= numCols_, std::out_of_range,
                             "SerialBandDenseMatrix<T>::checkIndex: "
                             "Col index " << colIndex << " out of range [0, "<< numCols_ << ")");

}

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

template<typename OrdinalType, typename ScalarType>
void SerialBandDenseMatrix<OrdinalType, ScalarType>::copyMat(
  ScalarType* inputMatrix, OrdinalType strideInput, OrdinalType numRows_in,
  OrdinalType numCols_in, ScalarType* outputMatrix, OrdinalType strideOutput, OrdinalType startCol, ScalarType alpha
  )
{
  OrdinalType i, j;
  ScalarType* ptr1 = 0;
  ScalarType* ptr2 = 0;

  for(j = 0; j < numCols_in; j++) {
    ptr1 = outputMatrix + (j * strideOutput) + TEUCHOS_MAX(0, ku_-j);
    ptr2 = inputMatrix + (j + startCol) * strideInput + TEUCHOS_MAX(0, ku_-j);
    if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
      for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_in-1,j+kl_); i++) {
        *ptr1++ += alpha*(*ptr2++);
      }
    } else {
      for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_in-1,j+kl_); i++) {
        *ptr1++ = *ptr2++;
      }
    }
  }
}

template<typename OrdinalType, typename ScalarType> 
struct SerialBandDenseMatrixPrinter {
public:
  const SerialBandDenseMatrix<OrdinalType,ScalarType> &obj; 
  SerialBandDenseMatrixPrinter(
        const SerialBandDenseMatrix<OrdinalType,ScalarType> &obj_in)
      : obj(obj_in) {}
};

template<typename OrdinalType, typename ScalarType> 
std::ostream&
operator<<(std::ostream &out, 
           const SerialBandDenseMatrixPrinter<OrdinalType,ScalarType> printer)
{
  printer.obj.print(out);
  return out;
}

template<typename OrdinalType, typename ScalarType> 
SerialBandDenseMatrixPrinter<OrdinalType,ScalarType>
printMat(const SerialBandDenseMatrix<OrdinalType,ScalarType> &obj) 
{ 
  return SerialBandDenseMatrixPrinter<OrdinalType,ScalarType>(obj); 
}


} // namespace Teuchos


#endif /* _TEUCHOS_SERIALBANDDENSEMATRIX_HPP_ */