Teuchos Package Browser (Single Doxygen Collection) Version of the Day
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 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00039 //
00040 // ***********************************************************************
00041 // @HEADER
00042 
00043 #ifndef _TEUCHOS_SERIALSYMDENSEMATRIX_HPP_
00044 #define _TEUCHOS_SERIALSYMDENSEMATRIX_HPP_
00045 
00049 #include "Teuchos_CompObject.hpp"
00050 #include "Teuchos_BLAS.hpp"
00051 #include "Teuchos_ScalarTraits.hpp"
00052 #include "Teuchos_DataAccess.hpp"
00053 #include "Teuchos_ConfigDefs.hpp"
00054 #include "Teuchos_TestForException.hpp"
00055 
00118 namespace Teuchos {
00119  
00120 template<typename OrdinalType, typename ScalarType>
00121 class SerialSymDenseMatrix : public CompObject, public Object, public BLAS<OrdinalType,ScalarType>
00122 {
00123  public:
00124 
00126   typedef OrdinalType ordinalType;
00128   typedef ScalarType scalarType;
00129 
00131 
00132 
00133 
00141   SerialSymDenseMatrix();
00142 
00144 
00154   SerialSymDenseMatrix(OrdinalType numRowsCols, bool zeroOut = true);
00155 
00157 
00168   SerialSymDenseMatrix(DataAccess CV, bool upper, ScalarType* values, OrdinalType stride, OrdinalType numRowsCols);
00169   
00171   SerialSymDenseMatrix(const SerialSymDenseMatrix<OrdinalType, ScalarType> &Source);
00172 
00174 
00183   SerialSymDenseMatrix(DataAccess CV, const SerialSymDenseMatrix<OrdinalType, ScalarType> &Source, OrdinalType numRowCols, OrdinalType startRowCol=0);
00184 
00186   virtual ~SerialSymDenseMatrix ();
00188 
00190 
00191 
00193 
00202   int shape(OrdinalType numRowsCols); 
00203 
00205 
00214   int shapeUninitialized(OrdinalType numRowsCols); 
00215 
00217 
00227   int reshape(OrdinalType numRowsCols);
00228 
00230 
00232   void setLower();
00233 
00235 
00237   void setUpper();
00238 
00240 
00242 
00243 
00245 
00251   SerialSymDenseMatrix<OrdinalType, ScalarType>& operator= (const SerialSymDenseMatrix<OrdinalType, ScalarType>& Source);
00252 
00254 
00258   SerialSymDenseMatrix<OrdinalType, ScalarType>& assign (const SerialSymDenseMatrix<OrdinalType, ScalarType>& Source);
00259 
00261 
00264   SerialSymDenseMatrix<OrdinalType, ScalarType>& operator= (const ScalarType value) { putScalar(value); return(*this); }
00265 
00267 
00272   int putScalar( const ScalarType value = Teuchos::ScalarTraits<ScalarType>::zero(), bool fullMatrix = false );
00273 
00275 
00277   int random( const ScalarType bias = 0.1*Teuchos::ScalarTraits<ScalarType>::one() );
00278 
00280 
00282 
00283 
00285 
00292   ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex);
00293 
00295 
00302   const ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex) const;
00303 
00305 
00307   ScalarType* values() const { return(values_); }
00308 
00310 
00312 
00313 
00315   bool upper() const {return(upper_);};
00316 
00318   char UPLO() const {return(UPLO_);};
00320 
00322 
00323 
00325 
00331   SerialSymDenseMatrix<OrdinalType, ScalarType>& operator*= (const ScalarType alpha);
00332 
00334 
00337   SerialSymDenseMatrix<OrdinalType, ScalarType>& operator+= (const SerialSymDenseMatrix<OrdinalType, ScalarType>& Source);
00338 
00340 
00343   SerialSymDenseMatrix<OrdinalType, ScalarType>& operator-= (const SerialSymDenseMatrix<OrdinalType, ScalarType>& Source);
00344 
00346 
00348 
00349 
00351 
00354   bool operator== (const SerialSymDenseMatrix<OrdinalType, ScalarType> &Operand) const;
00355 
00357 
00360   bool operator!= (const SerialSymDenseMatrix<OrdinalType, ScalarType> &Operand) const;
00361 
00363 
00365 
00366 
00368   OrdinalType numRows() const { return(numRowCols_); }
00369 
00371   OrdinalType numCols() const { return(numRowCols_); }
00372 
00374   OrdinalType stride() const { return(stride_); }
00375 
00377   bool empty() const { return(numRowCols_ == 0); }
00378 
00380 
00382 
00383 
00385   typename ScalarTraits<ScalarType>::magnitudeType normOne() const;
00386 
00388   typename ScalarTraits<ScalarType>::magnitudeType normInf() const;
00389 
00391   typename ScalarTraits<ScalarType>::magnitudeType normFrobenius() const; 
00393 
00395 
00396 
00397   virtual void print(std::ostream& os) const;
00398 
00400 
00401  protected:
00402 
00403   // In-place scaling of matrix.
00404   void scale( const ScalarType alpha );
00405 
00406   // Copy the values from one matrix to the other.  
00407   void copyMat(bool inputUpper, ScalarType* inputMatrix, OrdinalType inputStride,
00408          OrdinalType numRowCols, bool outputUpper, ScalarType* outputMatrix, 
00409          OrdinalType outputStride, OrdinalType startRowCol, 
00410          ScalarType alpha = ScalarTraits<ScalarType>::zero() );
00411 
00412   // Copy the values from the active triangle of the matrix to the other to make the matrix full symmetric.
00413   void copyUPLOMat(bool inputUpper, ScalarType* inputMatrix, 
00414                    OrdinalType inputStride, OrdinalType inputRows);
00415 
00416   void deleteArrays();
00417   void checkIndex( OrdinalType rowIndex, OrdinalType colIndex = 0 ) const;
00418   OrdinalType numRowCols_;
00419   OrdinalType stride_;
00420   bool valuesCopied_;
00421   ScalarType* values_;
00422   bool upper_;
00423   char UPLO_;
00424 
00425 
00426 };
00427 
00428 //----------------------------------------------------------------------------------------------------
00429 //  Constructors and Destructor
00430 //----------------------------------------------------------------------------------------------------
00431 
00432 template<typename OrdinalType, typename ScalarType>
00433 SerialSymDenseMatrix<OrdinalType, ScalarType>::SerialSymDenseMatrix()
00434   : CompObject(), Object("Teuchos::SerialSymDenseMatrix"), numRowCols_(0), stride_(0), valuesCopied_(false), values_(0), upper_(false), UPLO_('L')
00435 {}
00436 
00437 template<typename OrdinalType, typename ScalarType>
00438 SerialSymDenseMatrix<OrdinalType, ScalarType>::SerialSymDenseMatrix(OrdinalType numRowCols_in, bool zeroOut)
00439   : CompObject(), Object("Teuchos::SerialSymDenseMatrix"), numRowCols_(numRowCols_in), stride_(numRowCols_in), valuesCopied_(false), upper_(false), UPLO_('L')
00440 {
00441   values_ = new ScalarType[stride_*numRowCols_];
00442   valuesCopied_ = true;
00443   if (zeroOut == true)
00444     putScalar( Teuchos::ScalarTraits<ScalarType>::zero(), true );
00445 }
00446   
00447 template<typename OrdinalType, typename ScalarType>
00448 SerialSymDenseMatrix<OrdinalType, ScalarType>::SerialSymDenseMatrix(
00449   DataAccess CV, bool upper_in, ScalarType* values_in, OrdinalType stride_in, OrdinalType numRowCols_in
00450   )
00451   : CompObject(), Object("Teuchos::SerialSymDenseMatrix"), numRowCols_(numRowCols_in), stride_(stride_in), valuesCopied_(false), 
00452     values_(values_in), upper_(upper_in)
00453 {
00454   if (upper_)
00455     UPLO_ = 'U';
00456   else
00457     UPLO_ = 'L';
00458 
00459   if(CV == Copy)
00460   {
00461     stride_ = numRowCols_;
00462     values_ = new ScalarType[stride_*numRowCols_];
00463     copyMat(upper_in, values_in, stride_in, numRowCols_, upper_, values_, stride_, 0);
00464     valuesCopied_ = true;
00465   }
00466 }
00467   
00468 template<typename OrdinalType, typename ScalarType>
00469 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_) 
00470 {
00471   values_ = new ScalarType[stride_*numRowCols_];
00472   copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0);
00473   valuesCopied_ = true;
00474 }
00475   
00476 template<typename OrdinalType, typename ScalarType>
00477 SerialSymDenseMatrix<OrdinalType, ScalarType>::SerialSymDenseMatrix(
00478                     DataAccess CV, const SerialSymDenseMatrix<OrdinalType, 
00479                     ScalarType> &Source, OrdinalType numRowCols_in, OrdinalType startRowCol )
00480   : CompObject(), Object("Teuchos::SerialSymDenseMatrix"), numRowCols_(numRowCols_in), stride_(Source.stride_), valuesCopied_(false), upper_(Source.upper_), UPLO_(Source.UPLO_)
00481 {
00482   if(CV == Copy)
00483   {
00484     stride_ = numRowCols_in;
00485     values_ = new ScalarType[stride_ * numRowCols_in];
00486     copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_in, upper_, values_, stride_, startRowCol);
00487     valuesCopied_ = true;
00488   }
00489   else // CV == View
00490   {
00491     values_ = Source.values_ + (stride_ * startRowCol) + startRowCol;
00492   }
00493 }
00494     
00495 template<typename OrdinalType, typename ScalarType>
00496 SerialSymDenseMatrix<OrdinalType, ScalarType>::~SerialSymDenseMatrix()
00497 {
00498   deleteArrays();
00499 }
00500 
00501 //----------------------------------------------------------------------------------------------------
00502 //  Shape methods 
00503 //----------------------------------------------------------------------------------------------------
00504 
00505 template<typename OrdinalType, typename ScalarType>
00506 int SerialSymDenseMatrix<OrdinalType, ScalarType>::shape( OrdinalType numRowCols_in )
00507 {
00508   deleteArrays(); // Get rid of anything that might be already allocated
00509   numRowCols_ = numRowCols_in;
00510   stride_ = numRowCols_;
00511   values_ = new ScalarType[stride_*numRowCols_];
00512   putScalar( Teuchos::ScalarTraits<ScalarType>::zero(), true );
00513   valuesCopied_ = true;
00514   return(0);
00515 }
00516 
00517 template<typename OrdinalType, typename ScalarType>
00518 int SerialSymDenseMatrix<OrdinalType, ScalarType>::shapeUninitialized( OrdinalType numRowCols_in )
00519 {
00520   deleteArrays(); // Get rid of anything that might be already allocated
00521   numRowCols_ = numRowCols_in;
00522   stride_ = numRowCols_;
00523   values_ = new ScalarType[stride_*numRowCols_];
00524   valuesCopied_ = true;
00525   return(0);
00526 }
00527   
00528 template<typename OrdinalType, typename ScalarType>
00529 int SerialSymDenseMatrix<OrdinalType, ScalarType>::reshape( OrdinalType numRowCols_in )
00530 {
00531   // Allocate space for new matrix
00532   ScalarType* values_tmp = new ScalarType[numRowCols_in * numRowCols_in];
00533   ScalarType zero = ScalarTraits<ScalarType>::zero();
00534   for(OrdinalType k = 0; k < numRowCols_in * numRowCols_in; k++)
00535   {
00536     values_tmp[k] = zero;
00537   }
00538   OrdinalType numRowCols_tmp = TEUCHOS_MIN(numRowCols_, numRowCols_in);
00539   if(values_ != 0)
00540   {
00541     copyMat(upper_, values_, stride_, numRowCols_tmp, upper_, values_tmp, numRowCols_in, 0); // Copy principal submatrix of A to new A
00542   }
00543   deleteArrays(); // Get rid of anything that might be already allocated
00544   numRowCols_ = numRowCols_in;
00545   stride_ = numRowCols_;
00546   values_ = values_tmp; // Set pointer to new A
00547   valuesCopied_ = true;
00548   return(0);
00549 }
00550    
00551 //----------------------------------------------------------------------------------------------------
00552 //   Set methods 
00553 //----------------------------------------------------------------------------------------------------
00554 
00555 template<typename OrdinalType, typename ScalarType>
00556 void SerialSymDenseMatrix<OrdinalType, ScalarType>::setLower() 
00557 {
00558   // Do nothing if the matrix is already a lower triangular matrix
00559   if (upper_ != false) {
00560     copyUPLOMat( true, values_, stride_, numRowCols_ );
00561     upper_ = false;
00562     UPLO_ = 'L';
00563   }
00564 }
00565 
00566 template<typename OrdinalType, typename ScalarType>
00567 void SerialSymDenseMatrix<OrdinalType, ScalarType>::setUpper() 
00568 {
00569   // Do nothing if the matrix is already an upper triangular matrix
00570   if (upper_ == false) {
00571     copyUPLOMat( false, values_, stride_, numRowCols_ );
00572     upper_ = true;
00573     UPLO_ = 'U';
00574   }
00575 }
00576 
00577 template<typename OrdinalType, typename ScalarType>
00578 int SerialSymDenseMatrix<OrdinalType, ScalarType>::putScalar( const ScalarType value_in, bool fullMatrix )
00579 {
00580   // Set each value of the dense matrix to "value".
00581   if (fullMatrix == true) {
00582     for(OrdinalType j = 0; j < numRowCols_; j++) 
00583       {
00584   for(OrdinalType i = 0; i < numRowCols_; i++) 
00585     {
00586       values_[i + j*stride_] = value_in;
00587     }
00588       }
00589   }
00590   // Set the active upper or lower triangular part of the matrix to "value"
00591   else {
00592     if (upper_) {
00593       for(OrdinalType j = 0; j < numRowCols_; j++) {
00594   for(OrdinalType i = 0; i <= j; i++) {
00595     values_[i + j*stride_] = value_in;
00596   }
00597       }
00598     }
00599     else {
00600       for(OrdinalType j = 0; j < numRowCols_; j++) {
00601   for(OrdinalType i = j; i < numRowCols_; i++) {
00602     values_[i + j*stride_] = value_in;
00603   }
00604       }
00605     }
00606   }
00607   return 0;
00608 }    
00609     
00610 template<typename OrdinalType, typename ScalarType>
00611 int SerialSymDenseMatrix<OrdinalType, ScalarType>::random( const ScalarType bias )
00612 {
00613   typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MT;
00614 
00615   // Set each value of the dense matrix to a random value.
00616   std::vector<MT> diagSum( numRowCols_, Teuchos::ScalarTraits<MT>::zero() );
00617   if (upper_) {
00618     for(OrdinalType j = 0; j < numRowCols_; j++) {
00619       for(OrdinalType i = 0; i < j; i++) {
00620   values_[i + j*stride_] = ScalarTraits<ScalarType>::random();
00621   diagSum[i] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
00622   diagSum[j] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
00623       }
00624     }
00625   }
00626   else {
00627     for(OrdinalType j = 0; j < numRowCols_; j++) {
00628       for(OrdinalType i = j+1; i < numRowCols_; i++) {
00629   values_[i + j*stride_] = ScalarTraits<ScalarType>::random();
00630   diagSum[i] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
00631   diagSum[j] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
00632       }
00633     }
00634   }
00635   
00636   // Set the diagonal.
00637   for(OrdinalType i = 0; i < numRowCols_; i++) {
00638     values_[i + i*stride_] = diagSum[i] + bias;
00639   }
00640   return 0;
00641 }
00642 
00643 template<typename OrdinalType, typename ScalarType>
00644 SerialSymDenseMatrix<OrdinalType,ScalarType>&
00645 SerialSymDenseMatrix<OrdinalType, ScalarType>::operator=( const SerialSymDenseMatrix<OrdinalType,ScalarType>& Source )
00646 {
00647   if(this == &Source)
00648     return(*this); // Special case of source same as target
00649   if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_)) {
00650     upper_ = Source.upper_;  // Might have to change the active part of the matrix.
00651     return(*this); // Special case of both are views to same data.
00652   }
00653   
00654   // If the source is a view then we will return a view, else we will return a copy.
00655   if (!Source.valuesCopied_) {
00656     if(valuesCopied_)       { 
00657       // Clean up stored data if this was previously a copy.
00658       deleteArrays();
00659     }
00660     numRowCols_ = Source.numRowCols_; 
00661     stride_ = Source.stride_;
00662     values_ = Source.values_;
00663     upper_ = Source.upper_;
00664     UPLO_ = Source.UPLO_;
00665   }
00666   else {
00667     // If we were a view, we will now be a copy.
00668     if(!valuesCopied_) {
00669       numRowCols_ = Source.numRowCols_;
00670       stride_ = Source.numRowCols_;
00671       upper_ = Source.upper_;
00672       UPLO_ = Source.UPLO_;
00673       const OrdinalType newsize = stride_ * numRowCols_;
00674       if(newsize > 0) {
00675         values_ = new ScalarType[newsize];
00676         valuesCopied_ = true;
00677       }
00678       else {
00679         values_ = 0;
00680       }
00681     }
00682     // If we were a copy, we will stay a copy.
00683     else {
00684       if((Source.numRowCols_ <= stride_) && (Source.numRowCols_ == numRowCols_)) { // we don't need to reallocate
00685         numRowCols_ = Source.numRowCols_;
00686   upper_ = Source.upper_;
00687   UPLO_ = Source.UPLO_;
00688       }
00689       else { // we need to allocate more space (or less space)
00690         deleteArrays();
00691         numRowCols_ = Source.numRowCols_;
00692         stride_ = Source.numRowCols_;
00693   upper_ = Source.upper_;
00694   UPLO_ = Source.UPLO_;
00695         const OrdinalType newsize = stride_ * numRowCols_;
00696         if(newsize > 0) {
00697           values_ = new ScalarType[newsize];
00698           valuesCopied_ = true;
00699         }
00700       }
00701     }
00702     copyMat(Source.upper_, Source.values_, Source.stride_, Source.numRowCols_, upper_, values_, stride_, 0);
00703   } 
00704   return(*this);
00705 }
00706 
00707 template<typename OrdinalType, typename ScalarType>
00708 SerialSymDenseMatrix<OrdinalType, ScalarType>& SerialSymDenseMatrix<OrdinalType, ScalarType>::operator*= (const ScalarType alpha)
00709 {
00710   this->scale(alpha);
00711   return(*this);
00712 }
00713 
00714 template<typename OrdinalType, typename ScalarType>
00715 SerialSymDenseMatrix<OrdinalType, ScalarType>& SerialSymDenseMatrix<OrdinalType, ScalarType>::operator+= (const SerialSymDenseMatrix<OrdinalType,ScalarType>& Source )
00716 {
00717   // Check for compatible dimensions
00718   if ((numRowCols_ != Source.numRowCols_))
00719     {
00720       TEUCHOS_CHK_REF(*this); // Return *this without altering it.
00721     }
00722   copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0, ScalarTraits<ScalarType>::one());
00723   return(*this);
00724 }
00725 
00726 template<typename OrdinalType, typename ScalarType>
00727 SerialSymDenseMatrix<OrdinalType, ScalarType>& SerialSymDenseMatrix<OrdinalType, ScalarType>::operator-= (const SerialSymDenseMatrix<OrdinalType,ScalarType>& Source )
00728 {
00729   // Check for compatible dimensions
00730   if ((numRowCols_ != Source.numRowCols_))
00731   {
00732     TEUCHOS_CHK_REF(*this); // Return *this without altering it.
00733   }
00734   copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0, -ScalarTraits<ScalarType>::one());
00735   return(*this);
00736 }
00737 
00738 template<typename OrdinalType, typename ScalarType>
00739 SerialSymDenseMatrix<OrdinalType,ScalarType>& SerialSymDenseMatrix<OrdinalType, ScalarType>::assign (const SerialSymDenseMatrix<OrdinalType,ScalarType>& Source) {
00740   if(this == &Source)
00741     return(*this); // Special case of source same as target
00742   if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_)) {
00743     upper_ = Source.upper_; // We may have to change the active part of the matrix.
00744     return(*this); // Special case of both are views to same data.
00745   }
00746 
00747   // Check for compatible dimensions
00748   if ((numRowCols_ != Source.numRowCols_))
00749   {
00750     TEUCHOS_CHK_REF(*this); // Return *this without altering it.
00751   }
00752   copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0 );
00753   return(*this);
00754 }
00755 
00756 //----------------------------------------------------------------------------------------------------
00757 //   Accessor methods 
00758 //----------------------------------------------------------------------------------------------------
00759 
00760 template<typename OrdinalType, typename ScalarType>
00761 inline ScalarType& SerialSymDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex)
00762 {
00763 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
00764   checkIndex( rowIndex, colIndex );
00765 #endif
00766   if ( rowIndex <= colIndex ) {
00767     // Accessing upper triangular part of matrix
00768     if (upper_)
00769       return(values_[colIndex * stride_ + rowIndex]);
00770     else
00771       return(values_[rowIndex * stride_ + colIndex]);   
00772   }
00773   else {
00774     // Accessing lower triangular part of matrix
00775     if (upper_)
00776       return(values_[rowIndex * stride_ + colIndex]);
00777     else
00778       return(values_[colIndex * stride_ + rowIndex]);
00779   }
00780 }
00781   
00782 template<typename OrdinalType, typename ScalarType>
00783 inline const ScalarType& SerialSymDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex) const
00784 {
00785 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
00786   checkIndex( rowIndex, colIndex );
00787 #endif
00788   if ( rowIndex <= colIndex ) {
00789     // Accessing upper triangular part of matrix
00790     if (upper_)
00791       return(values_[colIndex * stride_ + rowIndex]);
00792     else
00793       return(values_[rowIndex * stride_ + colIndex]);   
00794   }
00795   else {
00796     // Accessing lower triangular part of matrix
00797     if (upper_)
00798       return(values_[rowIndex * stride_ + colIndex]);
00799     else
00800       return(values_[colIndex * stride_ + rowIndex]);
00801   }
00802 }  
00803 
00804 //----------------------------------------------------------------------------------------------------
00805 //   Norm methods 
00806 //----------------------------------------------------------------------------------------------------
00807   
00808 template<typename OrdinalType, typename ScalarType>
00809 typename ScalarTraits<ScalarType>::magnitudeType SerialSymDenseMatrix<OrdinalType, ScalarType>::normOne() const
00810 {
00811   return(normInf());
00812 }
00813   
00814 template<typename OrdinalType, typename ScalarType>
00815 typename ScalarTraits<ScalarType>::magnitudeType SerialSymDenseMatrix<OrdinalType, ScalarType>::normInf() const
00816 {
00817   typedef typename ScalarTraits<ScalarType>::magnitudeType MT;
00818 
00819   OrdinalType i, j;
00820   
00821   MT sum, anorm = ScalarTraits<MT>::zero();
00822   ScalarType* ptr;
00823   
00824   if (upper_) {
00825     for (j=0; j<numRowCols_; j++) {
00826       sum = ScalarTraits<MT>::zero();
00827       ptr = values_ + j*stride_;
00828       for (i=0; i<j; i++) {
00829         sum += ScalarTraits<ScalarType>::magnitude( *ptr++ );
00830       }
00831       ptr = values_ + j + j*stride_;
00832       for (i=j; i<numRowCols_; i++) {
00833         sum += ScalarTraits<ScalarType>::magnitude( *ptr );
00834         ptr += stride_;
00835       }
00836       anorm = TEUCHOS_MAX( anorm, sum );
00837     }
00838   }
00839   else {
00840     for (j=0; j<numRowCols_; j++) {
00841       sum = ScalarTraits<MT>::zero();
00842       ptr = values_ + j + j*stride_;
00843       for (i=j; i<numRowCols_; i++) {
00844         sum += ScalarTraits<ScalarType>::magnitude( *ptr++ );
00845       }
00846       ptr = values_ + j;
00847       for (i=0; i<j; i++) {
00848         sum += ScalarTraits<ScalarType>::magnitude( *ptr );
00849         ptr += stride_;
00850       }
00851       anorm = TEUCHOS_MAX( anorm, sum );
00852     }
00853   }
00854   return(anorm);
00855 }
00856   
00857 template<typename OrdinalType, typename ScalarType>
00858 typename ScalarTraits<ScalarType>::magnitudeType SerialSymDenseMatrix<OrdinalType, ScalarType>::normFrobenius() const
00859 {
00860   typedef typename ScalarTraits<ScalarType>::magnitudeType MT;
00861 
00862   OrdinalType i, j;
00863   
00864   MT sum = ScalarTraits<MT>::zero(), anorm = ScalarTraits<MT>::zero();
00865  
00866   if (upper_) { 
00867     for (j = 0; j < numRowCols_; j++) {
00868       for (i = 0; i < j; i++) {
00869         sum += ScalarTraits<ScalarType>::magnitude(2.0*values_[i+j*stride_]*values_[i+j*stride_]);
00870       }
00871       sum += ScalarTraits<ScalarType>::magnitude(values_[j + j*stride_]*values_[j + j*stride_]);
00872     }
00873   }
00874   else {
00875     for (j = 0; j < numRowCols_; j++) {
00876       sum += ScalarTraits<ScalarType>::magnitude(values_[j + j*stride_]*values_[j + j*stride_]);
00877       for (i = j+1; i < numRowCols_; i++) {
00878         sum += ScalarTraits<ScalarType>::magnitude(2.0*values_[i+j*stride_]*values_[i+j*stride_]);
00879       }
00880     }
00881   }       
00882   anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::squareroot(sum));
00883   return(anorm);
00884 }
00885   
00886 //----------------------------------------------------------------------------------------------------
00887 //   Comparison methods 
00888 //----------------------------------------------------------------------------------------------------
00889   
00890 template<typename OrdinalType, typename ScalarType>
00891 bool SerialSymDenseMatrix<OrdinalType, ScalarType>::operator== (const SerialSymDenseMatrix<OrdinalType, ScalarType> &Operand) const
00892 {
00893   bool result = 1;
00894   if((numRowCols_ != Operand.numRowCols_)) {
00895     result = 0;
00896   }
00897   else {
00898     OrdinalType i, j;
00899     for(i = 0; i < numRowCols_; i++) {
00900       for(j = 0; j < numRowCols_; j++) {
00901   if((*this)(i, j) != Operand(i, j)) {
00902     return 0;
00903   }
00904       }
00905     }
00906   }
00907   return result;
00908 }
00909 
00910 template<typename OrdinalType, typename ScalarType>
00911 bool SerialSymDenseMatrix<OrdinalType, ScalarType>::operator!= (const SerialSymDenseMatrix<OrdinalType, ScalarType> &Operand) const
00912 {
00913   return !((*this) == Operand);
00914 }
00915 
00916 //----------------------------------------------------------------------------------------------------
00917 //   Multiplication method 
00918 //----------------------------------------------------------------------------------------------------  
00919 
00920 template<typename OrdinalType, typename ScalarType>
00921 void SerialSymDenseMatrix<OrdinalType, ScalarType>::scale( const ScalarType alpha )
00922 {
00923   OrdinalType i, j;
00924   ScalarType* ptr;
00925   
00926   if (upper_) {
00927     for (j=0; j<numRowCols_; j++) {
00928       ptr = values_ + j*stride_;
00929       for (i=0; i<= j; i++) { *ptr = alpha * (*ptr); ptr++; }
00930     }
00931   }
00932   else {
00933     for (j=0; j<numRowCols_; j++) {
00934       ptr = values_ + j*stride_ + j;
00935       for (i=j; i<numRowCols_; i++) { *ptr = alpha * (*ptr); ptr++; }
00936     }
00937   }
00938 }
00939 
00940 /*
00941 template<typename OrdinalType, typename ScalarType>
00942 int SerialSymDenseMatrix<OrdinalType, ScalarType>::scale( const SerialSymDenseMatrix<OrdinalType,ScalarType>& A )
00943 {
00944   OrdinalType i, j;
00945   ScalarType* ptr;
00946     
00947   // Check for compatible dimensions
00948   if ((numRowCols_ != A.numRowCols_)) {
00949     TEUCHOS_CHK_ERR(-1); // Return error
00950   }    
00951 
00952   if (upper_) {
00953     for (j=0; j<numRowCols_; j++) {
00954       ptr = values_ + j*stride_;
00955       for (i=0; i<= j; i++) { *ptr = A(i,j) * (*ptr); ptr++; }
00956     }
00957   }
00958   else {
00959     for (j=0; j<numRowCols_; j++) {
00960       ptr = values_ + j*stride_;
00961       for (i=j; i<numRowCols_; i++) { *ptr = A(i,j) * (*ptr); ptr++; }
00962     }
00963   }
00964 
00965   return(0);
00966 }
00967 */
00968 
00969 template<typename OrdinalType, typename ScalarType>
00970 void SerialSymDenseMatrix<OrdinalType, ScalarType>::print(std::ostream& os) const
00971 {
00972   os << std::endl;
00973   if(valuesCopied_)
00974     os << "Values_copied : yes" << std::endl;
00975   else
00976     os << "Values_copied : no" << std::endl;
00977   os << "Rows / Columns : " << numRowCols_ << std::endl;
00978   os << "LDA : " << stride_ << std::endl;
00979   if (upper_) 
00980     os << "Storage: Upper" << std::endl;
00981   else
00982     os << "Storage: Lower" << std::endl;
00983   if(numRowCols_ == 0) {
00984     os << "(matrix is empty, no values to display)" << std::endl;
00985   } else {
00986     for(OrdinalType i = 0; i < numRowCols_; i++) {
00987       for(OrdinalType j = 0; j < numRowCols_; j++){
00988         os << (*this)(i,j) << " ";
00989       }
00990       os << std::endl;
00991     }
00992   }
00993 }
00994 
00995 //----------------------------------------------------------------------------------------------------
00996 //   Protected methods 
00997 //----------------------------------------------------------------------------------------------------  
00998 
00999 template<typename OrdinalType, typename ScalarType>
01000 inline void SerialSymDenseMatrix<OrdinalType, ScalarType>::checkIndex( OrdinalType rowIndex, OrdinalType colIndex ) const {
01001   TEST_FOR_EXCEPTION(rowIndex < 0 || rowIndex >= numRowCols_, std::out_of_range,
01002     "SerialSymDenseMatrix<T>::checkIndex: "
01003     "Row index " << rowIndex << " out of range [0, "<< numRowCols_ << ")");
01004   TEST_FOR_EXCEPTION(colIndex < 0 || colIndex >= numRowCols_, std::out_of_range,
01005     "SerialSymDenseMatrix<T>::checkIndex: "
01006     "Col index " << colIndex << " out of range [0, "<< numRowCols_ << ")");
01007 }
01008 
01009 template<typename OrdinalType, typename ScalarType>
01010 void SerialSymDenseMatrix<OrdinalType, ScalarType>::deleteArrays(void)
01011 {
01012   if (valuesCopied_)
01013   {
01014     delete [] values_;
01015     values_ = 0;
01016     valuesCopied_ = false;
01017   }
01018 }
01019 
01020 template<typename OrdinalType, typename ScalarType>
01021 void SerialSymDenseMatrix<OrdinalType, ScalarType>::copyMat(
01022                   bool inputUpper, ScalarType* inputMatrix, 
01023                   OrdinalType inputStride, OrdinalType numRowCols_in, 
01024                   bool outputUpper, ScalarType* outputMatrix, 
01025                   OrdinalType outputStride, OrdinalType startRowCol, 
01026                   ScalarType alpha
01027                   )
01028 {
01029   OrdinalType i, j;
01030   ScalarType* ptr1 = 0;
01031   ScalarType* ptr2 = 0;
01032 
01033   for (j = 0; j < numRowCols_in; j++) {
01034     if (inputUpper == true) {
01035       // The input matrix is upper triangular, start at the beginning of each column.
01036       ptr2 = inputMatrix + (j + startRowCol) * inputStride + startRowCol;
01037       if (outputUpper == true) {
01038   // The output matrix matches the same pattern as the input matrix.
01039   ptr1 = outputMatrix + j*outputStride;
01040   if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
01041     for(i = 0; i <= j; i++) {
01042       *ptr1++ += alpha*(*ptr2++);
01043     }
01044   } else {
01045     for(i = 0; i <= j; i++) {
01046       *ptr1++ = *ptr2++;
01047     }
01048   }
01049       }
01050       else {
01051   // The output matrix has the opposite pattern as the input matrix.
01052   // Copy down across rows of the output matrix, but down columns of the input matrix.
01053   ptr1 = outputMatrix + j;
01054   if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
01055     for(i = 0; i <= j; i++) {
01056       *ptr1 += alpha*(*ptr2++);
01057       ptr1 += outputStride;
01058     }
01059   } else {
01060     for(i = 0; i <= j; i++) {
01061       *ptr1 = *ptr2++;
01062       ptr1 += outputStride;
01063     }
01064   }
01065       }
01066     }
01067     else {
01068       // The input matrix is lower triangular, start at the diagonal of each row.
01069       ptr2 = inputMatrix + (startRowCol+j) * inputStride + startRowCol + j;
01070       if (outputUpper == true) {
01071   // The output matrix has the opposite pattern as the input matrix.
01072   // Copy across rows of the output matrix, but down columns of the input matrix.
01073   ptr1 = outputMatrix + j*outputStride + j;
01074   if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
01075     for(i = j; i < numRowCols_in; i++) {
01076       *ptr1 += alpha*(*ptr2++);
01077       ptr1 += outputStride;
01078     }
01079   } else {
01080     for(i = j; i < numRowCols_in; i++) {
01081       *ptr1 = *ptr2++;
01082       ptr1 += outputStride;
01083     }
01084   }
01085       }
01086       else {
01087   // The output matrix matches the same pattern as the input matrix.
01088   ptr1 = outputMatrix + j*outputStride + j;
01089   if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
01090     for(i = j; i < numRowCols_in; i++) {
01091       *ptr1++ += alpha*(*ptr2++);
01092     }
01093   } else {
01094     for(i = j; i < numRowCols_in; i++) {
01095       *ptr1++ = *ptr2++;
01096     }
01097   }
01098       }
01099     }
01100   }
01101 }
01102   
01103 template<typename OrdinalType, typename ScalarType>
01104 void SerialSymDenseMatrix<OrdinalType, ScalarType>::copyUPLOMat(
01105                 bool inputUpper, ScalarType* inputMatrix, 
01106                 OrdinalType inputStride, OrdinalType inputRows 
01107                 )
01108 {
01109   OrdinalType i, j;
01110   ScalarType * ptr1 = 0;
01111   ScalarType * ptr2 = 0;
01112 
01113   if (inputUpper) {
01114     for (j=1; j<inputRows; j++) {
01115       ptr1 = inputMatrix + j;
01116       ptr2 = inputMatrix + j*inputStride;
01117       for (i=0; i<j; i++) {
01118   *ptr1 = *ptr2++;
01119   ptr1+=inputStride;
01120       }
01121     }
01122   }
01123   else {
01124     for (i=1; i<inputRows; i++) {
01125       ptr1 = inputMatrix + i;
01126       ptr2 = inputMatrix + i*inputStride;
01127       for (j=0; j<i; j++) {
01128   *ptr2++ = *ptr1;
01129   ptr1+=inputStride;
01130       }
01131     }
01132   }
01133 }
01134   
01135 } // namespace Teuchos
01136 
01137 #endif /* _TEUCHOS_SERIALSYMDENSEMATRIX_HPP_ */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines