Teuchos - Trilinos Tools Package  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_;
439  bool valuesCopied_;
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  numRows_ (0),
452  numCols_ (0),
453  stride_ (0),
454  kl_ (0),
455  ku_ (0),
456  valuesCopied_ (false),
457  values_ (0)
458 {}
459 
460 template<typename OrdinalType, typename ScalarType>
462 SerialBandDenseMatrix (OrdinalType numRows_in,
463  OrdinalType numCols_in,
464  OrdinalType kl_in,
465  OrdinalType ku_in,
466  bool zeroOut)
467  : CompObject (),
468  numRows_ (numRows_in),
469  numCols_ (numCols_in),
470  stride_ (kl_in+ku_in+1),
471  kl_ (kl_in),
472  ku_ (ku_in),
473  valuesCopied_ (true),
474  values_ (NULL) // for safety, in case allocation fails below
475 {
476  values_ = new ScalarType[stride_ * numCols_];
477  if (zeroOut) {
478  putScalar ();
479  }
480 }
481 
482 template<typename OrdinalType, typename ScalarType>
485  ScalarType* values_in,
486  OrdinalType stride_in,
487  OrdinalType numRows_in,
488  OrdinalType numCols_in,
489  OrdinalType kl_in,
490  OrdinalType ku_in)
491  : CompObject (),
492  numRows_ (numRows_in),
493  numCols_ (numCols_in),
494  stride_ (stride_in),
495  kl_ (kl_in),
496  ku_ (ku_in),
497  valuesCopied_ (false),
498  values_ (values_in)
499 {
500  if (CV == Copy) {
501  stride_ = kl_+ku_+1;
502  values_ = new ScalarType[stride_*numCols_];
503  copyMat (values_in, stride_in, numRows_, numCols_, values_, stride_, 0);
504  valuesCopied_ = true;
505  }
506 }
507 
508 template<typename OrdinalType, typename ScalarType>
511  : CompObject (),
512  numRows_ (0),
513  numCols_ (0),
514  stride_ (0),
515  kl_ (0),
516  ku_ (0),
517  valuesCopied_ (true),
518  values_ (NULL)
519 {
520  if (trans == NO_TRANS) {
521  numRows_ = Source.numRows_;
522  numCols_ = Source.numCols_;
523  kl_ = Source.kl_;
524  ku_ = Source.ku_;
525  stride_ = kl_+ku_+1;
526  values_ = new ScalarType[stride_*numCols_];
527  copyMat (Source.values_, Source.stride_, numRows_, numCols_,
528  values_, stride_, 0);
529  }
530  else if (trans == CONJ_TRANS && ScalarTraits<ScalarType>::isComplex) {
531  numRows_ = Source.numCols_;
532  numCols_ = Source.numRows_;
533  kl_ = Source.ku_;
534  ku_ = Source.kl_;
535  stride_ = kl_+ku_+1;
536  values_ = new ScalarType[stride_*numCols_];
537  for (OrdinalType j = 0; j < numCols_; ++j) {
538  for (OrdinalType i = TEUCHOS_MAX(0,j-ku_);
539  i <= TEUCHOS_MIN(numRows_-1,j+kl_); ++i) {
540  values_[j*stride_ + (ku_+i-j)] =
541  ScalarTraits<ScalarType>::conjugate (Source.values_[i*Source.stride_ + (Source.ku_+j-i)]);
542  }
543  }
544  }
545  else {
546  numRows_ = Source.numCols_;
547  numCols_ = Source.numRows_;
548  kl_ = Source.ku_;
549  ku_ = Source.kl_;
550  stride_ = kl_+ku_+1;
551  values_ = new ScalarType[stride_*numCols_];
552  for (OrdinalType j=0; j<numCols_; j++) {
553  for (OrdinalType i = TEUCHOS_MAX(0,j-ku_);
554  i <= TEUCHOS_MIN(numRows_-1,j+kl_); ++i) {
555  values_[j*stride_ + (ku_+i-j)] = Source.values_[i*Source.stride_ + (Source.ku_+j-i)];
556  }
557  }
558  }
559 }
560 
561 template<typename OrdinalType, typename ScalarType>
565  OrdinalType numRows_in,
566  OrdinalType numCols_in,
567  OrdinalType startCol)
568  : CompObject (),
569  numRows_ (numRows_in),
570  numCols_ (numCols_in),
571  stride_ (Source.stride_),
572  kl_ (Source.kl_),
573  ku_ (Source.ku_),
574  valuesCopied_ (false),
575  values_ (Source.values_)
576 {
577  if (CV == Copy) {
578  values_ = new ScalarType[stride_ * numCols_in];
579  copyMat (Source.values_, Source.stride_, numRows_in, numCols_in,
580  values_, stride_, startCol);
581  valuesCopied_ = true;
582  } else { // CV = View
583  values_ = values_ + (stride_ * startCol);
584  }
585 }
586 
587 template<typename OrdinalType, typename ScalarType>
589 {
590  deleteArrays();
591 }
592 
593 //----------------------------------------------------------------------------------------------------
594 // Shape methods
595 //----------------------------------------------------------------------------------------------------
596 
597 template<typename OrdinalType, typename ScalarType>
599  OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType kl_in, OrdinalType ku_in
600  )
601 {
602  deleteArrays(); // Get rid of anything that might be already allocated
603  numRows_ = numRows_in;
604  numCols_ = numCols_in;
605  kl_ = kl_in;
606  ku_ = ku_in;
607  stride_ = kl_+ku_+1;
608  values_ = new ScalarType[stride_*numCols_];
609  putScalar();
610  valuesCopied_ = true;
611 
612  return(0);
613 }
614 
615 template<typename OrdinalType, typename ScalarType>
617  OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType kl_in, OrdinalType ku_in
618  )
619 {
620  deleteArrays(); // Get rid of anything that might be already allocated
621  numRows_ = numRows_in;
622  numCols_ = numCols_in;
623  kl_ = kl_in;
624  ku_ = ku_in;
625  stride_ = kl_+ku_+1;
626  values_ = new ScalarType[stride_*numCols_];
627  valuesCopied_ = true;
628 
629  return(0);
630 }
631 
632 template<typename OrdinalType, typename ScalarType>
634  OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType kl_in, OrdinalType ku_in
635  )
636 {
637 
638  // Allocate space for new matrix
639  ScalarType* values_tmp = new ScalarType[(kl_in+ku_in+1) * numCols_in];
640  ScalarType zero = ScalarTraits<ScalarType>::zero();
641  for(OrdinalType k = 0; k < (kl_in+ku_in+1) * numCols_in; k++) {
642  values_tmp[k] = zero;
643  }
644  OrdinalType numRows_tmp = TEUCHOS_MIN(numRows_, numRows_in);
645  OrdinalType numCols_tmp = TEUCHOS_MIN(numCols_, numCols_in);
646  if(values_ != 0) {
647  copyMat(values_, stride_, numRows_tmp, numCols_tmp, values_tmp,
648  kl_in+ku_in+1, 0); // Copy principal submatrix of A to new A
649  }
650  deleteArrays(); // Get rid of anything that might be already allocated
651  numRows_ = numRows_in;
652  numCols_ = numCols_in;
653  kl_ = kl_in;
654  ku_ = ku_in;
655  stride_ = kl_+ku_+1;
656  values_ = values_tmp; // Set pointer to new A
657  valuesCopied_ = true;
658 
659  return(0);
660 }
661 
662 //----------------------------------------------------------------------------------------------------
663 // Set methods
664 //----------------------------------------------------------------------------------------------------
665 
666 template<typename OrdinalType, typename ScalarType>
668 {
669 
670  // Set each value of the dense matrix to "value".
671  for(OrdinalType j = 0; j < numCols_; j++) {
672  for (OrdinalType i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
673  values_[(ku_+i-j) + j*stride_] = value_in;
674  }
675  }
676 
677  return 0;
678 }
679 
680 template<typename OrdinalType, typename ScalarType>
682 {
683 
684  // Set each value of the dense matrix to a random value.
685  for(OrdinalType j = 0; j < numCols_; j++) {
686  for (OrdinalType i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
687  values_[(ku_+i-j) + j*stride_] = ScalarTraits<ScalarType>::random();
688  }
689  }
690 
691  return 0;
692 }
693 
694 template<typename OrdinalType, typename ScalarType>
698  )
699 {
700 
701  if(this == &Source)
702  return(*this); // Special case of source same as target
703  if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
704  return(*this); // Special case of both are views to same data.
705 
706  // If the source is a view then we will return a view, else we will return a copy.
707  if (!Source.valuesCopied_) {
708  if(valuesCopied_) {
709  // Clean up stored data if this was previously a copy.
710  deleteArrays();
711  }
712  numRows_ = Source.numRows_;
713  numCols_ = Source.numCols_;
714  kl_ = Source.kl_;
715  ku_ = Source.ku_;
716  stride_ = Source.stride_;
717  values_ = Source.values_;
718  } else {
719  // If we were a view, we will now be a copy.
720  if(!valuesCopied_) {
721  numRows_ = Source.numRows_;
722  numCols_ = Source.numCols_;
723  kl_ = Source.kl_;
724  ku_ = Source.ku_;
725  stride_ = kl_+ku_+1;
726  const OrdinalType newsize = stride_ * numCols_;
727  if(newsize > 0) {
728  values_ = new ScalarType[newsize];
729  valuesCopied_ = true;
730  } else {
731  values_ = 0;
732  }
733  } else {
734  // If we were a copy, we will stay a copy.
735  if((Source.numRows_ <= stride_) && (Source.numCols_ == numCols_)) { // we don't need to reallocate
736  numRows_ = Source.numRows_;
737  numCols_ = Source.numCols_;
738  kl_ = Source.kl_;
739  ku_ = Source.ku_;
740  } else {
741  // we need to allocate more space (or less space)
742  deleteArrays();
743  numRows_ = Source.numRows_;
744  numCols_ = Source.numCols_;
745  kl_ = Source.kl_;
746  ku_ = Source.ku_;
747  stride_ = kl_+ku_+1;
748  const OrdinalType newsize = stride_ * numCols_;
749  if(newsize > 0) {
750  values_ = new ScalarType[newsize];
751  valuesCopied_ = true;
752  }
753  }
754  }
755  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0);
756  }
757  return(*this);
758 
759 }
760 
761 template<typename OrdinalType, typename ScalarType>
763 {
764 
765  // Check for compatible dimensions
766  if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_) || (kl_ != Source.kl_) || (ku_ != Source.ku_)) {
767  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
768  }
769  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, ScalarTraits<ScalarType>::one());
770  return(*this);
771 
772 }
773 
774 template<typename OrdinalType, typename ScalarType>
776 {
777 
778  // Check for compatible dimensions
779  if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_) || (kl_ != Source.kl_) || (ku_ != Source.ku_)) {
780  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
781  }
782  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, -ScalarTraits<ScalarType>::one());
783  return(*this);
784 
785 }
786 
787 template<typename OrdinalType, typename ScalarType>
789 
790  if(this == &Source)
791  return(*this); // Special case of source same as target
792  if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
793  return(*this); // Special case of both are views to same data.
794 
795  // Check for compatible dimensions
796  if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_) || (kl_ != Source.kl_) || (ku_ != Source.ku_)) {
797  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
798  }
799  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0);
800  return(*this);
801 
802 }
803 
804 //----------------------------------------------------------------------------------------------------
805 // Accessor methods
806 //----------------------------------------------------------------------------------------------------
807 
808 template<typename OrdinalType, typename ScalarType>
809 inline ScalarType& SerialBandDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex)
810 {
811 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
812  checkIndex( rowIndex, colIndex );
813 #endif
814  return(values_[colIndex * stride_ + ku_+rowIndex-colIndex]);
815 }
816 
817 template<typename OrdinalType, typename ScalarType>
818 inline const ScalarType& SerialBandDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex) const
819 {
820 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
821  checkIndex( rowIndex, colIndex );
822 #endif
823  return(values_[colIndex * stride_ + ku_+rowIndex-colIndex]);
824 }
825 
826 template<typename OrdinalType, typename ScalarType>
827 inline const ScalarType* SerialBandDenseMatrix<OrdinalType, ScalarType>::operator [] (OrdinalType colIndex) const
828 {
829 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
830  checkIndex( 0, colIndex );
831 #endif
832  return(values_ + colIndex * stride_);
833 }
834 
835 template<typename OrdinalType, typename ScalarType>
836 inline ScalarType* SerialBandDenseMatrix<OrdinalType, ScalarType>::operator [] (OrdinalType colIndex)
837 {
838 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
839  checkIndex( 0, colIndex );
840 #endif
841  return(values_ + colIndex * stride_);
842 }
843 
844 //----------------------------------------------------------------------------------------------------
845 // Norm methods
846 //----------------------------------------------------------------------------------------------------
847 
848 template<typename OrdinalType, typename ScalarType>
850 {
851  OrdinalType i, j;
854 
855  ScalarType* ptr;
856  for(j = 0; j < numCols_; j++) {
857  ScalarType sum = 0;
858  ptr = values_ + j * stride_ + TEUCHOS_MAX(0, ku_-j);
859  for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
861  }
863  if(absSum > anorm) {
864  anorm = absSum;
865  }
866  }
867  updateFlops((kl_+ku_+1) * numCols_);
868 
869  return(anorm);
870 }
871 
872 template<typename OrdinalType, typename ScalarType>
874 {
875  OrdinalType i, j;
877 
878  for (i = 0; i < numRows_; i++) {
880  for (j=TEUCHOS_MAX(0,i-kl_); j<=TEUCHOS_MIN(numCols_-1,i+ku_); j++) {
881  sum += ScalarTraits<ScalarType>::magnitude(*(values_+(ku_+i-j)+j*stride_));
882  }
883  anorm = TEUCHOS_MAX( anorm, sum );
884  }
885  updateFlops((kl_+ku_+1) * numCols_);
886 
887  return(anorm);
888 }
889 
890 template<typename OrdinalType, typename ScalarType>
892 {
893  OrdinalType i, j;
895 
896  for (j = 0; j < numCols_; j++) {
897  for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
898  anorm += ScalarTraits<ScalarType>::magnitude(values_[(ku_+i-j)+j*stride_]*values_[(ku_+i-j)+j*stride_]);
899  }
900  }
902  updateFlops((kl_+ku_+1) * numCols_);
903 
904  return(anorm);
905 }
906 
907 //----------------------------------------------------------------------------------------------------
908 // Comparison methods
909 //----------------------------------------------------------------------------------------------------
910 
911 template<typename OrdinalType, typename ScalarType>
913 {
914  bool result = 1;
915 
916  if((numRows_ != Operand.numRows_) || (numCols_ != Operand.numCols_) || (kl_ != Operand.kl_) || (ku_ != Operand.ku_)) {
917  result = 0;
918  } else {
919  OrdinalType i, j;
920  for(j = 0; j < numCols_; j++) {
921  for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
922  if((*this)(i, j) != Operand(i, j)) {
923  return 0;
924  }
925  }
926  }
927  }
928 
929  return result;
930 }
931 
932 template<typename OrdinalType, typename ScalarType>
934 {
935  return !((*this) == Operand);
936 }
937 
938 //----------------------------------------------------------------------------------------------------
939 // Multiplication method
940 //----------------------------------------------------------------------------------------------------
941 
942 template<typename OrdinalType, typename ScalarType>
944 {
945  this->scale( alpha );
946  return(*this);
947 }
948 
949 template<typename OrdinalType, typename ScalarType>
951 {
952 
953  OrdinalType i, j;
954  ScalarType* ptr;
955 
956  for (j=0; j<numCols_; j++) {
957  ptr = values_ + j*stride_ + TEUCHOS_MAX(0, ku_-j);
958  for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
959  *ptr = alpha * (*ptr); ptr++;
960  }
961  }
962  updateFlops( (kl_+ku_+1)*numCols_ );
963 
964  return(0);
965 }
966 
967 template<typename OrdinalType, typename ScalarType>
969 {
970 
971  OrdinalType i, j;
972  ScalarType* ptr;
973 
974  // Check for compatible dimensions
975  if ((numRows_ != A.numRows_) || (numCols_ != A.numCols_) || (kl_ != A.kl_) || (ku_ != A.ku_)) {
976  TEUCHOS_CHK_ERR(-1); // Return error
977  }
978  for (j=0; j<numCols_; j++) {
979  ptr = values_ + j*stride_ + TEUCHOS_MAX(0, ku_-j);
980  for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_-1,j+kl_); i++) {
981  *ptr = A(i,j) * (*ptr); ptr++;
982  }
983  }
984  updateFlops( (kl_+ku_+1)*numCols_ );
985 
986  return(0);
987 }
988 
989 template<typename OrdinalType, typename ScalarType>
991 {
992  os << std::endl;
993  if(valuesCopied_)
994  os << "Values_copied : yes" << std::endl;
995  else
996  os << "Values_copied : no" << std::endl;
997  os << "Rows : " << numRows_ << std::endl;
998  os << "Columns : " << numCols_ << std::endl;
999  os << "Lower Bandwidth : " << kl_ << std::endl;
1000  os << "Upper Bandwidth : " << ku_ << std::endl;
1001  os << "LDA : " << stride_ << std::endl;
1002  if(numRows_ == 0 || numCols_ == 0) {
1003  os << "(matrix is empty, no values to display)" << std::endl;
1004  } else {
1005 
1006  for(OrdinalType i = 0; i < numRows_; i++) {
1007  for (OrdinalType j=TEUCHOS_MAX(0,i-kl_); j<=TEUCHOS_MIN(numCols_-1,i+ku_); j++) {
1008  os << (*this)(i,j) << " ";
1009  }
1010  os << std::endl;
1011  }
1012  }
1013 }
1014 
1015 //----------------------------------------------------------------------------------------------------
1016 // Protected methods
1017 //----------------------------------------------------------------------------------------------------
1018 
1019 template<typename OrdinalType, typename ScalarType>
1020 inline void SerialBandDenseMatrix<OrdinalType, ScalarType>::checkIndex( OrdinalType rowIndex, OrdinalType colIndex ) const {
1021 
1022  TEUCHOS_TEST_FOR_EXCEPTION(rowIndex < 0 || rowIndex >= numRows_ ||
1023  rowIndex < TEUCHOS_MAX(0,colIndex-ku_) || rowIndex > TEUCHOS_MIN(numRows_-1,colIndex+kl_),
1024  std::out_of_range,
1025  "SerialBandDenseMatrix<T>::checkIndex: "
1026  "Row index " << rowIndex << " out of range [0, "<< numRows_ << ")");
1027  TEUCHOS_TEST_FOR_EXCEPTION(colIndex < 0 || colIndex >= numCols_, std::out_of_range,
1028  "SerialBandDenseMatrix<T>::checkIndex: "
1029  "Col index " << colIndex << " out of range [0, "<< numCols_ << ")");
1030 
1031 }
1032 
1033 template<typename OrdinalType, typename ScalarType>
1035 {
1036  if (valuesCopied_) {
1037  delete [] values_;
1038  values_ = 0;
1039  valuesCopied_ = false;
1040  }
1041 }
1042 
1043 template<typename OrdinalType, typename ScalarType>
1045  ScalarType* inputMatrix, OrdinalType strideInput, OrdinalType numRows_in,
1046  OrdinalType numCols_in, ScalarType* outputMatrix, OrdinalType strideOutput, OrdinalType startCol, ScalarType alpha
1047  )
1048 {
1049  OrdinalType i, j;
1050  ScalarType* ptr1 = 0;
1051  ScalarType* ptr2 = 0;
1052 
1053  for(j = 0; j < numCols_in; j++) {
1054  ptr1 = outputMatrix + (j * strideOutput) + TEUCHOS_MAX(0, ku_-j);
1055  ptr2 = inputMatrix + (j + startCol) * strideInput + TEUCHOS_MAX(0, ku_-j);
1056  if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
1057  for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_in-1,j+kl_); i++) {
1058  *ptr1++ += alpha*(*ptr2++);
1059  }
1060  } else {
1061  for (i=TEUCHOS_MAX(0,j-ku_); i<=TEUCHOS_MIN(numRows_in-1,j+kl_); i++) {
1062  *ptr1++ = *ptr2++;
1063  }
1064  }
1065  }
1066 }
1067 
1068 } // namespace Teuchos
1069 
1070 
1071 #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.
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.
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...
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.
virtual void print(std::ostream &os) const
Print method. Defines the behavior of the std::ostream << operator inherited from the Object class...
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.
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.
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.
The Teuchos namespace contains all of the classes, structs and enums used by Teuchos, as well as a number of utility routines.
void updateFlops(int addflops) const
Increment Flop count for this object.
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.