Teuchos_SerialDenseMatrix.hpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //                    Teuchos: Common Tools Package
00005 //                 Copyright (2004) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 // Kris
00030 // 06.18.03 -- Removed comments/documentation; file too hard to edit otherwise. Will replace later.
00031 //          -- Begin conversion from <ScalarType> template to <OrdinalType, ScalarType>
00032 // 06.23.03 -- Finished conversion from <ScalarType> to <OrdinalType, ScalarType>
00033 //          -- Tpetra_DenseMatrix.cpp is now obsolete
00034 //          -- Added new constructor to allow construction of a submatrix
00035 //          -- Altered copyMat to enable its use in new constructor
00036 //          -- Commented out broken print() function
00037 //          -- Fixed oneNorm() (uninitialized return variable was causing erroneous results)
00038 // 06.24.03 -- Minor formatting changes
00039 // 07.01.03 -- Added TempPrint() function to temporarily take the place of print() and operator<< while I figure out how to fix them
00040 // 07.02.03 -- Added operator== and operator!= to make testing programs easier to write/read. Implementation of == isn't the most
00041 //             efficient/robust, but it works. Will consider optimizing later.
00042 //          -- Warning! Constructor DenseMatrix(DataAccess, const DenseMatrix<OrdinalType, ScalarType> &, int, int, int, int) (the
00043 //             "submatrix grabber" constructor) does not work correctly when used with CV == View (always grabs submatrix from top
00044 //             left corner).
00045 // 07.07.03 -- Constructor bug detailed above (07.02) is now corrected (hopefully).
00046 // 07.08.03 -- Move into Teuchos package/namespace
00047 
00048 #ifndef _TEUCHOS_SERIALDENSEMATRIX_HPP_
00049 #define _TEUCHOS_SERIALDENSEMATRIX_HPP_
00050 
00054 #include "Teuchos_CompObject.hpp"
00055 #include "Teuchos_BLAS.hpp"
00056 #include "Teuchos_ScalarTraits.hpp"
00057 #include "Teuchos_DataAccess.hpp"
00058 #include "Teuchos_ConfigDefs.hpp"
00059 #include "Teuchos_TestForException.hpp"
00060 #include "Teuchos_SerialSymDenseMatrix.hpp"
00061 
00070 namespace Teuchos {  
00071 
00072 template<typename OrdinalType, typename ScalarType>
00073 class SerialDenseMatrix : public CompObject, public Object, public BLAS<OrdinalType, ScalarType>
00074 {
00075 public:
00076 
00078   typedef OrdinalType ordinalType;
00080   typedef ScalarType scalarType;
00081 
00083 
00084 
00086 
00089   SerialDenseMatrix();
00090 
00092 
00100   SerialDenseMatrix(OrdinalType numRows, OrdinalType numCols, bool zeroOut = true);
00101 
00103 
00111   SerialDenseMatrix(DataAccess CV, ScalarType* values, OrdinalType stride, OrdinalType numRows, OrdinalType numCols);
00112 
00114 
00119   SerialDenseMatrix(const SerialDenseMatrix<OrdinalType, ScalarType> &Source, ETransp trans = Teuchos::NO_TRANS);
00120 
00122 
00134   SerialDenseMatrix(DataAccess CV, const SerialDenseMatrix<OrdinalType, ScalarType> &Source, OrdinalType numRows, OrdinalType numCols, OrdinalType startRow=0, OrdinalType startCol=0);
00135 
00137   virtual ~SerialDenseMatrix();
00139 
00141 
00142 
00143 
00153   int shape(OrdinalType numRows, OrdinalType numCols);
00154 
00156   int shapeUninitialized(OrdinalType numRows, OrdinalType numCols);
00157 
00159 
00169   int reshape(OrdinalType numRows, OrdinalType numCols);
00170 
00171 
00173 
00175 
00176 
00178 
00184   SerialDenseMatrix<OrdinalType, ScalarType>& operator= (const SerialDenseMatrix<OrdinalType, ScalarType>& Source);
00185 
00187 
00192   SerialDenseMatrix<OrdinalType, ScalarType>& assign (const SerialDenseMatrix<OrdinalType, ScalarType>& Source);
00193 
00195 
00198   SerialDenseMatrix<OrdinalType, ScalarType>& operator= (const ScalarType value) { putScalar(value); return(*this); }
00199 
00201 
00205   int putScalar( const ScalarType value = Teuchos::ScalarTraits<ScalarType>::zero() );
00206 
00208   int random();
00209 
00211 
00213 
00214 
00216 
00223   ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex);
00224 
00226 
00233   const ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex) const;
00234 
00236 
00243   ScalarType* operator [] (OrdinalType colIndex);
00244 
00246 
00253   const ScalarType* operator [] (OrdinalType colIndex) const;
00254 
00256 
00257   ScalarType* values() const { return(values_); }
00258 
00260 
00262 
00263 
00265 
00268   SerialDenseMatrix<OrdinalType, ScalarType>& operator+= (const SerialDenseMatrix<OrdinalType, ScalarType>& Source);
00269 
00271 
00274   SerialDenseMatrix<OrdinalType, ScalarType>& operator-= (const SerialDenseMatrix<OrdinalType, ScalarType>& Source);
00275 
00277 
00280   SerialDenseMatrix<OrdinalType, ScalarType>& operator*= (const ScalarType alpha);
00281 
00283 
00287   int scale ( const ScalarType alpha );
00288 
00290 
00296   int scale ( const SerialDenseMatrix<OrdinalType, ScalarType>& A );
00297 
00299 
00313   int multiply (ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix<OrdinalType, ScalarType> &A, const SerialDenseMatrix<OrdinalType, ScalarType> &B, ScalarType beta);
00314 
00316 
00327   int multiply (ESide sideA, ScalarType alpha, const SerialSymDenseMatrix<OrdinalType, ScalarType> &A, const SerialDenseMatrix<OrdinalType, ScalarType> &B, ScalarType beta);
00328 
00330 
00332 
00333 
00335 
00338   bool operator== (const SerialDenseMatrix<OrdinalType, ScalarType> &Operand) const;
00339 
00341 
00344   bool operator!= (const SerialDenseMatrix<OrdinalType, ScalarType> &Operand) const;
00345 
00347 
00349 
00350 
00352   OrdinalType numRows() const { return(numRows_); }
00353 
00355   OrdinalType numCols() const { return(numCols_); }
00356 
00358   OrdinalType stride() const { return(stride_); }
00359 
00361   bool empty() const { return(numRows_ == 0 || numCols_ == 0); }
00363 
00365 
00366 
00368   typename ScalarTraits<ScalarType>::magnitudeType normOne() const;
00369 
00371   typename ScalarTraits<ScalarType>::magnitudeType normInf() const;
00372 
00374   typename ScalarTraits<ScalarType>::magnitudeType normFrobenius() const; 
00376 
00378 
00379 
00380   virtual void print(std::ostream& os) const;
00381 
00383 protected:
00384   void copyMat(ScalarType* inputMatrix, OrdinalType strideInput,
00385     OrdinalType numRows, OrdinalType numCols, ScalarType* outputMatrix,
00386     OrdinalType strideOutput, OrdinalType startRow, OrdinalType startCol,
00387     ScalarType alpha = ScalarTraits<ScalarType>::zero() );
00388   void deleteArrays();
00389   void checkIndex( OrdinalType rowIndex, OrdinalType colIndex = 0 ) const;
00390   OrdinalType numRows_;
00391   OrdinalType numCols_;
00392   OrdinalType stride_;
00393   bool valuesCopied_;
00394   ScalarType* values_;
00395 
00396 }; // class Teuchos_SerialDenseMatrix
00397 
00398 //----------------------------------------------------------------------------------------------------
00399 //  Constructors and Destructor
00400 //----------------------------------------------------------------------------------------------------
00401 
00402 template<typename OrdinalType, typename ScalarType>
00403 SerialDenseMatrix<OrdinalType, ScalarType>::SerialDenseMatrix()
00404   : CompObject(), Object("Teuchos::SerialDenseMatrix"), numRows_(0), numCols_(0), stride_(0), valuesCopied_(false), values_(0)
00405 {}
00406   
00407 template<typename OrdinalType, typename ScalarType>
00408 SerialDenseMatrix<OrdinalType, ScalarType>::SerialDenseMatrix(
00409   OrdinalType numRows_in, OrdinalType numCols_in, bool zeroOut
00410   )
00411   : CompObject(), Object("Teuchos::SerialDenseMatrix"), numRows_(numRows_in), numCols_(numCols_in), stride_(numRows_in)
00412 {
00413   values_ = new ScalarType[stride_*numCols_];
00414   valuesCopied_ = true;
00415   if (zeroOut == true)  
00416     putScalar();
00417 }
00418 
00419 template<typename OrdinalType, typename ScalarType>
00420 SerialDenseMatrix<OrdinalType, ScalarType>::SerialDenseMatrix(
00421   DataAccess CV, ScalarType* values_in, OrdinalType stride_in, OrdinalType numRows_in,
00422   OrdinalType numCols_in
00423   )
00424   : CompObject(), Object("Teuchos::SerialDenseMatrix"), numRows_(numRows_in), numCols_(numCols_in), stride_(stride_in),
00425     valuesCopied_(false), values_(values_in)
00426 {
00427   if(CV == Copy)
00428   {
00429     stride_ = numRows_;
00430     values_ = new ScalarType[stride_*numCols_];
00431     copyMat(values_in, stride_in, numRows_, numCols_, values_, stride_, 0, 0);
00432     valuesCopied_ = true;
00433   }
00434 }
00435   
00436 template<typename OrdinalType, typename ScalarType>
00437 SerialDenseMatrix<OrdinalType, ScalarType>::SerialDenseMatrix(const SerialDenseMatrix<OrdinalType, ScalarType> &Source, ETransp trans) : CompObject(), Object("Teuchos::SerialDenseMatrix"), numRows_(0), numCols_(0), stride_(0), valuesCopied_(true), values_(0)
00438 {
00439   if ( trans == Teuchos::NO_TRANS ) 
00440   {
00441     numRows_ = Source.numRows_;
00442     numCols_ = Source.numCols_;
00443     stride_ = numRows_;
00444     values_ = new ScalarType[stride_*numCols_];
00445     copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0);
00446   } 
00447   else if ( trans == Teuchos::CONJ_TRANS && ScalarTraits<ScalarType>::isComplex ) 
00448   {   
00449     numRows_ = Source.numCols_;
00450     numCols_ = Source.numRows_;
00451     stride_ = numRows_;
00452     values_ = new ScalarType[stride_*numCols_];
00453     for (OrdinalType j=0; j<numCols_; j++) {
00454       for (OrdinalType i=0; i<numRows_; i++) {
00455         values_[j*stride_ + i] = Teuchos::ScalarTraits<ScalarType>::conjugate(Source.values_[i*Source.stride_ + j]);
00456       }
00457     }
00458   } 
00459   else 
00460   {
00461     numRows_ = Source.numCols_;
00462     numCols_ = Source.numRows_;
00463     stride_ = numRows_;
00464     values_ = new ScalarType[stride_*numCols_];
00465     for (OrdinalType j=0; j<numCols_; j++) {
00466       for (OrdinalType i=0; i<numRows_; i++) {
00467         values_[j*stride_ + i] = Source.values_[i*Source.stride_ + j];
00468       }
00469     }
00470   }
00471 }
00472   
00473 template<typename OrdinalType, typename ScalarType>
00474 SerialDenseMatrix<OrdinalType, ScalarType>::SerialDenseMatrix(
00475   DataAccess CV, const SerialDenseMatrix<OrdinalType, ScalarType> &Source,
00476   OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType startRow,
00477   OrdinalType startCol
00478   )
00479   : CompObject(), Object("Teuchos::SerialDenseMatrix"), numRows_(numRows_in), numCols_(numCols_in), stride_(Source.stride_),
00480     valuesCopied_(false), values_(Source.values_)
00481 {
00482   if(CV == Copy)
00483   {
00484     stride_ = numRows_in;
00485     values_ = new ScalarType[stride_ * numCols_in];
00486     copyMat(Source.values_, Source.stride_, numRows_in, numCols_in, values_, stride_, startRow, startCol);
00487     valuesCopied_ = true;
00488   }
00489   else // CV == View
00490   {
00491     values_ = values_ + (stride_ * startCol) + startRow;
00492   }
00493 }
00494     
00495 template<typename OrdinalType, typename ScalarType>
00496 SerialDenseMatrix<OrdinalType, ScalarType>::~SerialDenseMatrix()
00497 {
00498   deleteArrays();
00499 }
00500   
00501 //----------------------------------------------------------------------------------------------------
00502 //  Shape methods 
00503 //----------------------------------------------------------------------------------------------------
00504 
00505 template<typename OrdinalType, typename ScalarType>
00506 int SerialDenseMatrix<OrdinalType, ScalarType>::shape(
00507   OrdinalType numRows_in, OrdinalType numCols_in
00508   )
00509 {
00510   deleteArrays(); // Get rid of anything that might be already allocated
00511   numRows_ = numRows_in;
00512   numCols_ = numCols_in;
00513   stride_ = numRows_;
00514   values_ = new ScalarType[stride_*numCols_];
00515   putScalar();
00516   valuesCopied_ = true;
00517   return(0);
00518 }
00519 
00520 template<typename OrdinalType, typename ScalarType>
00521 int SerialDenseMatrix<OrdinalType, ScalarType>::shapeUninitialized(
00522   OrdinalType numRows_in, OrdinalType numCols_in
00523   )
00524 {
00525   deleteArrays(); // Get rid of anything that might be already allocated
00526   numRows_ = numRows_in;
00527   numCols_ = numCols_in;
00528   stride_ = numRows_;
00529   values_ = new ScalarType[stride_*numCols_];
00530   valuesCopied_ = true;
00531   return(0);
00532 }
00533   
00534 template<typename OrdinalType, typename ScalarType>
00535 int SerialDenseMatrix<OrdinalType, ScalarType>::reshape(
00536   OrdinalType numRows_in, OrdinalType numCols_in
00537   )
00538 {
00539   // Allocate space for new matrix
00540   ScalarType* values_tmp = new ScalarType[numRows_in * numCols_in];
00541   ScalarType zero = ScalarTraits<ScalarType>::zero();
00542   for(OrdinalType k = 0; k < numRows_in * numCols_in; k++)
00543   {
00544     values_tmp[k] = zero;
00545   }
00546   OrdinalType numRows_tmp = TEUCHOS_MIN(numRows_, numRows_in);
00547   OrdinalType numCols_tmp = TEUCHOS_MIN(numCols_, numCols_in);
00548   if(values_ != 0)
00549   {
00550     copyMat(values_, stride_, numRows_tmp, numCols_tmp, values_tmp,
00551       numRows_in, 0, 0); // Copy principal submatrix of A to new A
00552   }
00553   deleteArrays(); // Get rid of anything that might be already allocated
00554   numRows_ = numRows_in;
00555   numCols_ = numCols_in;
00556   stride_ = numRows_;
00557   values_ = values_tmp; // Set pointer to new A
00558   valuesCopied_ = true;
00559   return(0);
00560 }
00561    
00562 //----------------------------------------------------------------------------------------------------
00563 //   Set methods 
00564 //----------------------------------------------------------------------------------------------------
00565 
00566 template<typename OrdinalType, typename ScalarType>
00567 int SerialDenseMatrix<OrdinalType, ScalarType>::putScalar( const ScalarType value_in )
00568 {
00569   // Set each value of the dense matrix to "value".
00570   for(OrdinalType j = 0; j < numCols_; j++) 
00571   {
00572     for(OrdinalType i = 0; i < numRows_; i++) 
00573     {
00574       values_[i + j*stride_] = value_in;
00575     }
00576   }
00577   return 0;
00578 }    
00579     
00580 template<typename OrdinalType, typename ScalarType>
00581 int SerialDenseMatrix<OrdinalType, ScalarType>::random()
00582 {
00583   // Set each value of the dense matrix to a random value.
00584   for(OrdinalType j = 0; j < numCols_; j++) 
00585   {
00586     for(OrdinalType i = 0; i < numRows_; i++) 
00587     {
00588       values_[i + j*stride_] = ScalarTraits<ScalarType>::random();
00589     }
00590   }
00591   return 0;
00592 }
00593 
00594 template<typename OrdinalType, typename ScalarType>
00595 SerialDenseMatrix<OrdinalType,ScalarType>&
00596 SerialDenseMatrix<OrdinalType, ScalarType>::operator=(
00597   const SerialDenseMatrix<OrdinalType,ScalarType>& Source
00598   )
00599 {
00600   if(this == &Source)
00601     return(*this); // Special case of source same as target
00602   if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
00603     return(*this); // Special case of both are views to same data.
00604 
00605   // If the source is a view then we will return a view, else we will return a copy.
00606   if (!Source.valuesCopied_) {
00607     if(valuesCopied_)       { 
00608       // Clean up stored data if this was previously a copy.
00609       deleteArrays();
00610     }
00611     numRows_ = Source.numRows_; 
00612     numCols_ = Source.numCols_;
00613     stride_ = Source.stride_;
00614     values_ = Source.values_;
00615   }
00616   else {
00617     // If we were a view, we will now be a copy.
00618     if(!valuesCopied_) {
00619       numRows_ = Source.numRows_;
00620       numCols_ = Source.numCols_;
00621       stride_ = Source.numRows_;
00622       const OrdinalType newsize = stride_ * numCols_;
00623       if(newsize > 0) {
00624         values_ = new ScalarType[newsize];
00625         valuesCopied_ = true;
00626       }
00627       else {
00628         values_ = 0;
00629       }
00630     }
00631     // If we were a copy, we will stay a copy.
00632     else {
00633       if((Source.numRows_ <= stride_) && (Source.numCols_ == numCols_)) { // we don't need to reallocate
00634         numRows_ = Source.numRows_;
00635         numCols_ = Source.numCols_;
00636       }
00637       else { // we need to allocate more space (or less space)
00638         deleteArrays();
00639         numRows_ = Source.numRows_;
00640         numCols_ = Source.numCols_;
00641         stride_ = Source.numRows_;
00642         const OrdinalType newsize = stride_ * numCols_;
00643         if(newsize > 0) {
00644           values_ = new ScalarType[newsize];
00645           valuesCopied_ = true;
00646         }
00647       }
00648     }
00649     copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0);
00650   } 
00651   return(*this);
00652 }
00653 
00654 template<typename OrdinalType, typename ScalarType>
00655 SerialDenseMatrix<OrdinalType, ScalarType>& SerialDenseMatrix<OrdinalType, ScalarType>::operator+= (const SerialDenseMatrix<OrdinalType,ScalarType>& Source )
00656 {
00657   // Check for compatible dimensions
00658   if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_))
00659   {
00660     TEUCHOS_CHK_REF(*this); // Return *this without altering it.
00661   }
00662   copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0, ScalarTraits<ScalarType>::one());
00663   return(*this);
00664 }
00665 
00666 template<typename OrdinalType, typename ScalarType>
00667 SerialDenseMatrix<OrdinalType, ScalarType>& SerialDenseMatrix<OrdinalType, ScalarType>::operator-= (const SerialDenseMatrix<OrdinalType,ScalarType>& Source )
00668 {
00669   // Check for compatible dimensions
00670   if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_))
00671   {
00672     TEUCHOS_CHK_REF(*this); // Return *this without altering it.
00673   }
00674   copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0, -ScalarTraits<ScalarType>::one());
00675   return(*this);
00676 }
00677 
00678 template<typename OrdinalType, typename ScalarType>
00679 SerialDenseMatrix<OrdinalType,ScalarType>& SerialDenseMatrix<OrdinalType, ScalarType>::assign (const SerialDenseMatrix<OrdinalType,ScalarType>& Source) {
00680   if(this == &Source)
00681     return(*this); // Special case of source same as target
00682   if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
00683     return(*this); // Special case of both are views to same data.
00684 
00685   // Check for compatible dimensions
00686   if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_))
00687   {
00688     TEUCHOS_CHK_REF(*this); // Return *this without altering it.
00689   }
00690   copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0);
00691   return(*this);
00692 }
00693 
00694 //----------------------------------------------------------------------------------------------------
00695 //   Accessor methods 
00696 //----------------------------------------------------------------------------------------------------
00697 
00698 template<typename OrdinalType, typename ScalarType>
00699 inline ScalarType& SerialDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex)
00700 {
00701 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
00702   checkIndex( rowIndex, colIndex );
00703 #endif
00704   return(values_[colIndex * stride_ + rowIndex]);
00705 }
00706   
00707 template<typename OrdinalType, typename ScalarType>
00708 inline const ScalarType& SerialDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex) const
00709 {
00710 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
00711   checkIndex( rowIndex, colIndex );
00712 #endif
00713   return(values_[colIndex * stride_ + rowIndex]);
00714 }
00715   
00716 template<typename OrdinalType, typename ScalarType>
00717 inline const ScalarType* SerialDenseMatrix<OrdinalType, ScalarType>::operator [] (OrdinalType colIndex) const
00718 {
00719 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
00720   checkIndex( 0, colIndex );
00721 #endif
00722   return(values_ + colIndex * stride_);
00723 }
00724   
00725 template<typename OrdinalType, typename ScalarType>
00726 inline ScalarType* SerialDenseMatrix<OrdinalType, ScalarType>::operator [] (OrdinalType colIndex)
00727 {
00728 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
00729   checkIndex( 0, colIndex );
00730 #endif
00731   return(values_ + colIndex * stride_);
00732 }
00733 
00734 //----------------------------------------------------------------------------------------------------
00735 //   Norm methods 
00736 //----------------------------------------------------------------------------------------------------
00737   
00738 template<typename OrdinalType, typename ScalarType>
00739 typename ScalarTraits<ScalarType>::magnitudeType SerialDenseMatrix<OrdinalType, ScalarType>::normOne() const
00740 {
00741   OrdinalType i, j;
00742   typename ScalarTraits<ScalarType>::magnitudeType anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00743   typename ScalarTraits<ScalarType>::magnitudeType absSum = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00744   ScalarType* ptr;
00745   for(j = 0; j < numCols_; j++)
00746   {
00747     ScalarType sum = 0;
00748     ptr = values_ + j * stride_;
00749     for(i = 0; i < numRows_; i++)
00750     {
00751       sum += ScalarTraits<ScalarType>::magnitude(*ptr++);
00752     }
00753     absSum = ScalarTraits<ScalarType>::magnitude(sum);
00754     if(absSum > anorm)
00755     {
00756       anorm = absSum;
00757     }
00758   }
00759   updateFlops(numRows_ * numCols_);
00760   return(anorm);
00761 }
00762   
00763 template<typename OrdinalType, typename ScalarType>
00764 typename ScalarTraits<ScalarType>::magnitudeType SerialDenseMatrix<OrdinalType, ScalarType>::normInf() const
00765 {
00766   OrdinalType i, j;
00767   typename ScalarTraits<ScalarType>::magnitudeType sum, anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00768     
00769   for (i = 0; i < numRows_; i++) {
00770     sum = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00771     for (j=0; j< numCols_; j++) {
00772       sum += ScalarTraits<ScalarType>::magnitude(*(values_+i+j*stride_));
00773     }
00774     anorm = TEUCHOS_MAX( anorm, sum );
00775   }
00776   updateFlops(numRows_ * numCols_);
00777   return(anorm);
00778 }
00779   
00780 template<typename OrdinalType, typename ScalarType>
00781 typename ScalarTraits<ScalarType>::magnitudeType SerialDenseMatrix<OrdinalType, ScalarType>::normFrobenius() const
00782 {
00783   OrdinalType i, j;
00784   typename ScalarTraits<ScalarType>::magnitudeType anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00785   for (j = 0; j < numCols_; j++) {
00786     for (i = 0; i < numRows_; i++) {
00787       anorm += ScalarTraits<ScalarType>::magnitude(values_[i+j*stride_]*values_[i+j*stride_]);
00788     }
00789   }
00790   anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::squareroot(anorm));
00791   updateFlops(numRows_ * numCols_);
00792   return(anorm);
00793 }
00794   
00795 //----------------------------------------------------------------------------------------------------
00796 //   Comparison methods 
00797 //----------------------------------------------------------------------------------------------------
00798   
00799 template<typename OrdinalType, typename ScalarType>
00800 bool SerialDenseMatrix<OrdinalType, ScalarType>::operator== (const SerialDenseMatrix<OrdinalType, ScalarType> &Operand) const
00801 {
00802   bool result = 1;
00803   if((numRows_ != Operand.numRows_) || (numCols_ != Operand.numCols_))
00804   {
00805     result = 0;
00806   }
00807   else
00808   {
00809     OrdinalType i, j;
00810     for(i = 0; i < numRows_; i++)
00811     {
00812       for(j = 0; j < numCols_; j++)
00813       {
00814         if((*this)(i, j) != Operand(i, j))
00815         {
00816           return 0;
00817         }
00818       }
00819     }
00820   }
00821   return result;
00822 }
00823   
00824 template<typename OrdinalType, typename ScalarType>
00825 bool SerialDenseMatrix<OrdinalType, ScalarType>::operator!= (const SerialDenseMatrix<OrdinalType, ScalarType> &Operand) const
00826 {
00827   return !((*this) == Operand);
00828 }
00829   
00830 //----------------------------------------------------------------------------------------------------
00831 //   Multiplication method 
00832 //----------------------------------------------------------------------------------------------------  
00833 
00834 template<typename OrdinalType, typename ScalarType>
00835 SerialDenseMatrix<OrdinalType, ScalarType>& SerialDenseMatrix<OrdinalType, ScalarType>::operator*= (const ScalarType alpha )
00836 {
00837   this->scale( alpha );
00838   return(*this);
00839 }
00840 
00841 template<typename OrdinalType, typename ScalarType>
00842 int SerialDenseMatrix<OrdinalType, ScalarType>::scale( const ScalarType alpha )
00843 {
00844   OrdinalType i, j;
00845   ScalarType* ptr;
00846     
00847   for (j=0; j<numCols_; j++) {
00848     ptr = values_ + j*stride_;
00849     for (i=0; i<numRows_; i++) { *ptr = alpha * (*ptr); ptr++; }
00850   }
00851   updateFlops( numRows_*numCols_ );
00852   return(0);
00853 }
00854 
00855 template<typename OrdinalType, typename ScalarType>
00856 int SerialDenseMatrix<OrdinalType, ScalarType>::scale( const SerialDenseMatrix<OrdinalType,ScalarType>& A )
00857 {
00858   OrdinalType i, j;
00859   ScalarType* ptr;
00860     
00861   // Check for compatible dimensions
00862   if ((numRows_ != A.numRows_) || (numCols_ != A.numCols_))
00863   {
00864     TEUCHOS_CHK_ERR(-1); // Return error
00865   }    
00866   for (j=0; j<numCols_; j++) {
00867     ptr = values_ + j*stride_;
00868     for (i=0; i<numRows_; i++) { *ptr = A(i,j) * (*ptr); ptr++; }
00869   }
00870   updateFlops( numRows_*numCols_ );
00871   return(0);
00872 }
00873 
00874 template<typename OrdinalType, typename ScalarType>
00875 int  SerialDenseMatrix<OrdinalType, ScalarType>::multiply(ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix<OrdinalType, ScalarType> &A, const SerialDenseMatrix<OrdinalType, ScalarType> &B, ScalarType beta)
00876 {
00877   // Check for compatible dimensions
00878   OrdinalType A_nrows = (ETranspChar[transa]!='N') ? A.numCols() : A.numRows();
00879   OrdinalType A_ncols = (ETranspChar[transa]!='N') ? A.numRows() : A.numCols();
00880   OrdinalType B_nrows = (ETranspChar[transb]!='N') ? B.numCols() : B.numRows();
00881   OrdinalType B_ncols = (ETranspChar[transb]!='N') ? B.numRows() : B.numCols();
00882   if ((numRows_ != A_nrows) || (A_ncols != B_nrows) || (numCols_ != B_ncols))
00883   {
00884     TEUCHOS_CHK_ERR(-1); // Return error
00885   }
00886   // Call GEMM function
00887   this->GEMM(transa, transb, numRows_, numCols_, A_ncols, alpha, A.values(), A.stride(), B.values(), B.stride(), beta, values_, stride_);
00888   double nflops = 2 * numRows_;
00889   nflops *= numCols_;
00890   nflops *= A_ncols;
00891   updateFlops(nflops);
00892   return(0);
00893 }
00894   
00895 template<typename OrdinalType, typename ScalarType>
00896 int SerialDenseMatrix<OrdinalType, ScalarType>::multiply (ESide sideA, ScalarType alpha, const SerialSymDenseMatrix<OrdinalType, ScalarType> &A, const SerialDenseMatrix<OrdinalType, ScalarType> &B, ScalarType beta)
00897 {
00898   // Check for compatible dimensions
00899   OrdinalType A_nrows = A.numRows(), A_ncols = A.numCols();
00900   OrdinalType B_nrows = B.numRows(), B_ncols = B.numCols();
00901   
00902   if (ESideChar[sideA]=='L') {
00903     if ((numRows_ != A_nrows) || (A_ncols != B_nrows) || (numCols_ != B_ncols)) {
00904       TEUCHOS_CHK_ERR(-1); // Return error
00905     }
00906   } else {
00907     if ((numRows_ != B_nrows) || (B_ncols != A_nrows) || (numCols_ != A_ncols)) {
00908       TEUCHOS_CHK_ERR(-1); // Return error
00909     }
00910   } 
00911   
00912   // Call SYMM function
00913   EUplo uplo = (A.upper() ? Teuchos::UPPER_TRI : Teuchos::LOWER_TRI);
00914   this->SYMM(sideA, uplo, numRows_, numCols_, alpha, A.values(), A.stride(), B.values(), B.stride(), beta, values_, stride_);
00915   double nflops = 2 * numRows_;
00916   nflops *= numCols_;
00917   nflops *= A_ncols;
00918   updateFlops(nflops);
00919   return(0);
00920 }
00921 
00922 template<typename OrdinalType, typename ScalarType>
00923 void SerialDenseMatrix<OrdinalType, ScalarType>::print(std::ostream& os) const
00924 {
00925   os << std::endl;
00926   if(valuesCopied_)
00927     os << "Values_copied : yes" << std::endl;
00928   else
00929     os << "Values_copied : no" << std::endl;
00930   os << "Rows : " << numRows_ << std::endl;
00931   os << "Columns : " << numCols_ << std::endl;
00932   os << "LDA : " << stride_ << std::endl;
00933   if(numRows_ == 0 || numCols_ == 0) {
00934     os << "(matrix is empty, no values to display)" << std::endl;
00935   } else {
00936     for(OrdinalType i = 0; i < numRows_; i++) {
00937       for(OrdinalType j = 0; j < numCols_; j++){
00938         os << (*this)(i,j) << " ";
00939       }
00940       os << std::endl;
00941     }
00942   }
00943 }
00944 
00945 //----------------------------------------------------------------------------------------------------
00946 //   Protected methods 
00947 //----------------------------------------------------------------------------------------------------  
00948 
00949 template<typename OrdinalType, typename ScalarType>
00950 inline void SerialDenseMatrix<OrdinalType, ScalarType>::checkIndex( OrdinalType rowIndex, OrdinalType colIndex ) const {
00951   TEST_FOR_EXCEPTION(rowIndex < 0 || rowIndex >= numRows_, std::out_of_range,
00952     "SerialDenseMatrix<T>::checkIndex: "
00953     "Row index " << rowIndex << " out of range [0, "<< numRows_ << ")");
00954   TEST_FOR_EXCEPTION(colIndex < 0 || colIndex >= numCols_, std::out_of_range,
00955     "SerialDenseMatrix<T>::checkIndex: "
00956     "Col index " << colIndex << " out of range [0, "<< numCols_ << ")");
00957 }
00958 
00959 template<typename OrdinalType, typename ScalarType>
00960 void SerialDenseMatrix<OrdinalType, ScalarType>::deleteArrays(void)
00961 {
00962   if (valuesCopied_)
00963   {
00964     delete [] values_;
00965     values_ = 0;
00966     valuesCopied_ = false;
00967   }
00968 }
00969   
00970 template<typename OrdinalType, typename ScalarType>
00971 void SerialDenseMatrix<OrdinalType, ScalarType>::copyMat(
00972   ScalarType* inputMatrix, OrdinalType strideInput, OrdinalType numRows_in,
00973   OrdinalType numCols_in, ScalarType* outputMatrix, OrdinalType strideOutput,
00974   OrdinalType startRow, OrdinalType startCol, ScalarType alpha
00975   )
00976 {
00977   OrdinalType i, j;
00978   ScalarType* ptr1 = 0;
00979   ScalarType* ptr2 = 0;
00980   for(j = 0; j < numCols_in; j++) {
00981     ptr1 = outputMatrix + (j * strideOutput);
00982     ptr2 = inputMatrix + (j + startCol) * strideInput + startRow;
00983     if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
00984       for(i = 0; i < numRows_in; i++)
00985       {
00986         *ptr1++ += alpha*(*ptr2++);
00987       }
00988     } else {
00989       for(i = 0; i < numRows_in; i++)
00990       {
00991         *ptr1++ = *ptr2++;
00992       }
00993     }
00994   }
00995 }
00996   
00997 } // namespace Teuchos
00998 
00999 
01000 #endif /* _TEUCHOS_SERIALDENSEMATRIX_HPP_ */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

Generated on Tue Oct 20 10:14:00 2009 for Teuchos Package Browser (Single Doxygen Collection) by  doxygen 1.6.1