Teuchos Package Browser (Single Doxygen Collection)  Version of the Day
Teuchos_SerialBandDenseMatrix.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Teuchos: Common Tools Package
5 // Copyright (2004) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #ifndef _TEUCHOS_SERIALBANDDENSEMATRIX_HPP_
43 #define _TEUCHOS_SERIALBANDDENSEMATRIX_HPP_
44 
48 #include "Teuchos_CompObject.hpp"
49 #include "Teuchos_BLAS.hpp"
50 #include "Teuchos_ScalarTraits.hpp"
51 #include "Teuchos_DataAccess.hpp"
52 #include "Teuchos_ConfigDefs.hpp"
53 #include "Teuchos_RCP.hpp"
54 #include "Teuchos_Assert.hpp"
55 
130 namespace Teuchos {
131 
132 template<typename OrdinalType, typename ScalarType>
133 class SerialBandDenseMatrix : public CompObject, public Object, public BLAS<OrdinalType, ScalarType>
134 {
135 public:
136 
138  typedef OrdinalType ordinalType;
140  typedef ScalarType scalarType;
141 
143 
144 
146 
150 
152 
162  SerialBandDenseMatrix(OrdinalType numRows, OrdinalType numCols, OrdinalType kl, OrdinalType ku, bool zeroOut = true);
163 
165 
175  SerialBandDenseMatrix(DataAccess CV, ScalarType* values, OrdinalType stride, OrdinalType numRows, OrdinalType numCols, OrdinalType kl, OrdinalType ku);
176 
178 
184 
186 
200  SerialBandDenseMatrix(DataAccess CV, const SerialBandDenseMatrix<OrdinalType, ScalarType> &Source, OrdinalType numRows, OrdinalType numCols, OrdinalType startCol=0);
201 
203  virtual ~SerialBandDenseMatrix();
205 
207 
208 
221  int shape(OrdinalType numRows, OrdinalType numCols, OrdinalType kl, OrdinalType ku);
222 
224  int shapeUninitialized(OrdinalType numRows, OrdinalType numCols, OrdinalType kl, OrdinalType ku);
225 
227 
239  int reshape(OrdinalType numRows, OrdinalType numCols, OrdinalType kl, OrdinalType ku);
240 
241 
243 
245 
246 
248 
255 
257 
263 
265 
268  SerialBandDenseMatrix<OrdinalType, ScalarType>& operator= (const ScalarType value) { putScalar(value); return(*this); }
269 
271 
275  int putScalar( const ScalarType value = Teuchos::ScalarTraits<ScalarType>::zero() );
276 
278  int random();
279 
281 
283 
284 
286 
293  ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex);
294 
296 
303  const ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex) const;
304 
306 
313  ScalarType* operator [] (OrdinalType colIndex);
314 
316 
323  const ScalarType* operator [] (OrdinalType colIndex) const;
324 
326 
327  ScalarType* values() const { return(values_); }
328 
330 
332 
333 
335 
339 
341 
345 
347 
351 
353 
357  int scale ( const ScalarType alpha );
358 
360 
367 
369 
371 
372 
374 
378 
380 
384 
386 
388 
389 
391  OrdinalType numRows() const { return(numRows_); }
392 
394  OrdinalType numCols() const { return(numCols_); }
395 
397  OrdinalType lowerBandwidth() const { return(kl_); }
398 
400  OrdinalType upperBandwidth() const { return(ku_); }
401 
403  OrdinalType stride() const { return(stride_); }
404 
406  bool empty() const { return(numRows_ == 0 || numCols_ == 0); }
408 
410 
411 
414 
417 
421 
423 
424  virtual void print(std::ostream& os) const;
426 
428 protected:
429  void copyMat(ScalarType* inputMatrix, OrdinalType strideInput,
430  OrdinalType numRows, OrdinalType numCols, ScalarType* outputMatrix,
431  OrdinalType strideOutput, OrdinalType startCol, ScalarType alpha = ScalarTraits<ScalarType>::zero() );
432  void deleteArrays();
433  void checkIndex( OrdinalType rowIndex, OrdinalType colIndex = 0 ) const;
434  OrdinalType numRows_;
435  OrdinalType numCols_;
436  OrdinalType stride_;
437  OrdinalType kl_;
438  OrdinalType ku_;
440  ScalarType* values_;
441 
442 }; // class Teuchos_SerialBandDenseMatrix
443 
444 //----------------------------------------------------------------------------------------------------
445 // Constructors and Destructor
446 //----------------------------------------------------------------------------------------------------
447 
448 template<typename OrdinalType, typename ScalarType>
450  : CompObject (),
451  Object(),
452  BLAS<OrdinalType,ScalarType>(),
453  numRows_ (0),
454  numCols_ (0),
455  stride_ (0),
456  kl_ (0),
457  ku_ (0),
458  valuesCopied_ (false),
459  values_ (0)
460 {}
461 
462 template<typename OrdinalType, typename ScalarType>
464 SerialBandDenseMatrix (OrdinalType numRows_in,
465  OrdinalType numCols_in,
466  OrdinalType kl_in,
467  OrdinalType ku_in,
468  bool zeroOut)
469  : CompObject (),
470  Object(),
471  BLAS<OrdinalType,ScalarType>(),
472  numRows_ (numRows_in),
473  numCols_ (numCols_in),
474  stride_ (kl_in+ku_in+1),
475  kl_ (kl_in),
476  ku_ (ku_in),
477  valuesCopied_ (true),
478  values_ (NULL) // for safety, in case allocation fails below
479 {
480  values_ = new ScalarType[stride_ * numCols_];
481  if (zeroOut) {
482  putScalar ();
483  }
484 }
485 
486 template<typename OrdinalType, typename ScalarType>
489  ScalarType* values_in,
490  OrdinalType stride_in,
491  OrdinalType numRows_in,
492  OrdinalType numCols_in,
493  OrdinalType kl_in,
494  OrdinalType ku_in)
495  : CompObject (),
496  Object(),
497  BLAS<OrdinalType,ScalarType>(),
498  numRows_ (numRows_in),
499  numCols_ (numCols_in),
500  stride_ (stride_in),
501  kl_ (kl_in),
502  ku_ (ku_in),
503  valuesCopied_ (false),
504  values_ (values_in)
505 {
506  if (CV == Copy) {
507  stride_ = kl_+ku_+1;
508  values_ = new ScalarType[stride_*numCols_];
509  copyMat (values_in, stride_in, numRows_, numCols_, values_, stride_, 0);
510  valuesCopied_ = true;
511  }
512 }
513 
514 template<typename OrdinalType, typename ScalarType>
517  : CompObject (),
518  Object(),
519  BLAS<OrdinalType,ScalarType>(),
520  numRows_ (0),
521  numCols_ (0),
522  stride_ (0),
523  kl_ (0),
524  ku_ (0),
525  valuesCopied_ (true),
526  values_ (NULL)
527 {
528  if (trans == NO_TRANS) {
529  numRows_ = Source.numRows_;
530  numCols_ = Source.numCols_;
531  kl_ = Source.kl_;
532  ku_ = Source.ku_;
533  stride_ = kl_+ku_+1;
534  values_ = new ScalarType[stride_*numCols_];
535  copyMat (Source.values_, Source.stride_, numRows_, numCols_,
536  values_, stride_, 0);
537  }
538  else if (trans == CONJ_TRANS && ScalarTraits<ScalarType>::isComplex) {
539  numRows_ = Source.numCols_;
540  numCols_ = Source.numRows_;
541  kl_ = Source.ku_;
542  ku_ = Source.kl_;
543  stride_ = kl_+ku_+1;
544  values_ = new ScalarType[stride_*numCols_];
545  for (OrdinalType j = 0; j < numCols_; ++j) {
546  for (OrdinalType i = TEUCHOS_MAX(0,j-ku_);
547  i <= TEUCHOS_MIN(numRows_-1,j+kl_); ++i) {
548  values_[j*stride_ + (ku_+i-j)] =
549  ScalarTraits<ScalarType>::conjugate (Source.values_[i*Source.stride_ + (Source.ku_+j-i)]);
550  }
551  }
552  }
553  else {
554  numRows_ = Source.numCols_;
555  numCols_ = Source.numRows_;
556  kl_ = Source.ku_;
557  ku_ = Source.kl_;
558  stride_ = kl_+ku_+1;
559  values_ = new ScalarType[stride_*numCols_];
560  for (OrdinalType j=0; j<numCols_; j++) {
561  for (OrdinalType i = TEUCHOS_MAX(0,j-ku_);
562  i <= TEUCHOS_MIN(numRows_-1,j+kl_); ++i) {
563  values_[j*stride_ + (ku_+i-j)] = Source.values_[i*Source.stride_ + (Source.ku_+j-i)];
564  }
565  }
566  }
567 }
568 
569 template<typename OrdinalType, typename ScalarType>
573  OrdinalType numRows_in,
574  OrdinalType numCols_in,
575  OrdinalType startCol)
576  : CompObject (),
577  Object(),
578  BLAS<OrdinalType,ScalarType>(),
579  numRows_ (numRows_in),
580  numCols_ (numCols_in),
581  stride_ (Source.stride_),
582  kl_ (Source.kl_),
583  ku_ (Source.ku_),
584  valuesCopied_ (false),
585  values_ (Source.values_)
586 {
587  if (CV == Copy) {
588  values_ = new ScalarType[stride_ * numCols_in];
589  copyMat (Source.values_, Source.stride_, numRows_in, numCols_in,
590  values_, stride_, startCol);
591  valuesCopied_ = true;
592  } else { // CV = View
593  values_ = values_ + (stride_ * startCol);
594  }
595 }
596 
597 template<typename OrdinalType, typename ScalarType>
599 {
600  deleteArrays();
601 }
602 
603 //----------------------------------------------------------------------------------------------------
604 // Shape methods
605 //----------------------------------------------------------------------------------------------------
606 
607 template<typename OrdinalType, typename ScalarType>
609  OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType kl_in, OrdinalType ku_in
610  )
611 {
612  deleteArrays(); // Get rid of anything that might be already allocated
613  numRows_ = numRows_in;
614  numCols_ = numCols_in;
615  kl_ = kl_in;
616  ku_ = ku_in;
617  stride_ = kl_+ku_+1;
618  values_ = new ScalarType[stride_*numCols_];
619  putScalar();
620  valuesCopied_ = true;
621 
622  return(0);
623 }
624 
625 template<typename OrdinalType, typename ScalarType>
627  OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType kl_in, OrdinalType ku_in
628  )
629 {
630  deleteArrays(); // Get rid of anything that might be already allocated
631  numRows_ = numRows_in;
632  numCols_ = numCols_in;
633  kl_ = kl_in;
634  ku_ = ku_in;
635  stride_ = kl_+ku_+1;
636  values_ = new ScalarType[stride_*numCols_];
637  valuesCopied_ = true;
638 
639  return(0);
640 }
641 
642 template<typename OrdinalType, typename ScalarType>
644  OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType kl_in, OrdinalType ku_in
645  )
646 {
647 
648  // Allocate space for new matrix
649  ScalarType* values_tmp = new ScalarType[(kl_in+ku_in+1) * numCols_in];
650  ScalarType zero = ScalarTraits<ScalarType>::zero();
651  for(OrdinalType k = 0; k < (kl_in+ku_in+1) * numCols_in; k++) {
652  values_tmp[k] = zero;
653  }
654  OrdinalType numRows_tmp = TEUCHOS_MIN(numRows_, numRows_in);
655  OrdinalType numCols_tmp = TEUCHOS_MIN(numCols_, numCols_in);
656  if(values_ != 0) {
657  copyMat(values_, stride_, numRows_tmp, numCols_tmp, values_tmp,
658  kl_in+ku_in+1, 0); // Copy principal submatrix of A to new A
659  }
660  deleteArrays(); // Get rid of anything that might be already allocated
661  numRows_ = numRows_in;
662  numCols_ = numCols_in;
663  kl_ = kl_in;
664  ku_ = ku_in;
665  stride_ = kl_+ku_+1;
666  values_ = values_tmp; // Set pointer to new A
667  valuesCopied_ = true;
668 
669  return(0);
670 }
671 
672 //----------------------------------------------------------------------------------------------------
673 // Set methods
674 //----------------------------------------------------------------------------------------------------
675 
676 template<typename OrdinalType, typename ScalarType>
678 {
679 
680  // Set each value of the dense matrix to "value".
681  for(OrdinalType j = 0; j < numCols_; j++) {
682  for (OrdinalType i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
683  values_[(ku_+i-j) + j*stride_] = value_in;
684  }
685  }
686 
687  return 0;
688 }
689 
690 template<typename OrdinalType, typename ScalarType>
692 {
693 
694  // Set each value of the dense matrix to a random value.
695  for(OrdinalType j = 0; j < numCols_; j++) {
696  for (OrdinalType i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
697  values_[(ku_+i-j) + j*stride_] = ScalarTraits<ScalarType>::random();
698  }
699  }
700 
701  return 0;
702 }
703 
704 template<typename OrdinalType, typename ScalarType>
708  )
709 {
710 
711  if(this == &Source)
712  return(*this); // Special case of source same as target
713  if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
714  return(*this); // Special case of both are views to same data.
715 
716  // If the source is a view then we will return a view, else we will return a copy.
717  if (!Source.valuesCopied_) {
718  if(valuesCopied_) {
719  // Clean up stored data if this was previously a copy.
720  deleteArrays();
721  }
722  numRows_ = Source.numRows_;
723  numCols_ = Source.numCols_;
724  kl_ = Source.kl_;
725  ku_ = Source.ku_;
726  stride_ = Source.stride_;
727  values_ = Source.values_;
728  } else {
729  // If we were a view, we will now be a copy.
730  if(!valuesCopied_) {
731  numRows_ = Source.numRows_;
732  numCols_ = Source.numCols_;
733  kl_ = Source.kl_;
734  ku_ = Source.ku_;
735  stride_ = kl_+ku_+1;
736  const OrdinalType newsize = stride_ * numCols_;
737  if(newsize > 0) {
738  values_ = new ScalarType[newsize];
739  valuesCopied_ = true;
740  } else {
741  values_ = 0;
742  }
743  } else {
744  // If we were a copy, we will stay a copy.
745  if((Source.numRows_ <= stride_) && (Source.numCols_ == numCols_)) { // we don't need to reallocate
746  numRows_ = Source.numRows_;
747  numCols_ = Source.numCols_;
748  kl_ = Source.kl_;
749  ku_ = Source.ku_;
750  } else {
751  // we need to allocate more space (or less space)
752  deleteArrays();
753  numRows_ = Source.numRows_;
754  numCols_ = Source.numCols_;
755  kl_ = Source.kl_;
756  ku_ = Source.ku_;
757  stride_ = kl_+ku_+1;
758  const OrdinalType newsize = stride_ * numCols_;
759  if(newsize > 0) {
760  values_ = new ScalarType[newsize];
761  valuesCopied_ = true;
762  }
763  }
764  }
765  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0);
766  }
767  return(*this);
768 
769 }
770 
771 template<typename OrdinalType, typename ScalarType>
773 {
774 
775  // Check for compatible dimensions
776  if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_) || (kl_ != Source.kl_) || (ku_ != Source.ku_)) {
777  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
778  }
779  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, ScalarTraits<ScalarType>::one());
780  return(*this);
781 
782 }
783 
784 template<typename OrdinalType, typename ScalarType>
786 {
787 
788  // Check for compatible dimensions
789  if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_) || (kl_ != Source.kl_) || (ku_ != Source.ku_)) {
790  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
791  }
792  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, -ScalarTraits<ScalarType>::one());
793  return(*this);
794 
795 }
796 
797 template<typename OrdinalType, typename ScalarType>
799 
800  if(this == &Source)
801  return(*this); // Special case of source same as target
802  if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
803  return(*this); // Special case of both are views to same data.
804 
805  // Check for compatible dimensions
806  if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_) || (kl_ != Source.kl_) || (ku_ != Source.ku_)) {
807  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
808  }
809  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0);
810  return(*this);
811 
812 }
813 
814 //----------------------------------------------------------------------------------------------------
815 // Accessor methods
816 //----------------------------------------------------------------------------------------------------
817 
818 template<typename OrdinalType, typename ScalarType>
819 inline ScalarType& SerialBandDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex)
820 {
821 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
822  checkIndex( rowIndex, colIndex );
823 #endif
824  return(values_[colIndex * stride_ + ku_+rowIndex-colIndex]);
825 }
826 
827 template<typename OrdinalType, typename ScalarType>
828 inline const ScalarType& SerialBandDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex) const
829 {
830 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
831  checkIndex( rowIndex, colIndex );
832 #endif
833  return(values_[colIndex * stride_ + ku_+rowIndex-colIndex]);
834 }
835 
836 template<typename OrdinalType, typename ScalarType>
837 inline const ScalarType* SerialBandDenseMatrix<OrdinalType, ScalarType>::operator [] (OrdinalType colIndex) const
838 {
839 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
840  checkIndex( 0, colIndex );
841 #endif
842  return(values_ + colIndex * stride_);
843 }
844 
845 template<typename OrdinalType, typename ScalarType>
846 inline ScalarType* SerialBandDenseMatrix<OrdinalType, ScalarType>::operator [] (OrdinalType colIndex)
847 {
848 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
849  checkIndex( 0, colIndex );
850 #endif
851  return(values_ + colIndex * stride_);
852 }
853 
854 //----------------------------------------------------------------------------------------------------
855 // Norm methods
856 //----------------------------------------------------------------------------------------------------
857 
858 template<typename OrdinalType, typename ScalarType>
860 {
861  OrdinalType i, j;
864 
865  ScalarType* ptr;
866  for(j = 0; j < numCols_; j++) {
867  ScalarType sum = 0;
868  ptr = values_ + j * stride_ + TEUCHOS_MAX(0, ku_-j);
869  for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
871  }
873  if(absSum > anorm) {
874  anorm = absSum;
875  }
876  }
877  updateFlops((kl_+ku_+1) * numCols_);
878 
879  return(anorm);
880 }
881 
882 template<typename OrdinalType, typename ScalarType>
884 {
885  OrdinalType i, j;
887 
888  for (i = 0; i < numRows_; i++) {
890  for (j=TEUCHOS_MAX(0,i-kl_); j<=TEUCHOS_MIN(numCols_-1,i+ku_); j++) {
891  sum += ScalarTraits<ScalarType>::magnitude(*(values_+(ku_+i-j)+j*stride_));
892  }
893  anorm = TEUCHOS_MAX( anorm, sum );
894  }
895  updateFlops((kl_+ku_+1) * numCols_);
896 
897  return(anorm);
898 }
899 
900 template<typename OrdinalType, typename ScalarType>
902 {
903  OrdinalType i, j;
905 
906  for (j = 0; j < numCols_; j++) {
907  for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
908  anorm += ScalarTraits<ScalarType>::magnitude(values_[(ku_+i-j)+j*stride_]*values_[(ku_+i-j)+j*stride_]);
909  }
910  }
912  updateFlops((kl_+ku_+1) * numCols_);
913 
914  return(anorm);
915 }
916 
917 //----------------------------------------------------------------------------------------------------
918 // Comparison methods
919 //----------------------------------------------------------------------------------------------------
920 
921 template<typename OrdinalType, typename ScalarType>
923 {
924  bool result = 1;
925 
926  if((numRows_ != Operand.numRows_) || (numCols_ != Operand.numCols_) || (kl_ != Operand.kl_) || (ku_ != Operand.ku_)) {
927  result = 0;
928  } else {
929  OrdinalType i, j;
930  for(j = 0; j < numCols_; j++) {
931  for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
932  if((*this)(i, j) != Operand(i, j)) {
933  return 0;
934  }
935  }
936  }
937  }
938 
939  return result;
940 }
941 
942 template<typename OrdinalType, typename ScalarType>
944 {
945  return !((*this) == Operand);
946 }
947 
948 //----------------------------------------------------------------------------------------------------
949 // Multiplication method
950 //----------------------------------------------------------------------------------------------------
951 
952 template<typename OrdinalType, typename ScalarType>
954 {
955  this->scale( alpha );
956  return(*this);
957 }
958 
959 template<typename OrdinalType, typename ScalarType>
961 {
962 
963  OrdinalType i, j;
964  ScalarType* ptr;
965 
966  for (j=0; j<numCols_; j++) {
967  ptr = values_ + j*stride_ + TEUCHOS_MAX(0, ku_-j);
968  for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
969  *ptr = alpha * (*ptr); ptr++;
970  }
971  }
972  updateFlops( (kl_+ku_+1)*numCols_ );
973 
974  return(0);
975 }
976 
977 template<typename OrdinalType, typename ScalarType>
979 {
980 
981  OrdinalType i, j;
982  ScalarType* ptr;
983 
984  // Check for compatible dimensions
985  if ((numRows_ != A.numRows_) || (numCols_ != A.numCols_) || (kl_ != A.kl_) || (ku_ != A.ku_)) {
986  TEUCHOS_CHK_ERR(-1); // Return error
987  }
988  for (j=0; j<numCols_; j++) {
989  ptr = values_ + j*stride_ + TEUCHOS_MAX(0, ku_-j);
990  for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
991  *ptr = A(i,j) * (*ptr); ptr++;
992  }
993  }
994  updateFlops( (kl_+ku_+1)*numCols_ );
995 
996  return(0);
997 }
998 
999 template<typename OrdinalType, typename ScalarType>
1001 {
1002  os << std::endl;
1003  if(valuesCopied_)
1004  os << "Values_copied : yes" << std::endl;
1005  else
1006  os << "Values_copied : no" << std::endl;
1007  os << "Rows : " << numRows_ << std::endl;
1008  os << "Columns : " << numCols_ << std::endl;
1009  os << "Lower Bandwidth : " << kl_ << std::endl;
1010  os << "Upper Bandwidth : " << ku_ << std::endl;
1011  os << "LDA : " << stride_ << std::endl;
1012  if(numRows_ == 0 || numCols_ == 0) {
1013  os << "(matrix is empty, no values to display)" << std::endl;
1014  } else {
1015 
1016  for(OrdinalType i = 0; i < numRows_; i++) {
1017  for (OrdinalType j=TEUCHOS_MAX(0,i-kl_); j<=TEUCHOS_MIN(numCols_-1,i+ku_); j++) {
1018  os << (*this)(i,j) << " ";
1019  }
1020  os << std::endl;
1021  }
1022  }
1023 }
1024 
1025 //----------------------------------------------------------------------------------------------------
1026 // Protected methods
1027 //----------------------------------------------------------------------------------------------------
1028 
1029 template<typename OrdinalType, typename ScalarType>
1030 inline void SerialBandDenseMatrix<OrdinalType, ScalarType>::checkIndex( OrdinalType rowIndex, OrdinalType colIndex ) const {
1031 
1032  TEUCHOS_TEST_FOR_EXCEPTION(rowIndex < 0 || rowIndex >= numRows_ ||
1033  rowIndex < TEUCHOS_MAX(0,colIndex-ku_) || rowIndex > TEUCHOS_MIN(numRows_-1,colIndex+kl_),
1034  std::out_of_range,
1035  "SerialBandDenseMatrix<T>::checkIndex: "
1036  "Row index " << rowIndex << " out of range [0, "<< numRows_ << ")");
1037  TEUCHOS_TEST_FOR_EXCEPTION(colIndex < 0 || colIndex >= numCols_, std::out_of_range,
1038  "SerialBandDenseMatrix<T>::checkIndex: "
1039  "Col index " << colIndex << " out of range [0, "<< numCols_ << ")");
1040 
1041 }
1042 
1043 template<typename OrdinalType, typename ScalarType>
1045 {
1046  if (valuesCopied_) {
1047  delete [] values_;
1048  values_ = 0;
1049  valuesCopied_ = false;
1050  }
1051 }
1052 
1053 template<typename OrdinalType, typename ScalarType>
1055  ScalarType* inputMatrix, OrdinalType strideInput, OrdinalType numRows_in,
1056  OrdinalType numCols_in, ScalarType* outputMatrix, OrdinalType strideOutput, OrdinalType startCol, ScalarType alpha
1057  )
1058 {
1059  OrdinalType i, j;
1060  ScalarType* ptr1 = 0;
1061  ScalarType* ptr2 = 0;
1062 
1063  for(j = 0; j < numCols_in; j++) {
1064  ptr1 = outputMatrix + (j * strideOutput) + TEUCHOS_MAX(0, ku_-j);
1065  ptr2 = inputMatrix + (j + startCol) * strideInput + TEUCHOS_MAX(0, ku_-j);
1066  if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
1067  for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_in-1,j+kl_); i++) {
1068  *ptr1++ += alpha*(*ptr2++);
1069  }
1070  } else {
1071  for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_in-1,j+kl_); i++) {
1072  *ptr1++ = *ptr2++;
1073  }
1074  }
1075  }
1076 }
1077 
1078 } // namespace Teuchos
1079 
1080 
1081 #endif /* _TEUCHOS_SERIALBANDDENSEMATRIX_HPP_ */
ScalarTraits< ScalarType >::magnitudeType normInf() const
Returns the Infinity-norm of the matrix.
bool operator==(const SerialBandDenseMatrix< OrdinalType, ScalarType > &Operand) const
Equality of two matrices.
SerialBandDenseMatrix< OrdinalType, ScalarType > & operator+=(const SerialBandDenseMatrix< OrdinalType, ScalarType > &Source)
Add another matrix to this matrix.
bool operator!=(const SerialBandDenseMatrix< OrdinalType, ScalarType > &Operand) const
Inequality of two matrices.
ScalarType & operator()(OrdinalType rowIndex, OrdinalType colIndex)
Element access method (non-const).
int scale(const ScalarType alpha)
Scale this matrix by alpha; *this = alpha**this.
void checkIndex(OrdinalType rowIndex, OrdinalType colIndex=0) const
int shape(OrdinalType numRows, OrdinalType numCols, OrdinalType kl, OrdinalType ku)
Shape method for changing the size of a SerialBandDenseMatrix, initializing entries to zero...
Templated interface class to BLAS routines.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Macro for throwing an exception with breakpointing to ease debugging.
int reshape(OrdinalType numRows, OrdinalType numCols, OrdinalType kl, OrdinalType ku)
Reshaping method for changing the size of a SerialBandDenseMatrix, keeping the entries.
OrdinalType upperBandwidth() const
Returns the upper bandwidth of this matrix.
Teuchos header file which uses auto-configuration information to include necessary C++ headers...
Templated BLAS wrapper.
SerialBandDenseMatrix< OrdinalType, ScalarType > & assign(const SerialBandDenseMatrix< OrdinalType, ScalarType > &Source)
Copies values from one matrix to another.
Object for storing data and providing functionality that is common to all computational classes...
Definition: PackageA.cpp:3
This structure defines some basic traits for a scalar field type.
#define TEUCHOS_CHK_ERR(a)
Functionality and data that is common to all computational classes.
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
Returns the Frobenius-norm of the matrix.
ScalarTraits< ScalarType >::magnitudeType normOne() const
Returns the 1-norm of the matrix.
The base Teuchos class.
This class creates and provides basic support for banded dense matrices of templated type...
SerialBandDenseMatrix< OrdinalType, ScalarType > & operator*=(const ScalarType alpha)
Scale this matrix by alpha; *this = alpha**this.
ScalarType * values() const
Data array access method.
#define TEUCHOS_MAX(x, y)
SerialBandDenseMatrix< OrdinalType, ScalarType > & operator-=(const SerialBandDenseMatrix< OrdinalType, ScalarType > &Source)
Subtract another matrix from this matrix.
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero())
Set all values in the matrix to a constant value.
int random()
Set all values in the matrix to be random numbers.
bool empty() const
Returns whether this matrix is empty.
OrdinalType numRows() const
Returns the row dimension of this matrix.
static magnitudeType magnitude(T a)
Returns the magnitudeType of the scalar type a.
OrdinalType ordinalType
Typedef for ordinal type.
SerialBandDenseMatrix< OrdinalType, ScalarType > & operator=(const SerialBandDenseMatrix< OrdinalType, ScalarType > &Source)
Copies values from one matrix to another.
void copyMat(ScalarType *inputMatrix, OrdinalType strideInput, OrdinalType numRows, OrdinalType numCols, ScalarType *outputMatrix, OrdinalType strideOutput, OrdinalType startCol, ScalarType alpha=ScalarTraits< ScalarType >::zero())
Defines basic traits for the scalar field type.
static T zero()
Returns representation of zero for this scalar type.
ScalarType * operator[](OrdinalType colIndex)
Column access method (non-const).
ScalarType scalarType
Typedef for scalar type.
#define TEUCHOS_CHK_REF(a)
OrdinalType numCols() const
Returns the column dimension of this matrix.
OrdinalType stride() const
Returns the stride between the columns of this matrix in memory.
#define TEUCHOS_MIN(x, y)
Reference-counted pointer class and non-member templated function implementations.
OrdinalType lowerBandwidth() const
Returns the lower bandwidth of this matrix.
virtual void print(std::ostream &os) const
Print method. Defines the behavior of the std::ostream << operator inherited from the Object class...
int shapeUninitialized(OrdinalType numRows, OrdinalType numCols, OrdinalType kl, OrdinalType ku)
Same as shape() except leaves uninitialized.
Teuchos::DataAccess Mode enumerable type.