Teuchos - Trilinos Tools Package Version of the Day
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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 // Kris
00043 // 06.18.03 -- Removed comments/documentation; file too hard to edit otherwise. Will replace later.
00044 //          -- Begin conversion from <ScalarType> template to <OrdinalType, ScalarType>
00045 // 06.23.03 -- Finished conversion from <ScalarType> to <OrdinalType, ScalarType>
00046 //          -- Tpetra_DenseMatrix.cpp is now obsolete
00047 //          -- Added new constructor to allow construction of a submatrix
00048 //          -- Altered copyMat to enable its use in new constructor
00049 //          -- Commented out broken print() function
00050 //          -- Fixed oneNorm() (uninitialized return variable was causing erroneous results)
00051 // 06.24.03 -- Minor formatting changes
00052 // 07.01.03 -- Added TempPrint() function to temporarily take the place of print() and operator<< while I figure out how to fix them
00053 // 07.02.03 -- Added operator== and operator!= to make testing programs easier to write/read. Implementation of == isn't the most
00054 //             efficient/robust, but it works. Will consider optimizing later.
00055 //          -- Warning! Constructor DenseMatrix(DataAccess, const DenseMatrix<OrdinalType, ScalarType> &, int, int, int, int) (the
00056 //             "submatrix grabber" constructor) does not work correctly when used with CV == View (always grabs submatrix from top
00057 //             left corner).
00058 // 07.07.03 -- Constructor bug detailed above (07.02) is now corrected (hopefully).
00059 // 07.08.03 -- Move into Teuchos package/namespace
00060 
00061 #ifndef _TEUCHOS_SERIALDENSEMATRIX_HPP_
00062 #define _TEUCHOS_SERIALDENSEMATRIX_HPP_
00063 
00067 #include "Teuchos_CompObject.hpp"
00068 #include "Teuchos_BLAS.hpp"
00069 #include "Teuchos_ScalarTraits.hpp"
00070 #include "Teuchos_DataAccess.hpp"
00071 #include "Teuchos_ConfigDefs.hpp"
00072 #include "Teuchos_TestForException.hpp"
00073 #include "Teuchos_SerialSymDenseMatrix.hpp"
00074 
00083 namespace Teuchos {  
00084 
00085 template<typename OrdinalType, typename ScalarType>
00086 class SerialDenseMatrix : public CompObject, public Object, public BLAS<OrdinalType, ScalarType>
00087 {
00088 public:
00089 
00091   typedef OrdinalType ordinalType;
00093   typedef ScalarType scalarType;
00094 
00096 
00097 
00099 
00102   SerialDenseMatrix();
00103 
00105 
00113   SerialDenseMatrix(OrdinalType numRows, OrdinalType numCols, bool zeroOut = true);
00114 
00116 
00124   SerialDenseMatrix(DataAccess CV, ScalarType* values, OrdinalType stride, OrdinalType numRows, OrdinalType numCols);
00125 
00127 
00132   SerialDenseMatrix(const SerialDenseMatrix<OrdinalType, ScalarType> &Source, ETransp trans = Teuchos::NO_TRANS);
00133 
00135 
00147   SerialDenseMatrix(DataAccess CV, const SerialDenseMatrix<OrdinalType, ScalarType> &Source, OrdinalType numRows, OrdinalType numCols, OrdinalType startRow=0, OrdinalType startCol=0);
00148 
00150   virtual ~SerialDenseMatrix();
00152 
00154 
00155 
00156 
00166   int shape(OrdinalType numRows, OrdinalType numCols);
00167 
00169   int shapeUninitialized(OrdinalType numRows, OrdinalType numCols);
00170 
00172 
00182   int reshape(OrdinalType numRows, OrdinalType numCols);
00183 
00184 
00186 
00188 
00189 
00191 
00197   SerialDenseMatrix<OrdinalType, ScalarType>& operator= (const SerialDenseMatrix<OrdinalType, ScalarType>& Source);
00198 
00200 
00205   SerialDenseMatrix<OrdinalType, ScalarType>& assign (const SerialDenseMatrix<OrdinalType, ScalarType>& Source);
00206 
00208 
00211   SerialDenseMatrix<OrdinalType, ScalarType>& operator= (const ScalarType value) { putScalar(value); return(*this); }
00212 
00214 
00218   int putScalar( const ScalarType value = Teuchos::ScalarTraits<ScalarType>::zero() );
00219 
00221   int random();
00222 
00224 
00226 
00227 
00229 
00236   ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex);
00237 
00239 
00246   const ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex) const;
00247 
00249 
00256   ScalarType* operator [] (OrdinalType colIndex);
00257 
00259 
00266   const ScalarType* operator [] (OrdinalType colIndex) const;
00267 
00269 
00270   ScalarType* values() const { return(values_); }
00271 
00273 
00275 
00276 
00278 
00281   SerialDenseMatrix<OrdinalType, ScalarType>& operator+= (const SerialDenseMatrix<OrdinalType, ScalarType>& Source);
00282 
00284 
00287   SerialDenseMatrix<OrdinalType, ScalarType>& operator-= (const SerialDenseMatrix<OrdinalType, ScalarType>& Source);
00288 
00290 
00293   SerialDenseMatrix<OrdinalType, ScalarType>& operator*= (const ScalarType alpha);
00294 
00296 
00300   int scale ( const ScalarType alpha );
00301 
00303 
00309   int scale ( const SerialDenseMatrix<OrdinalType, ScalarType>& A );
00310 
00312 
00326   int multiply (ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix<OrdinalType, ScalarType> &A, const SerialDenseMatrix<OrdinalType, ScalarType> &B, ScalarType beta);
00327 
00329 
00340   int multiply (ESide sideA, ScalarType alpha, const SerialSymDenseMatrix<OrdinalType, ScalarType> &A, const SerialDenseMatrix<OrdinalType, ScalarType> &B, ScalarType beta);
00341 
00343 
00345 
00346 
00348 
00351   bool operator== (const SerialDenseMatrix<OrdinalType, ScalarType> &Operand) const;
00352 
00354 
00357   bool operator!= (const SerialDenseMatrix<OrdinalType, ScalarType> &Operand) const;
00358 
00360 
00362 
00363 
00365   OrdinalType numRows() const { return(numRows_); }
00366 
00368   OrdinalType numCols() const { return(numCols_); }
00369 
00371   OrdinalType stride() const { return(stride_); }
00372 
00374   bool empty() const { return(numRows_ == 0 || numCols_ == 0); }
00376 
00378 
00379 
00381   typename ScalarTraits<ScalarType>::magnitudeType normOne() const;
00382 
00384   typename ScalarTraits<ScalarType>::magnitudeType normInf() const;
00385 
00387   typename ScalarTraits<ScalarType>::magnitudeType normFrobenius() const; 
00389 
00391 
00392 
00393   virtual void print(std::ostream& os) const;
00394 
00396 protected:
00397   void copyMat(ScalarType* inputMatrix, OrdinalType strideInput,
00398     OrdinalType numRows, OrdinalType numCols, ScalarType* outputMatrix,
00399     OrdinalType strideOutput, OrdinalType startRow, OrdinalType startCol,
00400     ScalarType alpha = ScalarTraits<ScalarType>::zero() );
00401   void deleteArrays();
00402   void checkIndex( OrdinalType rowIndex, OrdinalType colIndex = 0 ) const;
00403   OrdinalType numRows_;
00404   OrdinalType numCols_;
00405   OrdinalType stride_;
00406   bool valuesCopied_;
00407   ScalarType* values_;
00408 
00409 }; // class Teuchos_SerialDenseMatrix
00410 
00411 //----------------------------------------------------------------------------------------------------
00412 //  Constructors and Destructor
00413 //----------------------------------------------------------------------------------------------------
00414 
00415 template<typename OrdinalType, typename ScalarType>
00416 SerialDenseMatrix<OrdinalType, ScalarType>::SerialDenseMatrix()
00417   : CompObject(), Object("Teuchos::SerialDenseMatrix"), numRows_(0), numCols_(0), stride_(0), valuesCopied_(false), values_(0)
00418 {}
00419   
00420 template<typename OrdinalType, typename ScalarType>
00421 SerialDenseMatrix<OrdinalType, ScalarType>::SerialDenseMatrix(
00422   OrdinalType numRows_in, OrdinalType numCols_in, bool zeroOut
00423   )
00424   : CompObject(), Object("Teuchos::SerialDenseMatrix"), numRows_(numRows_in), numCols_(numCols_in), stride_(numRows_in)
00425 {
00426   values_ = new ScalarType[stride_*numCols_];
00427   valuesCopied_ = true;
00428   if (zeroOut == true)  
00429     putScalar();
00430 }
00431 
00432 template<typename OrdinalType, typename ScalarType>
00433 SerialDenseMatrix<OrdinalType, ScalarType>::SerialDenseMatrix(
00434   DataAccess CV, ScalarType* values_in, OrdinalType stride_in, OrdinalType numRows_in,
00435   OrdinalType numCols_in
00436   )
00437   : CompObject(), Object("Teuchos::SerialDenseMatrix"), numRows_(numRows_in), numCols_(numCols_in), stride_(stride_in),
00438     valuesCopied_(false), values_(values_in)
00439 {
00440   if(CV == Copy)
00441   {
00442     stride_ = numRows_;
00443     values_ = new ScalarType[stride_*numCols_];
00444     copyMat(values_in, stride_in, numRows_, numCols_, values_, stride_, 0, 0);
00445     valuesCopied_ = true;
00446   }
00447 }
00448   
00449 template<typename OrdinalType, typename ScalarType>
00450 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)
00451 {
00452   if ( trans == Teuchos::NO_TRANS ) 
00453   {
00454     numRows_ = Source.numRows_;
00455     numCols_ = Source.numCols_;
00456     stride_ = numRows_;
00457     values_ = new ScalarType[stride_*numCols_];
00458     copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0);
00459   } 
00460   else if ( trans == Teuchos::CONJ_TRANS && ScalarTraits<ScalarType>::isComplex ) 
00461   {   
00462     numRows_ = Source.numCols_;
00463     numCols_ = Source.numRows_;
00464     stride_ = numRows_;
00465     values_ = new ScalarType[stride_*numCols_];
00466     for (OrdinalType j=0; j<numCols_; j++) {
00467       for (OrdinalType i=0; i<numRows_; i++) {
00468         values_[j*stride_ + i] = Teuchos::ScalarTraits<ScalarType>::conjugate(Source.values_[i*Source.stride_ + j]);
00469       }
00470     }
00471   } 
00472   else 
00473   {
00474     numRows_ = Source.numCols_;
00475     numCols_ = Source.numRows_;
00476     stride_ = numRows_;
00477     values_ = new ScalarType[stride_*numCols_];
00478     for (OrdinalType j=0; j<numCols_; j++) {
00479       for (OrdinalType i=0; i<numRows_; i++) {
00480         values_[j*stride_ + i] = Source.values_[i*Source.stride_ + j];
00481       }
00482     }
00483   }
00484 }
00485   
00486 template<typename OrdinalType, typename ScalarType>
00487 SerialDenseMatrix<OrdinalType, ScalarType>::SerialDenseMatrix(
00488   DataAccess CV, const SerialDenseMatrix<OrdinalType, ScalarType> &Source,
00489   OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType startRow,
00490   OrdinalType startCol
00491   )
00492   : CompObject(), Object("Teuchos::SerialDenseMatrix"), numRows_(numRows_in), numCols_(numCols_in), stride_(Source.stride_),
00493     valuesCopied_(false), values_(Source.values_)
00494 {
00495   if(CV == Copy)
00496   {
00497     stride_ = numRows_in;
00498     values_ = new ScalarType[stride_ * numCols_in];
00499     copyMat(Source.values_, Source.stride_, numRows_in, numCols_in, values_, stride_, startRow, startCol);
00500     valuesCopied_ = true;
00501   }
00502   else // CV == View
00503   {
00504     values_ = values_ + (stride_ * startCol) + startRow;
00505   }
00506 }
00507     
00508 template<typename OrdinalType, typename ScalarType>
00509 SerialDenseMatrix<OrdinalType, ScalarType>::~SerialDenseMatrix()
00510 {
00511   deleteArrays();
00512 }
00513   
00514 //----------------------------------------------------------------------------------------------------
00515 //  Shape methods 
00516 //----------------------------------------------------------------------------------------------------
00517 
00518 template<typename OrdinalType, typename ScalarType>
00519 int SerialDenseMatrix<OrdinalType, ScalarType>::shape(
00520   OrdinalType numRows_in, OrdinalType numCols_in
00521   )
00522 {
00523   deleteArrays(); // Get rid of anything that might be already allocated
00524   numRows_ = numRows_in;
00525   numCols_ = numCols_in;
00526   stride_ = numRows_;
00527   values_ = new ScalarType[stride_*numCols_];
00528   putScalar();
00529   valuesCopied_ = true;
00530   return(0);
00531 }
00532 
00533 template<typename OrdinalType, typename ScalarType>
00534 int SerialDenseMatrix<OrdinalType, ScalarType>::shapeUninitialized(
00535   OrdinalType numRows_in, OrdinalType numCols_in
00536   )
00537 {
00538   deleteArrays(); // Get rid of anything that might be already allocated
00539   numRows_ = numRows_in;
00540   numCols_ = numCols_in;
00541   stride_ = numRows_;
00542   values_ = new ScalarType[stride_*numCols_];
00543   valuesCopied_ = true;
00544   return(0);
00545 }
00546   
00547 template<typename OrdinalType, typename ScalarType>
00548 int SerialDenseMatrix<OrdinalType, ScalarType>::reshape(
00549   OrdinalType numRows_in, OrdinalType numCols_in
00550   )
00551 {
00552   // Allocate space for new matrix
00553   ScalarType* values_tmp = new ScalarType[numRows_in * numCols_in];
00554   ScalarType zero = ScalarTraits<ScalarType>::zero();
00555   for(OrdinalType k = 0; k < numRows_in * numCols_in; k++)
00556   {
00557     values_tmp[k] = zero;
00558   }
00559   OrdinalType numRows_tmp = TEUCHOS_MIN(numRows_, numRows_in);
00560   OrdinalType numCols_tmp = TEUCHOS_MIN(numCols_, numCols_in);
00561   if(values_ != 0)
00562   {
00563     copyMat(values_, stride_, numRows_tmp, numCols_tmp, values_tmp,
00564       numRows_in, 0, 0); // Copy principal submatrix of A to new A
00565   }
00566   deleteArrays(); // Get rid of anything that might be already allocated
00567   numRows_ = numRows_in;
00568   numCols_ = numCols_in;
00569   stride_ = numRows_;
00570   values_ = values_tmp; // Set pointer to new A
00571   valuesCopied_ = true;
00572   return(0);
00573 }
00574    
00575 //----------------------------------------------------------------------------------------------------
00576 //   Set methods 
00577 //----------------------------------------------------------------------------------------------------
00578 
00579 template<typename OrdinalType, typename ScalarType>
00580 int SerialDenseMatrix<OrdinalType, ScalarType>::putScalar( const ScalarType value_in )
00581 {
00582   // Set each value of the dense matrix to "value".
00583   for(OrdinalType j = 0; j < numCols_; j++) 
00584   {
00585     for(OrdinalType i = 0; i < numRows_; i++) 
00586     {
00587       values_[i + j*stride_] = value_in;
00588     }
00589   }
00590   return 0;
00591 }    
00592     
00593 template<typename OrdinalType, typename ScalarType>
00594 int SerialDenseMatrix<OrdinalType, ScalarType>::random()
00595 {
00596   // Set each value of the dense matrix to a random value.
00597   for(OrdinalType j = 0; j < numCols_; j++) 
00598   {
00599     for(OrdinalType i = 0; i < numRows_; i++) 
00600     {
00601       values_[i + j*stride_] = ScalarTraits<ScalarType>::random();
00602     }
00603   }
00604   return 0;
00605 }
00606 
00607 template<typename OrdinalType, typename ScalarType>
00608 SerialDenseMatrix<OrdinalType,ScalarType>&
00609 SerialDenseMatrix<OrdinalType, ScalarType>::operator=(
00610   const SerialDenseMatrix<OrdinalType,ScalarType>& Source
00611   )
00612 {
00613   if(this == &Source)
00614     return(*this); // Special case of source same as target
00615   if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
00616     return(*this); // Special case of both are views to same data.
00617 
00618   // If the source is a view then we will return a view, else we will return a copy.
00619   if (!Source.valuesCopied_) {
00620     if(valuesCopied_)       { 
00621       // Clean up stored data if this was previously a copy.
00622       deleteArrays();
00623     }
00624     numRows_ = Source.numRows_; 
00625     numCols_ = Source.numCols_;
00626     stride_ = Source.stride_;
00627     values_ = Source.values_;
00628   }
00629   else {
00630     // If we were a view, we will now be a copy.
00631     if(!valuesCopied_) {
00632       numRows_ = Source.numRows_;
00633       numCols_ = Source.numCols_;
00634       stride_ = Source.numRows_;
00635       const OrdinalType newsize = stride_ * numCols_;
00636       if(newsize > 0) {
00637         values_ = new ScalarType[newsize];
00638         valuesCopied_ = true;
00639       }
00640       else {
00641         values_ = 0;
00642       }
00643     }
00644     // If we were a copy, we will stay a copy.
00645     else {
00646       if((Source.numRows_ <= stride_) && (Source.numCols_ == numCols_)) { // we don't need to reallocate
00647         numRows_ = Source.numRows_;
00648         numCols_ = Source.numCols_;
00649       }
00650       else { // we need to allocate more space (or less space)
00651         deleteArrays();
00652         numRows_ = Source.numRows_;
00653         numCols_ = Source.numCols_;
00654         stride_ = Source.numRows_;
00655         const OrdinalType newsize = stride_ * numCols_;
00656         if(newsize > 0) {
00657           values_ = new ScalarType[newsize];
00658           valuesCopied_ = true;
00659         }
00660       }
00661     }
00662     copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0);
00663   } 
00664   return(*this);
00665 }
00666 
00667 template<typename OrdinalType, typename ScalarType>
00668 SerialDenseMatrix<OrdinalType, ScalarType>& SerialDenseMatrix<OrdinalType, ScalarType>::operator+= (const SerialDenseMatrix<OrdinalType,ScalarType>& Source )
00669 {
00670   // Check for compatible dimensions
00671   if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_))
00672   {
00673     TEUCHOS_CHK_REF(*this); // Return *this without altering it.
00674   }
00675   copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0, ScalarTraits<ScalarType>::one());
00676   return(*this);
00677 }
00678 
00679 template<typename OrdinalType, typename ScalarType>
00680 SerialDenseMatrix<OrdinalType, ScalarType>& SerialDenseMatrix<OrdinalType, ScalarType>::operator-= (const SerialDenseMatrix<OrdinalType,ScalarType>& Source )
00681 {
00682   // Check for compatible dimensions
00683   if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_))
00684   {
00685     TEUCHOS_CHK_REF(*this); // Return *this without altering it.
00686   }
00687   copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0, -ScalarTraits<ScalarType>::one());
00688   return(*this);
00689 }
00690 
00691 template<typename OrdinalType, typename ScalarType>
00692 SerialDenseMatrix<OrdinalType,ScalarType>& SerialDenseMatrix<OrdinalType, ScalarType>::assign (const SerialDenseMatrix<OrdinalType,ScalarType>& Source) {
00693   if(this == &Source)
00694     return(*this); // Special case of source same as target
00695   if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
00696     return(*this); // Special case of both are views to same data.
00697 
00698   // Check for compatible dimensions
00699   if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_))
00700   {
00701     TEUCHOS_CHK_REF(*this); // Return *this without altering it.
00702   }
00703   copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0);
00704   return(*this);
00705 }
00706 
00707 //----------------------------------------------------------------------------------------------------
00708 //   Accessor methods 
00709 //----------------------------------------------------------------------------------------------------
00710 
00711 template<typename OrdinalType, typename ScalarType>
00712 inline ScalarType& SerialDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex)
00713 {
00714 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
00715   checkIndex( rowIndex, colIndex );
00716 #endif
00717   return(values_[colIndex * stride_ + rowIndex]);
00718 }
00719   
00720 template<typename OrdinalType, typename ScalarType>
00721 inline const ScalarType& SerialDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex) const
00722 {
00723 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
00724   checkIndex( rowIndex, colIndex );
00725 #endif
00726   return(values_[colIndex * stride_ + rowIndex]);
00727 }
00728   
00729 template<typename OrdinalType, typename ScalarType>
00730 inline const ScalarType* SerialDenseMatrix<OrdinalType, ScalarType>::operator [] (OrdinalType colIndex) const
00731 {
00732 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
00733   checkIndex( 0, colIndex );
00734 #endif
00735   return(values_ + colIndex * stride_);
00736 }
00737   
00738 template<typename OrdinalType, typename ScalarType>
00739 inline ScalarType* SerialDenseMatrix<OrdinalType, ScalarType>::operator [] (OrdinalType colIndex)
00740 {
00741 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
00742   checkIndex( 0, colIndex );
00743 #endif
00744   return(values_ + colIndex * stride_);
00745 }
00746 
00747 //----------------------------------------------------------------------------------------------------
00748 //   Norm methods 
00749 //----------------------------------------------------------------------------------------------------
00750   
00751 template<typename OrdinalType, typename ScalarType>
00752 typename ScalarTraits<ScalarType>::magnitudeType SerialDenseMatrix<OrdinalType, ScalarType>::normOne() const
00753 {
00754   OrdinalType i, j;
00755   typename ScalarTraits<ScalarType>::magnitudeType anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00756   typename ScalarTraits<ScalarType>::magnitudeType absSum = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00757   ScalarType* ptr;
00758   for(j = 0; j < numCols_; j++)
00759   {
00760     ScalarType sum = 0;
00761     ptr = values_ + j * stride_;
00762     for(i = 0; i < numRows_; i++)
00763     {
00764       sum += ScalarTraits<ScalarType>::magnitude(*ptr++);
00765     }
00766     absSum = ScalarTraits<ScalarType>::magnitude(sum);
00767     if(absSum > anorm)
00768     {
00769       anorm = absSum;
00770     }
00771   }
00772   updateFlops(numRows_ * numCols_);
00773   return(anorm);
00774 }
00775   
00776 template<typename OrdinalType, typename ScalarType>
00777 typename ScalarTraits<ScalarType>::magnitudeType SerialDenseMatrix<OrdinalType, ScalarType>::normInf() const
00778 {
00779   OrdinalType i, j;
00780   typename ScalarTraits<ScalarType>::magnitudeType sum, anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00781     
00782   for (i = 0; i < numRows_; i++) {
00783     sum = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00784     for (j=0; j< numCols_; j++) {
00785       sum += ScalarTraits<ScalarType>::magnitude(*(values_+i+j*stride_));
00786     }
00787     anorm = TEUCHOS_MAX( anorm, sum );
00788   }
00789   updateFlops(numRows_ * numCols_);
00790   return(anorm);
00791 }
00792   
00793 template<typename OrdinalType, typename ScalarType>
00794 typename ScalarTraits<ScalarType>::magnitudeType SerialDenseMatrix<OrdinalType, ScalarType>::normFrobenius() const
00795 {
00796   OrdinalType i, j;
00797   typename ScalarTraits<ScalarType>::magnitudeType anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00798   for (j = 0; j < numCols_; j++) {
00799     for (i = 0; i < numRows_; i++) {
00800       anorm += ScalarTraits<ScalarType>::magnitude(values_[i+j*stride_]*values_[i+j*stride_]);
00801     }
00802   }
00803   anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::squareroot(anorm));
00804   updateFlops(numRows_ * numCols_);
00805   return(anorm);
00806 }
00807   
00808 //----------------------------------------------------------------------------------------------------
00809 //   Comparison methods 
00810 //----------------------------------------------------------------------------------------------------
00811   
00812 template<typename OrdinalType, typename ScalarType>
00813 bool SerialDenseMatrix<OrdinalType, ScalarType>::operator== (const SerialDenseMatrix<OrdinalType, ScalarType> &Operand) const
00814 {
00815   bool result = 1;
00816   if((numRows_ != Operand.numRows_) || (numCols_ != Operand.numCols_))
00817   {
00818     result = 0;
00819   }
00820   else
00821   {
00822     OrdinalType i, j;
00823     for(i = 0; i < numRows_; i++)
00824     {
00825       for(j = 0; j < numCols_; j++)
00826       {
00827         if((*this)(i, j) != Operand(i, j))
00828         {
00829           return 0;
00830         }
00831       }
00832     }
00833   }
00834   return result;
00835 }
00836   
00837 template<typename OrdinalType, typename ScalarType>
00838 bool SerialDenseMatrix<OrdinalType, ScalarType>::operator!= (const SerialDenseMatrix<OrdinalType, ScalarType> &Operand) const
00839 {
00840   return !((*this) == Operand);
00841 }
00842   
00843 //----------------------------------------------------------------------------------------------------
00844 //   Multiplication method 
00845 //----------------------------------------------------------------------------------------------------  
00846 
00847 template<typename OrdinalType, typename ScalarType>
00848 SerialDenseMatrix<OrdinalType, ScalarType>& SerialDenseMatrix<OrdinalType, ScalarType>::operator*= (const ScalarType alpha )
00849 {
00850   this->scale( alpha );
00851   return(*this);
00852 }
00853 
00854 template<typename OrdinalType, typename ScalarType>
00855 int SerialDenseMatrix<OrdinalType, ScalarType>::scale( const ScalarType alpha )
00856 {
00857   OrdinalType i, j;
00858   ScalarType* ptr;
00859     
00860   for (j=0; j<numCols_; j++) {
00861     ptr = values_ + j*stride_;
00862     for (i=0; i<numRows_; i++) { *ptr = alpha * (*ptr); ptr++; }
00863   }
00864   updateFlops( numRows_*numCols_ );
00865   return(0);
00866 }
00867 
00868 template<typename OrdinalType, typename ScalarType>
00869 int SerialDenseMatrix<OrdinalType, ScalarType>::scale( const SerialDenseMatrix<OrdinalType,ScalarType>& A )
00870 {
00871   OrdinalType i, j;
00872   ScalarType* ptr;
00873     
00874   // Check for compatible dimensions
00875   if ((numRows_ != A.numRows_) || (numCols_ != A.numCols_))
00876   {
00877     TEUCHOS_CHK_ERR(-1); // Return error
00878   }    
00879   for (j=0; j<numCols_; j++) {
00880     ptr = values_ + j*stride_;
00881     for (i=0; i<numRows_; i++) { *ptr = A(i,j) * (*ptr); ptr++; }
00882   }
00883   updateFlops( numRows_*numCols_ );
00884   return(0);
00885 }
00886 
00887 template<typename OrdinalType, typename ScalarType>
00888 int  SerialDenseMatrix<OrdinalType, ScalarType>::multiply(ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix<OrdinalType, ScalarType> &A, const SerialDenseMatrix<OrdinalType, ScalarType> &B, ScalarType beta)
00889 {
00890   // Check for compatible dimensions
00891   OrdinalType A_nrows = (ETranspChar[transa]!='N') ? A.numCols() : A.numRows();
00892   OrdinalType A_ncols = (ETranspChar[transa]!='N') ? A.numRows() : A.numCols();
00893   OrdinalType B_nrows = (ETranspChar[transb]!='N') ? B.numCols() : B.numRows();
00894   OrdinalType B_ncols = (ETranspChar[transb]!='N') ? B.numRows() : B.numCols();
00895   if ((numRows_ != A_nrows) || (A_ncols != B_nrows) || (numCols_ != B_ncols))
00896   {
00897     TEUCHOS_CHK_ERR(-1); // Return error
00898   }
00899   // Call GEMM function
00900   this->GEMM(transa, transb, numRows_, numCols_, A_ncols, alpha, A.values(), A.stride(), B.values(), B.stride(), beta, values_, stride_);
00901   double nflops = 2 * numRows_;
00902   nflops *= numCols_;
00903   nflops *= A_ncols;
00904   updateFlops(nflops);
00905   return(0);
00906 }
00907   
00908 template<typename OrdinalType, typename ScalarType>
00909 int SerialDenseMatrix<OrdinalType, ScalarType>::multiply (ESide sideA, ScalarType alpha, const SerialSymDenseMatrix<OrdinalType, ScalarType> &A, const SerialDenseMatrix<OrdinalType, ScalarType> &B, ScalarType beta)
00910 {
00911   // Check for compatible dimensions
00912   OrdinalType A_nrows = A.numRows(), A_ncols = A.numCols();
00913   OrdinalType B_nrows = B.numRows(), B_ncols = B.numCols();
00914   
00915   if (ESideChar[sideA]=='L') {
00916     if ((numRows_ != A_nrows) || (A_ncols != B_nrows) || (numCols_ != B_ncols)) {
00917       TEUCHOS_CHK_ERR(-1); // Return error
00918     }
00919   } else {
00920     if ((numRows_ != B_nrows) || (B_ncols != A_nrows) || (numCols_ != A_ncols)) {
00921       TEUCHOS_CHK_ERR(-1); // Return error
00922     }
00923   } 
00924   
00925   // Call SYMM function
00926   EUplo uplo = (A.upper() ? Teuchos::UPPER_TRI : Teuchos::LOWER_TRI);
00927   this->SYMM(sideA, uplo, numRows_, numCols_, alpha, A.values(), A.stride(), B.values(), B.stride(), beta, values_, stride_);
00928   double nflops = 2 * numRows_;
00929   nflops *= numCols_;
00930   nflops *= A_ncols;
00931   updateFlops(nflops);
00932   return(0);
00933 }
00934 
00935 template<typename OrdinalType, typename ScalarType>
00936 void SerialDenseMatrix<OrdinalType, ScalarType>::print(std::ostream& os) const
00937 {
00938   os << std::endl;
00939   if(valuesCopied_)
00940     os << "Values_copied : yes" << std::endl;
00941   else
00942     os << "Values_copied : no" << std::endl;
00943   os << "Rows : " << numRows_ << std::endl;
00944   os << "Columns : " << numCols_ << std::endl;
00945   os << "LDA : " << stride_ << std::endl;
00946   if(numRows_ == 0 || numCols_ == 0) {
00947     os << "(matrix is empty, no values to display)" << std::endl;
00948   } else {
00949     for(OrdinalType i = 0; i < numRows_; i++) {
00950       for(OrdinalType j = 0; j < numCols_; j++){
00951         os << (*this)(i,j) << " ";
00952       }
00953       os << std::endl;
00954     }
00955   }
00956 }
00957 
00958 //----------------------------------------------------------------------------------------------------
00959 //   Protected methods 
00960 //----------------------------------------------------------------------------------------------------  
00961 
00962 template<typename OrdinalType, typename ScalarType>
00963 inline void SerialDenseMatrix<OrdinalType, ScalarType>::checkIndex( OrdinalType rowIndex, OrdinalType colIndex ) const {
00964   TEST_FOR_EXCEPTION(rowIndex < 0 || rowIndex >= numRows_, std::out_of_range,
00965     "SerialDenseMatrix<T>::checkIndex: "
00966     "Row index " << rowIndex << " out of range [0, "<< numRows_ << ")");
00967   TEST_FOR_EXCEPTION(colIndex < 0 || colIndex >= numCols_, std::out_of_range,
00968     "SerialDenseMatrix<T>::checkIndex: "
00969     "Col index " << colIndex << " out of range [0, "<< numCols_ << ")");
00970 }
00971 
00972 template<typename OrdinalType, typename ScalarType>
00973 void SerialDenseMatrix<OrdinalType, ScalarType>::deleteArrays(void)
00974 {
00975   if (valuesCopied_)
00976   {
00977     delete [] values_;
00978     values_ = 0;
00979     valuesCopied_ = false;
00980   }
00981 }
00982   
00983 template<typename OrdinalType, typename ScalarType>
00984 void SerialDenseMatrix<OrdinalType, ScalarType>::copyMat(
00985   ScalarType* inputMatrix, OrdinalType strideInput, OrdinalType numRows_in,
00986   OrdinalType numCols_in, ScalarType* outputMatrix, OrdinalType strideOutput,
00987   OrdinalType startRow, OrdinalType startCol, ScalarType alpha
00988   )
00989 {
00990   OrdinalType i, j;
00991   ScalarType* ptr1 = 0;
00992   ScalarType* ptr2 = 0;
00993   for(j = 0; j < numCols_in; j++) {
00994     ptr1 = outputMatrix + (j * strideOutput);
00995     ptr2 = inputMatrix + (j + startCol) * strideInput + startRow;
00996     if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
00997       for(i = 0; i < numRows_in; i++)
00998       {
00999         *ptr1++ += alpha*(*ptr2++);
01000       }
01001     } else {
01002       for(i = 0; i < numRows_in; i++)
01003       {
01004         *ptr1++ = *ptr2++;
01005       }
01006     }
01007   }
01008 }
01009   
01010 } // namespace Teuchos
01011 
01012 
01013 #endif /* _TEUCHOS_SERIALDENSEMATRIX_HPP_ */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines