Teuchos - Trilinos Tools Package  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Teuchos_SerialBandDenseMatrix.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Teuchos: Common Tools Package
4 //
5 // Copyright 2004 NTESS and the Teuchos contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef _TEUCHOS_SERIALBANDDENSEMATRIX_HPP_
11 #define _TEUCHOS_SERIALBANDDENSEMATRIX_HPP_
12 
16 #include "Teuchos_CompObject.hpp"
17 #include "Teuchos_BLAS.hpp"
18 #include "Teuchos_ScalarTraits.hpp"
19 #include "Teuchos_DataAccess.hpp"
20 #include "Teuchos_ConfigDefs.hpp"
21 #include "Teuchos_RCP.hpp"
22 #include "Teuchos_Assert.hpp"
23 
98 namespace Teuchos {
99 
100 template<typename OrdinalType, typename ScalarType>
101 class SerialBandDenseMatrix : public CompObject, public BLAS<OrdinalType, ScalarType>
102 {
103 public:
104 
106  typedef OrdinalType ordinalType;
108  typedef ScalarType scalarType;
109 
111 
112 
114 
118 
120 
130  SerialBandDenseMatrix(OrdinalType numRows, OrdinalType numCols, OrdinalType kl, OrdinalType ku, bool zeroOut = true);
131 
133 
143  SerialBandDenseMatrix(DataAccess CV, ScalarType* values, OrdinalType stride, OrdinalType numRows, OrdinalType numCols, OrdinalType kl, OrdinalType ku);
144 
146 
152 
154 
168  SerialBandDenseMatrix(DataAccess CV, const SerialBandDenseMatrix<OrdinalType, ScalarType> &Source, OrdinalType numRows, OrdinalType numCols, OrdinalType startCol=0);
169 
171  virtual ~SerialBandDenseMatrix();
173 
175 
176 
189  int shape(OrdinalType numRows, OrdinalType numCols, OrdinalType kl, OrdinalType ku);
190 
192  int shapeUninitialized(OrdinalType numRows, OrdinalType numCols, OrdinalType kl, OrdinalType ku);
193 
195 
207  int reshape(OrdinalType numRows, OrdinalType numCols, OrdinalType kl, OrdinalType ku);
208 
209 
211 
213 
214 
216 
223 
225 
231 
233 
236  SerialBandDenseMatrix<OrdinalType, ScalarType>& operator= (const ScalarType value) { putScalar(value); return(*this); }
237 
239 
243  int putScalar( const ScalarType value = Teuchos::ScalarTraits<ScalarType>::zero() );
244 
246  int random();
247 
249 
251 
252 
254 
261  ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex);
262 
264 
271  const ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex) const;
272 
274 
281  ScalarType* operator [] (OrdinalType colIndex);
282 
284 
291  const ScalarType* operator [] (OrdinalType colIndex) const;
292 
294 
295  ScalarType* values() const { return(values_); }
296 
298 
300 
301 
303 
307 
309 
313 
315 
319 
321 
325  int scale ( const ScalarType alpha );
326 
328 
335 
337 
339 
340 
342 
346 
348 
352 
354 
356 
357 
359  OrdinalType numRows() const { return(numRows_); }
360 
362  OrdinalType numCols() const { return(numCols_); }
363 
365  OrdinalType lowerBandwidth() const { return(kl_); }
366 
368  OrdinalType upperBandwidth() const { return(ku_); }
369 
371  OrdinalType stride() const { return(stride_); }
372 
374  bool empty() const { return(numRows_ == 0 || numCols_ == 0); }
376 
378 
379 
382 
385 
389 
391 
392 
394  virtual std::ostream& print(std::ostream& os) const;
395 
396 
398 protected:
399  void copyMat(ScalarType* inputMatrix, OrdinalType strideInput,
400  OrdinalType numRows, OrdinalType numCols, ScalarType* outputMatrix,
401  OrdinalType strideOutput, OrdinalType startCol, ScalarType alpha = ScalarTraits<ScalarType>::zero() );
402  void deleteArrays();
403  void checkIndex( OrdinalType rowIndex, OrdinalType colIndex = 0 ) const;
404  OrdinalType numRows_;
405  OrdinalType numCols_;
406  OrdinalType stride_;
407  OrdinalType kl_;
408  OrdinalType ku_;
409  bool valuesCopied_;
410  ScalarType* values_;
411 
412 }; // class Teuchos_SerialBandDenseMatrix
413 
414 //----------------------------------------------------------------------------------------------------
415 // Constructors and Destructor
416 //----------------------------------------------------------------------------------------------------
417 
418 template<typename OrdinalType, typename ScalarType>
420  : CompObject (),
421  BLAS<OrdinalType,ScalarType>(),
422  numRows_ (0),
423  numCols_ (0),
424  stride_ (0),
425  kl_ (0),
426  ku_ (0),
427  valuesCopied_ (false),
428  values_ (0)
429 {}
430 
431 template<typename OrdinalType, typename ScalarType>
433 SerialBandDenseMatrix (OrdinalType numRows_in,
434  OrdinalType numCols_in,
435  OrdinalType kl_in,
436  OrdinalType ku_in,
437  bool zeroOut)
438  : CompObject (),
439  BLAS<OrdinalType,ScalarType>(),
440  numRows_ (numRows_in),
441  numCols_ (numCols_in),
442  stride_ (kl_in+ku_in+1),
443  kl_ (kl_in),
444  ku_ (ku_in),
445  valuesCopied_ (true),
446  values_ (NULL) // for safety, in case allocation fails below
447 {
448  values_ = new ScalarType[stride_ * numCols_];
449  if (zeroOut) {
450  putScalar ();
451  }
452 }
453 
454 template<typename OrdinalType, typename ScalarType>
457  ScalarType* values_in,
458  OrdinalType stride_in,
459  OrdinalType numRows_in,
460  OrdinalType numCols_in,
461  OrdinalType kl_in,
462  OrdinalType ku_in)
463  : CompObject (),
464  BLAS<OrdinalType,ScalarType>(),
465  numRows_ (numRows_in),
466  numCols_ (numCols_in),
467  stride_ (stride_in),
468  kl_ (kl_in),
469  ku_ (ku_in),
470  valuesCopied_ (false),
471  values_ (values_in)
472 {
473  if (CV == Copy) {
474  stride_ = kl_+ku_+1;
475  values_ = new ScalarType[stride_*numCols_];
476  copyMat (values_in, stride_in, numRows_, numCols_, values_, stride_, 0);
477  valuesCopied_ = true;
478  }
479 }
480 
481 template<typename OrdinalType, typename ScalarType>
484  : CompObject (),
485  BLAS<OrdinalType,ScalarType>(),
486  numRows_ (0),
487  numCols_ (0),
488  stride_ (0),
489  kl_ (0),
490  ku_ (0),
491  valuesCopied_ (true),
492  values_ (NULL)
493 {
494  if (trans == NO_TRANS) {
495  numRows_ = Source.numRows_;
496  numCols_ = Source.numCols_;
497  kl_ = Source.kl_;
498  ku_ = Source.ku_;
499  stride_ = kl_+ku_+1;
500  values_ = new ScalarType[stride_*numCols_];
501  copyMat (Source.values_, Source.stride_, numRows_, numCols_,
502  values_, stride_, 0);
503  }
504  else if (trans == CONJ_TRANS && ScalarTraits<ScalarType>::isComplex) {
505  numRows_ = Source.numCols_;
506  numCols_ = Source.numRows_;
507  kl_ = Source.ku_;
508  ku_ = Source.kl_;
509  stride_ = kl_+ku_+1;
510  values_ = new ScalarType[stride_*numCols_];
511  for (OrdinalType j = 0; j < numCols_; ++j) {
512  for (OrdinalType i = TEUCHOS_MAX(0,j-ku_);
513  i <= TEUCHOS_MIN(numRows_-1,j+kl_); ++i) {
514  values_[j*stride_ + (ku_+i-j)] =
515  ScalarTraits<ScalarType>::conjugate (Source.values_[i*Source.stride_ + (Source.ku_+j-i)]);
516  }
517  }
518  }
519  else {
520  numRows_ = Source.numCols_;
521  numCols_ = Source.numRows_;
522  kl_ = Source.ku_;
523  ku_ = Source.kl_;
524  stride_ = kl_+ku_+1;
525  values_ = new ScalarType[stride_*numCols_];
526  for (OrdinalType j=0; j<numCols_; j++) {
527  for (OrdinalType i = TEUCHOS_MAX(0,j-ku_);
528  i <= TEUCHOS_MIN(numRows_-1,j+kl_); ++i) {
529  values_[j*stride_ + (ku_+i-j)] = Source.values_[i*Source.stride_ + (Source.ku_+j-i)];
530  }
531  }
532  }
533 }
534 
535 template<typename OrdinalType, typename ScalarType>
539  OrdinalType numRows_in,
540  OrdinalType numCols_in,
541  OrdinalType startCol)
542  : CompObject (),
543  BLAS<OrdinalType,ScalarType>(),
544  numRows_ (numRows_in),
545  numCols_ (numCols_in),
546  stride_ (Source.stride_),
547  kl_ (Source.kl_),
548  ku_ (Source.ku_),
549  valuesCopied_ (false),
550  values_ (Source.values_)
551 {
552  if (CV == Copy) {
553  values_ = new ScalarType[stride_ * numCols_in];
554  copyMat (Source.values_, Source.stride_, numRows_in, numCols_in,
555  values_, stride_, startCol);
556  valuesCopied_ = true;
557  } else { // CV = View
558  values_ = values_ + (stride_ * startCol);
559  }
560 }
561 
562 template<typename OrdinalType, typename ScalarType>
564 {
565  deleteArrays();
566 }
567 
568 //----------------------------------------------------------------------------------------------------
569 // Shape methods
570 //----------------------------------------------------------------------------------------------------
571 
572 template<typename OrdinalType, typename ScalarType>
574  OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType kl_in, OrdinalType ku_in
575  )
576 {
577  deleteArrays(); // Get rid of anything that might be already allocated
578  numRows_ = numRows_in;
579  numCols_ = numCols_in;
580  kl_ = kl_in;
581  ku_ = ku_in;
582  stride_ = kl_+ku_+1;
583  values_ = new ScalarType[stride_*numCols_];
584  putScalar();
585  valuesCopied_ = true;
586 
587  return(0);
588 }
589 
590 template<typename OrdinalType, typename ScalarType>
592  OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType kl_in, OrdinalType ku_in
593  )
594 {
595  deleteArrays(); // Get rid of anything that might be already allocated
596  numRows_ = numRows_in;
597  numCols_ = numCols_in;
598  kl_ = kl_in;
599  ku_ = ku_in;
600  stride_ = kl_+ku_+1;
601  values_ = new ScalarType[stride_*numCols_];
602  valuesCopied_ = true;
603 
604  return(0);
605 }
606 
607 template<typename OrdinalType, typename ScalarType>
609  OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType kl_in, OrdinalType ku_in
610  )
611 {
612 
613  // Allocate space for new matrix
614  ScalarType* values_tmp = new ScalarType[(kl_in+ku_in+1) * numCols_in];
615  ScalarType zero = ScalarTraits<ScalarType>::zero();
616  for(OrdinalType k = 0; k < (kl_in+ku_in+1) * numCols_in; k++) {
617  values_tmp[k] = zero;
618  }
619  OrdinalType numRows_tmp = TEUCHOS_MIN(numRows_, numRows_in);
620  OrdinalType numCols_tmp = TEUCHOS_MIN(numCols_, numCols_in);
621  if(values_ != 0) {
622  copyMat(values_, stride_, numRows_tmp, numCols_tmp, values_tmp,
623  kl_in+ku_in+1, 0); // Copy principal submatrix of A to new A
624  }
625  deleteArrays(); // Get rid of anything that might be already allocated
626  numRows_ = numRows_in;
627  numCols_ = numCols_in;
628  kl_ = kl_in;
629  ku_ = ku_in;
630  stride_ = kl_+ku_+1;
631  values_ = values_tmp; // Set pointer to new A
632  valuesCopied_ = true;
633 
634  return(0);
635 }
636 
637 //----------------------------------------------------------------------------------------------------
638 // Set methods
639 //----------------------------------------------------------------------------------------------------
640 
641 template<typename OrdinalType, typename ScalarType>
643 {
644 
645  // Set each value of the dense matrix to "value".
646  for(OrdinalType j = 0; j < numCols_; j++) {
647  for (OrdinalType i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
648  values_[(ku_+i-j) + j*stride_] = value_in;
649  }
650  }
651 
652  return 0;
653 }
654 
655 template<typename OrdinalType, typename ScalarType>
657 {
658 
659  // Set each value of the dense matrix to a random value.
660  for(OrdinalType j = 0; j < numCols_; j++) {
661  for (OrdinalType i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
662  values_[(ku_+i-j) + j*stride_] = ScalarTraits<ScalarType>::random();
663  }
664  }
665 
666  return 0;
667 }
668 
669 template<typename OrdinalType, typename ScalarType>
673  )
674 {
675 
676  if(this == &Source)
677  return(*this); // Special case of source same as target
678  if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
679  return(*this); // Special case of both are views to same data.
680 
681  // If the source is a view then we will return a view, else we will return a copy.
682  if (!Source.valuesCopied_) {
683  if(valuesCopied_) {
684  // Clean up stored data if this was previously a copy.
685  deleteArrays();
686  }
687  numRows_ = Source.numRows_;
688  numCols_ = Source.numCols_;
689  kl_ = Source.kl_;
690  ku_ = Source.ku_;
691  stride_ = Source.stride_;
692  values_ = Source.values_;
693  } else {
694  // If we were a view, we will now be a copy.
695  if(!valuesCopied_) {
696  numRows_ = Source.numRows_;
697  numCols_ = Source.numCols_;
698  kl_ = Source.kl_;
699  ku_ = Source.ku_;
700  stride_ = kl_+ku_+1;
701  const OrdinalType newsize = stride_ * numCols_;
702  if(newsize > 0) {
703  values_ = new ScalarType[newsize];
704  valuesCopied_ = true;
705  } else {
706  values_ = 0;
707  }
708  } else {
709  // If we were a copy, we will stay a copy.
710  if((Source.numRows_ <= stride_) && (Source.numCols_ == numCols_)) { // we don't need to reallocate
711  numRows_ = Source.numRows_;
712  numCols_ = Source.numCols_;
713  kl_ = Source.kl_;
714  ku_ = Source.ku_;
715  } else {
716  // we need to allocate more space (or less space)
717  deleteArrays();
718  numRows_ = Source.numRows_;
719  numCols_ = Source.numCols_;
720  kl_ = Source.kl_;
721  ku_ = Source.ku_;
722  stride_ = kl_+ku_+1;
723  const OrdinalType newsize = stride_ * numCols_;
724  if(newsize > 0) {
725  values_ = new ScalarType[newsize];
726  valuesCopied_ = true;
727  }
728  }
729  }
730  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0);
731  }
732  return(*this);
733 
734 }
735 
736 template<typename OrdinalType, typename ScalarType>
738 {
739 
740  // Check for compatible dimensions
741  if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_) || (kl_ != Source.kl_) || (ku_ != Source.ku_)) {
742  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
743  }
744  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, ScalarTraits<ScalarType>::one());
745  return(*this);
746 
747 }
748 
749 template<typename OrdinalType, typename ScalarType>
751 {
752 
753  // Check for compatible dimensions
754  if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_) || (kl_ != Source.kl_) || (ku_ != Source.ku_)) {
755  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
756  }
757  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, -ScalarTraits<ScalarType>::one());
758  return(*this);
759 
760 }
761 
762 template<typename OrdinalType, typename ScalarType>
764 
765  if(this == &Source)
766  return(*this); // Special case of source same as target
767  if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
768  return(*this); // Special case of both are views to same data.
769 
770  // Check for compatible dimensions
771  if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_) || (kl_ != Source.kl_) || (ku_ != Source.ku_)) {
772  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
773  }
774  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0);
775  return(*this);
776 
777 }
778 
779 //----------------------------------------------------------------------------------------------------
780 // Accessor methods
781 //----------------------------------------------------------------------------------------------------
782 
783 template<typename OrdinalType, typename ScalarType>
784 inline ScalarType& SerialBandDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex)
785 {
786 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
787  checkIndex( rowIndex, colIndex );
788 #endif
789  return(values_[colIndex * stride_ + ku_+rowIndex-colIndex]);
790 }
791 
792 template<typename OrdinalType, typename ScalarType>
793 inline const ScalarType& SerialBandDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex) const
794 {
795 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
796  checkIndex( rowIndex, colIndex );
797 #endif
798  return(values_[colIndex * stride_ + ku_+rowIndex-colIndex]);
799 }
800 
801 template<typename OrdinalType, typename ScalarType>
802 inline const ScalarType* SerialBandDenseMatrix<OrdinalType, ScalarType>::operator [] (OrdinalType colIndex) const
803 {
804 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
805  checkIndex( 0, colIndex );
806 #endif
807  return(values_ + colIndex * stride_);
808 }
809 
810 template<typename OrdinalType, typename ScalarType>
811 inline ScalarType* SerialBandDenseMatrix<OrdinalType, ScalarType>::operator [] (OrdinalType colIndex)
812 {
813 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
814  checkIndex( 0, colIndex );
815 #endif
816  return(values_ + colIndex * stride_);
817 }
818 
819 //----------------------------------------------------------------------------------------------------
820 // Norm methods
821 //----------------------------------------------------------------------------------------------------
822 
823 template<typename OrdinalType, typename ScalarType>
825 {
826  OrdinalType i, j;
829 
830  ScalarType* ptr;
831  for(j = 0; j < numCols_; j++) {
832  ScalarType sum = 0;
833  ptr = values_ + j * stride_ + TEUCHOS_MAX(0, ku_-j);
834  for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
836  }
838  if(absSum > anorm) {
839  anorm = absSum;
840  }
841  }
842  updateFlops((kl_+ku_+1) * numCols_);
843 
844  return(anorm);
845 }
846 
847 template<typename OrdinalType, typename ScalarType>
849 {
850  OrdinalType i, j;
852 
853  for (i = 0; i < numRows_; i++) {
855  for (j=TEUCHOS_MAX(0,i-kl_); j<=TEUCHOS_MIN(numCols_-1,i+ku_); j++) {
856  sum += ScalarTraits<ScalarType>::magnitude(*(values_+(ku_+i-j)+j*stride_));
857  }
858  anorm = TEUCHOS_MAX( anorm, sum );
859  }
860  updateFlops((kl_+ku_+1) * numCols_);
861 
862  return(anorm);
863 }
864 
865 template<typename OrdinalType, typename ScalarType>
867 {
868  OrdinalType i, j;
870 
871  for (j = 0; j < numCols_; j++) {
872  for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
873  anorm += ScalarTraits<ScalarType>::magnitude(values_[(ku_+i-j)+j*stride_]*values_[(ku_+i-j)+j*stride_]);
874  }
875  }
877  updateFlops((kl_+ku_+1) * numCols_);
878 
879  return(anorm);
880 }
881 
882 //----------------------------------------------------------------------------------------------------
883 // Comparison methods
884 //----------------------------------------------------------------------------------------------------
885 
886 template<typename OrdinalType, typename ScalarType>
888 {
889  bool result = 1;
890 
891  if((numRows_ != Operand.numRows_) || (numCols_ != Operand.numCols_) || (kl_ != Operand.kl_) || (ku_ != Operand.ku_)) {
892  result = 0;
893  } else {
894  OrdinalType i, j;
895  for(j = 0; j < numCols_; j++) {
896  for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
897  if((*this)(i, j) != Operand(i, j)) {
898  return 0;
899  }
900  }
901  }
902  }
903 
904  return result;
905 }
906 
907 template<typename OrdinalType, typename ScalarType>
909 {
910  return !((*this) == Operand);
911 }
912 
913 //----------------------------------------------------------------------------------------------------
914 // Multiplication method
915 //----------------------------------------------------------------------------------------------------
916 
917 template<typename OrdinalType, typename ScalarType>
919 {
920  this->scale( alpha );
921  return(*this);
922 }
923 
924 template<typename OrdinalType, typename ScalarType>
926 {
927 
928  OrdinalType i, j;
929  ScalarType* ptr;
930 
931  for (j=0; j<numCols_; j++) {
932  ptr = values_ + j*stride_ + TEUCHOS_MAX(0, ku_-j);
933  for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
934  *ptr = alpha * (*ptr); ptr++;
935  }
936  }
937  updateFlops( (kl_+ku_+1)*numCols_ );
938 
939  return(0);
940 }
941 
942 template<typename OrdinalType, typename ScalarType>
944 {
945 
946  OrdinalType i, j;
947  ScalarType* ptr;
948 
949  // Check for compatible dimensions
950  if ((numRows_ != A.numRows_) || (numCols_ != A.numCols_) || (kl_ != A.kl_) || (ku_ != A.ku_)) {
951  TEUCHOS_CHK_ERR(-1); // Return error
952  }
953  for (j=0; j<numCols_; j++) {
954  ptr = values_ + j*stride_ + TEUCHOS_MAX(0, ku_-j);
955  for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
956  *ptr = A(i,j) * (*ptr); ptr++;
957  }
958  }
959  updateFlops( (kl_+ku_+1)*numCols_ );
960 
961  return(0);
962 }
963 
964 template<typename OrdinalType, typename ScalarType>
965 std::ostream& SerialBandDenseMatrix<OrdinalType, ScalarType>::print(std::ostream& os) const
966 {
967  os << std::endl;
968  if(valuesCopied_)
969  os << "Values_copied : yes" << std::endl;
970  else
971  os << "Values_copied : no" << std::endl;
972  os << "Rows : " << numRows_ << std::endl;
973  os << "Columns : " << numCols_ << std::endl;
974  os << "Lower Bandwidth : " << kl_ << std::endl;
975  os << "Upper Bandwidth : " << ku_ << std::endl;
976  os << "LDA : " << stride_ << std::endl;
977  if(numRows_ == 0 || numCols_ == 0) {
978  os << "(matrix is empty, no values to display)" << std::endl;
979  } else {
980 
981  for(OrdinalType i = 0; i < numRows_; i++) {
982  for (OrdinalType j=TEUCHOS_MAX(0,i-kl_); j<=TEUCHOS_MIN(numCols_-1,i+ku_); j++) {
983  os << (*this)(i,j) << " ";
984  }
985  os << std::endl;
986  }
987  }
988  return os;
989 }
990 
991 //----------------------------------------------------------------------------------------------------
992 // Protected methods
993 //----------------------------------------------------------------------------------------------------
994 
995 template<typename OrdinalType, typename ScalarType>
996 inline void SerialBandDenseMatrix<OrdinalType, ScalarType>::checkIndex( OrdinalType rowIndex, OrdinalType colIndex ) const {
997 
998  TEUCHOS_TEST_FOR_EXCEPTION(rowIndex < 0 || rowIndex >= numRows_ ||
999  rowIndex < TEUCHOS_MAX(0,colIndex-ku_) || rowIndex > TEUCHOS_MIN(numRows_-1,colIndex+kl_),
1000  std::out_of_range,
1001  "SerialBandDenseMatrix<T>::checkIndex: "
1002  "Row index " << rowIndex << " out of range [0, "<< numRows_ << ")");
1003  TEUCHOS_TEST_FOR_EXCEPTION(colIndex < 0 || colIndex >= numCols_, std::out_of_range,
1004  "SerialBandDenseMatrix<T>::checkIndex: "
1005  "Col index " << colIndex << " out of range [0, "<< numCols_ << ")");
1006 
1007 }
1008 
1009 template<typename OrdinalType, typename ScalarType>
1010 void SerialBandDenseMatrix<OrdinalType, ScalarType>::deleteArrays(void)
1011 {
1012  if (valuesCopied_) {
1013  delete [] values_;
1014  values_ = 0;
1015  valuesCopied_ = false;
1016  }
1017 }
1018 
1019 template<typename OrdinalType, typename ScalarType>
1020 void SerialBandDenseMatrix<OrdinalType, ScalarType>::copyMat(
1021  ScalarType* inputMatrix, OrdinalType strideInput, OrdinalType numRows_in,
1022  OrdinalType numCols_in, ScalarType* outputMatrix, OrdinalType strideOutput, OrdinalType startCol, ScalarType alpha
1023  )
1024 {
1025  OrdinalType i, j;
1026  ScalarType* ptr1 = 0;
1027  ScalarType* ptr2 = 0;
1028 
1029  for(j = 0; j < numCols_in; j++) {
1030  ptr1 = outputMatrix + (j * strideOutput) + TEUCHOS_MAX(0, ku_-j);
1031  ptr2 = inputMatrix + (j + startCol) * strideInput + TEUCHOS_MAX(0, ku_-j);
1032  if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
1033  for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_in-1,j+kl_); i++) {
1034  *ptr1++ += alpha*(*ptr2++);
1035  }
1036  } else {
1037  for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_in-1,j+kl_); i++) {
1038  *ptr1++ = *ptr2++;
1039  }
1040  }
1041  }
1042 }
1043 
1045 template<typename OrdinalType, typename ScalarType>
1047 public:
1051  : obj(obj_in) {}
1052 };
1053 
1055 template<typename OrdinalType, typename ScalarType>
1056 std::ostream&
1057 operator<<(std::ostream &out,
1059 {
1060  printer.obj.print(out);
1061  return out;
1062 }
1063 
1065 template<typename OrdinalType, typename ScalarType>
1066 SerialBandDenseMatrixPrinter<OrdinalType,ScalarType>
1068 {
1070 }
1071 
1072 
1073 } // namespace Teuchos
1074 
1075 
1076 #endif /* _TEUCHOS_SERIALBANDDENSEMATRIX_HPP_ */
bool empty() const
Returns whether this matrix is empty.
OrdinalType numRows() const
Returns the row dimension of this matrix.
SerialBandDenseMatrix< OrdinalType, ScalarType > & operator+=(const SerialBandDenseMatrix< OrdinalType, ScalarType > &Source)
Add another matrix to this matrix.
ScalarType & operator()(OrdinalType rowIndex, OrdinalType colIndex)
Element access method (non-const).
int scale(const ScalarType alpha)
Scale this matrix by alpha; *this = alpha**this.
int shape(OrdinalType numRows, OrdinalType numCols, OrdinalType kl, OrdinalType ku)
Shape method for changing the size of a SerialBandDenseMatrix, initializing entries to zero...
ScalarTraits< ScalarType >::magnitudeType normInf() const
Returns the Infinity-norm of the matrix.
Ostream manipulator for SerialBandDenseMatrix.
Templated interface class to BLAS routines.
SerialBandDenseMatrixPrinter< OrdinalType, ScalarType > printMat(const SerialBandDenseMatrix< OrdinalType, ScalarType > &obj)
Return SerialBandDenseMatrix ostream manipulator Use as:
#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.
Teuchos header file which uses auto-configuration information to include necessary C++ headers...
ScalarTraits< ScalarType >::magnitudeType normOne() const
Returns the 1-norm of the matrix.
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...
virtual std::ostream & print(std::ostream &os) const
Print method. Defines the behavior of the std::ostream &lt;&lt; operator.
This structure defines some basic traits for a scalar field type.
OrdinalType stride() const
Returns the stride between the columns of this matrix in memory.
ScalarType * values() const
Data array access method.
Functionality and data that is common to all computational classes.
bool operator!=(const SerialBandDenseMatrix< OrdinalType, ScalarType > &Operand) const
Inequality of two matrices.
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
Returns the Frobenius-norm of the matrix.
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.
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 operator==(const SerialBandDenseMatrix< OrdinalType, ScalarType > &Operand) const
Equality of two matrices.
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.
OrdinalType upperBandwidth() const
Returns the upper bandwidth of this matrix.
Defines basic traits for the scalar field type.
static T zero()
Returns representation of zero for this scalar type.
OrdinalType lowerBandwidth() const
Returns the lower bandwidth of this matrix.
ScalarType * operator[](OrdinalType colIndex)
Column access method (non-const).
ScalarType scalarType
Typedef for scalar type.
OrdinalType numCols() const
Returns the column dimension of this matrix.
Reference-counted pointer class and non-member templated function implementations.
int shapeUninitialized(OrdinalType numRows, OrdinalType numCols, OrdinalType kl, OrdinalType ku)
Same as shape() except leaves uninitialized.
Teuchos::DataAccess Mode enumerable type.