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