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