Teuchos_SerialSymDenseMatrix.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.
//
// ***********************************************************************
// @HEADER

#ifndef _TEUCHOS_SERIALSYMDENSEMATRIX_HPP_
#define _TEUCHOS_SERIALSYMDENSEMATRIX_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"
#include <utility>
#include <vector>

namespace Teuchos {

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

  typedef OrdinalType ordinalType;
  typedef ScalarType scalarType;



  SerialSymDenseMatrix() = default;


  SerialSymDenseMatrix(OrdinalType numRowsCols, bool zeroOut = true);


  SerialSymDenseMatrix(DataAccess CV, bool upper, ScalarType* values, OrdinalType stride, OrdinalType numRowsCols);

  SerialSymDenseMatrix(const SerialSymDenseMatrix<OrdinalType, ScalarType> &Source);


  SerialSymDenseMatrix(DataAccess CV, const SerialSymDenseMatrix<OrdinalType, ScalarType> &Source, OrdinalType numRowCols, OrdinalType startRowCol=0);

  virtual ~SerialSymDenseMatrix ();




  int shape(OrdinalType numRowsCols);


  int shapeUninitialized(OrdinalType numRowsCols);


  int reshape(OrdinalType numRowsCols);


  void setLower();


  void setUpper();





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


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


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


  int putScalar( const ScalarType value = Teuchos::ScalarTraits<ScalarType>::zero(), bool fullMatrix = false );


  void swap (SerialSymDenseMatrix<OrdinalType, ScalarType> &B);


  int random( const ScalarType bias = 0.1*Teuchos::ScalarTraits<ScalarType>::one() );





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


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


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




  bool upper() const {return(upper_);};

  char UPLO() const {return(UPLO_);};




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


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


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





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


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




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

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

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

  bool empty() const { return(numRowCols_ == 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:

  // In-place scaling of matrix.
  void scale( const ScalarType alpha );

  // Copy the values from one matrix to the other.
  void copyMat(bool inputUpper, ScalarType* inputMatrix, OrdinalType inputStride,
               OrdinalType numRowCols, bool outputUpper, ScalarType* outputMatrix,
               OrdinalType outputStride, OrdinalType startRowCol,
               ScalarType alpha = ScalarTraits<ScalarType>::zero() );

  // Copy the values from the active triangle of the matrix to the other to make the matrix full symmetric.
  void copyUPLOMat(bool inputUpper, ScalarType* inputMatrix,
                   OrdinalType inputStride, OrdinalType inputRows);

  void deleteArrays();
  void checkIndex( OrdinalType rowIndex, OrdinalType colIndex = 0 ) const;

  static ScalarType*
  allocateValues(const OrdinalType numRows,
                 const OrdinalType numCols)
  {
    const size_t size = size_t(numRows) * size_t(numCols);
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wvla"
    return new ScalarType[size];
#pragma GCC diagnostic pop
  }

  OrdinalType numRowCols_ = 0;
  OrdinalType stride_ = 0;
  bool valuesCopied_ = false;
  ScalarType* values_ = nullptr;
  bool upper_ = false;
  char UPLO_ {'L'};
};

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

template<typename OrdinalType, typename ScalarType>
SerialSymDenseMatrix<OrdinalType, ScalarType>::SerialSymDenseMatrix(OrdinalType numRowCols_in, bool zeroOut)
  : numRowCols_(numRowCols_in), stride_(numRowCols_in), valuesCopied_(false), upper_(false), UPLO_('L')
{
  values_ = allocateValues(stride_, numRowCols_);
  valuesCopied_ = true;
  if (zeroOut) {
    const ScalarType ZERO = Teuchos::ScalarTraits<ScalarType>::zero();
    putScalar(ZERO, true);
  }
}

template<typename OrdinalType, typename ScalarType>
SerialSymDenseMatrix<OrdinalType, ScalarType>::SerialSymDenseMatrix(
  DataAccess CV, bool upper_in, ScalarType* values_in, OrdinalType stride_in, OrdinalType numRowCols_in
  )
  : numRowCols_(numRowCols_in), stride_(stride_in), valuesCopied_(false),
    values_(values_in), upper_(upper_in)
{
  if (upper_)
    UPLO_ = 'U';
  else
    UPLO_ = 'L';

  if(CV == Copy)
  {
    stride_ = numRowCols_;
    values_ = allocateValues(stride_, numRowCols_);
    copyMat(upper_in, values_in, stride_in, numRowCols_, upper_, values_, stride_, 0);
    valuesCopied_ = true;
  }
}

template<typename OrdinalType, typename ScalarType>
SerialSymDenseMatrix<OrdinalType, ScalarType>::SerialSymDenseMatrix(const SerialSymDenseMatrix<OrdinalType, ScalarType> &Source)
  : numRowCols_(Source.numRowCols_), stride_(0), valuesCopied_(true),
    values_(0), upper_(Source.upper_), UPLO_(Source.UPLO_)
{
  if (!Source.valuesCopied_)
  {
    stride_ = Source.stride_;
    values_ = Source.values_;
    valuesCopied_ = false;
  }
  else
  {
    stride_ = numRowCols_;
    if(stride_ > 0 && numRowCols_ > 0) {
      values_ = allocateValues(stride_, numRowCols_);
      copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0);
    }
    else {
      numRowCols_ = 0;
      stride_ = 0;
      valuesCopied_ = false;
    }
  }
}

template<typename OrdinalType, typename ScalarType>
SerialSymDenseMatrix<OrdinalType, ScalarType>::SerialSymDenseMatrix(
                                                                    DataAccess CV, const SerialSymDenseMatrix<OrdinalType,
                                                                    ScalarType> &Source, OrdinalType numRowCols_in, OrdinalType startRowCol )
  : CompObject(), BLAS<OrdinalType,ScalarType>(),
    numRowCols_(numRowCols_in), stride_(Source.stride_), valuesCopied_(false), upper_(Source.upper_), UPLO_(Source.UPLO_)
{
  if(CV == Copy)
  {
    stride_ = numRowCols_in;
    values_ = allocateValues(stride_, numRowCols_in);
    copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_in, upper_, values_, stride_, startRowCol);
    valuesCopied_ = true;
  }
  else // CV == View
  {
    values_ = Source.values_ + (stride_ * startRowCol) + startRowCol;
  }
}

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

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

template<typename OrdinalType, typename ScalarType>
int SerialSymDenseMatrix<OrdinalType, ScalarType>::shape( OrdinalType numRowCols_in )
{
  deleteArrays(); // Get rid of anything that might be already allocated
  numRowCols_ = numRowCols_in;
  stride_ = numRowCols_;
  values_ = allocateValues(stride_, numRowCols_);
  putScalar( Teuchos::ScalarTraits<ScalarType>::zero(), true );
  valuesCopied_ = true;
  return(0);
}

template<typename OrdinalType, typename ScalarType>
int SerialSymDenseMatrix<OrdinalType, ScalarType>::shapeUninitialized( OrdinalType numRowCols_in )
{
  deleteArrays(); // Get rid of anything that might be already allocated
  numRowCols_ = numRowCols_in;
  stride_ = numRowCols_;
  values_ = allocateValues(stride_, numRowCols_);
  valuesCopied_ = true;
  return(0);
}

template<typename OrdinalType, typename ScalarType>
int SerialSymDenseMatrix<OrdinalType, ScalarType>::reshape( OrdinalType numRowCols_in )
{
  // Allocate space for new matrix
  ScalarType* values_tmp = allocateValues(numRowCols_in, numRowCols_in);
  const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
  for(OrdinalType k = 0; k < numRowCols_in * numRowCols_in; k++)
  {
    values_tmp[k] = zero;
  }
  OrdinalType numRowCols_tmp = TEUCHOS_MIN(numRowCols_, numRowCols_in);
  if(values_ != 0)
  {
    copyMat(upper_, values_, stride_, numRowCols_tmp, upper_, values_tmp, numRowCols_in, 0); // Copy principal submatrix of A to new A
  }
  deleteArrays(); // Get rid of anything that might be already allocated
  numRowCols_ = numRowCols_in;
  stride_ = numRowCols_;
  values_ = values_tmp; // Set pointer to new A
  valuesCopied_ = true;
  return(0);
}

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

template<typename OrdinalType, typename ScalarType>
void SerialSymDenseMatrix<OrdinalType, ScalarType>::setLower()
{
  // Do nothing if the matrix is already a lower triangular matrix
  if (upper_ != false) {
    copyUPLOMat( true, values_, stride_, numRowCols_ );
    upper_ = false;
    UPLO_ = 'L';
  }
}

template<typename OrdinalType, typename ScalarType>
void SerialSymDenseMatrix<OrdinalType, ScalarType>::setUpper()
{
  // Do nothing if the matrix is already an upper triangular matrix
  if (upper_ == false) {
    copyUPLOMat( false, values_, stride_, numRowCols_ );
    upper_ = true;
    UPLO_ = 'U';
  }
}

template<typename OrdinalType, typename ScalarType>
int SerialSymDenseMatrix<OrdinalType, ScalarType>::putScalar( const ScalarType value_in, bool fullMatrix )
{
  // Set each value of the dense matrix to "value".
  if (fullMatrix == true) {
    for(OrdinalType j = 0; j < numRowCols_; j++)
      {
        for(OrdinalType i = 0; i < numRowCols_; i++)
          {
            values_[i + j*stride_] = value_in;
          }
      }
  }
  // Set the active upper or lower triangular part of the matrix to "value"
  else {
    if (upper_) {
      for(OrdinalType j = 0; j < numRowCols_; j++) {
        for(OrdinalType i = 0; i <= j; i++) {
          values_[i + j*stride_] = value_in;
        }
      }
    }
    else {
      for(OrdinalType j = 0; j < numRowCols_; j++) {
        for(OrdinalType i = j; i < numRowCols_; i++) {
          values_[i + j*stride_] = value_in;
        }
      }
    }
  }
  return 0;
}

template<typename OrdinalType, typename ScalarType> void
SerialSymDenseMatrix<OrdinalType, ScalarType>::swap(
  SerialSymDenseMatrix<OrdinalType, ScalarType> &B)
{
  std::swap(values_ ,      B.values_);
  std::swap(numRowCols_,   B.numRowCols_);
  std::swap(stride_,       B.stride_);
  std::swap(valuesCopied_, B.valuesCopied_);
  std::swap(upper_,        B.upper_);
  std::swap(UPLO_,         B.UPLO_);
}

template<typename OrdinalType, typename ScalarType>
int SerialSymDenseMatrix<OrdinalType, ScalarType>::random( const ScalarType bias )
{
  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MT;

  // Set each value of the dense matrix to a random value.
  std::vector<MT> diagSum( numRowCols_, Teuchos::ScalarTraits<MT>::zero() );
  if (upper_) {
    for(OrdinalType j = 0; j < numRowCols_; j++) {
      for(OrdinalType i = 0; i < j; i++) {
        values_[i + j*stride_] = ScalarTraits<ScalarType>::random();
        diagSum[i] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
        diagSum[j] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
      }
    }
  }
  else {
    for(OrdinalType j = 0; j < numRowCols_; j++) {
      for(OrdinalType i = j+1; i < numRowCols_; i++) {
        values_[i + j*stride_] = ScalarTraits<ScalarType>::random();
        diagSum[i] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
        diagSum[j] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
      }
    }
  }

  // Set the diagonal.
  for(OrdinalType i = 0; i < numRowCols_; i++) {
    values_[i + i*stride_] = diagSum[i] + bias;
  }
  return 0;
}

template<typename OrdinalType, typename ScalarType>
SerialSymDenseMatrix<OrdinalType,ScalarType>&
SerialSymDenseMatrix<OrdinalType, ScalarType>::operator=( const SerialSymDenseMatrix<OrdinalType,ScalarType>& Source )
{
  if(this == &Source)
    return(*this); // Special case of source same as target
  if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_)) {
    upper_ = Source.upper_;  // Might have to change the active part of the matrix.
    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();
    }
    numRowCols_ = Source.numRowCols_;
    stride_ = Source.stride_;
    values_ = Source.values_;
    upper_ = Source.upper_;
    UPLO_ = Source.UPLO_;
  }
  else {
    // If we were a view, we will now be a copy.
    if(!valuesCopied_) {
      numRowCols_ = Source.numRowCols_;
      stride_ = Source.numRowCols_;
      upper_ = Source.upper_;
      UPLO_ = Source.UPLO_;
      if(stride_ > 0 && numRowCols_ > 0) {
        values_ = allocateValues(stride_, numRowCols_);
        valuesCopied_ = true;
      }
      else {
        values_ = 0;
      }
    }
    // If we were a copy, we will stay a copy.
    else {
      if((Source.numRowCols_ <= stride_) && (Source.numRowCols_ == numRowCols_)) { // we don't need to reallocate
        numRowCols_ = Source.numRowCols_;
        upper_ = Source.upper_;
        UPLO_ = Source.UPLO_;
      }
      else { // we need to allocate more space (or less space)
        deleteArrays();
        numRowCols_ = Source.numRowCols_;
        stride_ = Source.numRowCols_;
        upper_ = Source.upper_;
        UPLO_ = Source.UPLO_;
        if(stride_ > 0 && numRowCols_ > 0) {
          values_ = allocateValues(stride_, numRowCols_);
          valuesCopied_ = true;
        }
      }
    }
    copyMat(Source.upper_, Source.values_, Source.stride_, Source.numRowCols_, upper_, values_, stride_, 0);
  }
  return(*this);
}

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

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

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

template<typename OrdinalType, typename ScalarType>
SerialSymDenseMatrix<OrdinalType,ScalarType>& SerialSymDenseMatrix<OrdinalType, ScalarType>::assign (const SerialSymDenseMatrix<OrdinalType,ScalarType>& Source) {
  if(this == &Source)
    return(*this); // Special case of source same as target
  if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_)) {
    upper_ = Source.upper_; // We may have to change the active part of the matrix.
    return(*this); // Special case of both are views to same data.
  }

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

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

template<typename OrdinalType, typename ScalarType>
inline ScalarType& SerialSymDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex)
{
#ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
  checkIndex( rowIndex, colIndex );
#endif
  if ( rowIndex <= colIndex ) {
    // Accessing upper triangular part of matrix
    if (upper_)
      return(values_[colIndex * stride_ + rowIndex]);
    else
      return(values_[rowIndex * stride_ + colIndex]);
  }
  else {
    // Accessing lower triangular part of matrix
    if (upper_)
      return(values_[rowIndex * stride_ + colIndex]);
    else
      return(values_[colIndex * stride_ + rowIndex]);
  }
}

template<typename OrdinalType, typename ScalarType>
inline const ScalarType& SerialSymDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex) const
{
#ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
  checkIndex( rowIndex, colIndex );
#endif
  if ( rowIndex <= colIndex ) {
    // Accessing upper triangular part of matrix
    if (upper_)
      return(values_[colIndex * stride_ + rowIndex]);
    else
      return(values_[rowIndex * stride_ + colIndex]);
  }
  else {
    // Accessing lower triangular part of matrix
    if (upper_)
      return(values_[rowIndex * stride_ + colIndex]);
    else
      return(values_[colIndex * stride_ + rowIndex]);
  }
}

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

template<typename OrdinalType, typename ScalarType>
typename ScalarTraits<ScalarType>::magnitudeType SerialSymDenseMatrix<OrdinalType, ScalarType>::normOne() const
{
  return(normInf());
}

template<typename OrdinalType, typename ScalarType>
typename ScalarTraits<ScalarType>::magnitudeType SerialSymDenseMatrix<OrdinalType, ScalarType>::normInf() const
{
  typedef typename ScalarTraits<ScalarType>::magnitudeType MT;

  OrdinalType i, j;

  MT sum, anorm = ScalarTraits<MT>::zero();
  ScalarType* ptr;

  if (upper_) {
    for (j=0; j<numRowCols_; j++) {
      sum = ScalarTraits<MT>::zero();
      ptr = values_ + j*stride_;
      for (i=0; i<j; i++) {
        sum += ScalarTraits<ScalarType>::magnitude( *ptr++ );
      }
      ptr = values_ + j + j*stride_;
      for (i=j; i<numRowCols_; i++) {
        sum += ScalarTraits<ScalarType>::magnitude( *ptr );
        ptr += stride_;
      }
      anorm = TEUCHOS_MAX( anorm, sum );
    }
  }
  else {
    for (j=0; j<numRowCols_; j++) {
      sum = ScalarTraits<MT>::zero();
      ptr = values_ + j + j*stride_;
      for (i=j; i<numRowCols_; i++) {
        sum += ScalarTraits<ScalarType>::magnitude( *ptr++ );
      }
      ptr = values_ + j;
      for (i=0; i<j; i++) {
        sum += ScalarTraits<ScalarType>::magnitude( *ptr );
        ptr += stride_;
      }
      anorm = TEUCHOS_MAX( anorm, sum );
    }
  }
  return(anorm);
}

template<typename OrdinalType, typename ScalarType>
typename ScalarTraits<ScalarType>::magnitudeType SerialSymDenseMatrix<OrdinalType, ScalarType>::normFrobenius() const
{
  typedef typename ScalarTraits<ScalarType>::magnitudeType MT;

  OrdinalType i, j;

  MT sum = ScalarTraits<MT>::zero(), anorm = ScalarTraits<MT>::zero();

  if (upper_) {
    for (j = 0; j < numRowCols_; j++) {
      for (i = 0; i < j; i++) {
        sum += ScalarTraits<ScalarType>::magnitude(2.0*values_[i+j*stride_]*values_[i+j*stride_]);
      }
      sum += ScalarTraits<ScalarType>::magnitude(values_[j + j*stride_]*values_[j + j*stride_]);
    }
  }
  else {
    for (j = 0; j < numRowCols_; j++) {
      sum += ScalarTraits<ScalarType>::magnitude(values_[j + j*stride_]*values_[j + j*stride_]);
      for (i = j+1; i < numRowCols_; i++) {
        sum += ScalarTraits<ScalarType>::magnitude(2.0*values_[i+j*stride_]*values_[i+j*stride_]);
      }
    }
  }
  anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::squareroot(sum));
  return(anorm);
}

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

template<typename OrdinalType, typename ScalarType>
bool SerialSymDenseMatrix<OrdinalType, ScalarType>::operator== (const SerialSymDenseMatrix<OrdinalType, ScalarType> &Operand) const
{
  bool result = 1;
  if((numRowCols_ != Operand.numRowCols_)) {
    result = 0;
  }
  else {
    OrdinalType i, j;
    for(i = 0; i < numRowCols_; i++) {
      for(j = 0; j < numRowCols_; j++) {
        if((*this)(i, j) != Operand(i, j)) {
          return 0;
        }
      }
    }
  }
  return result;
}

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

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

template<typename OrdinalType, typename ScalarType>
void SerialSymDenseMatrix<OrdinalType, ScalarType>::scale( const ScalarType alpha )
{
  OrdinalType i, j;
  ScalarType* ptr;

  if (upper_) {
    for (j=0; j<numRowCols_; j++) {
      ptr = values_ + j*stride_;
      for (i=0; i<= j; i++) { *ptr = alpha * (*ptr); ptr++; }
    }
  }
  else {
    for (j=0; j<numRowCols_; j++) {
      ptr = values_ + j*stride_ + j;
      for (i=j; i<numRowCols_; i++) { *ptr = alpha * (*ptr); ptr++; }
    }
  }
}

/*
template<typename OrdinalType, typename ScalarType>
int SerialSymDenseMatrix<OrdinalType, ScalarType>::scale( const SerialSymDenseMatrix<OrdinalType,ScalarType>& A )
{
  OrdinalType i, j;
  ScalarType* ptr;

  // Check for compatible dimensions
  if ((numRowCols_ != A.numRowCols_)) {
    TEUCHOS_CHK_ERR(-1); // Return error
  }

  if (upper_) {
    for (j=0; j<numRowCols_; j++) {
      ptr = values_ + j*stride_;
      for (i=0; i<= j; i++) { *ptr = A(i,j) * (*ptr); ptr++; }
    }
  }
  else {
    for (j=0; j<numRowCols_; j++) {
      ptr = values_ + j*stride_;
      for (i=j; i<numRowCols_; i++) { *ptr = A(i,j) * (*ptr); ptr++; }
    }
  }

  return(0);
}
*/

template<typename OrdinalType, typename ScalarType>
std::ostream& SerialSymDenseMatrix<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 / Columns : " << numRowCols_ << std::endl;
  os << "LDA : " << stride_ << std::endl;
  if (upper_)
    os << "Storage: Upper" << std::endl;
  else
    os << "Storage: Lower" << std::endl;
  if(numRowCols_ == 0) {
    os << "(matrix is empty, no values to display)" << std::endl;
  } else {
    for(OrdinalType i = 0; i < numRowCols_; i++) {
      for(OrdinalType j = 0; j < numRowCols_; j++){
        os << (*this)(i,j) << " ";
      }
      os << std::endl;
    }
  }
  return os;
}

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

template<typename OrdinalType, typename ScalarType>
inline void SerialSymDenseMatrix<OrdinalType, ScalarType>::checkIndex( OrdinalType rowIndex, OrdinalType colIndex ) const {
  TEUCHOS_TEST_FOR_EXCEPTION(rowIndex < 0 || rowIndex >= numRowCols_, std::out_of_range,
    "SerialSymDenseMatrix<T>::checkIndex: "
    "Row index " << rowIndex << " out of range [0, "<< numRowCols_ << ")");
  TEUCHOS_TEST_FOR_EXCEPTION(colIndex < 0 || colIndex >= numRowCols_, std::out_of_range,
    "SerialSymDenseMatrix<T>::checkIndex: "
    "Col index " << colIndex << " out of range [0, "<< numRowCols_ << ")");
}

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

template<typename OrdinalType, typename ScalarType>
void SerialSymDenseMatrix<OrdinalType, ScalarType>::copyMat(
                                                            bool inputUpper, ScalarType* inputMatrix,
                                                            OrdinalType inputStride, OrdinalType numRowCols_in,
                                                            bool outputUpper, ScalarType* outputMatrix,
                                                            OrdinalType outputStride, OrdinalType startRowCol,
                                                            ScalarType alpha
                                                            )
{
  OrdinalType i, j;
  ScalarType* ptr1 = 0;
  ScalarType* ptr2 = 0;

  for (j = 0; j < numRowCols_in; j++) {
    if (inputUpper == true) {
      // The input matrix is upper triangular, start at the beginning of each column.
      ptr2 = inputMatrix + (j + startRowCol) * inputStride + startRowCol;
      if (outputUpper == true) {
        // The output matrix matches the same pattern as the input matrix.
        ptr1 = outputMatrix + j*outputStride;
        if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
          for(i = 0; i <= j; i++) {
            *ptr1++ += alpha*(*ptr2++);
          }
        } else {
          for(i = 0; i <= j; i++) {
            *ptr1++ = *ptr2++;
          }
        }
      }
      else {
        // The output matrix has the opposite pattern as the input matrix.
        // Copy down across rows of the output matrix, but down columns of the input matrix.
        ptr1 = outputMatrix + j;
        if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
          for(i = 0; i <= j; i++) {
            *ptr1 += alpha*(*ptr2++);
            ptr1 += outputStride;
          }
        } else {
          for(i = 0; i <= j; i++) {
            *ptr1 = *ptr2++;
            ptr1 += outputStride;
          }
        }
      }
    }
    else {
      // The input matrix is lower triangular, start at the diagonal of each row.
      ptr2 = inputMatrix + (startRowCol+j) * inputStride + startRowCol + j;
      if (outputUpper == true) {
        // The output matrix has the opposite pattern as the input matrix.
        // Copy across rows of the output matrix, but down columns of the input matrix.
        ptr1 = outputMatrix + j*outputStride + j;
        if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
          for(i = j; i < numRowCols_in; i++) {
            *ptr1 += alpha*(*ptr2++);
            ptr1 += outputStride;
          }
        } else {
          for(i = j; i < numRowCols_in; i++) {
            *ptr1 = *ptr2++;
            ptr1 += outputStride;
          }
        }
      }
      else {
        // The output matrix matches the same pattern as the input matrix.
        ptr1 = outputMatrix + j*outputStride + j;
        if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
          for(i = j; i < numRowCols_in; i++) {
            *ptr1++ += alpha*(*ptr2++);
          }
        } else {
          for(i = j; i < numRowCols_in; i++) {
            *ptr1++ = *ptr2++;
          }
        }
      }
    }
  }
}

template<typename OrdinalType, typename ScalarType>
void SerialSymDenseMatrix<OrdinalType, ScalarType>::copyUPLOMat(
                                                                bool inputUpper, ScalarType* inputMatrix,
                                                                OrdinalType inputStride, OrdinalType inputRows
                                                                )
{
  OrdinalType i, j;
  ScalarType * ptr1 = 0;
  ScalarType * ptr2 = 0;

  if (inputUpper) {
    for (j=1; j<inputRows; j++) {
      ptr1 = inputMatrix + j;
      ptr2 = inputMatrix + j*inputStride;
      for (i=0; i<j; i++) {
        *ptr1 = *ptr2++;
        ptr1+=inputStride;
      }
    }
  }
  else {
    for (i=1; i<inputRows; i++) {
      ptr1 = inputMatrix + i;
      ptr2 = inputMatrix + i*inputStride;
      for (j=0; j<i; j++) {
        *ptr2++ = *ptr1;
        ptr1+=inputStride;
      }
    }
  }
}

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

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

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


} // namespace Teuchos

#endif /* _TEUCHOS_SERIALSYMDENSEMATRIX_HPP_ */