Teuchos_SerialSymDenseMatrix.hpp

Go to the documentation of this file.
00001 
00002 //@HEADER
00003 // ***********************************************************************
00004 // 
00005 //                    Teuchos: Common Tools Package
00006 //                 Copyright (2004) Sandia Corporation
00007 // 
00008 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00009 // license for use of this work by or on behalf of the U.S. Government.
00010 // 
00011 // This library is free software; you can redistribute it and/or modify
00012 // it under the terms of the GNU Lesser General Public License as
00013 // published by the Free Software Foundation; either version 2.1 of the
00014 // License, or (at your option) any later version.
00015 //  
00016 // This library is distributed in the hope that it will be useful, but
00017 // WITHOUT ANY WARRANTY; without even the implied warranty of
00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00019 // Lesser General Public License for more details.
00020 //  
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License along with this library; if not, write to the Free Software
00023 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00024 // USA
00025 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00026 // 
00027 // ***********************************************************************
00028 // @HEADER
00029 
00030 #ifndef _TEUCHOS_SERIALSYMDENSEMATRIX_HPP_
00031 #define _TEUCHOS_SERIALSYMDENSEMATRIX_HPP_
00032 
00036 #include "Teuchos_CompObject.hpp"
00037 #include "Teuchos_BLAS.hpp"
00038 #include "Teuchos_ScalarTraits.hpp"
00039 #include "Teuchos_DataAccess.hpp"
00040 #include "Teuchos_ConfigDefs.hpp"
00041 #include "Teuchos_TestForException.hpp"
00042 
00104 namespace Teuchos {
00105  
00106 template<typename OrdinalType, typename ScalarType>
00107 class SerialSymDenseMatrix : public CompObject, public Object, public BLAS<OrdinalType,ScalarType>
00108 {
00109  public:
00110 
00112   typedef OrdinalType ordinalType;
00114   typedef ScalarType scalarType;
00115 
00117 
00118 
00119 
00127   SerialSymDenseMatrix();
00128 
00130 
00140   SerialSymDenseMatrix(OrdinalType numRowsCols, bool zeroOut = true);
00141 
00143 
00154   SerialSymDenseMatrix(DataAccess CV, bool upper, ScalarType* values, OrdinalType stride, OrdinalType numRowsCols);
00155   
00157   SerialSymDenseMatrix(const SerialSymDenseMatrix<OrdinalType, ScalarType> &Source);
00158 
00160 
00169   SerialSymDenseMatrix(DataAccess CV, const SerialSymDenseMatrix<OrdinalType, ScalarType> &Source, OrdinalType numRowCols, OrdinalType startRowCol=0);
00170 
00172   virtual ~SerialSymDenseMatrix ();
00174 
00176 
00177 
00179 
00188   int shape(OrdinalType numRowsCols); 
00189 
00191 
00200   int shapeUninitialized(OrdinalType numRowsCols); 
00201 
00203 
00213   int reshape(OrdinalType numRowsCols);
00214 
00216 
00218   void setLower();
00219 
00221 
00223   void setUpper();
00224 
00226 
00228 
00229 
00231 
00237   SerialSymDenseMatrix<OrdinalType, ScalarType>& operator= (const SerialSymDenseMatrix<OrdinalType, ScalarType>& Source);
00238 
00240 
00244   SerialSymDenseMatrix<OrdinalType, ScalarType>& assign (const SerialSymDenseMatrix<OrdinalType, ScalarType>& Source);
00245 
00247 
00250   SerialSymDenseMatrix<OrdinalType, ScalarType>& operator= (const ScalarType value) { putScalar(value); return(*this); }
00251 
00253 
00258   int putScalar( const ScalarType value = Teuchos::ScalarTraits<ScalarType>::zero(), bool fullMatrix = false );
00259 
00261 
00263   int random( const ScalarType bias = 0.1*Teuchos::ScalarTraits<ScalarType>::one() );
00264 
00266 
00268 
00269 
00271 
00278   ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex);
00279 
00281 
00288   const ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex) const;
00289 
00291 
00293   ScalarType* values() const { return(values_); }
00294 
00296 
00298 
00299 
00301   bool upper() const {return(upper_);};
00302 
00304   char UPLO() const {return(UPLO_);};
00306 
00308 
00309 
00311 
00317   SerialSymDenseMatrix<OrdinalType, ScalarType>& operator*= (const ScalarType alpha);
00318 
00320 
00323   SerialSymDenseMatrix<OrdinalType, ScalarType>& operator+= (const SerialSymDenseMatrix<OrdinalType, ScalarType>& Source);
00324 
00326 
00329   SerialSymDenseMatrix<OrdinalType, ScalarType>& operator-= (const SerialSymDenseMatrix<OrdinalType, ScalarType>& Source);
00330 
00332 
00334 
00335 
00337 
00340   bool operator== (const SerialSymDenseMatrix<OrdinalType, ScalarType> &Operand) const;
00341 
00343 
00346   bool operator!= (const SerialSymDenseMatrix<OrdinalType, ScalarType> &Operand) const;
00347 
00349 
00351 
00352 
00354   OrdinalType numRows() const { return(numRowCols_); }
00355 
00357   OrdinalType numCols() const { return(numRowCols_); }
00358 
00360   OrdinalType stride() const { return(stride_); }
00361 
00363   bool empty() const { return(numRowCols_ == 0); }
00364 
00366 
00368 
00369 
00371   typename ScalarTraits<ScalarType>::magnitudeType normOne() const;
00372 
00374   typename ScalarTraits<ScalarType>::magnitudeType normInf() const;
00375 
00377   typename ScalarTraits<ScalarType>::magnitudeType normFrobenius() const; 
00379 
00381 
00382 
00383   virtual void print(std::ostream& os) const;
00384 
00386 
00387  protected:
00388 
00389   // In-place scaling of matrix.
00390   void scale( const ScalarType alpha );
00391 
00392   // Copy the values from one matrix to the other.  
00393   void copyMat(bool inputUpper, ScalarType* inputMatrix, OrdinalType inputStride,
00394          OrdinalType numRowCols, bool outputUpper, ScalarType* outputMatrix, 
00395          OrdinalType outputStride, OrdinalType startRowCol, 
00396          ScalarType alpha = ScalarTraits<ScalarType>::zero() );
00397 
00398   // Copy the values from the active triangle of the matrix to the other to make the matrix full symmetric.
00399   void copyUPLOMat(bool inputUpper, ScalarType* inputMatrix, 
00400                    OrdinalType inputStride, OrdinalType inputRows);
00401 
00402   void deleteArrays();
00403   void checkIndex( OrdinalType rowIndex, OrdinalType colIndex = 0 ) const;
00404   OrdinalType numRowCols_;
00405   OrdinalType stride_;
00406   bool valuesCopied_;
00407   ScalarType* values_;
00408   bool upper_;
00409   char UPLO_;
00410 
00411 
00412 };
00413 
00414 //----------------------------------------------------------------------------------------------------
00415 //  Constructors and Destructor
00416 //----------------------------------------------------------------------------------------------------
00417 
00418 template<typename OrdinalType, typename ScalarType>
00419 SerialSymDenseMatrix<OrdinalType, ScalarType>::SerialSymDenseMatrix()
00420   : CompObject(), Object("Teuchos::SerialSymDenseMatrix"), numRowCols_(0), stride_(0), valuesCopied_(false), values_(0), upper_(false), UPLO_('L')
00421 {}
00422 
00423 template<typename OrdinalType, typename ScalarType>
00424 SerialSymDenseMatrix<OrdinalType, ScalarType>::SerialSymDenseMatrix(OrdinalType numRowCols_in, bool zeroOut)
00425   : CompObject(), Object("Teuchos::SerialSymDenseMatrix"), numRowCols_(numRowCols_in), stride_(numRowCols_in), valuesCopied_(false), upper_(false), UPLO_('L')
00426 {
00427   values_ = new ScalarType[stride_*numRowCols_];
00428   valuesCopied_ = true;
00429   if (zeroOut == true)
00430     putScalar( Teuchos::ScalarTraits<ScalarType>::zero(), true );
00431 }
00432   
00433 template<typename OrdinalType, typename ScalarType>
00434 SerialSymDenseMatrix<OrdinalType, ScalarType>::SerialSymDenseMatrix(
00435   DataAccess CV, bool upper_in, ScalarType* values_in, OrdinalType stride_in, OrdinalType numRowCols_in
00436   )
00437   : CompObject(), Object("Teuchos::SerialSymDenseMatrix"), numRowCols_(numRowCols_in), stride_(stride_in), valuesCopied_(false), 
00438     values_(values_in), upper_(upper_in)
00439 {
00440   if (upper_)
00441     UPLO_ = 'U';
00442   else
00443     UPLO_ = 'L';
00444 
00445   if(CV == Copy)
00446   {
00447     stride_ = numRowCols_;
00448     values_ = new ScalarType[stride_*numRowCols_];
00449     copyMat(upper_in, values_in, stride_in, numRowCols_, upper_, values_, stride_, 0);
00450     valuesCopied_ = true;
00451   }
00452 }
00453   
00454 template<typename OrdinalType, typename ScalarType>
00455 SerialSymDenseMatrix<OrdinalType, ScalarType>::SerialSymDenseMatrix(const SerialSymDenseMatrix<OrdinalType, ScalarType> &Source) : CompObject(), Object("Teuchos::SerialSymDenseMatrix"), numRowCols_(Source.numRowCols_), stride_(Source.numRowCols_), valuesCopied_(true), values_(0), upper_(Source.upper_), UPLO_(Source.UPLO_) 
00456 {
00457   values_ = new ScalarType[stride_*numRowCols_];
00458   copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0);
00459   valuesCopied_ = true;
00460 }
00461   
00462 template<typename OrdinalType, typename ScalarType>
00463 SerialSymDenseMatrix<OrdinalType, ScalarType>::SerialSymDenseMatrix(
00464                     DataAccess CV, const SerialSymDenseMatrix<OrdinalType, 
00465                     ScalarType> &Source, OrdinalType numRowCols_in, OrdinalType startRowCol )
00466   : CompObject(), Object("Teuchos::SerialSymDenseMatrix"), numRowCols_(numRowCols_in), stride_(Source.stride_), valuesCopied_(false), upper_(Source.upper_), UPLO_(Source.UPLO_)
00467 {
00468   if(CV == Copy)
00469   {
00470     stride_ = numRowCols_in;
00471     values_ = new ScalarType[stride_ * numRowCols_in];
00472     copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_in, upper_, values_, stride_, startRowCol);
00473     valuesCopied_ = true;
00474   }
00475   else // CV == View
00476   {
00477     values_ = Source.values_ + (stride_ * startRowCol) + startRowCol;
00478   }
00479 }
00480     
00481 template<typename OrdinalType, typename ScalarType>
00482 SerialSymDenseMatrix<OrdinalType, ScalarType>::~SerialSymDenseMatrix()
00483 {
00484   deleteArrays();
00485 }
00486 
00487 //----------------------------------------------------------------------------------------------------
00488 //  Shape methods 
00489 //----------------------------------------------------------------------------------------------------
00490 
00491 template<typename OrdinalType, typename ScalarType>
00492 int SerialSymDenseMatrix<OrdinalType, ScalarType>::shape( OrdinalType numRowCols_in )
00493 {
00494   deleteArrays(); // Get rid of anything that might be already allocated
00495   numRowCols_ = numRowCols_in;
00496   stride_ = numRowCols_;
00497   values_ = new ScalarType[stride_*numRowCols_];
00498   putScalar( Teuchos::ScalarTraits<ScalarType>::zero(), true );
00499   valuesCopied_ = true;
00500   return(0);
00501 }
00502 
00503 template<typename OrdinalType, typename ScalarType>
00504 int SerialSymDenseMatrix<OrdinalType, ScalarType>::shapeUninitialized( OrdinalType numRowCols_in )
00505 {
00506   deleteArrays(); // Get rid of anything that might be already allocated
00507   numRowCols_ = numRowCols_in;
00508   stride_ = numRowCols_;
00509   values_ = new ScalarType[stride_*numRowCols_];
00510   valuesCopied_ = true;
00511   return(0);
00512 }
00513   
00514 template<typename OrdinalType, typename ScalarType>
00515 int SerialSymDenseMatrix<OrdinalType, ScalarType>::reshape( OrdinalType numRowCols_in )
00516 {
00517   // Allocate space for new matrix
00518   ScalarType* values_tmp = new ScalarType[numRowCols_in * numRowCols_in];
00519   ScalarType zero = ScalarTraits<ScalarType>::zero();
00520   for(OrdinalType k = 0; k < numRowCols_in * numRowCols_in; k++)
00521   {
00522     values_tmp[k] = zero;
00523   }
00524   OrdinalType numRowCols_tmp = TEUCHOS_MIN(numRowCols_, numRowCols_in);
00525   if(values_ != 0)
00526   {
00527     copyMat(upper_, values_, stride_, numRowCols_tmp, upper_, values_tmp, numRowCols_in, 0); // Copy principal submatrix of A to new A
00528   }
00529   deleteArrays(); // Get rid of anything that might be already allocated
00530   numRowCols_ = numRowCols_in;
00531   stride_ = numRowCols_;
00532   values_ = values_tmp; // Set pointer to new A
00533   valuesCopied_ = true;
00534   return(0);
00535 }
00536    
00537 //----------------------------------------------------------------------------------------------------
00538 //   Set methods 
00539 //----------------------------------------------------------------------------------------------------
00540 
00541 template<typename OrdinalType, typename ScalarType>
00542 void SerialSymDenseMatrix<OrdinalType, ScalarType>::setLower() 
00543 {
00544   // Do nothing if the matrix is already a lower triangular matrix
00545   if (upper_ != false) {
00546     copyUPLOMat( true, values_, stride_, numRowCols_ );
00547     upper_ = false;
00548     UPLO_ = 'L';
00549   }
00550 }
00551 
00552 template<typename OrdinalType, typename ScalarType>
00553 void SerialSymDenseMatrix<OrdinalType, ScalarType>::setUpper() 
00554 {
00555   // Do nothing if the matrix is already an upper triangular matrix
00556   if (upper_ == false) {
00557     copyUPLOMat( false, values_, stride_, numRowCols_ );
00558     upper_ = true;
00559     UPLO_ = 'U';
00560   }
00561 }
00562 
00563 template<typename OrdinalType, typename ScalarType>
00564 int SerialSymDenseMatrix<OrdinalType, ScalarType>::putScalar( const ScalarType value_in, bool fullMatrix )
00565 {
00566   // Set each value of the dense matrix to "value".
00567   if (fullMatrix == true) {
00568     for(OrdinalType j = 0; j < numRowCols_; j++) 
00569       {
00570   for(OrdinalType i = 0; i < numRowCols_; i++) 
00571     {
00572       values_[i + j*stride_] = value_in;
00573     }
00574       }
00575   }
00576   // Set the active upper or lower triangular part of the matrix to "value"
00577   else {
00578     if (upper_) {
00579       for(OrdinalType j = 0; j < numRowCols_; j++) {
00580   for(OrdinalType i = 0; i <= j; i++) {
00581     values_[i + j*stride_] = value_in;
00582   }
00583       }
00584     }
00585     else {
00586       for(OrdinalType j = 0; j < numRowCols_; j++) {
00587   for(OrdinalType i = j; i < numRowCols_; i++) {
00588     values_[i + j*stride_] = value_in;
00589   }
00590       }
00591     }
00592   }
00593   return 0;
00594 }    
00595     
00596 template<typename OrdinalType, typename ScalarType>
00597 int SerialSymDenseMatrix<OrdinalType, ScalarType>::random( const ScalarType bias )
00598 {
00599   typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MT;
00600 
00601   // Set each value of the dense matrix to a random value.
00602   std::vector<MT> diagSum( numRowCols_, Teuchos::ScalarTraits<MT>::zero() );
00603   if (upper_) {
00604     for(OrdinalType j = 0; j < numRowCols_; j++) {
00605       for(OrdinalType i = 0; i < j; i++) {
00606   values_[i + j*stride_] = ScalarTraits<ScalarType>::random();
00607   diagSum[i] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
00608   diagSum[j] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
00609       }
00610     }
00611   }
00612   else {
00613     for(OrdinalType j = 0; j < numRowCols_; j++) {
00614       for(OrdinalType i = j+1; i < numRowCols_; i++) {
00615   values_[i + j*stride_] = ScalarTraits<ScalarType>::random();
00616   diagSum[i] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
00617   diagSum[j] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
00618       }
00619     }
00620   }
00621   
00622   // Set the diagonal.
00623   for(OrdinalType i = 0; i < numRowCols_; i++) {
00624     values_[i + i*stride_] = diagSum[i] + bias;
00625   }
00626   return 0;
00627 }
00628 
00629 template<typename OrdinalType, typename ScalarType>
00630 SerialSymDenseMatrix<OrdinalType,ScalarType>&
00631 SerialSymDenseMatrix<OrdinalType, ScalarType>::operator=( const SerialSymDenseMatrix<OrdinalType,ScalarType>& Source )
00632 {
00633   if(this == &Source)
00634     return(*this); // Special case of source same as target
00635   if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_)) {
00636     upper_ = Source.upper_;  // Might have to change the active part of the matrix.
00637     return(*this); // Special case of both are views to same data.
00638   }
00639   
00640   // If the source is a view then we will return a view, else we will return a copy.
00641   if (!Source.valuesCopied_) {
00642     if(valuesCopied_)       { 
00643       // Clean up stored data if this was previously a copy.
00644       deleteArrays();
00645     }
00646     numRowCols_ = Source.numRowCols_; 
00647     stride_ = Source.stride_;
00648     values_ = Source.values_;
00649     upper_ = Source.upper_;
00650     UPLO_ = Source.UPLO_;
00651   }
00652   else {
00653     // If we were a view, we will now be a copy.
00654     if(!valuesCopied_) {
00655       numRowCols_ = Source.numRowCols_;
00656       stride_ = Source.numRowCols_;
00657       upper_ = Source.upper_;
00658       UPLO_ = Source.UPLO_;
00659       const OrdinalType newsize = stride_ * numRowCols_;
00660       if(newsize > 0) {
00661         values_ = new ScalarType[newsize];
00662         valuesCopied_ = true;
00663       }
00664       else {
00665         values_ = 0;
00666       }
00667     }
00668     // If we were a copy, we will stay a copy.
00669     else {
00670       if((Source.numRowCols_ <= stride_) && (Source.numRowCols_ == numRowCols_)) { // we don't need to reallocate
00671         numRowCols_ = Source.numRowCols_;
00672   upper_ = Source.upper_;
00673   UPLO_ = Source.UPLO_;
00674       }
00675       else { // we need to allocate more space (or less space)
00676         deleteArrays();
00677         numRowCols_ = Source.numRowCols_;
00678         stride_ = Source.numRowCols_;
00679   upper_ = Source.upper_;
00680   UPLO_ = Source.UPLO_;
00681         const OrdinalType newsize = stride_ * numRowCols_;
00682         if(newsize > 0) {
00683           values_ = new ScalarType[newsize];
00684           valuesCopied_ = true;
00685         }
00686       }
00687     }
00688     copyMat(Source.upper_, Source.values_, Source.stride_, Source.numRowCols_, upper_, values_, stride_, 0);
00689   } 
00690   return(*this);
00691 }
00692 
00693 template<typename OrdinalType, typename ScalarType>
00694 SerialSymDenseMatrix<OrdinalType, ScalarType>& SerialSymDenseMatrix<OrdinalType, ScalarType>::operator*= (const ScalarType alpha)
00695 {
00696   this->scale(alpha);
00697   return(*this);
00698 }
00699 
00700 template<typename OrdinalType, typename ScalarType>
00701 SerialSymDenseMatrix<OrdinalType, ScalarType>& SerialSymDenseMatrix<OrdinalType, ScalarType>::operator+= (const SerialSymDenseMatrix<OrdinalType,ScalarType>& Source )
00702 {
00703   // Check for compatible dimensions
00704   if ((numRowCols_ != Source.numRowCols_))
00705     {
00706       TEUCHOS_CHK_REF(*this); // Return *this without altering it.
00707     }
00708   copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0, ScalarTraits<ScalarType>::one());
00709   return(*this);
00710 }
00711 
00712 template<typename OrdinalType, typename ScalarType>
00713 SerialSymDenseMatrix<OrdinalType, ScalarType>& SerialSymDenseMatrix<OrdinalType, ScalarType>::operator-= (const SerialSymDenseMatrix<OrdinalType,ScalarType>& Source )
00714 {
00715   // Check for compatible dimensions
00716   if ((numRowCols_ != Source.numRowCols_))
00717   {
00718     TEUCHOS_CHK_REF(*this); // Return *this without altering it.
00719   }
00720   copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0, -ScalarTraits<ScalarType>::one());
00721   return(*this);
00722 }
00723 
00724 template<typename OrdinalType, typename ScalarType>
00725 SerialSymDenseMatrix<OrdinalType,ScalarType>& SerialSymDenseMatrix<OrdinalType, ScalarType>::assign (const SerialSymDenseMatrix<OrdinalType,ScalarType>& Source) {
00726   if(this == &Source)
00727     return(*this); // Special case of source same as target
00728   if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_)) {
00729     upper_ = Source.upper_; // We may have to change the active part of the matrix.
00730     return(*this); // Special case of both are views to same data.
00731   }
00732 
00733   // Check for compatible dimensions
00734   if ((numRowCols_ != Source.numRowCols_))
00735   {
00736     TEUCHOS_CHK_REF(*this); // Return *this without altering it.
00737   }
00738   copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0 );
00739   return(*this);
00740 }
00741 
00742 //----------------------------------------------------------------------------------------------------
00743 //   Accessor methods 
00744 //----------------------------------------------------------------------------------------------------
00745 
00746 template<typename OrdinalType, typename ScalarType>
00747 inline ScalarType& SerialSymDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex)
00748 {
00749 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
00750   checkIndex( rowIndex, colIndex );
00751 #endif
00752   if ( rowIndex <= colIndex ) {
00753     // Accessing upper triangular part of matrix
00754     if (upper_)
00755       return(values_[colIndex * stride_ + rowIndex]);
00756     else
00757       return(values_[rowIndex * stride_ + colIndex]);   
00758   }
00759   else {
00760     // Accessing lower triangular part of matrix
00761     if (upper_)
00762       return(values_[rowIndex * stride_ + colIndex]);
00763     else
00764       return(values_[colIndex * stride_ + rowIndex]);
00765   }
00766 }
00767   
00768 template<typename OrdinalType, typename ScalarType>
00769 inline const ScalarType& SerialSymDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex) const
00770 {
00771 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
00772   checkIndex( rowIndex, colIndex );
00773 #endif
00774   if ( rowIndex <= colIndex ) {
00775     // Accessing upper triangular part of matrix
00776     if (upper_)
00777       return(values_[colIndex * stride_ + rowIndex]);
00778     else
00779       return(values_[rowIndex * stride_ + colIndex]);   
00780   }
00781   else {
00782     // Accessing lower triangular part of matrix
00783     if (upper_)
00784       return(values_[rowIndex * stride_ + colIndex]);
00785     else
00786       return(values_[colIndex * stride_ + rowIndex]);
00787   }
00788 }  
00789 
00790 //----------------------------------------------------------------------------------------------------
00791 //   Norm methods 
00792 //----------------------------------------------------------------------------------------------------
00793   
00794 template<typename OrdinalType, typename ScalarType>
00795 typename ScalarTraits<ScalarType>::magnitudeType SerialSymDenseMatrix<OrdinalType, ScalarType>::normOne() const
00796 {
00797   return(normInf());
00798 }
00799   
00800 template<typename OrdinalType, typename ScalarType>
00801 typename ScalarTraits<ScalarType>::magnitudeType SerialSymDenseMatrix<OrdinalType, ScalarType>::normInf() const
00802 {
00803   typedef typename ScalarTraits<ScalarType>::magnitudeType MT;
00804 
00805   OrdinalType i, j;
00806   
00807   MT sum, anorm = ScalarTraits<MT>::zero();
00808   ScalarType* ptr;
00809   
00810   if (upper_) {
00811     for (j=0; j<numRowCols_; j++) {
00812       sum = ScalarTraits<MT>::zero();
00813       ptr = values_ + j*stride_;
00814       for (i=0; i<j; i++) {
00815         sum += ScalarTraits<ScalarType>::magnitude( *ptr++ );
00816       }
00817       ptr = values_ + j + j*stride_;
00818       for (i=j; i<numRowCols_; i++) {
00819         sum += ScalarTraits<ScalarType>::magnitude( *ptr );
00820         ptr += stride_;
00821       }
00822       anorm = TEUCHOS_MAX( anorm, sum );
00823     }
00824   }
00825   else {
00826     for (j=0; j<numRowCols_; j++) {
00827       sum = ScalarTraits<MT>::zero();
00828       ptr = values_ + j + j*stride_;
00829       for (i=j; i<numRowCols_; i++) {
00830         sum += ScalarTraits<ScalarType>::magnitude( *ptr++ );
00831       }
00832       ptr = values_ + j;
00833       for (i=0; i<j; i++) {
00834         sum += ScalarTraits<ScalarType>::magnitude( *ptr );
00835         ptr += stride_;
00836       }
00837       anorm = TEUCHOS_MAX( anorm, sum );
00838     }
00839   }
00840   return(anorm);
00841 }
00842   
00843 template<typename OrdinalType, typename ScalarType>
00844 typename ScalarTraits<ScalarType>::magnitudeType SerialSymDenseMatrix<OrdinalType, ScalarType>::normFrobenius() const
00845 {
00846   typedef typename ScalarTraits<ScalarType>::magnitudeType MT;
00847 
00848   OrdinalType i, j;
00849   
00850   MT sum = ScalarTraits<MT>::zero(), anorm = ScalarTraits<MT>::zero();
00851  
00852   if (upper_) { 
00853     for (j = 0; j < numRowCols_; j++) {
00854       for (i = 0; i < j; i++) {
00855         sum += ScalarTraits<ScalarType>::magnitude(2.0*values_[i+j*stride_]*values_[i+j*stride_]);
00856       }
00857       sum += ScalarTraits<ScalarType>::magnitude(values_[j + j*stride_]*values_[j + j*stride_]);
00858     }
00859   }
00860   else {
00861     for (j = 0; j < numRowCols_; j++) {
00862       sum += ScalarTraits<ScalarType>::magnitude(values_[j + j*stride_]*values_[j + j*stride_]);
00863       for (i = j+1; i < numRowCols_; i++) {
00864         sum += ScalarTraits<ScalarType>::magnitude(2.0*values_[i+j*stride_]*values_[i+j*stride_]);
00865       }
00866     }
00867   }       
00868   anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::squareroot(sum));
00869   return(anorm);
00870 }
00871   
00872 //----------------------------------------------------------------------------------------------------
00873 //   Comparison methods 
00874 //----------------------------------------------------------------------------------------------------
00875   
00876 template<typename OrdinalType, typename ScalarType>
00877 bool SerialSymDenseMatrix<OrdinalType, ScalarType>::operator== (const SerialSymDenseMatrix<OrdinalType, ScalarType> &Operand) const
00878 {
00879   bool result = 1;
00880   if((numRowCols_ != Operand.numRowCols_)) {
00881     result = 0;
00882   }
00883   else {
00884     OrdinalType i, j;
00885     for(i = 0; i < numRowCols_; i++) {
00886       for(j = 0; j < numRowCols_; j++) {
00887   if((*this)(i, j) != Operand(i, j)) {
00888     return 0;
00889   }
00890       }
00891     }
00892   }
00893   return result;
00894 }
00895 
00896 template<typename OrdinalType, typename ScalarType>
00897 bool SerialSymDenseMatrix<OrdinalType, ScalarType>::operator!= (const SerialSymDenseMatrix<OrdinalType, ScalarType> &Operand) const
00898 {
00899   return !((*this) == Operand);
00900 }
00901 
00902 //----------------------------------------------------------------------------------------------------
00903 //   Multiplication method 
00904 //----------------------------------------------------------------------------------------------------  
00905 
00906 template<typename OrdinalType, typename ScalarType>
00907 void SerialSymDenseMatrix<OrdinalType, ScalarType>::scale( const ScalarType alpha )
00908 {
00909   OrdinalType i, j;
00910   ScalarType* ptr;
00911   
00912   if (upper_) {
00913     for (j=0; j<numRowCols_; j++) {
00914       ptr = values_ + j*stride_;
00915       for (i=0; i<= j; i++) { *ptr = alpha * (*ptr); ptr++; }
00916     }
00917   }
00918   else {
00919     for (j=0; j<numRowCols_; j++) {
00920       ptr = values_ + j*stride_;
00921       for (i=j; i<numRowCols_; i++) { *ptr = alpha * (*ptr); ptr++; }
00922     }
00923   }
00924 }
00925 
00926 /*
00927 template<typename OrdinalType, typename ScalarType>
00928 int SerialSymDenseMatrix<OrdinalType, ScalarType>::scale( const SerialSymDenseMatrix<OrdinalType,ScalarType>& A )
00929 {
00930   OrdinalType i, j;
00931   ScalarType* ptr;
00932     
00933   // Check for compatible dimensions
00934   if ((numRowCols_ != A.numRowCols_)) {
00935     TEUCHOS_CHK_ERR(-1); // Return error
00936   }    
00937 
00938   if (upper_) {
00939     for (j=0; j<numRowCols_; j++) {
00940       ptr = values_ + j*stride_;
00941       for (i=0; i<= j; i++) { *ptr = A(i,j) * (*ptr); ptr++; }
00942     }
00943   }
00944   else {
00945     for (j=0; j<numRowCols_; j++) {
00946       ptr = values_ + j*stride_;
00947       for (i=j; i<numRowCols_; i++) { *ptr = A(i,j) * (*ptr); ptr++; }
00948     }
00949   }
00950 
00951   return(0);
00952 }
00953 */
00954 
00955 template<typename OrdinalType, typename ScalarType>
00956 void SerialSymDenseMatrix<OrdinalType, ScalarType>::print(std::ostream& os) const
00957 {
00958   os << std::endl;
00959   if(valuesCopied_)
00960     os << "Values_copied : yes" << std::endl;
00961   else
00962     os << "Values_copied : no" << std::endl;
00963   os << "Rows / Columns : " << numRowCols_ << std::endl;
00964   os << "LDA : " << stride_ << std::endl;
00965   if (upper_) 
00966     os << "Storage: Upper" << std::endl;
00967   else
00968     os << "Storage: Lower" << std::endl;
00969   if(numRowCols_ == 0) {
00970     os << "(matrix is empty, no values to display)" << std::endl;
00971   } else {
00972     for(OrdinalType i = 0; i < numRowCols_; i++) {
00973       for(OrdinalType j = 0; j < numRowCols_; j++){
00974         os << (*this)(i,j) << " ";
00975       }
00976       os << std::endl;
00977     }
00978   }
00979 }
00980 
00981 //----------------------------------------------------------------------------------------------------
00982 //   Protected methods 
00983 //----------------------------------------------------------------------------------------------------  
00984 
00985 template<typename OrdinalType, typename ScalarType>
00986 inline void SerialSymDenseMatrix<OrdinalType, ScalarType>::checkIndex( OrdinalType rowIndex, OrdinalType colIndex ) const {
00987   TEST_FOR_EXCEPTION(rowIndex < 0 || rowIndex >= numRowCols_, std::out_of_range,
00988     "SerialSymDenseMatrix<T>::checkIndex: "
00989     "Row index " << rowIndex << " out of range [0, "<< numRowCols_ << ")");
00990   TEST_FOR_EXCEPTION(colIndex < 0 || colIndex >= numRowCols_, std::out_of_range,
00991     "SerialSymDenseMatrix<T>::checkIndex: "
00992     "Col index " << colIndex << " out of range [0, "<< numRowCols_ << ")");
00993 }
00994 
00995 template<typename OrdinalType, typename ScalarType>
00996 void SerialSymDenseMatrix<OrdinalType, ScalarType>::deleteArrays(void)
00997 {
00998   if (valuesCopied_)
00999   {
01000     delete [] values_;
01001     values_ = 0;
01002     valuesCopied_ = false;
01003   }
01004 }
01005 
01006 template<typename OrdinalType, typename ScalarType>
01007 void SerialSymDenseMatrix<OrdinalType, ScalarType>::copyMat(
01008                   bool inputUpper, ScalarType* inputMatrix, 
01009                   OrdinalType inputStride, OrdinalType numRowCols_in, 
01010                   bool outputUpper, ScalarType* outputMatrix, 
01011                   OrdinalType outputStride, OrdinalType startRowCol, 
01012                   ScalarType alpha
01013                   )
01014 {
01015   OrdinalType i, j;
01016   ScalarType* ptr1 = 0;
01017   ScalarType* ptr2 = 0;
01018 
01019   for (j = 0; j < numRowCols_in; j++) {
01020     if (inputUpper == true) {
01021       // The input matrix is upper triangular, start at the beginning of each column.
01022       ptr2 = inputMatrix + (j + startRowCol) * inputStride + startRowCol;
01023       if (outputUpper == true) {
01024   // The output matrix matches the same pattern as the input matrix.
01025   ptr1 = outputMatrix + j*outputStride;
01026   if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
01027     for(i = 0; i <= j; i++) {
01028       *ptr1++ += alpha*(*ptr2++);
01029     }
01030   } else {
01031     for(i = 0; i <= j; i++) {
01032       *ptr1++ = *ptr2++;
01033     }
01034   }
01035       }
01036       else {
01037   // The output matrix has the opposite pattern as the input matrix.
01038   // Copy down across rows of the output matrix, but down columns of the input matrix.
01039   ptr1 = outputMatrix + j;
01040   if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
01041     for(i = 0; i <= j; i++) {
01042       *ptr1 += alpha*(*ptr2++);
01043       ptr1 += outputStride;
01044     }
01045   } else {
01046     for(i = 0; i <= j; i++) {
01047       *ptr1 = *ptr2++;
01048       ptr1 += outputStride;
01049     }
01050   }
01051       }
01052     }
01053     else {
01054       // The input matrix is lower triangular, start at the diagonal of each row.
01055       ptr2 = inputMatrix + (startRowCol+j) * inputStride + startRowCol + j;
01056       if (outputUpper == true) {
01057   // The output matrix has the opposite pattern as the input matrix.
01058   // Copy across rows of the output matrix, but down columns of the input matrix.
01059   ptr1 = outputMatrix + j*outputStride + j;
01060   if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
01061     for(i = j; i < numRowCols_in; i++) {
01062       *ptr1 += alpha*(*ptr2++);
01063       ptr1 += outputStride;
01064     }
01065   } else {
01066     for(i = j; i < numRowCols_in; i++) {
01067       *ptr1 = *ptr2++;
01068       ptr1 += outputStride;
01069     }
01070   }
01071       }
01072       else {
01073   // The output matrix matches the same pattern as the input matrix.
01074   ptr1 = outputMatrix + j*outputStride + j;
01075   if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
01076     for(i = j; i < numRowCols_in; i++) {
01077       *ptr1++ += alpha*(*ptr2++);
01078     }
01079   } else {
01080     for(i = j; i < numRowCols_in; i++) {
01081       *ptr1++ = *ptr2++;
01082     }
01083   }
01084       }
01085     }
01086   }
01087 }
01088   
01089 template<typename OrdinalType, typename ScalarType>
01090 void SerialSymDenseMatrix<OrdinalType, ScalarType>::copyUPLOMat(
01091                 bool inputUpper, ScalarType* inputMatrix, 
01092                 OrdinalType inputStride, OrdinalType inputRows 
01093                 )
01094 {
01095   OrdinalType i, j;
01096   ScalarType * ptr1 = 0;
01097   ScalarType * ptr2 = 0;
01098 
01099   if (inputUpper) {
01100     for (j=1; j<inputRows; j++) {
01101       ptr1 = inputMatrix + j;
01102       ptr2 = inputMatrix + j*inputStride;
01103       for (i=0; i<j; i++) {
01104   *ptr1 = *ptr2++;
01105   ptr1+=inputStride;
01106       }
01107     }
01108   }
01109   else {
01110     for (i=1; i<inputRows; i++) {
01111       ptr1 = inputMatrix + i;
01112       ptr2 = inputMatrix + i*inputStride;
01113       for (j=0; j<i; j++) {
01114   *ptr2++ = *ptr1;
01115   ptr1+=inputStride;
01116       }
01117     }
01118   }
01119 }
01120   
01121 } // namespace Teuchos
01122 
01123 #endif /* _TEUCHOS_SERIALSYMDENSEMATRIX_HPP_ */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 09:57:30 2011 for Teuchos Package Browser (Single Doxygen Collection) by  doxygen 1.6.3