Teuchos - Trilinos Tools Package: Teuchos_SerialDenseMatrix.hpp Source 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 
00069 namespace Teuchos {  
00070 
00071   template<typename OrdinalType, typename ScalarType>
00072   class SerialDenseMatrix : public CompObject, public Object, public BLAS<OrdinalType, ScalarType>
00073   {
00074   public:
00075 
00077 
00078 
00080 
00083     SerialDenseMatrix();
00084 
00086 
00093     SerialDenseMatrix(int numRows, int numCols);
00094 
00096 
00104     SerialDenseMatrix(DataAccess CV, ScalarType* values, int stride, int numRows, int numCols);
00105 
00107 
00112     SerialDenseMatrix(const SerialDenseMatrix<OrdinalType, ScalarType> &Source, ETransp trans = Teuchos::NO_TRANS);
00113 
00115 
00127     SerialDenseMatrix(DataAccess CV, const SerialDenseMatrix<OrdinalType, ScalarType> &Source, int numRows, int numCols, int startRow=0, int startCol=0);
00128 
00130     virtual ~SerialDenseMatrix();
00132 
00134 
00135 
00136 
00146     int shape(int numRows, int numCols);
00147 
00149     int shapeUninitialized(int numRows, int numCols);
00150 
00152 
00162     int reshape(int numRows, int numCols);
00163 
00164 
00166 
00168 
00169 
00171 
00177     SerialDenseMatrix<OrdinalType, ScalarType>& operator= (const SerialDenseMatrix<OrdinalType, ScalarType>& Source);
00178 
00180 
00185     SerialDenseMatrix<OrdinalType, ScalarType>& assign (const SerialDenseMatrix<OrdinalType, ScalarType>& Source);
00186 
00188 
00192     int putScalar( const ScalarType value = Teuchos::ScalarTraits<ScalarType>::zero() );
00193 
00195     int random();
00196 
00198 
00200 
00201 
00203 
00210     ScalarType& operator () (int rowIndex, int colIndex);
00211 
00213 
00220     const ScalarType& operator () (int rowIndex, int colIndex) const;
00221 
00223 
00230     ScalarType* operator [] (int colIndex);
00231 
00233 
00240     const ScalarType* operator [] (int colIndex) const;
00241 
00243 
00244     ScalarType* values() const { return(values_); };
00245 
00247 
00249 
00250 
00252 
00255     SerialDenseMatrix<OrdinalType, ScalarType>& operator+= (const SerialDenseMatrix<OrdinalType, ScalarType>& Source);
00256 
00258 
00261     SerialDenseMatrix<OrdinalType, ScalarType>& operator-= (const SerialDenseMatrix<OrdinalType, ScalarType>& Source);
00262 
00264 
00268     int scale ( const ScalarType alpha );
00269 
00271 
00277     int scale ( const SerialDenseMatrix<OrdinalType, ScalarType>& A );
00278 
00280 
00294     int multiply (ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix<OrdinalType, ScalarType> &A, const SerialDenseMatrix<OrdinalType, ScalarType> &B, ScalarType beta);
00296 
00298 
00299 
00301 
00304     bool operator== (const SerialDenseMatrix<OrdinalType, ScalarType> &Operand);
00305 
00307 
00310     bool operator!= (const SerialDenseMatrix<OrdinalType, ScalarType> &Operand);
00311 
00313 
00315 
00316 
00318     int numRows() const { return(numRows_); };
00319 
00321     int numCols() const { return(numCols_); };
00322 
00324     int stride() const { return(stride_); };
00326 
00328 
00329 
00331     typename ScalarTraits<ScalarType>::magnitudeType normOne() const;
00332 
00334     typename ScalarTraits<ScalarType>::magnitudeType normInf() const;
00335 
00337     typename ScalarTraits<ScalarType>::magnitudeType normFrobenius() const; 
00339 
00341 
00342 
00343     virtual void print(ostream& os) const;
00344 
00346   protected:
00347     void copyMat(ScalarType* inputMatrix, int strideInput, int numRows, int numCols, ScalarType* outputMatrix, int strideOutput, int startRow, int startCol, ScalarType alpha = ScalarTraits<ScalarType>::zero() );
00348     void deleteArrays();
00349     void checkIndex( int rowIndex, int colIndex = 0 ) const;
00350     int numRows_;
00351     int numCols_;
00352     int stride_;
00353     bool valuesCopied_;
00354     ScalarType* values_;
00355 
00356   }; // class Teuchos_SerialDenseMatrix
00357 
00358   //----------------------------------------------------------------------------------------------------
00359   //  Constructors and Destructor
00360   //----------------------------------------------------------------------------------------------------
00361 
00362   template<typename OrdinalType, typename ScalarType>
00363   SerialDenseMatrix<OrdinalType, ScalarType>::SerialDenseMatrix() : CompObject(), numRows_(0), numCols_(0), stride_(0), valuesCopied_(false), values_(0) {}
00364   
00365   template<typename OrdinalType, typename ScalarType>
00366   SerialDenseMatrix<OrdinalType, ScalarType>::SerialDenseMatrix( int numRows, int numCols ) : CompObject(), numRows_(numRows), numCols_(numCols), stride_(numRows)
00367   {
00368     values_ = new ScalarType[stride_*numCols_];
00369     putScalar();
00370     valuesCopied_ = true;
00371   }
00372 
00373   template<typename OrdinalType, typename ScalarType>
00374   SerialDenseMatrix<OrdinalType, ScalarType>::SerialDenseMatrix(DataAccess CV, ScalarType* values, int stride, int numRows, int numCols) : CompObject(), numRows_(numRows), numCols_(numCols), stride_(stride), valuesCopied_(false), values_(values)
00375   {
00376     if(CV == Copy)
00377       {
00378   stride_ = numRows_;
00379   values_ = new ScalarType[stride_*numCols_];
00380   copyMat(values, stride, numRows_, numCols_, values_, stride_, 0, 0, false);
00381   valuesCopied_ = true;
00382       }
00383   }
00384   
00385   template<typename OrdinalType, typename ScalarType>
00386   SerialDenseMatrix<OrdinalType, ScalarType>::SerialDenseMatrix(const SerialDenseMatrix<OrdinalType, ScalarType> &Source, ETransp trans) : CompObject(), numRows_(0), numCols_(0), stride_(0), valuesCopied_(true), values_(0)
00387   {
00388     if ( trans == Teuchos::NO_TRANS ) 
00389       {
00390   numRows_ = Source.numRows_;
00391   numCols_ = Source.numCols_;
00392   stride_ = numRows_;
00393   values_ = new ScalarType[stride_*numCols_];
00394   copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0, false);
00395       } 
00396     else if ( trans == Teuchos::CONJ_TRANS && ScalarTraits<ScalarType>::isComplex ) 
00397       {   
00398   numRows_ = Source.numCols_;
00399   numCols_ = Source.numRows_;
00400   stride_ = numRows_;
00401   values_ = new ScalarType[stride_*numCols_];
00402   for (int j=0; j<numCols_; j++) {
00403     for (int i=0; i<numRows_; i++) {
00404       values_[j*stride_ + i] = Teuchos::ScalarTraits<ScalarType>::conjugate(Source.values_[i*Source.stride_ + j]);
00405     }
00406   }
00407       } 
00408     else 
00409       {
00410   numRows_ = Source.numCols_;
00411   numCols_ = Source.numRows_;
00412   stride_ = numRows_;
00413   values_ = new ScalarType[stride_*numCols_];
00414   for (int j=0; j<numCols_; j++) {
00415     for (int i=0; i<numRows_; i++) {
00416       values_[j*stride_ + i] = Source.values_[i*Source.stride_ + j];
00417     }
00418   }
00419       }
00420   }
00421   
00422   template<typename OrdinalType, typename ScalarType>
00423   SerialDenseMatrix<OrdinalType, ScalarType>::SerialDenseMatrix(DataAccess CV, const SerialDenseMatrix<OrdinalType, ScalarType> &Source, int numRows, int numCols, int startRow, int startCol) : CompObject(), numRows_(numRows), numCols_(numCols), stride_(Source.stride_), valuesCopied_(false), values_(Source.values_)
00424   {
00425     if(CV == Copy)
00426       {
00427   stride_ = numRows;
00428   values_ = new ScalarType[stride_ * numCols];
00429   copyMat(Source.values_, Source.stride_, numRows, numCols, values_, stride_, startRow, startCol, false);
00430   valuesCopied_ = true;
00431       }
00432     else // CV == View
00433       {
00434   values_ = values_ + (stride_ * startCol) + startRow;
00435       }
00436   }
00437     
00438   template<typename OrdinalType, typename ScalarType>
00439   SerialDenseMatrix<OrdinalType, ScalarType>::~SerialDenseMatrix()
00440   {
00441     deleteArrays();
00442   }
00443   
00444   //----------------------------------------------------------------------------------------------------
00445   //  Shape methods 
00446   //----------------------------------------------------------------------------------------------------
00447 
00448   template<typename OrdinalType, typename ScalarType>
00449   int SerialDenseMatrix<OrdinalType, ScalarType>::shape(int numRows, int numCols)
00450   {
00451     deleteArrays(); // Get rid of anything that might be already allocated
00452     numRows_ = numRows;
00453     numCols_ = numCols;
00454     stride_ = numRows_;
00455     values_ = new ScalarType[stride_*numCols_];
00456     putScalar();
00457     valuesCopied_ = true;
00458     return(0);
00459   }
00460 
00461   template<typename OrdinalType, typename ScalarType>
00462   int SerialDenseMatrix<OrdinalType, ScalarType>::shapeUninitialized(int numRows, int numCols)
00463   {
00464     deleteArrays(); // Get rid of anything that might be already allocated
00465     numRows_ = numRows;
00466     numCols_ = numCols;
00467     stride_ = numRows_;
00468     values_ = new ScalarType[stride_*numCols_];
00469     valuesCopied_ = true;
00470     return(0);
00471   }
00472   
00473   template<typename OrdinalType, typename ScalarType>
00474   int SerialDenseMatrix<OrdinalType, ScalarType>::reshape(int numRows, int numCols)
00475   {
00476     // Allocate space for new matrix
00477     ScalarType* values_tmp = new ScalarType[numRows * numCols];
00478     ScalarType zero = ScalarTraits<ScalarType>::zero();
00479     for(int k = 0; k < numRows * numCols; k++)
00480       {
00481   values_tmp[k] = zero;
00482       }
00483     int numRows_tmp = TEUCHOS_MIN(numRows_, numRows);
00484     int numCols_tmp = TEUCHOS_MIN(numCols_, numCols);
00485     if(values_ != 0)
00486       {
00487   copyMat(values_, stride_, numRows_tmp, numCols_tmp, values_tmp, numRows, 0, 0, false); // Copy principal submatrix of A to new A
00488       }
00489     deleteArrays(); // Get rid of anything that might be already allocated
00490     numRows_ = numRows;
00491     numCols_ = numCols;
00492     stride_ = numRows_;
00493     values_ = values_tmp; // Set pointer to new A
00494     valuesCopied_ = true;
00495     return(0);
00496   }
00497    
00498   //----------------------------------------------------------------------------------------------------
00499   //   Set methods 
00500   //----------------------------------------------------------------------------------------------------
00501 
00502   template<typename OrdinalType, typename ScalarType>
00503   int SerialDenseMatrix<OrdinalType, ScalarType>::putScalar( const ScalarType value )
00504   {
00505     // Set each value of the dense matrix to "value".
00506     for(int j = 0; j < numCols_; j++) 
00507       {
00508   for(int i = 0; i < numRows_; i++) 
00509     {
00510       values_[i + j*stride_] = value;
00511     }
00512       }
00513     return 0;
00514   }    
00515     
00516   template<typename OrdinalType, typename ScalarType>
00517   int SerialDenseMatrix<OrdinalType, ScalarType>::random()
00518   {
00519     // Set each value of the dense matrix to a random value.
00520     for(int j = 0; j < numCols_; j++) 
00521       {
00522   for(int i = 0; i < numRows_; i++) 
00523     {
00524       values_[i + j*stride_] = ScalarTraits<ScalarType>::random();
00525     }
00526       }
00527     return 0;
00528   }
00529 
00530   template<typename OrdinalType, typename ScalarType>
00531   SerialDenseMatrix<OrdinalType,ScalarType>& SerialDenseMatrix<OrdinalType, ScalarType>::operator= (const SerialDenseMatrix<OrdinalType,ScalarType>& Source) {
00532     if(this == &Source)
00533       return(*this); // Special case of source same as target
00534     if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
00535       return(*this); // Special case of both are views to same data.
00536 
00537     // If the source is a view then we will return a view, else we will return a copy.
00538     if (!Source.valuesCopied_) {
00539       if(valuesCopied_)       { 
00540   // Clean up stored data if this was previously a copy.
00541   deleteArrays();
00542       }
00543       numRows_ = Source.numRows_; 
00544       numCols_ = Source.numCols_;
00545       stride_ = Source.stride_;
00546       values_ = Source.values_;
00547     }
00548     else {
00549       // If we were a view, we will now be a copy.
00550       if(!valuesCopied_) {
00551   numRows_ = Source.numRows_;
00552   numCols_ = Source.numCols_;
00553   stride_ = Source.numRows_;
00554   const int newsize = stride_ * numCols_;
00555   if(newsize > 0) {
00556     values_ = new ScalarType[newsize];
00557     valuesCopied_ = true;
00558   }
00559   else {
00560     values_ = 0;
00561   }
00562       }
00563       // If we were a copy, we will stay a copy.
00564       else {
00565   if((Source.numRows_ <= stride_) && (Source.numCols_ == numCols_)) { // we don't need to reallocate
00566     numRows_ = Source.numRows_;
00567     numCols_ = Source.numCols_;
00568   }
00569   else { // we need to allocate more space (or less space)
00570     deleteArrays();
00571     numRows_ = Source.numRows_;
00572     numCols_ = Source.numCols_;
00573     stride_ = Source.numRows_;
00574     const int newsize = stride_ * numCols_;
00575     if(newsize > 0) {
00576       values_ = new ScalarType[newsize];
00577       valuesCopied_ = true;
00578     }
00579   }
00580       }
00581       copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0, false);
00582     } 
00583     return(*this);
00584   }
00585 
00586   template<typename OrdinalType, typename ScalarType>
00587   SerialDenseMatrix<OrdinalType, ScalarType>& SerialDenseMatrix<OrdinalType, ScalarType>::operator+= (const SerialDenseMatrix<OrdinalType,ScalarType>& Source )
00588   {
00589     // Check for compatible dimensions
00590     if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_))
00591       {
00592   TEUCHOS_CHK_REF(*this); // Return *this without altering it.
00593       }
00594     copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0, 1.0);
00595     return(*this);
00596   }
00597 
00598   template<typename OrdinalType, typename ScalarType>
00599   SerialDenseMatrix<OrdinalType, ScalarType>& SerialDenseMatrix<OrdinalType, ScalarType>::operator-= (const SerialDenseMatrix<OrdinalType,ScalarType>& Source )
00600   {
00601     // Check for compatible dimensions
00602     if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_))
00603       {
00604   TEUCHOS_CHK_REF(*this); // Return *this without altering it.
00605       }
00606     copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0, -1.0);
00607     return(*this);
00608   }
00609 
00610   template<typename OrdinalType, typename ScalarType>
00611   SerialDenseMatrix<OrdinalType,ScalarType>& SerialDenseMatrix<OrdinalType, ScalarType>::assign (const SerialDenseMatrix<OrdinalType,ScalarType>& Source) {
00612     if(this == &Source)
00613       return(*this); // Special case of source same as target
00614     if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
00615       return(*this); // Special case of both are views to same data.
00616 
00617     // Check for compatible dimensions
00618     if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_))
00619       {
00620   TEUCHOS_CHK_REF(*this); // Return *this without altering it.
00621       }
00622     copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0 );
00623     return(*this);
00624   }
00625 
00626   //----------------------------------------------------------------------------------------------------
00627   //   Accessor methods 
00628   //----------------------------------------------------------------------------------------------------
00629 
00630   template<typename OrdinalType, typename ScalarType>
00631   inline ScalarType& SerialDenseMatrix<OrdinalType, ScalarType>::operator () (int rowIndex, int colIndex)
00632   {
00633 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
00634     checkIndex( rowIndex, colIndex );
00635 #endif
00636     return(values_[colIndex * stride_ + rowIndex]);
00637   }
00638   
00639   template<typename OrdinalType, typename ScalarType>
00640   inline const ScalarType& SerialDenseMatrix<OrdinalType, ScalarType>::operator () (int rowIndex, int colIndex) const
00641   {
00642 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
00643     checkIndex( rowIndex, colIndex );
00644 #endif
00645     return(values_[colIndex * stride_ + rowIndex]);
00646   }
00647   
00648   template<typename OrdinalType, typename ScalarType>
00649   inline const ScalarType* SerialDenseMatrix<OrdinalType, ScalarType>::operator [] (int colIndex) const
00650   {
00651 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
00652     checkIndex( 0, colIndex );
00653 #endif
00654     return(values_ + colIndex * stride_);
00655   }
00656   
00657   template<typename OrdinalType, typename ScalarType>
00658   inline ScalarType* SerialDenseMatrix<OrdinalType, ScalarType>::operator [] (int colIndex)
00659   {
00660 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
00661     checkIndex( 0, colIndex );
00662 #endif
00663   return(values_ + colIndex * stride_);
00664   }
00665 
00666   //----------------------------------------------------------------------------------------------------
00667   //   Norm methods 
00668   //----------------------------------------------------------------------------------------------------
00669   
00670   template<typename OrdinalType, typename ScalarType>
00671   typename ScalarTraits<ScalarType>::magnitudeType SerialDenseMatrix<OrdinalType, ScalarType>::normOne() const
00672   {
00673     int i, j;
00674     typename ScalarTraits<ScalarType>::magnitudeType anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00675     typename ScalarTraits<ScalarType>::magnitudeType absSum = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00676     ScalarType* ptr;
00677     for(j = 0; j < numCols_; j++)
00678       {
00679   ScalarType sum = 0;
00680   ptr = values_ + j * stride_;
00681   for(i = 0; i < numRows_; i++)
00682     {
00683       sum += ScalarTraits<ScalarType>::magnitude(*ptr++);
00684     }
00685   absSum = ScalarTraits<ScalarType>::magnitude(sum);
00686   if(absSum > anorm)
00687     {
00688       anorm = absSum;
00689     }
00690       }
00691     updateFlops(numRows_ * numCols_);
00692     return(anorm);
00693   }
00694   
00695   template<typename OrdinalType, typename ScalarType>
00696   typename ScalarTraits<ScalarType>::magnitudeType SerialDenseMatrix<OrdinalType, ScalarType>::normInf() const
00697   {
00698     int i, j;
00699     typename ScalarTraits<ScalarType>::magnitudeType sum, anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00700     
00701     for (i = 0; i < numRows_; i++) {
00702       sum = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00703       for (j=0; j< numCols_; j++) {
00704   sum += ScalarTraits<ScalarType>::magnitude(*(values_+i+j*stride_));
00705       }
00706       anorm = TEUCHOS_MAX( anorm, sum );
00707     }
00708     updateFlops(numRows_ * numCols_);
00709     return(anorm);
00710   }
00711   
00712   template<typename OrdinalType, typename ScalarType>
00713   typename ScalarTraits<ScalarType>::magnitudeType SerialDenseMatrix<OrdinalType, ScalarType>::normFrobenius() const
00714   {
00715     int i, j;
00716     typename ScalarTraits<ScalarType>::magnitudeType anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00717     for (j = 0; j < numCols_; j++) {
00718       for (i = 0; i < numRows_; i++) {
00719   anorm += ScalarTraits<ScalarType>::magnitude(values_[i+j*stride_]*values_[i+j*stride_]);
00720       }
00721     }
00722     anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::squareroot(anorm));
00723     updateFlops(numRows_ * numCols_);
00724     return(anorm);
00725   }
00726   
00727   //----------------------------------------------------------------------------------------------------
00728   //   Comparison methods 
00729   //----------------------------------------------------------------------------------------------------
00730   
00731   template<typename OrdinalType, typename ScalarType>
00732   bool SerialDenseMatrix<OrdinalType, ScalarType>::operator== (const SerialDenseMatrix<OrdinalType, ScalarType> &Operand)
00733   {
00734     bool result = 1;
00735     if((numRows_ != Operand.numRows_) || (numCols_ != Operand.numCols_))
00736       {
00737   result = 0;
00738       }
00739     else
00740       {
00741   int i, j;
00742   for(i = 0; i < numRows_; i++)
00743     {
00744       for(j = 0; j < numCols_; j++)
00745         {
00746     if((*this)(i, j) != Operand(i, j))
00747       {
00748         return 0;
00749       }
00750         }
00751     }
00752       }
00753     return result;
00754   }
00755   
00756   template<typename OrdinalType, typename ScalarType>
00757   bool SerialDenseMatrix<OrdinalType, ScalarType>::operator!= (const SerialDenseMatrix<OrdinalType, ScalarType> &Operand)
00758   {
00759     return !((*this) == Operand);
00760   }
00761   
00762   //----------------------------------------------------------------------------------------------------
00763   //   Multiplication method 
00764   //----------------------------------------------------------------------------------------------------  
00765 
00766   template<typename OrdinalType, typename ScalarType>
00767   int SerialDenseMatrix<OrdinalType, ScalarType>::scale( const ScalarType alpha )
00768   {
00769     int i, j;
00770     ScalarType* ptr;
00771     
00772     for (j=0; j<numCols_; j++) {
00773       ptr = values_ + j*stride_;
00774       for (i=0; i<numRows_; i++) { *ptr = alpha * (*ptr); ptr++; }
00775     }
00776     updateFlops( numRows_*numCols_ );
00777     return(0);
00778   }
00779 
00780   template<typename OrdinalType, typename ScalarType>
00781   int SerialDenseMatrix<OrdinalType, ScalarType>::scale( const SerialDenseMatrix<OrdinalType,ScalarType>& A )
00782   {
00783     int i, j;
00784     ScalarType* ptr;
00785     
00786     // Check for compatible dimensions
00787     if ((numRows_ != A.numRows_) || (numCols_ != A.numCols_))
00788       {
00789         TEUCHOS_CHK_ERR(-1); // Return error
00790       }    
00791     for (j=0; j<numCols_; j++) {
00792       ptr = values_ + j*stride_;
00793       for (i=0; i<numRows_; i++) { *ptr = A(i,j) * (*ptr); ptr++; }
00794     }
00795     updateFlops( numRows_*numCols_ );
00796     return(0);
00797   }
00798 
00799   template<typename OrdinalType, typename ScalarType>
00800   int  SerialDenseMatrix<OrdinalType, ScalarType>::multiply(ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix<OrdinalType, ScalarType> &A, const SerialDenseMatrix<OrdinalType, ScalarType> &B, ScalarType beta)
00801   {
00802     // Check for compatible dimensions
00803     int A_nrows = (ETranspChar[transa]!='N') ? A.numCols() : A.numRows();
00804     int A_ncols = (ETranspChar[transa]!='N') ? A.numRows() : A.numCols();
00805     int B_nrows = (ETranspChar[transb]!='N') ? B.numCols() : B.numRows();
00806     int B_ncols = (ETranspChar[transb]!='N') ? B.numRows() : B.numCols();
00807     if ((numRows_ != A_nrows) || (A_ncols != B_nrows) || (numCols_ != B_ncols))
00808       {
00809         TEUCHOS_CHK_ERR(-1); // Return error
00810       }
00811     // Call GEMM function
00812     this->GEMM(transa, transb, numRows_, numCols_, A_ncols, alpha, A.values(), A.stride(), B.values(), B.stride(), beta, values_, stride_);
00813     double nflops = 2 * numRows_;
00814     nflops *= numCols_;
00815     nflops *= A_ncols;
00816     updateFlops(nflops);
00817     return(0);
00818   }
00819   
00820   
00821   template<typename OrdinalType, typename ScalarType>
00822   void SerialDenseMatrix<OrdinalType, ScalarType>::print(ostream& os) const
00823   {
00824     os << endl;
00825     if(valuesCopied_)
00826       os << "Values_copied : yes" << endl;
00827     else
00828       os << "Values_copied : no" << endl;
00829       os << "Rows : " << numRows_ << endl;
00830       os << "Columns : " << numCols_ << endl;
00831       os << "LDA : " << stride_ << endl;
00832     if(numRows_ == 0 || numCols_ == 0) {
00833       os << "(matrix is empty, no values to display)" << endl;
00834     } else {
00835       for(int i = 0; i < numRows_; i++) {
00836   for(int j = 0; j < numCols_; j++){
00837     os << (*this)(i,j) << " ";
00838   }
00839   os << endl;
00840       }
00841     }
00842   }
00843 
00844   //----------------------------------------------------------------------------------------------------
00845   //   Protected methods 
00846   //----------------------------------------------------------------------------------------------------  
00847 
00848   template<typename OrdinalType, typename ScalarType>
00849   inline void SerialDenseMatrix<OrdinalType, ScalarType>::checkIndex( int rowIndex, int colIndex ) const {
00850     TEST_FOR_EXCEPTION(rowIndex < 0 || rowIndex >= numRows_, std::out_of_range,
00851                        "SerialDenseMatrix<T>::checkIndex: "
00852                        "Row index " << rowIndex << " out of range [0, "<< numRows_ << ")");
00853     TEST_FOR_EXCEPTION(colIndex < 0 || colIndex >= numCols_, std::out_of_range,
00854                        "SerialDenseMatrix<T>::checkIndex: "
00855                        "Col index " << colIndex << " out of range [0, "<< numCols_ << ")");
00856   }
00857 
00858   template<typename OrdinalType, typename ScalarType>
00859   void SerialDenseMatrix<OrdinalType, ScalarType>::deleteArrays(void)
00860   {
00861     if (valuesCopied_)
00862       {
00863   delete [] values_;
00864   values_ = 0;
00865   valuesCopied_ = false;
00866       }
00867   }
00868   
00869   template<typename OrdinalType, typename ScalarType>
00870   void SerialDenseMatrix<OrdinalType, ScalarType>::copyMat(ScalarType* inputMatrix, int strideInput, int numRows, int numCols, ScalarType* outputMatrix, int strideOutput, int startRow, int startCol, ScalarType alpha)
00871   {
00872     int i, j;
00873     ScalarType* ptr1 = 0;
00874     ScalarType* ptr2 = 0;
00875     for(j = 0; j < numCols; j++) {
00876   ptr1 = outputMatrix + (j * strideOutput);
00877   ptr2 = inputMatrix + (j + startCol) * strideInput + startRow;
00878   if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
00879     for(i = 0; i < numRows; i++)
00880       {
00881         *ptr1++ += alpha*(*ptr2++);
00882       }
00883   } else {
00884     for(i = 0; i < numRows; i++)
00885       {
00886         *ptr1++ = *ptr2++;
00887       }
00888   }
00889     }
00890   }
00891   
00892 } // namespace Teuchos
00893 
00894 
00895 #endif /* _TEUCHOS_SERIALDENSEMATRIX_HPP_ */