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 
00079 
00081 
00084   SerialDenseMatrix();
00085 
00087 
00095   SerialDenseMatrix(OrdinalType numRows, OrdinalType numCols, bool zeroOut = true);
00096 
00098 
00106   SerialDenseMatrix(DataAccess CV, ScalarType* values, OrdinalType stride, OrdinalType numRows, OrdinalType numCols);
00107 
00109 
00114   SerialDenseMatrix(const SerialDenseMatrix<OrdinalType, ScalarType> &Source, ETransp trans = Teuchos::NO_TRANS);
00115 
00117 
00129   SerialDenseMatrix(DataAccess CV, const SerialDenseMatrix<OrdinalType, ScalarType> &Source, OrdinalType numRows, OrdinalType numCols, OrdinalType startRow=0, OrdinalType startCol=0);
00130 
00132   virtual ~SerialDenseMatrix();
00134 
00136 
00137 
00138 
00148   int shape(OrdinalType numRows, OrdinalType numCols);
00149 
00151   int shapeUninitialized(OrdinalType numRows, OrdinalType numCols);
00152 
00154 
00164   int reshape(OrdinalType numRows, OrdinalType numCols);
00165 
00166 
00168 
00170 
00171 
00173 
00179   SerialDenseMatrix<OrdinalType, ScalarType>& operator= (const SerialDenseMatrix<OrdinalType, ScalarType>& Source);
00180 
00182 
00187   SerialDenseMatrix<OrdinalType, ScalarType>& assign (const SerialDenseMatrix<OrdinalType, ScalarType>& Source);
00188 
00190 
00194   int putScalar( const ScalarType value = Teuchos::ScalarTraits<ScalarType>::zero() );
00195 
00197   int random();
00198 
00200 
00202 
00203 
00205 
00212   ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex);
00213 
00215 
00222   const ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex) const;
00223 
00225 
00232   ScalarType* operator [] (OrdinalType colIndex);
00233 
00235 
00242   const ScalarType* operator [] (OrdinalType colIndex) const;
00243 
00245 
00246   ScalarType* values() const { return(values_); }
00247 
00249 
00251 
00252 
00254 
00257   SerialDenseMatrix<OrdinalType, ScalarType>& operator+= (const SerialDenseMatrix<OrdinalType, ScalarType>& Source);
00258 
00260 
00263   SerialDenseMatrix<OrdinalType, ScalarType>& operator-= (const SerialDenseMatrix<OrdinalType, ScalarType>& Source);
00264 
00266 
00269   SerialDenseMatrix<OrdinalType, ScalarType>& operator*= (const ScalarType alpha);
00270 
00272 
00276   int scale ( const ScalarType alpha );
00277 
00279 
00285   int scale ( const SerialDenseMatrix<OrdinalType, ScalarType>& A );
00286 
00288 
00302   int multiply (ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix<OrdinalType, ScalarType> &A, const SerialDenseMatrix<OrdinalType, ScalarType> &B, ScalarType beta);
00303 
00305 
00316   int multiply (ESide sideA, ScalarType alpha, const SerialSymDenseMatrix<OrdinalType, ScalarType> &A, const SerialDenseMatrix<OrdinalType, ScalarType> &B, ScalarType beta);
00317 
00319 
00321 
00322 
00324 
00327   bool operator== (const SerialDenseMatrix<OrdinalType, ScalarType> &Operand);
00328 
00330 
00333   bool operator!= (const SerialDenseMatrix<OrdinalType, ScalarType> &Operand);
00334 
00336 
00338 
00339 
00341   OrdinalType numRows() const { return(numRows_); }
00342 
00344   OrdinalType numCols() const { return(numCols_); }
00345 
00347   OrdinalType stride() const { return(stride_); }
00348 
00350   bool empty() const { return(numRows_ == 0 || numCols_ == 0); }
00352 
00354 
00355 
00357   typename ScalarTraits<ScalarType>::magnitudeType normOne() const;
00358 
00360   typename ScalarTraits<ScalarType>::magnitudeType normInf() const;
00361 
00363   typename ScalarTraits<ScalarType>::magnitudeType normFrobenius() const; 
00365 
00367 
00368 
00369   virtual void print(std::ostream& os) const;
00370 
00372 protected:
00373   void copyMat(ScalarType* inputMatrix, OrdinalType strideInput,
00374     OrdinalType numRows, OrdinalType numCols, ScalarType* outputMatrix,
00375     OrdinalType strideOutput, OrdinalType startRow, OrdinalType startCol,
00376     ScalarType alpha = ScalarTraits<ScalarType>::zero() );
00377   void deleteArrays();
00378   void checkIndex( OrdinalType rowIndex, OrdinalType colIndex = 0 ) const;
00379   OrdinalType numRows_;
00380   OrdinalType numCols_;
00381   OrdinalType stride_;
00382   bool valuesCopied_;
00383   ScalarType* values_;
00384 
00385 }; // class Teuchos_SerialDenseMatrix
00386 
00387 //----------------------------------------------------------------------------------------------------
00388 //  Constructors and Destructor
00389 //----------------------------------------------------------------------------------------------------
00390 
00391 template<typename OrdinalType, typename ScalarType>
00392 SerialDenseMatrix<OrdinalType, ScalarType>::SerialDenseMatrix()
00393   : CompObject(), Object("Teuchos::SerialDenseMatrix"), numRows_(0), numCols_(0), stride_(0), valuesCopied_(false), values_(0)
00394 {}
00395   
00396 template<typename OrdinalType, typename ScalarType>
00397 SerialDenseMatrix<OrdinalType, ScalarType>::SerialDenseMatrix(
00398   OrdinalType numRows_in, OrdinalType numCols_in, bool zeroOut
00399   )
00400   : CompObject(), Object("Teuchos::SerialDenseMatrix"), numRows_(numRows_in), numCols_(numCols_in), stride_(numRows_in)
00401 {
00402   values_ = new ScalarType[stride_*numCols_];
00403   valuesCopied_ = true;
00404   if (zeroOut == true)  
00405     putScalar();
00406 }
00407 
00408 template<typename OrdinalType, typename ScalarType>
00409 SerialDenseMatrix<OrdinalType, ScalarType>::SerialDenseMatrix(
00410   DataAccess CV, ScalarType* values_in, OrdinalType stride_in, OrdinalType numRows_in,
00411   OrdinalType numCols_in
00412   )
00413   : CompObject(), Object("Teuchos::SerialDenseMatrix"), numRows_(numRows_in), numCols_(numCols_in), stride_(stride_in),
00414     valuesCopied_(false), values_(values_in)
00415 {
00416   if(CV == Copy)
00417   {
00418     stride_ = numRows_;
00419     values_ = new ScalarType[stride_*numCols_];
00420     copyMat(values_in, stride_in, numRows_, numCols_, values_, stride_, 0, 0, false);
00421     valuesCopied_ = true;
00422   }
00423 }
00424   
00425 template<typename OrdinalType, typename ScalarType>
00426 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)
00427 {
00428   if ( trans == Teuchos::NO_TRANS ) 
00429   {
00430     numRows_ = Source.numRows_;
00431     numCols_ = Source.numCols_;
00432     stride_ = numRows_;
00433     values_ = new ScalarType[stride_*numCols_];
00434     copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0, false);
00435   } 
00436   else if ( trans == Teuchos::CONJ_TRANS && ScalarTraits<ScalarType>::isComplex ) 
00437   {   
00438     numRows_ = Source.numCols_;
00439     numCols_ = Source.numRows_;
00440     stride_ = numRows_;
00441     values_ = new ScalarType[stride_*numCols_];
00442     for (OrdinalType j=0; j<numCols_; j++) {
00443       for (OrdinalType i=0; i<numRows_; i++) {
00444         values_[j*stride_ + i] = Teuchos::ScalarTraits<ScalarType>::conjugate(Source.values_[i*Source.stride_ + j]);
00445       }
00446     }
00447   } 
00448   else 
00449   {
00450     numRows_ = Source.numCols_;
00451     numCols_ = Source.numRows_;
00452     stride_ = numRows_;
00453     values_ = new ScalarType[stride_*numCols_];
00454     for (OrdinalType j=0; j<numCols_; j++) {
00455       for (OrdinalType i=0; i<numRows_; i++) {
00456         values_[j*stride_ + i] = Source.values_[i*Source.stride_ + j];
00457       }
00458     }
00459   }
00460 }
00461   
00462 template<typename OrdinalType, typename ScalarType>
00463 SerialDenseMatrix<OrdinalType, ScalarType>::SerialDenseMatrix(
00464   DataAccess CV, const SerialDenseMatrix<OrdinalType, ScalarType> &Source,
00465   OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType startRow,
00466   OrdinalType startCol
00467   )
00468   : CompObject(), Object("Teuchos::SerialDenseMatrix"), numRows_(numRows_in), numCols_(numCols_in), stride_(Source.stride_),
00469     valuesCopied_(false), values_(Source.values_)
00470 {
00471   if(CV == Copy)
00472   {
00473     stride_ = numRows_in;
00474     values_ = new ScalarType[stride_ * numCols_in];
00475     copyMat(Source.values_, Source.stride_, numRows_in, numCols_in, values_, stride_, startRow, startCol, false);
00476     valuesCopied_ = true;
00477   }
00478   else // CV == View
00479   {
00480     values_ = values_ + (stride_ * startCol) + startRow;
00481   }
00482 }
00483     
00484 template<typename OrdinalType, typename ScalarType>
00485 SerialDenseMatrix<OrdinalType, ScalarType>::~SerialDenseMatrix()
00486 {
00487   deleteArrays();
00488 }
00489   
00490 //----------------------------------------------------------------------------------------------------
00491 //  Shape methods 
00492 //----------------------------------------------------------------------------------------------------
00493 
00494 template<typename OrdinalType, typename ScalarType>
00495 int SerialDenseMatrix<OrdinalType, ScalarType>::shape(
00496   OrdinalType numRows_in, OrdinalType numCols_in
00497   )
00498 {
00499   deleteArrays(); // Get rid of anything that might be already allocated
00500   numRows_ = numRows_in;
00501   numCols_ = numCols_in;
00502   stride_ = numRows_;
00503   values_ = new ScalarType[stride_*numCols_];
00504   putScalar();
00505   valuesCopied_ = true;
00506   return(0);
00507 }
00508 
00509 template<typename OrdinalType, typename ScalarType>
00510 int SerialDenseMatrix<OrdinalType, ScalarType>::shapeUninitialized(
00511   OrdinalType numRows_in, OrdinalType numCols_in
00512   )
00513 {
00514   deleteArrays(); // Get rid of anything that might be already allocated
00515   numRows_ = numRows_in;
00516   numCols_ = numCols_in;
00517   stride_ = numRows_;
00518   values_ = new ScalarType[stride_*numCols_];
00519   valuesCopied_ = true;
00520   return(0);
00521 }
00522   
00523 template<typename OrdinalType, typename ScalarType>
00524 int SerialDenseMatrix<OrdinalType, ScalarType>::reshape(
00525   OrdinalType numRows_in, OrdinalType numCols_in
00526   )
00527 {
00528   // Allocate space for new matrix
00529   ScalarType* values_tmp = new ScalarType[numRows_in * numCols_in];
00530   ScalarType zero = ScalarTraits<ScalarType>::zero();
00531   for(OrdinalType k = 0; k < numRows_in * numCols_in; k++)
00532   {
00533     values_tmp[k] = zero;
00534   }
00535   OrdinalType numRows_tmp = TEUCHOS_MIN(numRows_, numRows_in);
00536   OrdinalType numCols_tmp = TEUCHOS_MIN(numCols_, numCols_in);
00537   if(values_ != 0)
00538   {
00539     copyMat(values_, stride_, numRows_tmp, numCols_tmp, values_tmp,
00540       numRows_in, 0, 0, false); // Copy principal submatrix of A to new A
00541   }
00542   deleteArrays(); // Get rid of anything that might be already allocated
00543   numRows_ = numRows_in;
00544   numCols_ = numCols_in;
00545   stride_ = numRows_;
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 int SerialDenseMatrix<OrdinalType, ScalarType>::putScalar( const ScalarType value_in )
00557 {
00558   // Set each value of the dense matrix to "value".
00559   for(OrdinalType j = 0; j < numCols_; j++) 
00560   {
00561     for(OrdinalType i = 0; i < numRows_; i++) 
00562     {
00563       values_[i + j*stride_] = value_in;
00564     }
00565   }
00566   return 0;
00567 }    
00568     
00569 template<typename OrdinalType, typename ScalarType>
00570 int SerialDenseMatrix<OrdinalType, ScalarType>::random()
00571 {
00572   // Set each value of the dense matrix to a random value.
00573   for(OrdinalType j = 0; j < numCols_; j++) 
00574   {
00575     for(OrdinalType i = 0; i < numRows_; i++) 
00576     {
00577       values_[i + j*stride_] = ScalarTraits<ScalarType>::random();
00578     }
00579   }
00580   return 0;
00581 }
00582 
00583 template<typename OrdinalType, typename ScalarType>
00584 SerialDenseMatrix<OrdinalType,ScalarType>&
00585 SerialDenseMatrix<OrdinalType, ScalarType>::operator=(
00586   const SerialDenseMatrix<OrdinalType,ScalarType>& Source
00587   )
00588 {
00589   if(this == &Source)
00590     return(*this); // Special case of source same as target
00591   if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
00592     return(*this); // Special case of both are views to same data.
00593 
00594   // If the source is a view then we will return a view, else we will return a copy.
00595   if (!Source.valuesCopied_) {
00596     if(valuesCopied_)       { 
00597       // Clean up stored data if this was previously a copy.
00598       deleteArrays();
00599     }
00600     numRows_ = Source.numRows_; 
00601     numCols_ = Source.numCols_;
00602     stride_ = Source.stride_;
00603     values_ = Source.values_;
00604   }
00605   else {
00606     // If we were a view, we will now be a copy.
00607     if(!valuesCopied_) {
00608       numRows_ = Source.numRows_;
00609       numCols_ = Source.numCols_;
00610       stride_ = Source.numRows_;
00611       const OrdinalType newsize = stride_ * numCols_;
00612       if(newsize > 0) {
00613         values_ = new ScalarType[newsize];
00614         valuesCopied_ = true;
00615       }
00616       else {
00617         values_ = 0;
00618       }
00619     }
00620     // If we were a copy, we will stay a copy.
00621     else {
00622       if((Source.numRows_ <= stride_) && (Source.numCols_ == numCols_)) { // we don't need to reallocate
00623         numRows_ = Source.numRows_;
00624         numCols_ = Source.numCols_;
00625       }
00626       else { // we need to allocate more space (or less space)
00627         deleteArrays();
00628         numRows_ = Source.numRows_;
00629         numCols_ = Source.numCols_;
00630         stride_ = Source.numRows_;
00631         const OrdinalType newsize = stride_ * numCols_;
00632         if(newsize > 0) {
00633           values_ = new ScalarType[newsize];
00634           valuesCopied_ = true;
00635         }
00636       }
00637     }
00638     copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0, false);
00639   } 
00640   return(*this);
00641 }
00642 
00643 template<typename OrdinalType, typename ScalarType>
00644 SerialDenseMatrix<OrdinalType, ScalarType>& SerialDenseMatrix<OrdinalType, ScalarType>::operator+= (const SerialDenseMatrix<OrdinalType,ScalarType>& Source )
00645 {
00646   // Check for compatible dimensions
00647   if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_))
00648   {
00649     TEUCHOS_CHK_REF(*this); // Return *this without altering it.
00650   }
00651   copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0, 1.0);
00652   return(*this);
00653 }
00654 
00655 template<typename OrdinalType, typename ScalarType>
00656 SerialDenseMatrix<OrdinalType, ScalarType>& SerialDenseMatrix<OrdinalType, ScalarType>::operator-= (const SerialDenseMatrix<OrdinalType,ScalarType>& Source )
00657 {
00658   // Check for compatible dimensions
00659   if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_))
00660   {
00661     TEUCHOS_CHK_REF(*this); // Return *this without altering it.
00662   }
00663   copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0, -1.0);
00664   return(*this);
00665 }
00666 
00667 template<typename OrdinalType, typename ScalarType>
00668 SerialDenseMatrix<OrdinalType,ScalarType>& SerialDenseMatrix<OrdinalType, ScalarType>::assign (const SerialDenseMatrix<OrdinalType,ScalarType>& Source) {
00669   if(this == &Source)
00670     return(*this); // Special case of source same as target
00671   if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
00672     return(*this); // Special case of both are views to same data.
00673 
00674   // Check for compatible dimensions
00675   if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_))
00676   {
00677     TEUCHOS_CHK_REF(*this); // Return *this without altering it.
00678   }
00679   copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0 );
00680   return(*this);
00681 }
00682 
00683 //----------------------------------------------------------------------------------------------------
00684 //   Accessor methods 
00685 //----------------------------------------------------------------------------------------------------
00686 
00687 template<typename OrdinalType, typename ScalarType>
00688 inline ScalarType& SerialDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex)
00689 {
00690 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
00691   checkIndex( rowIndex, colIndex );
00692 #endif
00693   return(values_[colIndex * stride_ + rowIndex]);
00694 }
00695   
00696 template<typename OrdinalType, typename ScalarType>
00697 inline const ScalarType& SerialDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex) const
00698 {
00699 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
00700   checkIndex( rowIndex, colIndex );
00701 #endif
00702   return(values_[colIndex * stride_ + rowIndex]);
00703 }
00704   
00705 template<typename OrdinalType, typename ScalarType>
00706 inline const ScalarType* SerialDenseMatrix<OrdinalType, ScalarType>::operator [] (OrdinalType colIndex) const
00707 {
00708 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
00709   checkIndex( 0, colIndex );
00710 #endif
00711   return(values_ + colIndex * stride_);
00712 }
00713   
00714 template<typename OrdinalType, typename ScalarType>
00715 inline ScalarType* SerialDenseMatrix<OrdinalType, ScalarType>::operator [] (OrdinalType colIndex)
00716 {
00717 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
00718   checkIndex( 0, colIndex );
00719 #endif
00720   return(values_ + colIndex * stride_);
00721 }
00722 
00723 //----------------------------------------------------------------------------------------------------
00724 //   Norm methods 
00725 //----------------------------------------------------------------------------------------------------
00726   
00727 template<typename OrdinalType, typename ScalarType>
00728 typename ScalarTraits<ScalarType>::magnitudeType SerialDenseMatrix<OrdinalType, ScalarType>::normOne() const
00729 {
00730   OrdinalType i, j;
00731   typename ScalarTraits<ScalarType>::magnitudeType anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00732   typename ScalarTraits<ScalarType>::magnitudeType absSum = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00733   ScalarType* ptr;
00734   for(j = 0; j < numCols_; j++)
00735   {
00736     ScalarType sum = 0;
00737     ptr = values_ + j * stride_;
00738     for(i = 0; i < numRows_; i++)
00739     {
00740       sum += ScalarTraits<ScalarType>::magnitude(*ptr++);
00741     }
00742     absSum = ScalarTraits<ScalarType>::magnitude(sum);
00743     if(absSum > anorm)
00744     {
00745       anorm = absSum;
00746     }
00747   }
00748   updateFlops(numRows_ * numCols_);
00749   return(anorm);
00750 }
00751   
00752 template<typename OrdinalType, typename ScalarType>
00753 typename ScalarTraits<ScalarType>::magnitudeType SerialDenseMatrix<OrdinalType, ScalarType>::normInf() const
00754 {
00755   OrdinalType i, j;
00756   typename ScalarTraits<ScalarType>::magnitudeType sum, anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00757     
00758   for (i = 0; i < numRows_; i++) {
00759     sum = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00760     for (j=0; j< numCols_; j++) {
00761       sum += ScalarTraits<ScalarType>::magnitude(*(values_+i+j*stride_));
00762     }
00763     anorm = TEUCHOS_MAX( anorm, sum );
00764   }
00765   updateFlops(numRows_ * numCols_);
00766   return(anorm);
00767 }
00768   
00769 template<typename OrdinalType, typename ScalarType>
00770 typename ScalarTraits<ScalarType>::magnitudeType SerialDenseMatrix<OrdinalType, ScalarType>::normFrobenius() const
00771 {
00772   OrdinalType i, j;
00773   typename ScalarTraits<ScalarType>::magnitudeType anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00774   for (j = 0; j < numCols_; j++) {
00775     for (i = 0; i < numRows_; i++) {
00776       anorm += ScalarTraits<ScalarType>::magnitude(values_[i+j*stride_]*values_[i+j*stride_]);
00777     }
00778   }
00779   anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::squareroot(anorm));
00780   updateFlops(numRows_ * numCols_);
00781   return(anorm);
00782 }
00783   
00784 //----------------------------------------------------------------------------------------------------
00785 //   Comparison methods 
00786 //----------------------------------------------------------------------------------------------------
00787   
00788 template<typename OrdinalType, typename ScalarType>
00789 bool SerialDenseMatrix<OrdinalType, ScalarType>::operator== (const SerialDenseMatrix<OrdinalType, ScalarType> &Operand)
00790 {
00791   bool result = 1;
00792   if((numRows_ != Operand.numRows_) || (numCols_ != Operand.numCols_))
00793   {
00794     result = 0;
00795   }
00796   else
00797   {
00798     OrdinalType i, j;
00799     for(i = 0; i < numRows_; i++)
00800     {
00801       for(j = 0; j < numCols_; j++)
00802       {
00803         if((*this)(i, j) != Operand(i, j))
00804         {
00805           return 0;
00806         }
00807       }
00808     }
00809   }
00810   return result;
00811 }
00812   
00813 template<typename OrdinalType, typename ScalarType>
00814 bool SerialDenseMatrix<OrdinalType, ScalarType>::operator!= (const SerialDenseMatrix<OrdinalType, ScalarType> &Operand)
00815 {
00816   return !((*this) == Operand);
00817 }
00818   
00819 //----------------------------------------------------------------------------------------------------
00820 //   Multiplication method 
00821 //----------------------------------------------------------------------------------------------------  
00822 
00823 template<typename OrdinalType, typename ScalarType>
00824 SerialDenseMatrix<OrdinalType, ScalarType>& SerialDenseMatrix<OrdinalType, ScalarType>::operator*= (const ScalarType alpha )
00825 {
00826   this->scale( alpha );
00827   return(*this);
00828 }
00829 
00830 template<typename OrdinalType, typename ScalarType>
00831 int SerialDenseMatrix<OrdinalType, ScalarType>::scale( const ScalarType alpha )
00832 {
00833   OrdinalType i, j;
00834   ScalarType* ptr;
00835     
00836   for (j=0; j<numCols_; j++) {
00837     ptr = values_ + j*stride_;
00838     for (i=0; i<numRows_; i++) { *ptr = alpha * (*ptr); ptr++; }
00839   }
00840   updateFlops( numRows_*numCols_ );
00841   return(0);
00842 }
00843 
00844 template<typename OrdinalType, typename ScalarType>
00845 int SerialDenseMatrix<OrdinalType, ScalarType>::scale( const SerialDenseMatrix<OrdinalType,ScalarType>& A )
00846 {
00847   OrdinalType i, j;
00848   ScalarType* ptr;
00849     
00850   // Check for compatible dimensions
00851   if ((numRows_ != A.numRows_) || (numCols_ != A.numCols_))
00852   {
00853     TEUCHOS_CHK_ERR(-1); // Return error
00854   }    
00855   for (j=0; j<numCols_; j++) {
00856     ptr = values_ + j*stride_;
00857     for (i=0; i<numRows_; i++) { *ptr = A(i,j) * (*ptr); ptr++; }
00858   }
00859   updateFlops( numRows_*numCols_ );
00860   return(0);
00861 }
00862 
00863 template<typename OrdinalType, typename ScalarType>
00864 int  SerialDenseMatrix<OrdinalType, ScalarType>::multiply(ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix<OrdinalType, ScalarType> &A, const SerialDenseMatrix<OrdinalType, ScalarType> &B, ScalarType beta)
00865 {
00866   // Check for compatible dimensions
00867   OrdinalType A_nrows = (ETranspChar[transa]!='N') ? A.numCols() : A.numRows();
00868   OrdinalType A_ncols = (ETranspChar[transa]!='N') ? A.numRows() : A.numCols();
00869   OrdinalType B_nrows = (ETranspChar[transb]!='N') ? B.numCols() : B.numRows();
00870   OrdinalType B_ncols = (ETranspChar[transb]!='N') ? B.numRows() : B.numCols();
00871   if ((numRows_ != A_nrows) || (A_ncols != B_nrows) || (numCols_ != B_ncols))
00872   {
00873     TEUCHOS_CHK_ERR(-1); // Return error
00874   }
00875   // Call GEMM function
00876   this->GEMM(transa, transb, numRows_, numCols_, A_ncols, alpha, A.values(), A.stride(), B.values(), B.stride(), beta, values_, stride_);
00877   double nflops = 2 * numRows_;
00878   nflops *= numCols_;
00879   nflops *= A_ncols;
00880   updateFlops(nflops);
00881   return(0);
00882 }
00883   
00884 template<typename OrdinalType, typename ScalarType>
00885 int SerialDenseMatrix<OrdinalType, ScalarType>::multiply (ESide sideA, ScalarType alpha, const SerialSymDenseMatrix<OrdinalType, ScalarType> &A, const SerialDenseMatrix<OrdinalType, ScalarType> &B, ScalarType beta)
00886 {
00887   // Check for compatible dimensions
00888   OrdinalType A_nrows = A.numRows(), A_ncols = A.numCols();
00889   OrdinalType B_nrows = B.numRows(), B_ncols = B.numCols();
00890   
00891   if (ESideChar[sideA]=='L') {
00892     if ((numRows_ != A_nrows) || (A_ncols != B_nrows) || (numCols_ != B_ncols)) {
00893       TEUCHOS_CHK_ERR(-1); // Return error
00894     }
00895   } else {
00896     if ((numRows_ != B_nrows) || (B_ncols != A_nrows) || (numCols_ != A_ncols)) {
00897       TEUCHOS_CHK_ERR(-1); // Return error
00898     }
00899   } 
00900   
00901   // Call SYMM function
00902   EUplo uplo = (A.upper() ? Teuchos::UPPER_TRI : Teuchos::LOWER_TRI);
00903   this->SYMM(sideA, uplo, numRows_, numCols_, alpha, A.values(), A.stride(), B.values(), B.stride(), beta, values_, stride_);
00904   double nflops = 2 * numRows_;
00905   nflops *= numCols_;
00906   nflops *= A_ncols;
00907   updateFlops(nflops);
00908   return(0);
00909 }
00910 
00911 template<typename OrdinalType, typename ScalarType>
00912 void SerialDenseMatrix<OrdinalType, ScalarType>::print(std::ostream& os) const
00913 {
00914   os << std::endl;
00915   if(valuesCopied_)
00916     os << "Values_copied : yes" << std::endl;
00917   else
00918     os << "Values_copied : no" << std::endl;
00919   os << "Rows : " << numRows_ << std::endl;
00920   os << "Columns : " << numCols_ << std::endl;
00921   os << "LDA : " << stride_ << std::endl;
00922   if(numRows_ == 0 || numCols_ == 0) {
00923     os << "(matrix is empty, no values to display)" << std::endl;
00924   } else {
00925     for(OrdinalType i = 0; i < numRows_; i++) {
00926       for(OrdinalType j = 0; j < numCols_; j++){
00927         os << (*this)(i,j) << " ";
00928       }
00929       os << std::endl;
00930     }
00931   }
00932 }
00933 
00934 //----------------------------------------------------------------------------------------------------
00935 //   Protected methods 
00936 //----------------------------------------------------------------------------------------------------  
00937 
00938 template<typename OrdinalType, typename ScalarType>
00939 inline void SerialDenseMatrix<OrdinalType, ScalarType>::checkIndex( OrdinalType rowIndex, OrdinalType colIndex ) const {
00940   TEST_FOR_EXCEPTION(rowIndex < 0 || rowIndex >= numRows_, std::out_of_range,
00941     "SerialDenseMatrix<T>::checkIndex: "
00942     "Row index " << rowIndex << " out of range [0, "<< numRows_ << ")");
00943   TEST_FOR_EXCEPTION(colIndex < 0 || colIndex >= numCols_, std::out_of_range,
00944     "SerialDenseMatrix<T>::checkIndex: "
00945     "Col index " << colIndex << " out of range [0, "<< numCols_ << ")");
00946 }
00947 
00948 template<typename OrdinalType, typename ScalarType>
00949 void SerialDenseMatrix<OrdinalType, ScalarType>::deleteArrays(void)
00950 {
00951   if (valuesCopied_)
00952   {
00953     delete [] values_;
00954     values_ = 0;
00955     valuesCopied_ = false;
00956   }
00957 }
00958   
00959 template<typename OrdinalType, typename ScalarType>
00960 void SerialDenseMatrix<OrdinalType, ScalarType>::copyMat(
00961   ScalarType* inputMatrix, OrdinalType strideInput, OrdinalType numRows_in,
00962   OrdinalType numCols_in, ScalarType* outputMatrix, OrdinalType strideOutput,
00963   OrdinalType startRow, OrdinalType startCol, ScalarType alpha
00964   )
00965 {
00966   OrdinalType i, j;
00967   ScalarType* ptr1 = 0;
00968   ScalarType* ptr2 = 0;
00969   for(j = 0; j < numCols_in; j++) {
00970     ptr1 = outputMatrix + (j * strideOutput);
00971     ptr2 = inputMatrix + (j + startCol) * strideInput + startRow;
00972     if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
00973       for(i = 0; i < numRows_in; i++)
00974       {
00975         *ptr1++ += alpha*(*ptr2++);
00976       }
00977     } else {
00978       for(i = 0; i < numRows_in; i++)
00979       {
00980         *ptr1++ = *ptr2++;
00981       }
00982     }
00983   }
00984 }
00985   
00986 } // namespace Teuchos
00987 
00988 
00989 #endif /* _TEUCHOS_SERIALDENSEMATRIX_HPP_ */

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