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