Teuchos - Trilinos Tools Package Version of the Day
Loading...
Searching...
No Matches
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_
15
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
93
97
98namespace Teuchos {
99
100template<typename OrdinalType, typename ScalarType>
101class SerialBandDenseMatrix : public CompObject, public BLAS<OrdinalType, ScalarType>
102{
103public:
104
106 typedef OrdinalType ordinalType;
108 typedef ScalarType scalarType;
109
111
112
114
117 SerialBandDenseMatrix();
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
151 SerialBandDenseMatrix(const SerialBandDenseMatrix<OrdinalType, ScalarType> &Source, ETransp trans = Teuchos::NO_TRANS);
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
222 SerialBandDenseMatrix<OrdinalType, ScalarType>& operator= (const SerialBandDenseMatrix<OrdinalType, ScalarType>& Source);
223
225
230 SerialBandDenseMatrix<OrdinalType, ScalarType>& assign (const SerialBandDenseMatrix<OrdinalType, ScalarType>& Source);
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
306 SerialBandDenseMatrix<OrdinalType, ScalarType>& operator+= (const SerialBandDenseMatrix<OrdinalType, ScalarType>& Source);
307
309
312 SerialBandDenseMatrix<OrdinalType, ScalarType>& operator-= (const SerialBandDenseMatrix<OrdinalType, ScalarType>& Source);
313
315
318 SerialBandDenseMatrix<OrdinalType, ScalarType>& operator*= (const ScalarType alpha);
319
321
325 int scale ( const ScalarType alpha );
326
328
334 int scale ( const SerialBandDenseMatrix<OrdinalType, ScalarType>& A );
335
337
339
340
342
345 bool operator== (const SerialBandDenseMatrix<OrdinalType, ScalarType> &Operand) const;
346
348
351 bool operator!= (const SerialBandDenseMatrix<OrdinalType, ScalarType> &Operand) const;
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
381 typename ScalarTraits<ScalarType>::magnitudeType normOne() const;
382
384 typename ScalarTraits<ScalarType>::magnitudeType normInf() const;
385
387 typename ScalarTraits<ScalarType>::magnitudeType normFrobenius() const;
389
391
392
393
394 virtual std::ostream& print(std::ostream& os) const;
395
396
398protected:
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
418template<typename OrdinalType, typename ScalarType>
419SerialBandDenseMatrix<OrdinalType, ScalarType>::SerialBandDenseMatrix ()
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
431template<typename OrdinalType, typename ScalarType>
432SerialBandDenseMatrix<OrdinalType, ScalarType>::
433SerialBandDenseMatrix (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
454template<typename OrdinalType, typename ScalarType>
455SerialBandDenseMatrix<OrdinalType, ScalarType>::
456SerialBandDenseMatrix (DataAccess CV,
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
481template<typename OrdinalType, typename ScalarType>
482SerialBandDenseMatrix<OrdinalType, ScalarType>::
483SerialBandDenseMatrix (const SerialBandDenseMatrix<OrdinalType, ScalarType> &Source, ETransp trans)
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
535template<typename OrdinalType, typename ScalarType>
536SerialBandDenseMatrix<OrdinalType, ScalarType>::
537SerialBandDenseMatrix (DataAccess CV,
538 const SerialBandDenseMatrix<OrdinalType, ScalarType> &Source,
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
562template<typename OrdinalType, typename ScalarType>
567
568//----------------------------------------------------------------------------------------------------
569// Shape methods
570//----------------------------------------------------------------------------------------------------
571
572template<typename OrdinalType, typename ScalarType>
573int SerialBandDenseMatrix<OrdinalType, ScalarType>::shape(
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
590template<typename OrdinalType, typename ScalarType>
591int SerialBandDenseMatrix<OrdinalType, ScalarType>::shapeUninitialized(
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
607template<typename OrdinalType, typename ScalarType>
608int SerialBandDenseMatrix<OrdinalType, ScalarType>::reshape(
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
641template<typename OrdinalType, typename ScalarType>
642int SerialBandDenseMatrix<OrdinalType, ScalarType>::putScalar( const ScalarType value_in )
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
655template<typename OrdinalType, typename ScalarType>
656int SerialBandDenseMatrix<OrdinalType, ScalarType>::random()
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
669template<typename OrdinalType, typename ScalarType>
670SerialBandDenseMatrix<OrdinalType,ScalarType>&
671SerialBandDenseMatrix<OrdinalType, ScalarType>::operator=(
672 const SerialBandDenseMatrix<OrdinalType,ScalarType>& Source
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
736template<typename OrdinalType, typename ScalarType>
737SerialBandDenseMatrix<OrdinalType, ScalarType>& SerialBandDenseMatrix<OrdinalType, ScalarType>::operator+= (const SerialBandDenseMatrix<OrdinalType,ScalarType>& Source )
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
749template<typename OrdinalType, typename ScalarType>
750SerialBandDenseMatrix<OrdinalType, ScalarType>& SerialBandDenseMatrix<OrdinalType, ScalarType>::operator-= (const SerialBandDenseMatrix<OrdinalType,ScalarType>& Source )
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
762template<typename OrdinalType, typename ScalarType>
763SerialBandDenseMatrix<OrdinalType,ScalarType>& SerialBandDenseMatrix<OrdinalType, ScalarType>::assign (const SerialBandDenseMatrix<OrdinalType,ScalarType>& Source) {
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
783template<typename OrdinalType, typename ScalarType>
784inline 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
792template<typename OrdinalType, typename ScalarType>
793inline 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
801template<typename OrdinalType, typename ScalarType>
802inline 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
810template<typename OrdinalType, typename ScalarType>
811inline 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
823template<typename OrdinalType, typename ScalarType>
824typename ScalarTraits<ScalarType>::magnitudeType SerialBandDenseMatrix<OrdinalType, ScalarType>::normOne() const
825{
826 OrdinalType i, j;
827 typename ScalarTraits<ScalarType>::magnitudeType anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
828 typename ScalarTraits<ScalarType>::magnitudeType absSum = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
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++) {
835 sum += ScalarTraits<ScalarType>::magnitude(*ptr++);
836 }
837 absSum = ScalarTraits<ScalarType>::magnitude(sum);
838 if(absSum > anorm) {
839 anorm = absSum;
840 }
841 }
842 updateFlops((kl_+ku_+1) * numCols_);
843
844 return(anorm);
845}
846
847template<typename OrdinalType, typename ScalarType>
848typename ScalarTraits<ScalarType>::magnitudeType SerialBandDenseMatrix<OrdinalType, ScalarType>::normInf() const
849{
850 OrdinalType i, j;
851 typename ScalarTraits<ScalarType>::magnitudeType sum, anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
852
853 for (i = 0; i < numRows_; i++) {
854 sum = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
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
865template<typename OrdinalType, typename ScalarType>
866typename ScalarTraits<ScalarType>::magnitudeType SerialBandDenseMatrix<OrdinalType, ScalarType>::normFrobenius() const
867{
868 OrdinalType i, j;
869 typename ScalarTraits<ScalarType>::magnitudeType anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
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 }
876 anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::squareroot(anorm));
877 updateFlops((kl_+ku_+1) * numCols_);
878
879 return(anorm);
880}
881
882//----------------------------------------------------------------------------------------------------
883// Comparison methods
884//----------------------------------------------------------------------------------------------------
885
886template<typename OrdinalType, typename ScalarType>
887bool SerialBandDenseMatrix<OrdinalType, ScalarType>::operator== (const SerialBandDenseMatrix<OrdinalType, ScalarType> &Operand) const
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
907template<typename OrdinalType, typename ScalarType>
908bool SerialBandDenseMatrix<OrdinalType, ScalarType>::operator!= (const SerialBandDenseMatrix<OrdinalType, ScalarType> &Operand) const
909{
910 return !((*this) == Operand);
911}
912
913//----------------------------------------------------------------------------------------------------
914// Multiplication method
915//----------------------------------------------------------------------------------------------------
916
917template<typename OrdinalType, typename ScalarType>
918SerialBandDenseMatrix<OrdinalType, ScalarType>& SerialBandDenseMatrix<OrdinalType, ScalarType>::operator*= (const ScalarType alpha )
919{
920 this->scale( alpha );
921 return(*this);
922}
923
924template<typename OrdinalType, typename ScalarType>
925int SerialBandDenseMatrix<OrdinalType, ScalarType>::scale( const ScalarType alpha )
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
942template<typename OrdinalType, typename ScalarType>
943int SerialBandDenseMatrix<OrdinalType, ScalarType>::scale( const SerialBandDenseMatrix<OrdinalType,ScalarType>& A )
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
964template<typename OrdinalType, typename ScalarType>
965std::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
995template<typename OrdinalType, typename ScalarType>
996inline 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
1009template<typename OrdinalType, typename ScalarType>
1010void SerialBandDenseMatrix<OrdinalType, ScalarType>::deleteArrays(void)
1011{
1012 if (valuesCopied_) {
1013 delete [] values_;
1014 values_ = 0;
1015 valuesCopied_ = false;
1016 }
1017}
1018
1019template<typename OrdinalType, typename ScalarType>
1020void 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
1045template<typename OrdinalType, typename ScalarType>
1046struct SerialBandDenseMatrixPrinter {
1047public:
1048 const SerialBandDenseMatrix<OrdinalType,ScalarType> &obj;
1049 SerialBandDenseMatrixPrinter(
1050 const SerialBandDenseMatrix<OrdinalType,ScalarType> &obj_in)
1051 : obj(obj_in) {}
1052};
1053
1055template<typename OrdinalType, typename ScalarType>
1056std::ostream&
1057operator<<(std::ostream &out,
1058 const SerialBandDenseMatrixPrinter<OrdinalType,ScalarType> printer)
1059{
1060 printer.obj.print(out);
1061 return out;
1062}
1063
1065template<typename OrdinalType, typename ScalarType>
1066SerialBandDenseMatrixPrinter<OrdinalType,ScalarType>
1071
1072
1073} // namespace Teuchos
1074
1075
1076#endif /* _TEUCHOS_SERIALBANDDENSEMATRIX_HPP_ */
Teuchos_BLAS.hpp
Templated interface class to BLAS routines.
Teuchos_CompObject.hpp
Object for storing data and providing functionality that is common to all computational classes.
Teuchos_ConfigDefs.hpp
Teuchos header file which uses auto-configuration information to include necessary C++ headers.
Teuchos_DataAccess.hpp
Teuchos::DataAccess Mode enumerable type.
Teuchos_RCP.hpp
Reference-counted pointer class and non-member templated function implementations.
Teuchos_ScalarTraits.hpp
Defines basic traits for the scalar field type.
Teuchos::BLAS::BLAS
BLAS(void)
Default constructor.
Definition Teuchos_BLAS.hpp:222
Teuchos::CompObject::updateFlops
void updateFlops(int addflops) const
Increment Flop count for this object.
Definition Teuchos_CompObject.hpp:78
Teuchos::CompObject::CompObject
CompObject()
Default constructor.
Definition Teuchos_CompObject.cpp:20
Teuchos::Ptr::ptr
Ptr< T > ptr(T *p)
Create a pointer to an object from a raw pointer.
Definition Teuchos_PtrDecl.hpp:296
Teuchos::SerialBandDenseMatrix
This class creates and provides basic support for banded dense matrices of templated type.
Definition Teuchos_SerialBandDenseMatrix.hpp:102
Teuchos::SerialBandDenseMatrix::operator+=
SerialBandDenseMatrix< OrdinalType, ScalarType > & operator+=(const SerialBandDenseMatrix< OrdinalType, ScalarType > &Source)
Add another matrix to this matrix.
Definition Teuchos_SerialBandDenseMatrix.hpp:737
Teuchos::SerialBandDenseMatrix::shapeUninitialized
int shapeUninitialized(OrdinalType numRows, OrdinalType numCols, OrdinalType kl, OrdinalType ku)
Same as shape() except leaves uninitialized.
Definition Teuchos_SerialBandDenseMatrix.hpp:591
Teuchos::SerialBandDenseMatrix::random
int random()
Set all values in the matrix to be random numbers.
Definition Teuchos_SerialBandDenseMatrix.hpp:656
Teuchos::SerialBandDenseMatrix::numCols
OrdinalType numCols() const
Returns the column dimension of this matrix.
Definition Teuchos_SerialBandDenseMatrix.hpp:362
Teuchos::SerialBandDenseMatrix::print
virtual std::ostream & print(std::ostream &os) const
Print method. Defines the behavior of the std::ostream << operator.
Definition Teuchos_SerialBandDenseMatrix.hpp:965
Teuchos::SerialBandDenseMatrix::normFrobenius
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
Returns the Frobenius-norm of the matrix.
Definition Teuchos_SerialBandDenseMatrix.hpp:866
Teuchos::SerialBandDenseMatrix::scalarType
ScalarType scalarType
Typedef for scalar type.
Definition Teuchos_SerialBandDenseMatrix.hpp:108
Teuchos::SerialBandDenseMatrix::operator=
SerialBandDenseMatrix< OrdinalType, ScalarType > & operator=(const SerialBandDenseMatrix< OrdinalType, ScalarType > &Source)
Copies values from one matrix to another.
Definition Teuchos_SerialBandDenseMatrix.hpp:671
Teuchos::SerialBandDenseMatrix::empty
bool empty() const
Returns whether this matrix is empty.
Definition Teuchos_SerialBandDenseMatrix.hpp:374
Teuchos::SerialBandDenseMatrix::ordinalType
OrdinalType ordinalType
Typedef for ordinal type.
Definition Teuchos_SerialBandDenseMatrix.hpp:106
Teuchos::SerialBandDenseMatrix::values
ScalarType * values() const
Data array access method.
Definition Teuchos_SerialBandDenseMatrix.hpp:295
Teuchos::SerialBandDenseMatrix::numRows
OrdinalType numRows() const
Returns the row dimension of this matrix.
Definition Teuchos_SerialBandDenseMatrix.hpp:359
Teuchos::SerialBandDenseMatrix::scale
int scale(const ScalarType alpha)
Scale this matrix by alpha; *this = alpha**this.
Definition Teuchos_SerialBandDenseMatrix.hpp:925
Teuchos::SerialBandDenseMatrix::operator-=
SerialBandDenseMatrix< OrdinalType, ScalarType > & operator-=(const SerialBandDenseMatrix< OrdinalType, ScalarType > &Source)
Subtract another matrix from this matrix.
Definition Teuchos_SerialBandDenseMatrix.hpp:750
Teuchos::SerialBandDenseMatrix::operator[]
ScalarType * operator[](OrdinalType colIndex)
Column access method (non-const).
Definition Teuchos_SerialBandDenseMatrix.hpp:811
Teuchos::SerialBandDenseMatrix::putScalar
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero())
Set all values in the matrix to a constant value.
Definition Teuchos_SerialBandDenseMatrix.hpp:642
Teuchos::SerialBandDenseMatrix::operator*=
SerialBandDenseMatrix< OrdinalType, ScalarType > & operator*=(const ScalarType alpha)
Scale this matrix by alpha; *this = alpha**this.
Definition Teuchos_SerialBandDenseMatrix.hpp:918
Teuchos::SerialBandDenseMatrix::lowerBandwidth
OrdinalType lowerBandwidth() const
Returns the lower bandwidth of this matrix.
Definition Teuchos_SerialBandDenseMatrix.hpp:365
Teuchos::SerialBandDenseMatrix::operator==
bool operator==(const SerialBandDenseMatrix< OrdinalType, ScalarType > &Operand) const
Equality of two matrices.
Definition Teuchos_SerialBandDenseMatrix.hpp:887
Teuchos::SerialBandDenseMatrix::reshape
int reshape(OrdinalType numRows, OrdinalType numCols, OrdinalType kl, OrdinalType ku)
Reshaping method for changing the size of a SerialBandDenseMatrix, keeping the entries.
Definition Teuchos_SerialBandDenseMatrix.hpp:608
Teuchos::SerialBandDenseMatrix::stride
OrdinalType stride() const
Returns the stride between the columns of this matrix in memory.
Definition Teuchos_SerialBandDenseMatrix.hpp:371
Teuchos::SerialBandDenseMatrix::assign
SerialBandDenseMatrix< OrdinalType, ScalarType > & assign(const SerialBandDenseMatrix< OrdinalType, ScalarType > &Source)
Copies values from one matrix to another.
Definition Teuchos_SerialBandDenseMatrix.hpp:763
Teuchos::SerialBandDenseMatrix::operator!=
bool operator!=(const SerialBandDenseMatrix< OrdinalType, ScalarType > &Operand) const
Inequality of two matrices.
Definition Teuchos_SerialBandDenseMatrix.hpp:908
Teuchos::SerialBandDenseMatrix::operator()
ScalarType & operator()(OrdinalType rowIndex, OrdinalType colIndex)
Element access method (non-const).
Definition Teuchos_SerialBandDenseMatrix.hpp:784
Teuchos::SerialBandDenseMatrix::shape
int shape(OrdinalType numRows, OrdinalType numCols, OrdinalType kl, OrdinalType ku)
Shape method for changing the size of a SerialBandDenseMatrix, initializing entries to zero.
Definition Teuchos_SerialBandDenseMatrix.hpp:573
Teuchos::SerialBandDenseMatrix::upperBandwidth
OrdinalType upperBandwidth() const
Returns the upper bandwidth of this matrix.
Definition Teuchos_SerialBandDenseMatrix.hpp:368
Teuchos::SerialBandDenseMatrix::normInf
ScalarTraits< ScalarType >::magnitudeType normInf() const
Returns the Infinity-norm of the matrix.
Definition Teuchos_SerialBandDenseMatrix.hpp:848
Teuchos::SerialBandDenseMatrix::normOne
ScalarTraits< ScalarType >::magnitudeType normOne() const
Returns the 1-norm of the matrix.
Definition Teuchos_SerialBandDenseMatrix.hpp:824
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Macro for throwing an exception with breakpointing to ease debugging.
Definition Teuchos_TestForException.hpp:146
Teuchos
The Teuchos namespace contains all of the classes, structs and enums used by Teuchos,...
Teuchos::printMat
SerialBandDenseMatrixPrinter< OrdinalType, ScalarType > printMat(const SerialBandDenseMatrix< OrdinalType, ScalarType > &obj)
Return SerialBandDenseMatrix ostream manipulator Use as:
Definition Teuchos_SerialBandDenseMatrix.hpp:1067
Teuchos::ScalarTraits::magnitude
static magnitudeType magnitude(T a)
Returns the magnitudeType of the scalar type a.
Definition Teuchos_ScalarTraitsDecl.hpp:100
Teuchos::ScalarTraits::one
static T one()
Returns representation of one for this scalar type.
Definition Teuchos_ScalarTraitsDecl.hpp:104
Teuchos::ScalarTraits::magnitudeType
T magnitudeType
Mandatory typedef for result of magnitude.
Definition Teuchos_ScalarTraitsDecl.hpp:61
Teuchos::ScalarTraits::zero
static T zero()
Returns representation of zero for this scalar type.
Definition Teuchos_ScalarTraitsDecl.hpp:102
Teuchos::ScalarTraits::conjugate
static T conjugate(T a)
Returns the conjugate of the scalar type a.
Definition Teuchos_ScalarTraitsDecl.hpp:110
Teuchos::ScalarTraits::random
static T random()
Returns a random number (between -one() and +one()) of this scalar type.
Definition Teuchos_ScalarTraitsDecl.hpp:118
Teuchos::ScalarTraits::squareroot
static T squareroot(T x)
Returns a number of magnitudeType that is the square root of this scalar type x.
Definition Teuchos_ScalarTraitsDecl.hpp:122
Teuchos::ScalarTraits::isComplex
static const bool isComplex
Determines if scalar type is std::complex.
Definition Teuchos_ScalarTraitsDecl.hpp:69
Teuchos::SerialBandDenseMatrixPrinter
Ostream manipulator for SerialBandDenseMatrix.
Definition Teuchos_SerialBandDenseMatrix.hpp:1046
Generated on for Teuchos - Trilinos Tools Package by doxygen 1.15.0