Teuchos Package Browser (Single Doxygen Collection) Version of the Day
Teuchos_SerialBandDenseMatrix.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                    Teuchos: Common Tools Package
00005 //                 Copyright (2004) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #ifndef _TEUCHOS_SERIALBANDDENSEMATRIX_HPP_
00043 #define _TEUCHOS_SERIALBANDDENSEMATRIX_HPP_
00044 
00048 #include "Teuchos_CompObject.hpp"
00049 #include "Teuchos_BLAS.hpp"
00050 #include "Teuchos_ScalarTraits.hpp"
00051 #include "Teuchos_DataAccess.hpp"
00052 #include "Teuchos_ConfigDefs.hpp"
00053 #include "Teuchos_RCP.hpp"
00054 #include "Teuchos_Assert.hpp"
00055 
00130 namespace Teuchos {
00131 
00132 template<typename OrdinalType, typename ScalarType>
00133 class SerialBandDenseMatrix : public CompObject, public Object, public BLAS<OrdinalType, ScalarType>
00134 {
00135 public:
00136 
00138   typedef OrdinalType ordinalType;
00140   typedef ScalarType scalarType;
00141 
00143 
00144 
00146 
00149   SerialBandDenseMatrix();
00150 
00152 
00162   SerialBandDenseMatrix(OrdinalType numRows, OrdinalType numCols, OrdinalType kl, OrdinalType ku, bool zeroOut = true);
00163 
00165 
00175   SerialBandDenseMatrix(DataAccess CV, ScalarType* values, OrdinalType stride, OrdinalType numRows, OrdinalType numCols, OrdinalType kl, OrdinalType ku);
00176 
00178 
00183   SerialBandDenseMatrix(const SerialBandDenseMatrix<OrdinalType, ScalarType> &Source, ETransp trans = Teuchos::NO_TRANS);
00184 
00186 
00200   SerialBandDenseMatrix(DataAccess CV, const SerialBandDenseMatrix<OrdinalType, ScalarType> &Source, OrdinalType numRows, OrdinalType numCols, OrdinalType startCol=0);
00201 
00203   virtual ~SerialBandDenseMatrix();
00205 
00207 
00208 
00209 
00221   int shape(OrdinalType numRows, OrdinalType numCols, OrdinalType kl, OrdinalType ku);
00222 
00224   int shapeUninitialized(OrdinalType numRows, OrdinalType numCols, OrdinalType kl, OrdinalType ku);
00225 
00227 
00239   int reshape(OrdinalType numRows, OrdinalType numCols, OrdinalType kl, OrdinalType ku);
00240 
00241 
00243 
00245 
00246 
00248 
00254   SerialBandDenseMatrix<OrdinalType, ScalarType>& operator= (const SerialBandDenseMatrix<OrdinalType, ScalarType>& Source);
00255 
00257 
00262   SerialBandDenseMatrix<OrdinalType, ScalarType>& assign (const SerialBandDenseMatrix<OrdinalType, ScalarType>& Source);
00263 
00265 
00268   SerialBandDenseMatrix<OrdinalType, ScalarType>& operator= (const ScalarType value) { putScalar(value); return(*this); }
00269 
00271 
00275   int putScalar( const ScalarType value = Teuchos::ScalarTraits<ScalarType>::zero() );
00276 
00278   int random();
00279 
00281 
00283 
00284 
00286 
00293   ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex);
00294 
00296 
00303   const ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex) const;
00304 
00306 
00313   ScalarType* operator [] (OrdinalType colIndex);
00314 
00316 
00323   const ScalarType* operator [] (OrdinalType colIndex) const;
00324 
00326 
00327   ScalarType* values() const { return(values_); }
00328 
00330 
00332 
00333 
00335 
00338   SerialBandDenseMatrix<OrdinalType, ScalarType>& operator+= (const SerialBandDenseMatrix<OrdinalType, ScalarType>& Source);
00339 
00341 
00344   SerialBandDenseMatrix<OrdinalType, ScalarType>& operator-= (const SerialBandDenseMatrix<OrdinalType, ScalarType>& Source);
00345 
00347 
00350   SerialBandDenseMatrix<OrdinalType, ScalarType>& operator*= (const ScalarType alpha);
00351 
00353 
00357   int scale ( const ScalarType alpha );
00358 
00360 
00366   int scale ( const SerialBandDenseMatrix<OrdinalType, ScalarType>& A );
00367 
00369 
00371 
00372 
00374 
00377   bool operator== (const SerialBandDenseMatrix<OrdinalType, ScalarType> &Operand) const;
00378 
00380 
00383   bool operator!= (const SerialBandDenseMatrix<OrdinalType, ScalarType> &Operand) const;
00384 
00386 
00388 
00389 
00391   OrdinalType numRows() const { return(numRows_); }
00392 
00394   OrdinalType numCols() const { return(numCols_); }
00395 
00397   OrdinalType lowerBandwidth() const { return(kl_); }
00398 
00400   OrdinalType upperBandwidth() const { return(ku_); }
00401 
00403   OrdinalType stride() const { return(stride_); }
00404 
00406   bool empty() const { return(numRows_ == 0 || numCols_ == 0); }
00408 
00410 
00411 
00413   typename ScalarTraits<ScalarType>::magnitudeType normOne() const;
00414 
00416   typename ScalarTraits<ScalarType>::magnitudeType normInf() const;
00417 
00419   typename ScalarTraits<ScalarType>::magnitudeType normFrobenius() const;
00421 
00423 
00424 
00425   virtual void print(std::ostream& os) const;
00426 
00428 protected:
00429   void copyMat(ScalarType* inputMatrix, OrdinalType strideInput,
00430     OrdinalType numRows, OrdinalType numCols, ScalarType* outputMatrix,
00431     OrdinalType strideOutput, OrdinalType startCol, ScalarType alpha = ScalarTraits<ScalarType>::zero() );
00432   void deleteArrays();
00433   void checkIndex( OrdinalType rowIndex, OrdinalType colIndex = 0 ) const;
00434   OrdinalType numRows_;
00435   OrdinalType numCols_;
00436   OrdinalType stride_;
00437   OrdinalType kl_;
00438   OrdinalType ku_;
00439   bool valuesCopied_;
00440   ScalarType* values_;
00441 
00442 }; // class Teuchos_SerialBandDenseMatrix
00443 
00444 //----------------------------------------------------------------------------------------------------
00445 //  Constructors and Destructor
00446 //----------------------------------------------------------------------------------------------------
00447 
00448 template<typename OrdinalType, typename ScalarType>
00449 SerialBandDenseMatrix<OrdinalType, ScalarType>::SerialBandDenseMatrix()
00450   : CompObject(), Object("Teuchos::SerialBandDenseMatrix"), numRows_(0), numCols_(0), kl_(0), ku_(0), stride_(0), valuesCopied_(false), values_(0)
00451 {}
00452 
00453 template<typename OrdinalType, typename ScalarType>
00454 SerialBandDenseMatrix<OrdinalType, ScalarType>::SerialBandDenseMatrix(
00455   OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType kl_in, OrdinalType ku_in, bool zeroOut
00456   )
00457   : CompObject(), Object("Teuchos::SerialBandDenseMatrix"), numRows_(numRows_in), numCols_(numCols_in), kl_(kl_in), ku_(ku_in), stride_(kl_in+ku_in+1)
00458 {
00459   values_ = new ScalarType[stride_*numCols_];
00460   if (zeroOut == true)
00461     putScalar();
00462   valuesCopied_ = true;
00463 }
00464 
00465 template<typename OrdinalType, typename ScalarType>
00466 SerialBandDenseMatrix<OrdinalType, ScalarType>::SerialBandDenseMatrix(
00467   DataAccess CV, ScalarType* values_in, OrdinalType stride_in, OrdinalType numRows_in,
00468   OrdinalType numCols_in, OrdinalType kl_in, OrdinalType ku_in
00469   )
00470   : CompObject(), Object("Teuchos::SerialBandDenseMatrix"), numRows_(numRows_in), numCols_(numCols_in), kl_(kl_in), ku_(ku_in), stride_(stride_in), valuesCopied_(false), values_(values_in)
00471 {
00472   if(CV == Copy) {
00473     stride_ = kl_+ku_+1;
00474     values_ = new ScalarType[stride_*numCols_];
00475     copyMat(values_in, stride_in, numRows_, numCols_, values_, stride_, 0 );
00476     valuesCopied_ = true;
00477   }
00478 }
00479 
00480 template<typename OrdinalType, typename ScalarType>
00481 SerialBandDenseMatrix<OrdinalType, ScalarType>::SerialBandDenseMatrix(const SerialBandDenseMatrix<OrdinalType, ScalarType> &Source, ETransp trans) : CompObject(), Object("Teuchos::SerialBandDenseMatrix"), numRows_(0), numCols_(0), kl_(0), ku_(0), stride_(0), valuesCopied_(true), values_(0)
00482 {
00483 
00484   if ( trans == Teuchos::NO_TRANS ) {
00485 
00486     numRows_ = Source.numRows_;
00487     numCols_ = Source.numCols_;
00488     kl_ = Source.kl_;
00489     ku_ = Source.ku_;
00490     stride_ = kl_+ku_+1;
00491     values_ = new ScalarType[stride_*numCols_];
00492     copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0);
00493 
00494   } else if ( trans == Teuchos::CONJ_TRANS && ScalarTraits<ScalarType>::isComplex ) {
00495 
00496     numRows_ = Source.numCols_;
00497     numCols_ = Source.numRows_;
00498     kl_ = Source.ku_;
00499     ku_ = Source.kl_;
00500     stride_ = kl_+ku_+1;
00501     values_ = new ScalarType[stride_*numCols_];
00502     for (OrdinalType j=0; j<numCols_; j++) {
00503       for (OrdinalType i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
00504   values_[j*stride_ + (ku_+i-j)] =
00505     Teuchos::ScalarTraits<ScalarType>::conjugate(Source.values_[i*Source.stride_ + (Source.ku_+j-i)]);
00506       }
00507     }
00508 
00509   } else {
00510 
00511     numRows_ = Source.numCols_;
00512     numCols_ = Source.numRows_;
00513     kl_ = Source.ku_;
00514     ku_ = Source.kl_;
00515     stride_ = kl_+ku_+1;
00516     values_ = new ScalarType[stride_*numCols_];
00517     for (OrdinalType j=0; j<numCols_; j++) {
00518       for (OrdinalType i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
00519   values_[j*stride_ + (ku_+i-j)] = Source.values_[i*Source.stride_ + (Source.ku_+j-i)];
00520       }
00521     }
00522 
00523   }
00524 
00525 }
00526 
00527 template<typename OrdinalType, typename ScalarType>
00528 SerialBandDenseMatrix<OrdinalType, ScalarType>::SerialBandDenseMatrix(
00529   DataAccess CV, const SerialBandDenseMatrix<OrdinalType, ScalarType> &Source,
00530   OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType startCol )
00531   : CompObject(), Object("Teuchos::SerialBandDenseMatrix"), numRows_(numRows_in), numCols_(numCols_in), kl_(Source.kl_), ku_(Source.ku_), stride_(Source.stride_), valuesCopied_(false), values_(Source.values_)
00532 {
00533   if(CV == Copy) {
00534       values_ = new ScalarType[stride_ * numCols_in];
00535       copyMat(Source.values_, Source.stride_, numRows_in, numCols_in, values_, stride_, startCol);
00536       valuesCopied_ = true;
00537   } else { // CV = View
00538     values_ = values_ + (stride_ * startCol);
00539   }
00540 }
00541 
00542 template<typename OrdinalType, typename ScalarType>
00543 SerialBandDenseMatrix<OrdinalType, ScalarType>::~SerialBandDenseMatrix()
00544 {
00545   deleteArrays();
00546 }
00547 
00548 //----------------------------------------------------------------------------------------------------
00549 //  Shape methods
00550 //----------------------------------------------------------------------------------------------------
00551 
00552 template<typename OrdinalType, typename ScalarType>
00553 int SerialBandDenseMatrix<OrdinalType, ScalarType>::shape(
00554   OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType kl_in, OrdinalType ku_in
00555   )
00556 {
00557   deleteArrays(); // Get rid of anything that might be already allocated
00558   numRows_ = numRows_in;
00559   numCols_ = numCols_in;
00560   kl_ = kl_in;
00561   ku_ = ku_in;
00562   stride_ = kl_+ku_+1;
00563   values_ = new ScalarType[stride_*numCols_];
00564   putScalar();
00565   valuesCopied_ = true;
00566 
00567   return(0);
00568 }
00569 
00570 template<typename OrdinalType, typename ScalarType>
00571 int SerialBandDenseMatrix<OrdinalType, ScalarType>::shapeUninitialized(
00572   OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType kl_in, OrdinalType ku_in
00573   )
00574 {
00575   deleteArrays(); // Get rid of anything that might be already allocated
00576   numRows_ = numRows_in;
00577   numCols_ = numCols_in;
00578   kl_ = kl_in;
00579   ku_ = ku_in;
00580   stride_ = kl_+ku_+1;
00581   values_ = new ScalarType[stride_*numCols_];
00582   valuesCopied_ = true;
00583 
00584   return(0);
00585 }
00586 
00587 template<typename OrdinalType, typename ScalarType>
00588 int SerialBandDenseMatrix<OrdinalType, ScalarType>::reshape(
00589   OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType kl_in, OrdinalType ku_in
00590   )
00591 {
00592 
00593   // Allocate space for new matrix
00594   ScalarType* values_tmp = new ScalarType[(kl_in+ku_in+1) * numCols_in];
00595   ScalarType zero = ScalarTraits<ScalarType>::zero();
00596   for(OrdinalType k = 0; k < (kl_in+ku_in+1) * numCols_in; k++) {
00597     values_tmp[k] = zero;
00598   }
00599   OrdinalType numRows_tmp = TEUCHOS_MIN(numRows_, numRows_in);
00600   OrdinalType numCols_tmp = TEUCHOS_MIN(numCols_, numCols_in);
00601   if(values_ != 0) {
00602     copyMat(values_, stride_, numRows_tmp, numCols_tmp, values_tmp,
00603       kl_in+ku_in+1, 0); // Copy principal submatrix of A to new A
00604   }
00605   deleteArrays(); // Get rid of anything that might be already allocated
00606   numRows_ = numRows_in;
00607   numCols_ = numCols_in;
00608   kl_ = kl_in;
00609   ku_ = ku_in;
00610   stride_ = kl_+ku_+1;
00611   values_ = values_tmp; // Set pointer to new A
00612   valuesCopied_ = true;
00613 
00614   return(0);
00615 }
00616 
00617 //----------------------------------------------------------------------------------------------------
00618 //   Set methods
00619 //----------------------------------------------------------------------------------------------------
00620 
00621 template<typename OrdinalType, typename ScalarType>
00622 int SerialBandDenseMatrix<OrdinalType, ScalarType>::putScalar( const ScalarType value_in )
00623 {
00624 
00625   // Set each value of the dense matrix to "value".
00626   for(OrdinalType j = 0; j < numCols_; j++) {
00627     for (OrdinalType i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
00628       values_[(ku_+i-j) + j*stride_] = value_in;
00629     }
00630   }
00631 
00632   return 0;
00633 }
00634 
00635 template<typename OrdinalType, typename ScalarType>
00636 int SerialBandDenseMatrix<OrdinalType, ScalarType>::random()
00637 {
00638 
00639   // Set each value of the dense matrix to a random value.
00640   for(OrdinalType j = 0; j < numCols_; j++) {
00641     for (OrdinalType i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
00642       values_[(ku_+i-j) + j*stride_] = ScalarTraits<ScalarType>::random();
00643     }
00644   }
00645 
00646   return 0;
00647 }
00648 
00649 template<typename OrdinalType, typename ScalarType>
00650 SerialBandDenseMatrix<OrdinalType,ScalarType>&
00651 SerialBandDenseMatrix<OrdinalType, ScalarType>::operator=(
00652   const SerialBandDenseMatrix<OrdinalType,ScalarType>& Source
00653   )
00654 {
00655 
00656   if(this == &Source)
00657     return(*this); // Special case of source same as target
00658   if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
00659     return(*this); // Special case of both are views to same data.
00660 
00661   // If the source is a view then we will return a view, else we will return a copy.
00662   if (!Source.valuesCopied_) {
00663     if(valuesCopied_) {
00664       // Clean up stored data if this was previously a copy.
00665       deleteArrays();
00666     }
00667     numRows_ = Source.numRows_;
00668     numCols_ = Source.numCols_;
00669     kl_ = Source.kl_;
00670     ku_ = Source.ku_;
00671     stride_ = Source.stride_;
00672     values_ = Source.values_;
00673   } else {
00674     // If we were a view, we will now be a copy.
00675     if(!valuesCopied_) {
00676       numRows_ = Source.numRows_;
00677       numCols_ = Source.numCols_;
00678       kl_ = Source.kl_;
00679       ku_ = Source.ku_;
00680       stride_ = kl_+ku_+1;
00681       const OrdinalType newsize = stride_ * numCols_;
00682       if(newsize > 0) {
00683   values_ = new ScalarType[newsize];
00684   valuesCopied_ = true;
00685       } else {
00686   values_ = 0;
00687       }
00688     } else {
00689       // If we were a copy, we will stay a copy.
00690       if((Source.numRows_ <= stride_) && (Source.numCols_ == numCols_)) { // we don't need to reallocate
00691   numRows_ = Source.numRows_;
00692   numCols_ = Source.numCols_;
00693   kl_ = Source.kl_;
00694   ku_ = Source.ku_;
00695       } else {
00696   // we need to allocate more space (or less space)
00697   deleteArrays();
00698   numRows_ = Source.numRows_;
00699   numCols_ = Source.numCols_;
00700   kl_ = Source.kl_;
00701   ku_ = Source.ku_;
00702   stride_ = kl_+ku_+1;
00703   const OrdinalType newsize = stride_ * numCols_;
00704   if(newsize > 0) {
00705     values_ = new ScalarType[newsize];
00706     valuesCopied_ = true;
00707   }
00708       }
00709     }
00710     copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0);
00711   }
00712   return(*this);
00713 
00714 }
00715 
00716 template<typename OrdinalType, typename ScalarType>
00717 SerialBandDenseMatrix<OrdinalType, ScalarType>& SerialBandDenseMatrix<OrdinalType, ScalarType>::operator+= (const SerialBandDenseMatrix<OrdinalType,ScalarType>& Source )
00718 {
00719 
00720   // Check for compatible dimensions
00721   if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_) || (kl_ != Source.kl_) || (ku_ != Source.ku_)) {
00722     TEUCHOS_CHK_REF(*this); // Return *this without altering it.
00723   }
00724   copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, ScalarTraits<ScalarType>::one());
00725   return(*this);
00726 
00727 }
00728 
00729 template<typename OrdinalType, typename ScalarType>
00730 SerialBandDenseMatrix<OrdinalType, ScalarType>& SerialBandDenseMatrix<OrdinalType, ScalarType>::operator-= (const SerialBandDenseMatrix<OrdinalType,ScalarType>& Source )
00731 {
00732 
00733   // Check for compatible dimensions
00734   if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_) || (kl_ != Source.kl_) || (ku_ != Source.ku_)) {
00735     TEUCHOS_CHK_REF(*this); // Return *this without altering it.
00736   }
00737   copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, -ScalarTraits<ScalarType>::one());
00738   return(*this);
00739 
00740 }
00741 
00742 template<typename OrdinalType, typename ScalarType>
00743 SerialBandDenseMatrix<OrdinalType,ScalarType>& SerialBandDenseMatrix<OrdinalType, ScalarType>::assign (const SerialBandDenseMatrix<OrdinalType,ScalarType>& Source) {
00744 
00745   if(this == &Source)
00746     return(*this); // Special case of source same as target
00747   if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
00748     return(*this); // Special case of both are views to same data.
00749 
00750   // Check for compatible dimensions
00751   if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_) || (kl_ != Source.kl_) || (ku_ != Source.ku_)) {
00752     TEUCHOS_CHK_REF(*this); // Return *this without altering it.
00753   }
00754   copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0);
00755   return(*this);
00756 
00757 }
00758 
00759 //----------------------------------------------------------------------------------------------------
00760 //   Accessor methods
00761 //----------------------------------------------------------------------------------------------------
00762 
00763 template<typename OrdinalType, typename ScalarType>
00764 inline ScalarType& SerialBandDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex)
00765 {
00766 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
00767   checkIndex( rowIndex, colIndex );
00768 #endif
00769   return(values_[colIndex * stride_ + ku_+rowIndex-colIndex]);
00770 }
00771 
00772 template<typename OrdinalType, typename ScalarType>
00773 inline const ScalarType& SerialBandDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex) const
00774 {
00775 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
00776   checkIndex( rowIndex, colIndex );
00777 #endif
00778   return(values_[colIndex * stride_ + ku_+rowIndex-colIndex]);
00779 }
00780 
00781 template<typename OrdinalType, typename ScalarType>
00782 inline const ScalarType* SerialBandDenseMatrix<OrdinalType, ScalarType>::operator [] (OrdinalType colIndex) const
00783 {
00784 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
00785   checkIndex( 0, colIndex );
00786 #endif
00787   return(values_ + colIndex * stride_);
00788 }
00789 
00790 template<typename OrdinalType, typename ScalarType>
00791 inline ScalarType* SerialBandDenseMatrix<OrdinalType, ScalarType>::operator [] (OrdinalType colIndex)
00792 {
00793 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
00794   checkIndex( 0, colIndex );
00795 #endif
00796   return(values_ + colIndex * stride_);
00797 }
00798 
00799 //----------------------------------------------------------------------------------------------------
00800 //   Norm methods
00801 //----------------------------------------------------------------------------------------------------
00802 
00803 template<typename OrdinalType, typename ScalarType>
00804 typename ScalarTraits<ScalarType>::magnitudeType SerialBandDenseMatrix<OrdinalType, ScalarType>::normOne() const
00805 {
00806   OrdinalType i, j;
00807   typename ScalarTraits<ScalarType>::magnitudeType anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00808   typename ScalarTraits<ScalarType>::magnitudeType absSum = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00809 
00810   ScalarType* ptr;
00811   for(j = 0; j < numCols_; j++) {
00812     ScalarType sum = 0;
00813     ptr = values_ + j * stride_ + TEUCHOS_MAX(0, ku_-j);
00814     for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
00815       sum += ScalarTraits<ScalarType>::magnitude(*ptr++);
00816     }
00817     absSum = ScalarTraits<ScalarType>::magnitude(sum);
00818     if(absSum > anorm) {
00819       anorm = absSum;
00820     }
00821   }
00822   updateFlops((kl_+ku_+1) * numCols_);
00823 
00824   return(anorm);
00825 }
00826 
00827 template<typename OrdinalType, typename ScalarType>
00828 typename ScalarTraits<ScalarType>::magnitudeType SerialBandDenseMatrix<OrdinalType, ScalarType>::normInf() const
00829 {
00830   OrdinalType i, j;
00831   typename ScalarTraits<ScalarType>::magnitudeType sum, anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00832 
00833   for (i = 0; i < numRows_; i++) {
00834     sum = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00835     for (j=TEUCHOS_MAX(0,i-kl_); j<=TEUCHOS_MIN(numCols_-1,i+ku_); j++) {
00836       sum += ScalarTraits<ScalarType>::magnitude(*(values_+(ku_+i-j)+j*stride_));
00837     }
00838     anorm = TEUCHOS_MAX( anorm, sum );
00839   }
00840   updateFlops((kl_+ku_+1) * numCols_);
00841 
00842   return(anorm);
00843 }
00844 
00845 template<typename OrdinalType, typename ScalarType>
00846 typename ScalarTraits<ScalarType>::magnitudeType SerialBandDenseMatrix<OrdinalType, ScalarType>::normFrobenius() const
00847 {
00848   OrdinalType i, j;
00849   typename ScalarTraits<ScalarType>::magnitudeType anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00850 
00851   for (j = 0; j < numCols_; j++) {
00852     for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
00853       anorm += ScalarTraits<ScalarType>::magnitude(values_[(ku_+i-j)+j*stride_]*values_[(ku_+i-j)+j*stride_]);
00854     }
00855   }
00856   anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::squareroot(anorm));
00857   updateFlops((kl_+ku_+1) * numCols_);
00858 
00859   return(anorm);
00860 }
00861 
00862 //----------------------------------------------------------------------------------------------------
00863 //   Comparison methods
00864 //----------------------------------------------------------------------------------------------------
00865 
00866 template<typename OrdinalType, typename ScalarType>
00867 bool SerialBandDenseMatrix<OrdinalType, ScalarType>::operator== (const SerialBandDenseMatrix<OrdinalType, ScalarType> &Operand) const
00868 {
00869   bool result = 1;
00870 
00871   if((numRows_ != Operand.numRows_) || (numCols_ != Operand.numCols_) || (kl_ != Operand.kl_) || (ku_ != Operand.ku_)) {
00872     result = 0;
00873   } else {
00874     OrdinalType i, j;
00875     for(j = 0; j < numCols_; j++) {
00876       for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
00877   if((*this)(i, j) != Operand(i, j)) {
00878     return 0;
00879   }
00880       }
00881     }
00882   }
00883 
00884   return result;
00885 }
00886 
00887 template<typename OrdinalType, typename ScalarType>
00888 bool SerialBandDenseMatrix<OrdinalType, ScalarType>::operator!= (const SerialBandDenseMatrix<OrdinalType, ScalarType> &Operand) const
00889 {
00890   return !((*this) == Operand);
00891 }
00892 
00893 //----------------------------------------------------------------------------------------------------
00894 //   Multiplication method
00895 //----------------------------------------------------------------------------------------------------
00896 
00897 template<typename OrdinalType, typename ScalarType>
00898 SerialBandDenseMatrix<OrdinalType, ScalarType>& SerialBandDenseMatrix<OrdinalType, ScalarType>::operator*= (const ScalarType alpha )
00899 {
00900   this->scale( alpha );
00901   return(*this);
00902 }
00903 
00904 template<typename OrdinalType, typename ScalarType>
00905 int SerialBandDenseMatrix<OrdinalType, ScalarType>::scale( const ScalarType alpha )
00906 {
00907 
00908   OrdinalType i, j;
00909   ScalarType* ptr;
00910 
00911   for (j=0; j<numCols_; j++) {
00912     ptr = values_ + j*stride_ + TEUCHOS_MAX(0, ku_-j);
00913     for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
00914       *ptr = alpha * (*ptr); ptr++;
00915     }
00916   }
00917   updateFlops( (kl_+ku_+1)*numCols_ );
00918 
00919   return(0);
00920 }
00921 
00922 template<typename OrdinalType, typename ScalarType>
00923 int SerialBandDenseMatrix<OrdinalType, ScalarType>::scale( const SerialBandDenseMatrix<OrdinalType,ScalarType>& A )
00924 {
00925 
00926   OrdinalType i, j;
00927   ScalarType* ptr;
00928 
00929   // Check for compatible dimensions
00930   if ((numRows_ != A.numRows_) || (numCols_ != A.numCols_) || (kl_ != A.kl_) || (ku_ != A.ku_)) {
00931     TEUCHOS_CHK_ERR(-1); // Return error
00932   }
00933   for (j=0; j<numCols_; j++) {
00934     ptr = values_ + j*stride_ + TEUCHOS_MAX(0, ku_-j);
00935     for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
00936       *ptr = A(i,j) * (*ptr); ptr++;
00937     }
00938   }
00939   updateFlops( (kl_+ku_+1)*numCols_ );
00940 
00941   return(0);
00942 }
00943 
00944 template<typename OrdinalType, typename ScalarType>
00945 void SerialBandDenseMatrix<OrdinalType, ScalarType>::print(std::ostream& os) const
00946 {
00947   os << std::endl;
00948   if(valuesCopied_)
00949     os << "Values_copied : yes" << std::endl;
00950   else
00951     os << "Values_copied : no" << std::endl;
00952   os << "Rows : " << numRows_ << std::endl;
00953   os << "Columns : " << numCols_ << std::endl;
00954   os << "Lower Bandwidth : " << kl_ << std::endl;
00955   os << "Upper Bandwidth : " << ku_ << std::endl;
00956   os << "LDA : " << stride_ << std::endl;
00957   if(numRows_ == 0 || numCols_ == 0) {
00958     os << "(matrix is empty, no values to display)" << std::endl;
00959   } else {
00960 
00961     for(OrdinalType i = 0; i < numRows_; i++) {
00962       for (OrdinalType j=TEUCHOS_MAX(0,i-kl_); j<=TEUCHOS_MIN(numCols_-1,i+ku_); j++) {
00963   os << (*this)(i,j) << " ";
00964       }
00965       os << std::endl;
00966     }
00967   }
00968 }
00969 
00970 //----------------------------------------------------------------------------------------------------
00971 //   Protected methods
00972 //----------------------------------------------------------------------------------------------------
00973 
00974 template<typename OrdinalType, typename ScalarType>
00975 inline void SerialBandDenseMatrix<OrdinalType, ScalarType>::checkIndex( OrdinalType rowIndex, OrdinalType colIndex ) const {
00976 
00977   TEUCHOS_TEST_FOR_EXCEPTION(rowIndex < 0 || rowIndex >= numRows_ ||
00978            rowIndex < TEUCHOS_MAX(0,colIndex-ku_) || rowIndex > TEUCHOS_MIN(numRows_-1,colIndex+kl_),
00979            std::out_of_range,
00980            "SerialBandDenseMatrix<T>::checkIndex: "
00981            "Row index " << rowIndex << " out of range [0, "<< numRows_ << ")");
00982   TEUCHOS_TEST_FOR_EXCEPTION(colIndex < 0 || colIndex >= numCols_, std::out_of_range,
00983            "SerialBandDenseMatrix<T>::checkIndex: "
00984            "Col index " << colIndex << " out of range [0, "<< numCols_ << ")");
00985 
00986 }
00987 
00988 template<typename OrdinalType, typename ScalarType>
00989 void SerialBandDenseMatrix<OrdinalType, ScalarType>::deleteArrays(void)
00990 {
00991   if (valuesCopied_) {
00992     delete [] values_;
00993     values_ = 0;
00994     valuesCopied_ = false;
00995   }
00996 }
00997 
00998 template<typename OrdinalType, typename ScalarType>
00999 void SerialBandDenseMatrix<OrdinalType, ScalarType>::copyMat(
01000   ScalarType* inputMatrix, OrdinalType strideInput, OrdinalType numRows_in,
01001   OrdinalType numCols_in, ScalarType* outputMatrix, OrdinalType strideOutput, OrdinalType startCol, ScalarType alpha
01002   )
01003 {
01004   OrdinalType i, j;
01005   ScalarType* ptr1 = 0;
01006   ScalarType* ptr2 = 0;
01007 
01008   for(j = 0; j < numCols_in; j++) {
01009     ptr1 = outputMatrix + (j * strideOutput) + TEUCHOS_MAX(0, ku_-j);
01010     ptr2 = inputMatrix + (j + startCol) * strideInput + TEUCHOS_MAX(0, ku_-j);
01011     if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
01012       for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_in-1,j+kl_); i++) {
01013   *ptr1++ += alpha*(*ptr2++);
01014       }
01015     } else {
01016       for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_in-1,j+kl_); i++) {
01017   *ptr1++ = *ptr2++;
01018       }
01019     }
01020   }
01021 }
01022 
01023 } // namespace Teuchos
01024 
01025 
01026 #endif /* _TEUCHOS_SERIALBANDDENSEMATRIX_HPP_ */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines