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 (),
00451     Object ("Teuchos::SerialBandDenseMatrix"),
00452     numRows_ (0),
00453     numCols_ (0),
00454     stride_ (0),
00455     kl_ (0),
00456     ku_ (0),
00457     valuesCopied_ (false),
00458     values_ (0)
00459 {}
00460 
00461 template<typename OrdinalType, typename ScalarType>
00462 SerialBandDenseMatrix<OrdinalType, ScalarType>::
00463 SerialBandDenseMatrix (OrdinalType numRows_in,
00464                        OrdinalType numCols_in,
00465                        OrdinalType kl_in,
00466                        OrdinalType ku_in,
00467                        bool zeroOut)
00468   : CompObject (),
00469     Object ("Teuchos::SerialBandDenseMatrix"),
00470     numRows_ (numRows_in),
00471     numCols_ (numCols_in),
00472     stride_ (kl_in+ku_in+1),
00473     kl_ (kl_in),
00474     ku_ (ku_in),
00475     valuesCopied_ (true),
00476     values_ (NULL) // for safety, in case allocation fails below
00477 {
00478   values_ = new ScalarType[stride_ * numCols_];
00479   if (zeroOut) {
00480     putScalar ();
00481   }
00482 }
00483 
00484 template<typename OrdinalType, typename ScalarType>
00485 SerialBandDenseMatrix<OrdinalType, ScalarType>::
00486 SerialBandDenseMatrix (DataAccess CV,
00487                        ScalarType* values_in,
00488                        OrdinalType stride_in,
00489                        OrdinalType numRows_in,
00490                        OrdinalType numCols_in,
00491                        OrdinalType kl_in,
00492                        OrdinalType ku_in)
00493   : CompObject (),
00494     Object ("Teuchos::SerialBandDenseMatrix"),
00495     numRows_ (numRows_in),
00496     numCols_ (numCols_in),
00497     stride_ (stride_in),
00498     kl_ (kl_in),
00499     ku_ (ku_in),
00500     valuesCopied_ (false),
00501     values_ (values_in)
00502 {
00503   if (CV == Copy) {
00504     stride_ = kl_+ku_+1;
00505     values_ = new ScalarType[stride_*numCols_];
00506     copyMat (values_in, stride_in, numRows_, numCols_, values_, stride_, 0);
00507     valuesCopied_ = true;
00508   }
00509 }
00510 
00511 template<typename OrdinalType, typename ScalarType>
00512 SerialBandDenseMatrix<OrdinalType, ScalarType>::
00513 SerialBandDenseMatrix (const SerialBandDenseMatrix<OrdinalType, ScalarType> &Source, ETransp trans)
00514   : CompObject (),
00515     Object ("Teuchos::SerialBandDenseMatrix"),
00516     numRows_ (0),
00517     numCols_ (0),
00518     stride_ (0),
00519     kl_ (0),
00520     ku_ (0),
00521     valuesCopied_ (true),
00522     values_ (NULL)
00523 {
00524   if (trans == NO_TRANS) {
00525     numRows_ = Source.numRows_;
00526     numCols_ = Source.numCols_;
00527     kl_ = Source.kl_;
00528     ku_ = Source.ku_;
00529     stride_ = kl_+ku_+1;
00530     values_ = new ScalarType[stride_*numCols_];
00531     copyMat (Source.values_, Source.stride_, numRows_, numCols_,
00532              values_, stride_, 0);
00533   }
00534   else if (trans == CONJ_TRANS && ScalarTraits<ScalarType>::isComplex) {
00535     numRows_ = Source.numCols_;
00536     numCols_ = Source.numRows_;
00537     kl_ = Source.ku_;
00538     ku_ = Source.kl_;
00539     stride_ = kl_+ku_+1;
00540     values_ = new ScalarType[stride_*numCols_];
00541     for (OrdinalType j = 0; j < numCols_; ++j) {
00542       for (OrdinalType i = TEUCHOS_MAX(0,j-ku_);
00543            i <= TEUCHOS_MIN(numRows_-1,j+kl_); ++i) {
00544         values_[j*stride_ + (ku_+i-j)] =
00545           ScalarTraits<ScalarType>::conjugate (Source.values_[i*Source.stride_ + (Source.ku_+j-i)]);
00546       }
00547     }
00548   }
00549   else {
00550     numRows_ = Source.numCols_;
00551     numCols_ = Source.numRows_;
00552     kl_ = Source.ku_;
00553     ku_ = Source.kl_;
00554     stride_ = kl_+ku_+1;
00555     values_ = new ScalarType[stride_*numCols_];
00556     for (OrdinalType j=0; j<numCols_; j++) {
00557       for (OrdinalType i = TEUCHOS_MAX(0,j-ku_);
00558            i <= TEUCHOS_MIN(numRows_-1,j+kl_); ++i) {
00559         values_[j*stride_ + (ku_+i-j)] = Source.values_[i*Source.stride_ + (Source.ku_+j-i)];
00560       }
00561     }
00562   }
00563 }
00564 
00565 template<typename OrdinalType, typename ScalarType>
00566 SerialBandDenseMatrix<OrdinalType, ScalarType>::
00567 SerialBandDenseMatrix (DataAccess CV,
00568                        const SerialBandDenseMatrix<OrdinalType, ScalarType> &Source,
00569                        OrdinalType numRows_in,
00570                        OrdinalType numCols_in,
00571                        OrdinalType startCol)
00572   : CompObject (),
00573     Object ("Teuchos::SerialBandDenseMatrix"),
00574     numRows_ (numRows_in),
00575     numCols_ (numCols_in),
00576     stride_ (Source.stride_),
00577     kl_ (Source.kl_),
00578     ku_ (Source.ku_),
00579     valuesCopied_ (false),
00580     values_ (Source.values_)
00581 {
00582   if (CV == Copy) {
00583     values_ = new ScalarType[stride_ * numCols_in];
00584     copyMat (Source.values_, Source.stride_, numRows_in, numCols_in,
00585              values_, stride_, startCol);
00586     valuesCopied_ = true;
00587   } else { // CV = View
00588     values_ = values_ + (stride_ * startCol);
00589   }
00590 }
00591 
00592 template<typename OrdinalType, typename ScalarType>
00593 SerialBandDenseMatrix<OrdinalType, ScalarType>::~SerialBandDenseMatrix()
00594 {
00595   deleteArrays();
00596 }
00597 
00598 //----------------------------------------------------------------------------------------------------
00599 //  Shape methods
00600 //----------------------------------------------------------------------------------------------------
00601 
00602 template<typename OrdinalType, typename ScalarType>
00603 int SerialBandDenseMatrix<OrdinalType, ScalarType>::shape(
00604   OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType kl_in, OrdinalType ku_in
00605   )
00606 {
00607   deleteArrays(); // Get rid of anything that might be already allocated
00608   numRows_ = numRows_in;
00609   numCols_ = numCols_in;
00610   kl_ = kl_in;
00611   ku_ = ku_in;
00612   stride_ = kl_+ku_+1;
00613   values_ = new ScalarType[stride_*numCols_];
00614   putScalar();
00615   valuesCopied_ = true;
00616 
00617   return(0);
00618 }
00619 
00620 template<typename OrdinalType, typename ScalarType>
00621 int SerialBandDenseMatrix<OrdinalType, ScalarType>::shapeUninitialized(
00622   OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType kl_in, OrdinalType ku_in
00623   )
00624 {
00625   deleteArrays(); // Get rid of anything that might be already allocated
00626   numRows_ = numRows_in;
00627   numCols_ = numCols_in;
00628   kl_ = kl_in;
00629   ku_ = ku_in;
00630   stride_ = kl_+ku_+1;
00631   values_ = new ScalarType[stride_*numCols_];
00632   valuesCopied_ = true;
00633 
00634   return(0);
00635 }
00636 
00637 template<typename OrdinalType, typename ScalarType>
00638 int SerialBandDenseMatrix<OrdinalType, ScalarType>::reshape(
00639   OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType kl_in, OrdinalType ku_in
00640   )
00641 {
00642 
00643   // Allocate space for new matrix
00644   ScalarType* values_tmp = new ScalarType[(kl_in+ku_in+1) * numCols_in];
00645   ScalarType zero = ScalarTraits<ScalarType>::zero();
00646   for(OrdinalType k = 0; k < (kl_in+ku_in+1) * numCols_in; k++) {
00647     values_tmp[k] = zero;
00648   }
00649   OrdinalType numRows_tmp = TEUCHOS_MIN(numRows_, numRows_in);
00650   OrdinalType numCols_tmp = TEUCHOS_MIN(numCols_, numCols_in);
00651   if(values_ != 0) {
00652     copyMat(values_, stride_, numRows_tmp, numCols_tmp, values_tmp,
00653             kl_in+ku_in+1, 0); // Copy principal submatrix of A to new A
00654   }
00655   deleteArrays(); // Get rid of anything that might be already allocated
00656   numRows_ = numRows_in;
00657   numCols_ = numCols_in;
00658   kl_ = kl_in;
00659   ku_ = ku_in;
00660   stride_ = kl_+ku_+1;
00661   values_ = values_tmp; // Set pointer to new A
00662   valuesCopied_ = true;
00663 
00664   return(0);
00665 }
00666 
00667 //----------------------------------------------------------------------------------------------------
00668 //   Set methods
00669 //----------------------------------------------------------------------------------------------------
00670 
00671 template<typename OrdinalType, typename ScalarType>
00672 int SerialBandDenseMatrix<OrdinalType, ScalarType>::putScalar( const ScalarType value_in )
00673 {
00674 
00675   // Set each value of the dense matrix to "value".
00676   for(OrdinalType j = 0; j < numCols_; j++) {
00677     for (OrdinalType i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
00678       values_[(ku_+i-j) + j*stride_] = value_in;
00679     }
00680   }
00681 
00682   return 0;
00683 }
00684 
00685 template<typename OrdinalType, typename ScalarType>
00686 int SerialBandDenseMatrix<OrdinalType, ScalarType>::random()
00687 {
00688 
00689   // Set each value of the dense matrix to a random value.
00690   for(OrdinalType j = 0; j < numCols_; j++) {
00691     for (OrdinalType i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
00692       values_[(ku_+i-j) + j*stride_] = ScalarTraits<ScalarType>::random();
00693     }
00694   }
00695 
00696   return 0;
00697 }
00698 
00699 template<typename OrdinalType, typename ScalarType>
00700 SerialBandDenseMatrix<OrdinalType,ScalarType>&
00701 SerialBandDenseMatrix<OrdinalType, ScalarType>::operator=(
00702   const SerialBandDenseMatrix<OrdinalType,ScalarType>& Source
00703   )
00704 {
00705 
00706   if(this == &Source)
00707     return(*this); // Special case of source same as target
00708   if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
00709     return(*this); // Special case of both are views to same data.
00710 
00711   // If the source is a view then we will return a view, else we will return a copy.
00712   if (!Source.valuesCopied_) {
00713     if(valuesCopied_) {
00714       // Clean up stored data if this was previously a copy.
00715       deleteArrays();
00716     }
00717     numRows_ = Source.numRows_;
00718     numCols_ = Source.numCols_;
00719     kl_ = Source.kl_;
00720     ku_ = Source.ku_;
00721     stride_ = Source.stride_;
00722     values_ = Source.values_;
00723   } else {
00724     // If we were a view, we will now be a copy.
00725     if(!valuesCopied_) {
00726       numRows_ = Source.numRows_;
00727       numCols_ = Source.numCols_;
00728       kl_ = Source.kl_;
00729       ku_ = Source.ku_;
00730       stride_ = kl_+ku_+1;
00731       const OrdinalType newsize = stride_ * numCols_;
00732       if(newsize > 0) {
00733         values_ = new ScalarType[newsize];
00734         valuesCopied_ = true;
00735       } else {
00736         values_ = 0;
00737       }
00738     } else {
00739       // If we were a copy, we will stay a copy.
00740       if((Source.numRows_ <= stride_) && (Source.numCols_ == numCols_)) { // we don't need to reallocate
00741         numRows_ = Source.numRows_;
00742         numCols_ = Source.numCols_;
00743         kl_ = Source.kl_;
00744         ku_ = Source.ku_;
00745       } else {
00746         // we need to allocate more space (or less space)
00747         deleteArrays();
00748         numRows_ = Source.numRows_;
00749         numCols_ = Source.numCols_;
00750         kl_ = Source.kl_;
00751         ku_ = Source.ku_;
00752         stride_ = kl_+ku_+1;
00753         const OrdinalType newsize = stride_ * numCols_;
00754         if(newsize > 0) {
00755           values_ = new ScalarType[newsize];
00756           valuesCopied_ = true;
00757         }
00758       }
00759     }
00760     copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0);
00761   }
00762   return(*this);
00763 
00764 }
00765 
00766 template<typename OrdinalType, typename ScalarType>
00767 SerialBandDenseMatrix<OrdinalType, ScalarType>& SerialBandDenseMatrix<OrdinalType, ScalarType>::operator+= (const SerialBandDenseMatrix<OrdinalType,ScalarType>& Source )
00768 {
00769 
00770   // Check for compatible dimensions
00771   if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_) || (kl_ != Source.kl_) || (ku_ != Source.ku_)) {
00772     TEUCHOS_CHK_REF(*this); // Return *this without altering it.
00773   }
00774   copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, ScalarTraits<ScalarType>::one());
00775   return(*this);
00776 
00777 }
00778 
00779 template<typename OrdinalType, typename ScalarType>
00780 SerialBandDenseMatrix<OrdinalType, ScalarType>& SerialBandDenseMatrix<OrdinalType, ScalarType>::operator-= (const SerialBandDenseMatrix<OrdinalType,ScalarType>& Source )
00781 {
00782 
00783   // Check for compatible dimensions
00784   if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_) || (kl_ != Source.kl_) || (ku_ != Source.ku_)) {
00785     TEUCHOS_CHK_REF(*this); // Return *this without altering it.
00786   }
00787   copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, -ScalarTraits<ScalarType>::one());
00788   return(*this);
00789 
00790 }
00791 
00792 template<typename OrdinalType, typename ScalarType>
00793 SerialBandDenseMatrix<OrdinalType,ScalarType>& SerialBandDenseMatrix<OrdinalType, ScalarType>::assign (const SerialBandDenseMatrix<OrdinalType,ScalarType>& Source) {
00794 
00795   if(this == &Source)
00796     return(*this); // Special case of source same as target
00797   if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
00798     return(*this); // Special case of both are views to same data.
00799 
00800   // Check for compatible dimensions
00801   if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_) || (kl_ != Source.kl_) || (ku_ != Source.ku_)) {
00802     TEUCHOS_CHK_REF(*this); // Return *this without altering it.
00803   }
00804   copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0);
00805   return(*this);
00806 
00807 }
00808 
00809 //----------------------------------------------------------------------------------------------------
00810 //   Accessor methods
00811 //----------------------------------------------------------------------------------------------------
00812 
00813 template<typename OrdinalType, typename ScalarType>
00814 inline ScalarType& SerialBandDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex)
00815 {
00816 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
00817   checkIndex( rowIndex, colIndex );
00818 #endif
00819   return(values_[colIndex * stride_ + ku_+rowIndex-colIndex]);
00820 }
00821 
00822 template<typename OrdinalType, typename ScalarType>
00823 inline const ScalarType& SerialBandDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex) const
00824 {
00825 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
00826   checkIndex( rowIndex, colIndex );
00827 #endif
00828   return(values_[colIndex * stride_ + ku_+rowIndex-colIndex]);
00829 }
00830 
00831 template<typename OrdinalType, typename ScalarType>
00832 inline const ScalarType* SerialBandDenseMatrix<OrdinalType, ScalarType>::operator [] (OrdinalType colIndex) const
00833 {
00834 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
00835   checkIndex( 0, colIndex );
00836 #endif
00837   return(values_ + colIndex * stride_);
00838 }
00839 
00840 template<typename OrdinalType, typename ScalarType>
00841 inline ScalarType* SerialBandDenseMatrix<OrdinalType, ScalarType>::operator [] (OrdinalType colIndex)
00842 {
00843 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
00844   checkIndex( 0, colIndex );
00845 #endif
00846   return(values_ + colIndex * stride_);
00847 }
00848 
00849 //----------------------------------------------------------------------------------------------------
00850 //   Norm methods
00851 //----------------------------------------------------------------------------------------------------
00852 
00853 template<typename OrdinalType, typename ScalarType>
00854 typename ScalarTraits<ScalarType>::magnitudeType SerialBandDenseMatrix<OrdinalType, ScalarType>::normOne() const
00855 {
00856   OrdinalType i, j;
00857   typename ScalarTraits<ScalarType>::magnitudeType anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00858   typename ScalarTraits<ScalarType>::magnitudeType absSum = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00859 
00860   ScalarType* ptr;
00861   for(j = 0; j < numCols_; j++) {
00862     ScalarType sum = 0;
00863     ptr = values_ + j * stride_ + TEUCHOS_MAX(0, ku_-j);
00864     for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
00865       sum += ScalarTraits<ScalarType>::magnitude(*ptr++);
00866     }
00867     absSum = ScalarTraits<ScalarType>::magnitude(sum);
00868     if(absSum > anorm) {
00869       anorm = absSum;
00870     }
00871   }
00872   updateFlops((kl_+ku_+1) * numCols_);
00873 
00874   return(anorm);
00875 }
00876 
00877 template<typename OrdinalType, typename ScalarType>
00878 typename ScalarTraits<ScalarType>::magnitudeType SerialBandDenseMatrix<OrdinalType, ScalarType>::normInf() const
00879 {
00880   OrdinalType i, j;
00881   typename ScalarTraits<ScalarType>::magnitudeType sum, anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00882 
00883   for (i = 0; i < numRows_; i++) {
00884     sum = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00885     for (j=TEUCHOS_MAX(0,i-kl_); j<=TEUCHOS_MIN(numCols_-1,i+ku_); j++) {
00886       sum += ScalarTraits<ScalarType>::magnitude(*(values_+(ku_+i-j)+j*stride_));
00887     }
00888     anorm = TEUCHOS_MAX( anorm, sum );
00889   }
00890   updateFlops((kl_+ku_+1) * numCols_);
00891 
00892   return(anorm);
00893 }
00894 
00895 template<typename OrdinalType, typename ScalarType>
00896 typename ScalarTraits<ScalarType>::magnitudeType SerialBandDenseMatrix<OrdinalType, ScalarType>::normFrobenius() const
00897 {
00898   OrdinalType i, j;
00899   typename ScalarTraits<ScalarType>::magnitudeType anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00900 
00901   for (j = 0; j < numCols_; j++) {
00902     for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
00903       anorm += ScalarTraits<ScalarType>::magnitude(values_[(ku_+i-j)+j*stride_]*values_[(ku_+i-j)+j*stride_]);
00904     }
00905   }
00906   anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::squareroot(anorm));
00907   updateFlops((kl_+ku_+1) * numCols_);
00908 
00909   return(anorm);
00910 }
00911 
00912 //----------------------------------------------------------------------------------------------------
00913 //   Comparison methods
00914 //----------------------------------------------------------------------------------------------------
00915 
00916 template<typename OrdinalType, typename ScalarType>
00917 bool SerialBandDenseMatrix<OrdinalType, ScalarType>::operator== (const SerialBandDenseMatrix<OrdinalType, ScalarType> &Operand) const
00918 {
00919   bool result = 1;
00920 
00921   if((numRows_ != Operand.numRows_) || (numCols_ != Operand.numCols_) || (kl_ != Operand.kl_) || (ku_ != Operand.ku_)) {
00922     result = 0;
00923   } else {
00924     OrdinalType i, j;
00925     for(j = 0; j < numCols_; j++) {
00926       for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
00927         if((*this)(i, j) != Operand(i, j)) {
00928           return 0;
00929         }
00930       }
00931     }
00932   }
00933 
00934   return result;
00935 }
00936 
00937 template<typename OrdinalType, typename ScalarType>
00938 bool SerialBandDenseMatrix<OrdinalType, ScalarType>::operator!= (const SerialBandDenseMatrix<OrdinalType, ScalarType> &Operand) const
00939 {
00940   return !((*this) == Operand);
00941 }
00942 
00943 //----------------------------------------------------------------------------------------------------
00944 //   Multiplication method
00945 //----------------------------------------------------------------------------------------------------
00946 
00947 template<typename OrdinalType, typename ScalarType>
00948 SerialBandDenseMatrix<OrdinalType, ScalarType>& SerialBandDenseMatrix<OrdinalType, ScalarType>::operator*= (const ScalarType alpha )
00949 {
00950   this->scale( alpha );
00951   return(*this);
00952 }
00953 
00954 template<typename OrdinalType, typename ScalarType>
00955 int SerialBandDenseMatrix<OrdinalType, ScalarType>::scale( const ScalarType alpha )
00956 {
00957 
00958   OrdinalType i, j;
00959   ScalarType* ptr;
00960 
00961   for (j=0; j<numCols_; j++) {
00962     ptr = values_ + j*stride_ + TEUCHOS_MAX(0, ku_-j);
00963     for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
00964       *ptr = alpha * (*ptr); ptr++;
00965     }
00966   }
00967   updateFlops( (kl_+ku_+1)*numCols_ );
00968 
00969   return(0);
00970 }
00971 
00972 template<typename OrdinalType, typename ScalarType>
00973 int SerialBandDenseMatrix<OrdinalType, ScalarType>::scale( const SerialBandDenseMatrix<OrdinalType,ScalarType>& A )
00974 {
00975 
00976   OrdinalType i, j;
00977   ScalarType* ptr;
00978 
00979   // Check for compatible dimensions
00980   if ((numRows_ != A.numRows_) || (numCols_ != A.numCols_) || (kl_ != A.kl_) || (ku_ != A.ku_)) {
00981     TEUCHOS_CHK_ERR(-1); // Return error
00982   }
00983   for (j=0; j<numCols_; j++) {
00984     ptr = values_ + j*stride_ + TEUCHOS_MAX(0, ku_-j);
00985     for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
00986       *ptr = A(i,j) * (*ptr); ptr++;
00987     }
00988   }
00989   updateFlops( (kl_+ku_+1)*numCols_ );
00990 
00991   return(0);
00992 }
00993 
00994 template<typename OrdinalType, typename ScalarType>
00995 void SerialBandDenseMatrix<OrdinalType, ScalarType>::print(std::ostream& os) const
00996 {
00997   os << std::endl;
00998   if(valuesCopied_)
00999     os << "Values_copied : yes" << std::endl;
01000   else
01001     os << "Values_copied : no" << std::endl;
01002   os << "Rows : " << numRows_ << std::endl;
01003   os << "Columns : " << numCols_ << std::endl;
01004   os << "Lower Bandwidth : " << kl_ << std::endl;
01005   os << "Upper Bandwidth : " << ku_ << std::endl;
01006   os << "LDA : " << stride_ << std::endl;
01007   if(numRows_ == 0 || numCols_ == 0) {
01008     os << "(matrix is empty, no values to display)" << std::endl;
01009   } else {
01010 
01011     for(OrdinalType i = 0; i < numRows_; i++) {
01012       for (OrdinalType j=TEUCHOS_MAX(0,i-kl_); j<=TEUCHOS_MIN(numCols_-1,i+ku_); j++) {
01013         os << (*this)(i,j) << " ";
01014       }
01015       os << std::endl;
01016     }
01017   }
01018 }
01019 
01020 //----------------------------------------------------------------------------------------------------
01021 //   Protected methods
01022 //----------------------------------------------------------------------------------------------------
01023 
01024 template<typename OrdinalType, typename ScalarType>
01025 inline void SerialBandDenseMatrix<OrdinalType, ScalarType>::checkIndex( OrdinalType rowIndex, OrdinalType colIndex ) const {
01026 
01027   TEUCHOS_TEST_FOR_EXCEPTION(rowIndex < 0 || rowIndex >= numRows_ ||
01028                              rowIndex < TEUCHOS_MAX(0,colIndex-ku_) || rowIndex > TEUCHOS_MIN(numRows_-1,colIndex+kl_),
01029                              std::out_of_range,
01030                              "SerialBandDenseMatrix<T>::checkIndex: "
01031                              "Row index " << rowIndex << " out of range [0, "<< numRows_ << ")");
01032   TEUCHOS_TEST_FOR_EXCEPTION(colIndex < 0 || colIndex >= numCols_, std::out_of_range,
01033                              "SerialBandDenseMatrix<T>::checkIndex: "
01034                              "Col index " << colIndex << " out of range [0, "<< numCols_ << ")");
01035 
01036 }
01037 
01038 template<typename OrdinalType, typename ScalarType>
01039 void SerialBandDenseMatrix<OrdinalType, ScalarType>::deleteArrays(void)
01040 {
01041   if (valuesCopied_) {
01042     delete [] values_;
01043     values_ = 0;
01044     valuesCopied_ = false;
01045   }
01046 }
01047 
01048 template<typename OrdinalType, typename ScalarType>
01049 void SerialBandDenseMatrix<OrdinalType, ScalarType>::copyMat(
01050   ScalarType* inputMatrix, OrdinalType strideInput, OrdinalType numRows_in,
01051   OrdinalType numCols_in, ScalarType* outputMatrix, OrdinalType strideOutput, OrdinalType startCol, ScalarType alpha
01052   )
01053 {
01054   OrdinalType i, j;
01055   ScalarType* ptr1 = 0;
01056   ScalarType* ptr2 = 0;
01057 
01058   for(j = 0; j < numCols_in; j++) {
01059     ptr1 = outputMatrix + (j * strideOutput) + TEUCHOS_MAX(0, ku_-j);
01060     ptr2 = inputMatrix + (j + startCol) * strideInput + TEUCHOS_MAX(0, ku_-j);
01061     if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
01062       for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_in-1,j+kl_); i++) {
01063         *ptr1++ += alpha*(*ptr2++);
01064       }
01065     } else {
01066       for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_in-1,j+kl_); i++) {
01067         *ptr1++ = *ptr2++;
01068       }
01069     }
01070   }
01071 }
01072 
01073 } // namespace Teuchos
01074 
01075 
01076 #endif /* _TEUCHOS_SERIALBANDDENSEMATRIX_HPP_ */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines