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