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 
00113 
00114 
00122   SerialSymDenseMatrix();
00123 
00125 
00135   SerialSymDenseMatrix(OrdinalType numRowsCols, bool zeroOut = true);
00136 
00138 
00149   SerialSymDenseMatrix(DataAccess CV, bool upper, ScalarType* values, OrdinalType stride, OrdinalType numRowsCols);
00150   
00152   SerialSymDenseMatrix(const SerialSymDenseMatrix<OrdinalType, ScalarType> &Source);
00153 
00155 
00164   SerialSymDenseMatrix(DataAccess CV, const SerialSymDenseMatrix<OrdinalType, ScalarType> &Source, OrdinalType numRowCols, OrdinalType startRowCol=0);
00165 
00167   virtual ~SerialSymDenseMatrix ();
00169 
00171 
00172 
00174 
00183   int shape(OrdinalType numRowsCols); 
00184 
00186 
00195   int shapeUninitialized(OrdinalType numRowsCols); 
00196 
00198 
00208   int reshape(OrdinalType numRowsCols);
00209 
00211 
00213   void setLower();
00214 
00216 
00218   void setUpper();
00219 
00221 
00223 
00224 
00226 
00232   SerialSymDenseMatrix<OrdinalType, ScalarType>& operator= (const SerialSymDenseMatrix<OrdinalType, ScalarType>& Source);
00233 
00235 
00239   SerialSymDenseMatrix<OrdinalType, ScalarType>& assign (const SerialSymDenseMatrix<OrdinalType, ScalarType>& Source);
00240 
00242 
00247   int putScalar( const ScalarType value = Teuchos::ScalarTraits<ScalarType>::zero(), bool fullMatrix = false );
00248 
00250 
00252   int random( const ScalarType bias = 0.1*Teuchos::ScalarTraits<ScalarType>::one() );
00253 
00255 
00257 
00258 
00260 
00267   ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex);
00268 
00270 
00277   const ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex) const;
00278 
00280 
00282   ScalarType* values() const { return(values_); }
00283 
00285 
00287 
00288 
00290   bool upper() const {return(upper_);};
00291 
00293   char UPLO() const {return(UPLO_);};
00295 
00297 
00298 
00300 
00306   SerialSymDenseMatrix<OrdinalType, ScalarType>& operator*= (const ScalarType alpha);
00307 
00309 
00312   SerialSymDenseMatrix<OrdinalType, ScalarType>& operator+= (const SerialSymDenseMatrix<OrdinalType, ScalarType>& Source);
00313 
00315 
00318   SerialSymDenseMatrix<OrdinalType, ScalarType>& operator-= (const SerialSymDenseMatrix<OrdinalType, ScalarType>& Source);
00319 
00321 
00323 
00324 
00326 
00329   bool operator== (const SerialSymDenseMatrix<OrdinalType, ScalarType> &Operand);
00330 
00332 
00335   bool operator!= (const SerialSymDenseMatrix<OrdinalType, ScalarType> &Operand);
00336 
00338 
00340 
00341 
00343   OrdinalType numRows() const { return(numRowCols_); }
00344 
00346   OrdinalType numCols() const { return(numRowCols_); }
00347 
00349   OrdinalType stride() const { return(stride_); }
00350 
00352   bool empty() const { return(numRowCols_ == 0); }
00353 
00355 
00357 
00358 
00360   typename ScalarTraits<ScalarType>::magnitudeType normOne() const;
00361 
00363   typename ScalarTraits<ScalarType>::magnitudeType normInf() const;
00364 
00366   typename ScalarTraits<ScalarType>::magnitudeType normFrobenius() const; 
00368 
00370 
00371 
00372   virtual void print(std::ostream& os) const;
00373 
00375 
00376  protected:
00377 
00378   // In-place scaling of matrix.
00379   void scale( const ScalarType alpha );
00380 
00381   // Copy the values from one matrix to the other.  
00382   void copyMat(bool inputUpper, ScalarType* inputMatrix, OrdinalType inputStride,
00383          OrdinalType numRowCols, bool outputUpper, ScalarType* outputMatrix, 
00384          OrdinalType outputStride, OrdinalType startRowCol, 
00385          ScalarType alpha = ScalarTraits<ScalarType>::zero() );
00386 
00387   // Copy the values from the active triangle of the matrix to the other to make the matrix full symmetric.
00388   void copyUPLOMat(bool inputUpper, ScalarType* inputMatrix, 
00389                    OrdinalType inputStride, OrdinalType inputRows);
00390 
00391   void deleteArrays();
00392   void checkIndex( OrdinalType rowIndex, OrdinalType colIndex = 0 ) const;
00393   OrdinalType numRowCols_;
00394   OrdinalType stride_;
00395   bool valuesCopied_;
00396   ScalarType* values_;
00397   bool upper_;
00398   char UPLO_;
00399 
00400 
00401 };
00402 
00403 //----------------------------------------------------------------------------------------------------
00404 //  Constructors and Destructor
00405 //----------------------------------------------------------------------------------------------------
00406 
00407 template<typename OrdinalType, typename ScalarType>
00408 SerialSymDenseMatrix<OrdinalType, ScalarType>::SerialSymDenseMatrix()
00409   : CompObject(), Object("Teuchos::SerialSymDenseMatrix"), numRowCols_(0), stride_(0), valuesCopied_(false), values_(0), upper_(false), UPLO_('L')
00410 {}
00411 
00412 template<typename OrdinalType, typename ScalarType>
00413 SerialSymDenseMatrix<OrdinalType, ScalarType>::SerialSymDenseMatrix(OrdinalType numRowCols_in, bool zeroOut)
00414   : CompObject(), Object("Teuchos::SerialSymDenseMatrix"), numRowCols_(numRowCols_in), stride_(numRowCols_in), valuesCopied_(false), upper_(false), UPLO_('L')
00415 {
00416   values_ = new ScalarType[stride_*numRowCols_];
00417   valuesCopied_ = true;
00418   if (zeroOut == true)
00419     putScalar( Teuchos::ScalarTraits<ScalarType>::zero(), true );
00420 }
00421   
00422 template<typename OrdinalType, typename ScalarType>
00423 SerialSymDenseMatrix<OrdinalType, ScalarType>::SerialSymDenseMatrix(
00424   DataAccess CV, bool upper_in, ScalarType* values_in, OrdinalType stride_in, OrdinalType numRowCols_in
00425   )
00426   : CompObject(), Object("Teuchos::SerialSymDenseMatrix"), numRowCols_(numRowCols_in), stride_(stride_in), valuesCopied_(false), 
00427     values_(values_in), upper_(upper_in)
00428 {
00429   if (upper_)
00430     UPLO_ = 'U';
00431   else
00432     UPLO_ = 'L';
00433 
00434   if(CV == Copy)
00435   {
00436     stride_ = numRowCols_;
00437     values_ = new ScalarType[stride_*numRowCols_];
00438     copyMat(upper_in, values_in, stride_in, numRowCols_, upper_, values_, stride_, 0);
00439     valuesCopied_ = true;
00440   }
00441 }
00442   
00443 template<typename OrdinalType, typename ScalarType>
00444 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_) 
00445 {
00446   values_ = new ScalarType[stride_*numRowCols_];
00447   copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0);
00448   valuesCopied_ = true;
00449 }
00450   
00451 template<typename OrdinalType, typename ScalarType>
00452 SerialSymDenseMatrix<OrdinalType, ScalarType>::SerialSymDenseMatrix(
00453                     DataAccess CV, const SerialSymDenseMatrix<OrdinalType, 
00454                     ScalarType> &Source, OrdinalType numRowCols_in, OrdinalType startRowCol )
00455   : CompObject(), Object("Teuchos::SerialSymDenseMatrix"), numRowCols_(numRowCols_in), stride_(Source.stride_), valuesCopied_(false), upper_(Source.upper_), UPLO_(Source.UPLO_)
00456 {
00457   if(CV == Copy)
00458   {
00459     stride_ = numRowCols_in;
00460     values_ = new ScalarType[stride_ * numRowCols_in];
00461     copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_in, upper_, values_, stride_, startRowCol);
00462     valuesCopied_ = true;
00463   }
00464   else // CV == View
00465   {
00466     values_ = Source.values_ + (stride_ * startRowCol) + startRowCol;
00467   }
00468 }
00469     
00470 template<typename OrdinalType, typename ScalarType>
00471 SerialSymDenseMatrix<OrdinalType, ScalarType>::~SerialSymDenseMatrix()
00472 {
00473   deleteArrays();
00474 }
00475 
00476 //----------------------------------------------------------------------------------------------------
00477 //  Shape methods 
00478 //----------------------------------------------------------------------------------------------------
00479 
00480 template<typename OrdinalType, typename ScalarType>
00481 int SerialSymDenseMatrix<OrdinalType, ScalarType>::shape( OrdinalType numRowCols_in )
00482 {
00483   deleteArrays(); // Get rid of anything that might be already allocated
00484   numRowCols_ = numRowCols_in;
00485   stride_ = numRowCols_;
00486   values_ = new ScalarType[stride_*numRowCols_];
00487   putScalar( Teuchos::ScalarTraits<ScalarType>::zero(), true );
00488   valuesCopied_ = true;
00489   return(0);
00490 }
00491 
00492 template<typename OrdinalType, typename ScalarType>
00493 int SerialSymDenseMatrix<OrdinalType, ScalarType>::shapeUninitialized( OrdinalType numRowCols_in )
00494 {
00495   deleteArrays(); // Get rid of anything that might be already allocated
00496   numRowCols_ = numRowCols_in;
00497   stride_ = numRowCols_;
00498   values_ = new ScalarType[stride_*numRowCols_];
00499   valuesCopied_ = true;
00500   return(0);
00501 }
00502   
00503 template<typename OrdinalType, typename ScalarType>
00504 int SerialSymDenseMatrix<OrdinalType, ScalarType>::reshape( OrdinalType numRowCols_in )
00505 {
00506   // Allocate space for new matrix
00507   ScalarType* values_tmp = new ScalarType[numRowCols_in * numRowCols_in];
00508   ScalarType zero = ScalarTraits<ScalarType>::zero();
00509   for(OrdinalType k = 0; k < numRowCols_in * numRowCols_in; k++)
00510   {
00511     values_tmp[k] = zero;
00512   }
00513   OrdinalType numRowCols_tmp = TEUCHOS_MIN(numRowCols_, numRowCols_in);
00514   if(values_ != 0)
00515   {
00516     copyMat(upper_, values_, stride_, numRowCols_tmp, upper_, values_tmp, numRowCols_in, 0); // Copy principal submatrix of A to new A
00517   }
00518   deleteArrays(); // Get rid of anything that might be already allocated
00519   numRowCols_ = numRowCols_in;
00520   stride_ = numRowCols_;
00521   values_ = values_tmp; // Set pointer to new A
00522   valuesCopied_ = true;
00523   return(0);
00524 }
00525    
00526 //----------------------------------------------------------------------------------------------------
00527 //   Set methods 
00528 //----------------------------------------------------------------------------------------------------
00529 
00530 template<typename OrdinalType, typename ScalarType>
00531 void SerialSymDenseMatrix<OrdinalType, ScalarType>::setLower() 
00532 {
00533   // Do nothing if the matrix is already a lower triangular matrix
00534   if (upper_ != false) {
00535     copyUPLOMat( true, values_, stride_, numRowCols_ );
00536     upper_ = false;
00537     UPLO_ = 'L';
00538   }
00539 }
00540 
00541 template<typename OrdinalType, typename ScalarType>
00542 void SerialSymDenseMatrix<OrdinalType, ScalarType>::setUpper() 
00543 {
00544   // Do nothing if the matrix is already an upper triangular matrix
00545   if (upper_ == false) {
00546     copyUPLOMat( false, values_, stride_, numRowCols_ );
00547     upper_ = true;
00548     UPLO_ = 'U';
00549   }
00550 }
00551 
00552 template<typename OrdinalType, typename ScalarType>
00553 int SerialSymDenseMatrix<OrdinalType, ScalarType>::putScalar( const ScalarType value_in, bool fullMatrix )
00554 {
00555   // Set each value of the dense matrix to "value".
00556   if (fullMatrix == true) {
00557     for(OrdinalType j = 0; j < numRowCols_; j++) 
00558       {
00559   for(OrdinalType i = 0; i < numRowCols_; i++) 
00560     {
00561       values_[i + j*stride_] = value_in;
00562     }
00563       }
00564   }
00565   // Set the active upper or lower triangular part of the matrix to "value"
00566   else {
00567     if (upper_) {
00568       for(OrdinalType j = 0; j < numRowCols_; j++) {
00569   for(OrdinalType i = 0; i <= j; i++) {
00570     values_[i + j*stride_] = value_in;
00571   }
00572       }
00573     }
00574     else {
00575       for(OrdinalType j = 0; j < numRowCols_; j++) {
00576   for(OrdinalType i = j; i < numRowCols_; i++) {
00577     values_[i + j*stride_] = value_in;
00578   }
00579       }
00580     }
00581   }
00582   return 0;
00583 }    
00584     
00585 template<typename OrdinalType, typename ScalarType>
00586 int SerialSymDenseMatrix<OrdinalType, ScalarType>::random( const ScalarType bias )
00587 {
00588   typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MT;
00589 
00590   // Set each value of the dense matrix to a random value.
00591   std::vector<MT> diagSum( numRowCols_, Teuchos::ScalarTraits<MT>::zero() );
00592   if (upper_) {
00593     for(OrdinalType j = 0; j < numRowCols_; j++) {
00594       for(OrdinalType i = 0; i < j; i++) {
00595   values_[i + j*stride_] = ScalarTraits<ScalarType>::random();
00596   diagSum[i] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
00597   diagSum[j] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
00598       }
00599     }
00600   }
00601   else {
00602     for(OrdinalType j = 0; j < numRowCols_; j++) {
00603       for(OrdinalType i = j+1; i < numRowCols_; i++) {
00604   values_[i + j*stride_] = ScalarTraits<ScalarType>::random();
00605   diagSum[i] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
00606   diagSum[j] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
00607       }
00608     }
00609   }
00610   
00611   // Set the diagonal.
00612   for(OrdinalType i = 0; i < numRowCols_; i++) {
00613     values_[i + i*stride_] = diagSum[i] + bias;
00614   }
00615   return 0;
00616 }
00617 
00618 template<typename OrdinalType, typename ScalarType>
00619 SerialSymDenseMatrix<OrdinalType,ScalarType>&
00620 SerialSymDenseMatrix<OrdinalType, ScalarType>::operator=( const SerialSymDenseMatrix<OrdinalType,ScalarType>& Source )
00621 {
00622   if(this == &Source)
00623     return(*this); // Special case of source same as target
00624   if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_)) {
00625     upper_ = Source.upper_;  // Might have to change the active part of the matrix.
00626     return(*this); // Special case of both are views to same data.
00627   }
00628   
00629   // If the source is a view then we will return a view, else we will return a copy.
00630   if (!Source.valuesCopied_) {
00631     if(valuesCopied_)       { 
00632       // Clean up stored data if this was previously a copy.
00633       deleteArrays();
00634     }
00635     numRowCols_ = Source.numRowCols_; 
00636     stride_ = Source.stride_;
00637     values_ = Source.values_;
00638     upper_ = Source.upper_;
00639     UPLO_ = Source.UPLO_;
00640   }
00641   else {
00642     // If we were a view, we will now be a copy.
00643     if(!valuesCopied_) {
00644       numRowCols_ = Source.numRowCols_;
00645       stride_ = Source.numRowCols_;
00646       upper_ = Source.upper_;
00647       UPLO_ = Source.UPLO_;
00648       const OrdinalType newsize = stride_ * numRowCols_;
00649       if(newsize > 0) {
00650         values_ = new ScalarType[newsize];
00651         valuesCopied_ = true;
00652       }
00653       else {
00654         values_ = 0;
00655       }
00656     }
00657     // If we were a copy, we will stay a copy.
00658     else {
00659       if((Source.numRowCols_ <= stride_) && (Source.numRowCols_ == numRowCols_)) { // we don't need to reallocate
00660         numRowCols_ = Source.numRowCols_;
00661   upper_ = Source.upper_;
00662   UPLO_ = Source.UPLO_;
00663       }
00664       else { // we need to allocate more space (or less space)
00665         deleteArrays();
00666         numRowCols_ = Source.numRowCols_;
00667         stride_ = Source.numRowCols_;
00668   upper_ = Source.upper_;
00669   UPLO_ = Source.UPLO_;
00670         const OrdinalType newsize = stride_ * numRowCols_;
00671         if(newsize > 0) {
00672           values_ = new ScalarType[newsize];
00673           valuesCopied_ = true;
00674         }
00675       }
00676     }
00677     copyMat(Source.upper_, Source.values_, Source.stride_, Source.numRowCols_, upper_, values_, stride_, 0);
00678   } 
00679   return(*this);
00680 }
00681 
00682 template<typename OrdinalType, typename ScalarType>
00683 SerialSymDenseMatrix<OrdinalType, ScalarType>& SerialSymDenseMatrix<OrdinalType, ScalarType>::operator*= (const ScalarType alpha)
00684 {
00685   this->scale(alpha);
00686   return(*this);
00687 }
00688 
00689 template<typename OrdinalType, typename ScalarType>
00690 SerialSymDenseMatrix<OrdinalType, ScalarType>& SerialSymDenseMatrix<OrdinalType, ScalarType>::operator+= (const SerialSymDenseMatrix<OrdinalType,ScalarType>& Source )
00691 {
00692   // Check for compatible dimensions
00693   if ((numRowCols_ != Source.numRowCols_))
00694     {
00695       TEUCHOS_CHK_REF(*this); // Return *this without altering it.
00696     }
00697   copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0, 1.0);
00698   return(*this);
00699 }
00700 
00701 template<typename OrdinalType, typename ScalarType>
00702 SerialSymDenseMatrix<OrdinalType, ScalarType>& SerialSymDenseMatrix<OrdinalType, ScalarType>::operator-= (const SerialSymDenseMatrix<OrdinalType,ScalarType>& Source )
00703 {
00704   // Check for compatible dimensions
00705   if ((numRowCols_ != Source.numRowCols_))
00706   {
00707     TEUCHOS_CHK_REF(*this); // Return *this without altering it.
00708   }
00709   copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0, -1.0);
00710   return(*this);
00711 }
00712 
00713 template<typename OrdinalType, typename ScalarType>
00714 SerialSymDenseMatrix<OrdinalType,ScalarType>& SerialSymDenseMatrix<OrdinalType, ScalarType>::assign (const SerialSymDenseMatrix<OrdinalType,ScalarType>& Source) {
00715   if(this == &Source)
00716     return(*this); // Special case of source same as target
00717   if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_)) {
00718     upper_ = Source.upper_; // We may have to change the active part of the matrix.
00719     return(*this); // Special case of both are views to same data.
00720   }
00721 
00722   // Check for compatible dimensions
00723   if ((numRowCols_ != Source.numRowCols_))
00724   {
00725     TEUCHOS_CHK_REF(*this); // Return *this without altering it.
00726   }
00727   copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0 );
00728   return(*this);
00729 }
00730 
00731 //----------------------------------------------------------------------------------------------------
00732 //   Accessor methods 
00733 //----------------------------------------------------------------------------------------------------
00734 
00735 template<typename OrdinalType, typename ScalarType>
00736 inline ScalarType& SerialSymDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex)
00737 {
00738 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
00739   checkIndex( rowIndex, colIndex );
00740 #endif
00741   if ( rowIndex <= colIndex ) {
00742     // Accessing upper triangular part of matrix
00743     if (upper_)
00744       return(values_[colIndex * stride_ + rowIndex]);
00745     else
00746       return(values_[rowIndex * stride_ + colIndex]);   
00747   }
00748   else {
00749     // Accessing lower triangular part of matrix
00750     if (upper_)
00751       return(values_[rowIndex * stride_ + colIndex]);
00752     else
00753       return(values_[colIndex * stride_ + rowIndex]);
00754   }
00755 }
00756   
00757 template<typename OrdinalType, typename ScalarType>
00758 inline const ScalarType& SerialSymDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex) const
00759 {
00760 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
00761   checkIndex( rowIndex, colIndex );
00762 #endif
00763   if ( rowIndex <= colIndex ) {
00764     // Accessing upper triangular part of matrix
00765     if (upper_)
00766       return(values_[colIndex * stride_ + rowIndex]);
00767     else
00768       return(values_[rowIndex * stride_ + colIndex]);   
00769   }
00770   else {
00771     // Accessing lower triangular part of matrix
00772     if (upper_)
00773       return(values_[rowIndex * stride_ + colIndex]);
00774     else
00775       return(values_[colIndex * stride_ + rowIndex]);
00776   }
00777 }  
00778 
00779 //----------------------------------------------------------------------------------------------------
00780 //   Norm methods 
00781 //----------------------------------------------------------------------------------------------------
00782   
00783 template<typename OrdinalType, typename ScalarType>
00784 typename ScalarTraits<ScalarType>::magnitudeType SerialSymDenseMatrix<OrdinalType, ScalarType>::normOne() const
00785 {
00786   return(normInf());
00787 }
00788   
00789 template<typename OrdinalType, typename ScalarType>
00790 typename ScalarTraits<ScalarType>::magnitudeType SerialSymDenseMatrix<OrdinalType, ScalarType>::normInf() const
00791 {
00792   typedef typename ScalarTraits<ScalarType>::magnitudeType MT;
00793 
00794   OrdinalType i, j;
00795   
00796   MT sum, anorm = ScalarTraits<MT>::zero();
00797   ScalarType* ptr;
00798   
00799   if (upper_) {
00800     for (j=0; j<numRowCols_; j++) {
00801       sum = ScalarTraits<MT>::zero();
00802       ptr = values_ + j*stride_;
00803       for (i=0; i<j; i++) {
00804         sum += ScalarTraits<ScalarType>::magnitude( *ptr++ );
00805       }
00806       ptr = values_ + j + j*stride_;
00807       for (i=j; i<numRowCols_; i++) {
00808         sum += ScalarTraits<ScalarType>::magnitude( *ptr );
00809         ptr += stride_;
00810       }
00811       anorm = TEUCHOS_MAX( anorm, sum );
00812     }
00813   }
00814   else {
00815     for (j=0; j<numRowCols_; j++) {
00816       sum = ScalarTraits<MT>::zero();
00817       ptr = values_ + j + j*stride_;
00818       for (i=j; i<numRowCols_; i++) {
00819         sum += ScalarTraits<ScalarType>::magnitude( *ptr++ );
00820       }
00821       ptr = values_ + j;
00822       for (i=0; i<j; i++) {
00823         sum += ScalarTraits<ScalarType>::magnitude( *ptr );
00824         ptr += stride_;
00825       }
00826       anorm = TEUCHOS_MAX( anorm, sum );
00827     }
00828   }
00829   return(anorm);
00830 }
00831   
00832 template<typename OrdinalType, typename ScalarType>
00833 typename ScalarTraits<ScalarType>::magnitudeType SerialSymDenseMatrix<OrdinalType, ScalarType>::normFrobenius() const
00834 {
00835   typedef typename ScalarTraits<ScalarType>::magnitudeType MT;
00836 
00837   OrdinalType i, j;
00838   
00839   MT sum = ScalarTraits<MT>::zero(), anorm = ScalarTraits<MT>::zero();
00840  
00841   if (upper_) { 
00842     for (j = 0; j < numRowCols_; j++) {
00843       for (i = 0; i < j; i++) {
00844         sum += ScalarTraits<ScalarType>::magnitude(2.0*values_[i+j*stride_]*values_[i+j*stride_]);
00845       }
00846       sum += ScalarTraits<ScalarType>::magnitude(values_[j + j*stride_]*values_[j + j*stride_]);
00847     }
00848   }
00849   else {
00850     for (j = 0; j < numRowCols_; j++) {
00851       sum += ScalarTraits<ScalarType>::magnitude(values_[j + j*stride_]*values_[j + j*stride_]);
00852       for (i = j+1; i < numRowCols_; i++) {
00853         sum += ScalarTraits<ScalarType>::magnitude(2.0*values_[i+j*stride_]*values_[i+j*stride_]);
00854       }
00855     }
00856   }       
00857   anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::squareroot(sum));
00858   return(anorm);
00859 }
00860   
00861 //----------------------------------------------------------------------------------------------------
00862 //   Comparison methods 
00863 //----------------------------------------------------------------------------------------------------
00864   
00865 template<typename OrdinalType, typename ScalarType>
00866 bool SerialSymDenseMatrix<OrdinalType, ScalarType>::operator== (const SerialSymDenseMatrix<OrdinalType, ScalarType> &Operand)
00867 {
00868   bool result = 1;
00869   if((numRowCols_ != Operand.numRowCols_)) {
00870     result = 0;
00871   }
00872   else {
00873     OrdinalType i, j;
00874     for(i = 0; i < numRowCols_; i++) {
00875       for(j = 0; j < numRowCols_; j++) {
00876   if((*this)(i, j) != Operand(i, j)) {
00877     return 0;
00878   }
00879       }
00880     }
00881   }
00882   return result;
00883 }
00884 
00885 template<typename OrdinalType, typename ScalarType>
00886 bool SerialSymDenseMatrix<OrdinalType, ScalarType>::operator!= (const SerialSymDenseMatrix<OrdinalType, ScalarType> &Operand)
00887 {
00888   return !((*this) == Operand);
00889 }
00890 
00891 //----------------------------------------------------------------------------------------------------
00892 //   Multiplication method 
00893 //----------------------------------------------------------------------------------------------------  
00894 
00895 template<typename OrdinalType, typename ScalarType>
00896 void SerialSymDenseMatrix<OrdinalType, ScalarType>::scale( const ScalarType alpha )
00897 {
00898   OrdinalType i, j;
00899   ScalarType* ptr;
00900   
00901   if (upper_) {
00902     for (j=0; j<numRowCols_; j++) {
00903       ptr = values_ + j*stride_;
00904       for (i=0; i<= j; i++) { *ptr = alpha * (*ptr); ptr++; }
00905     }
00906   }
00907   else {
00908     for (j=0; j<numRowCols_; j++) {
00909       ptr = values_ + j*stride_;
00910       for (i=j; i<numRowCols_; i++) { *ptr = alpha * (*ptr); ptr++; }
00911     }
00912   }
00913 }
00914 
00915 /*
00916 template<typename OrdinalType, typename ScalarType>
00917 int SerialSymDenseMatrix<OrdinalType, ScalarType>::scale( const SerialSymDenseMatrix<OrdinalType,ScalarType>& A )
00918 {
00919   OrdinalType i, j;
00920   ScalarType* ptr;
00921     
00922   // Check for compatible dimensions
00923   if ((numRowCols_ != A.numRowCols_)) {
00924     TEUCHOS_CHK_ERR(-1); // Return error
00925   }    
00926 
00927   if (upper_) {
00928     for (j=0; j<numRowCols_; j++) {
00929       ptr = values_ + j*stride_;
00930       for (i=0; i<= j; i++) { *ptr = A(i,j) * (*ptr); ptr++; }
00931     }
00932   }
00933   else {
00934     for (j=0; j<numRowCols_; j++) {
00935       ptr = values_ + j*stride_;
00936       for (i=j; i<numRowCols_; i++) { *ptr = A(i,j) * (*ptr); ptr++; }
00937     }
00938   }
00939 
00940   return(0);
00941 }
00942 */
00943 
00944 template<typename OrdinalType, typename ScalarType>
00945 void SerialSymDenseMatrix<OrdinalType, ScalarType>::print(std::ostream& os) const
00946 {
00947   os << std::endl;
00948   if(valuesCopied_)
00949     os << "Values_copied : yes" << std::endl;
00950   else
00951     os << "Values_copied : no" << std::endl;
00952   os << "Rows / Columns : " << numRowCols_ << std::endl;
00953   os << "LDA : " << stride_ << std::endl;
00954   if (upper_) 
00955     os << "Storage: Upper" << std::endl;
00956   else
00957     os << "Storage: Lower" << std::endl;
00958   if(numRowCols_ == 0) {
00959     os << "(matrix is empty, no values to display)" << std::endl;
00960   } else {
00961     for(OrdinalType i = 0; i < numRowCols_; i++) {
00962       for(OrdinalType j = 0; j < numRowCols_; j++){
00963         os << (*this)(i,j) << " ";
00964       }
00965       os << std::endl;
00966     }
00967   }
00968 }
00969 
00970 //----------------------------------------------------------------------------------------------------
00971 //   Protected methods 
00972 //----------------------------------------------------------------------------------------------------  
00973 
00974 template<typename OrdinalType, typename ScalarType>
00975 inline void SerialSymDenseMatrix<OrdinalType, ScalarType>::checkIndex( OrdinalType rowIndex, OrdinalType colIndex ) const {
00976   TEST_FOR_EXCEPTION(rowIndex < 0 || rowIndex >= numRowCols_, std::out_of_range,
00977     "SerialSymDenseMatrix<T>::checkIndex: "
00978     "Row index " << rowIndex << " out of range [0, "<< numRowCols_ << ")");
00979   TEST_FOR_EXCEPTION(colIndex < 0 || colIndex >= numRowCols_, std::out_of_range,
00980     "SerialSymDenseMatrix<T>::checkIndex: "
00981     "Col index " << colIndex << " out of range [0, "<< numRowCols_ << ")");
00982 }
00983 
00984 template<typename OrdinalType, typename ScalarType>
00985 void SerialSymDenseMatrix<OrdinalType, ScalarType>::deleteArrays(void)
00986 {
00987   if (valuesCopied_)
00988   {
00989     delete [] values_;
00990     values_ = 0;
00991     valuesCopied_ = false;
00992   }
00993 }
00994 
00995 template<typename OrdinalType, typename ScalarType>
00996 void SerialSymDenseMatrix<OrdinalType, ScalarType>::copyMat(
00997                   bool inputUpper, ScalarType* inputMatrix, 
00998                   OrdinalType inputStride, OrdinalType numRowCols_in, 
00999                   bool outputUpper, ScalarType* outputMatrix, 
01000                   OrdinalType outputStride, OrdinalType startRowCol, 
01001                   ScalarType alpha
01002                   )
01003 {
01004   OrdinalType i, j;
01005   ScalarType* ptr1 = 0;
01006   ScalarType* ptr2 = 0;
01007 
01008   for (j = 0; j < numRowCols_in; j++) {
01009     if (inputUpper == true) {
01010       // The input matrix is upper triangular, start at the beginning of each column.
01011       ptr2 = inputMatrix + (j + startRowCol) * inputStride + startRowCol;
01012       if (outputUpper == true) {
01013   // The output matrix matches the same pattern as the input matrix.
01014   ptr1 = outputMatrix + j*outputStride;
01015   if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
01016     for(i = 0; i <= j; i++) {
01017       *ptr1++ += alpha*(*ptr2++);
01018     }
01019   } else {
01020     for(i = 0; i <= j; i++) {
01021       *ptr1++ = *ptr2++;
01022     }
01023   }
01024       }
01025       else {
01026   // The output matrix has the opposite pattern as the input matrix.
01027   // Copy down across rows of the output matrix, but down columns of the input matrix.
01028   ptr1 = outputMatrix + j;
01029   if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
01030     for(i = 0; i <= j; i++) {
01031       *ptr1 += alpha*(*ptr2++);
01032       ptr1 += outputStride;
01033     }
01034   } else {
01035     for(i = 0; i <= j; i++) {
01036       *ptr1 = *ptr2++;
01037       ptr1 += outputStride;
01038     }
01039   }
01040       }
01041     }
01042     else {
01043       // The input matrix is lower triangular, start at the diagonal of each row.
01044       ptr2 = inputMatrix + (startRowCol+j) * inputStride + startRowCol + j;
01045       if (outputUpper == true) {
01046   // The output matrix has the opposite pattern as the input matrix.
01047   // Copy across rows of the output matrix, but down columns of the input matrix.
01048   ptr1 = outputMatrix + j*outputStride + j;
01049   if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
01050     for(i = j; i < numRowCols_in; i++) {
01051       *ptr1 += alpha*(*ptr2++);
01052       ptr1 += outputStride;
01053     }
01054   } else {
01055     for(i = j; i < numRowCols_in; i++) {
01056       *ptr1 = *ptr2++;
01057       ptr1 += outputStride;
01058     }
01059   }
01060       }
01061       else {
01062   // The output matrix matches the same pattern as the input matrix.
01063   ptr1 = outputMatrix + j*outputStride + j;
01064   if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
01065     for(i = j; i < numRowCols_in; i++) {
01066       *ptr1++ += alpha*(*ptr2++);
01067     }
01068   } else {
01069     for(i = j; i < numRowCols_in; i++) {
01070       *ptr1++ = *ptr2++;
01071     }
01072   }
01073       }
01074     }
01075   }
01076 }
01077   
01078 template<typename OrdinalType, typename ScalarType>
01079 void SerialSymDenseMatrix<OrdinalType, ScalarType>::copyUPLOMat(
01080                 bool inputUpper, ScalarType* inputMatrix, 
01081                 OrdinalType inputStride, OrdinalType inputRows 
01082                 )
01083 {
01084   OrdinalType i, j;
01085   ScalarType * ptr1 = 0;
01086   ScalarType * ptr2 = 0;
01087 
01088   if (inputUpper) {
01089     for (j=1; j<inputRows; j++) {
01090       ptr1 = inputMatrix + j;
01091       ptr2 = inputMatrix + j*inputStride;
01092       for (i=0; i<j; i++) {
01093   *ptr1 = *ptr2++;
01094   ptr1+=inputStride;
01095       }
01096     }
01097   }
01098   else {
01099     for (i=1; i<inputRows; i++) {
01100       ptr1 = inputMatrix + i;
01101       ptr2 = inputMatrix + i*inputStride;
01102       for (j=0; j<i; j++) {
01103   *ptr2++ = *ptr1;
01104   ptr1+=inputStride;
01105       }
01106     }
01107   }
01108 }
01109   
01110 } // namespace Teuchos
01111 
01112 #endif /* _TEUCHOS_SERIALSYMDENSEMATRIX_HPP_ */

Generated on Wed May 12 21:40:32 2010 for Teuchos - Trilinos Tools Package by  doxygen 1.4.7