Teuchos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Teuchos_SerialTriDiMatrix.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
43#ifndef _TEUCHOS_SERIALTRIDIMATRIX_HPP_
44#define _TEUCHOS_SERIALTRIDIMATRIX_HPP_
50#include "Teuchos_BLAS.hpp"
54#include "Teuchos_Assert.hpp"
55
64namespace Teuchos {
65
66template<typename OrdinalType, typename ScalarType>
67class SerialTriDiMatrix : public CompObject, public Object, public BLAS<OrdinalType, ScalarType > {
68public:
69
71 typedef OrdinalType ordinalType;
73 typedef ScalarType scalarType;
74
76
77
79
83
85
93 SerialTriDiMatrix(OrdinalType numRows, OrdinalType numCols, bool zeroOut = true);
94
96
101 SerialTriDiMatrix(DataAccess CV, ScalarType* values, OrdinalType numRowsCols);
102
104
110
112
123 SerialTriDiMatrix(DataAccess CV, const SerialTriDiMatrix<OrdinalType, ScalarType> &Source, OrdinalType numRowsCols, OrdinalType startRowCols=0);
124
126 virtual ~SerialTriDiMatrix();
128
130
131
141 int shape(OrdinalType numRows);
142
144 int shapeUninitialized(OrdinalType numRows);
145
147
156 int reshape(OrdinalType numRowsCols);
157
159
161
162
164
171
173
179
181
184 SerialTriDiMatrix<OrdinalType, ScalarType>& operator= (const ScalarType value) { putScalar(value); return(*this); }
185
187
191 int putScalar( const ScalarType value = Teuchos::ScalarTraits<ScalarType>::zero() );
192
194 // int random();
195
197
199
200
202
209 ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex);
210
212
219 const ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex) const;
220
222
229 // ScalarType* operator [] (OrdinalType colIndex);
230
232
239 // const ScalarType* operator [] (OrdinalType colIndex) const;
240
242
243 ScalarType* values() const { return(values_); }
244
245 ScalarType* D() const { return D_;}
246 ScalarType* DL() const { return DL_;}
247 ScalarType* DU() const { return DU_;}
248 ScalarType* DU2() const { return DU2_;}
249
251
253
254
256
260
262
266
268
272
274
278 int scale ( const ScalarType alpha );
279
281
288
290
304 //int multiply (ETransp transa, ETransp transb, ScalarType alpha, const SerialTriDiMatrix<OrdinalType, ScalarType> &A, const SerialTriDiMatrix<OrdinalType, ScalarType> &B, ScalarType beta);
305
307
318 //int multiply (ESide sideA, ScalarType alpha, const SerialSymTriDiMatrix<OrdinalType, ScalarType> &A, const SerialTriDiMatrix<OrdinalType, ScalarType> &B, ScalarType beta);
319
321
323
324
326
330
332
336
338
340
341
343 OrdinalType numRowsCols() const { return(numRowsCols_); }
344
346 // OrdinalType numCols() const { return(numRowsCols_); }
347
349 // OrdinalType stride() const { return(stride_); }
350
352 bool empty() const { return(numRowsCols_ == 0); }
354
356
357
360
363
367
369
370
371 virtual void print(std::ostream& os) const;
372
374protected:
376 OrdinalType startCol,
377 ScalarType alpha = ScalarTraits<ScalarType>::zero() );
378 void deleteArrays();
379 void checkIndex( OrdinalType rowIndex, OrdinalType colIndex = 0 ) const;
380 OrdinalType numRowsCols_;
381
383 ScalarType* values_;
384 ScalarType* DL_;
385 ScalarType* D_;
386 ScalarType* DU_;
387 ScalarType* DU2_;
388
389}; // class Teuchos_SerialTriDiMatrix
390
391//----------------------------------------------------------------------------------------------------
392// Constructors and Destructor
393//----------------------------------------------------------------------------------------------------
394
395template<typename OrdinalType, typename ScalarType>
397 :
398 CompObject(),
399 numRowsCols_(0),
400 valuesCopied_(false),
401 values_(0),
402 DL_(NULL),
403 D_(NULL),
404 DU_(NULL),
405 DU2_(NULL)
406{}
407
408template<typename OrdinalType, typename ScalarType>
409SerialTriDiMatrix<OrdinalType, ScalarType>::SerialTriDiMatrix( OrdinalType numRowsCols_in, OrdinalType /* numCols_in */, bool zeroOut)
410 : CompObject(), numRowsCols_(numRowsCols_in) {
411
412 OrdinalType numvals = (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
413 values_ = new ScalarType [numvals];
414 DL_ = values_;
415 D_ = DL_ + (numRowsCols_-1);
416 DU_ = D_ + numRowsCols_;
417 DU2_ = DU_ + (numRowsCols_-1);
418
419 valuesCopied_ = true;
420 if (zeroOut == true)
421 putScalar();
422}
423
424template<typename OrdinalType, typename ScalarType>
425SerialTriDiMatrix<OrdinalType, ScalarType>::SerialTriDiMatrix(DataAccess CV, ScalarType* values_in, OrdinalType numRowsCols_in )
426 : CompObject(), numRowsCols_(numRowsCols_in),
427 valuesCopied_(false), values_(values_in)
428{
429 const OrdinalType numvals = (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
430 if(CV == Copy) {
431 values_ = new ScalarType[numvals];
432 valuesCopied_ = true;
433 }
434 else //CV == View
435 {
436 values_ = values_in;
437 valuesCopied_ = false;
438 }
439 DL_ = values_;
440 D_ = DL_ + (numRowsCols_-1);
441 DU_ = D_ + numRowsCols_;
442 DU2_ = DU_ + (numRowsCols_-1);
443 if(CV == Copy) {
444 for(OrdinalType i = 0 ; i < numRowsCols_ ; ++i )
445 values_[i] = values_in[i];
446 }
447}
448
449template<typename OrdinalType, typename ScalarType>
450SerialTriDiMatrix<OrdinalType, ScalarType>::SerialTriDiMatrix(const SerialTriDiMatrix<OrdinalType, ScalarType> &Source, ETransp trans) : CompObject(), BLAS<OrdinalType,ScalarType>(), numRowsCols_(0), valuesCopied_(true), values_(0)
451{
452 if ( trans == Teuchos::NO_TRANS ) {
453 numRowsCols_ = Source.numRowsCols_;
454
455 const OrdinalType numvals = (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
456 values_ = new ScalarType[numvals];
457 DL_ = values_;
458 D_ = DL_+ (numRowsCols_-1);
459 DU_ = D_ + numRowsCols_;
460 DU2_ = DU_ + (numRowsCols_-1);
461
462 copyMat(Source, 0, 0);
463 }
465 {
466 numRowsCols_ = Source.numRowsCols_;
467 const OrdinalType numvals = (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
468 values_ = new ScalarType[numvals];
469 DL_ = values_;
470 D_ = DL_+(numRowsCols_-1);
471 DU_ = D_ + numRowsCols_;
472 DU2_ = DU_ + (numRowsCols_-1);
473
474 OrdinalType min = numRowsCols_;
475 if(min > Source.numRowsCols_) min = Source.numRowsCols_;
476
477 for(OrdinalType i = 0 ; i< min ; ++i) {
479 if(i < (min-1)) {
482 }
483 if(i < (min-2)) {
485 }
486 }
487 }
488 else
489 {
490 numRowsCols_ = Source.numRowsCols_;
491 const OrdinalType numvals = (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
492 values_ = new ScalarType[numvals];
493 OrdinalType min = numRowsCols_;
494 if(min > Source.numRowsCols_) min = Source.numRowsCols_;
495 for(OrdinalType i = 0 ; i< min ; ++i) {
496 D_[i] = Source.D_[i];
497 if(i < (min-1)) {
498 DL_[i] = Source.DL_[i];
499 DU_[i] = Source.DU_[i];
500 }
501 if(i < (min-2)) {
502 DU2_[i] = Source.DU2_[i];
503 }
504 }
505 }
506}
507
508template<typename OrdinalType, typename ScalarType>
511 OrdinalType numRowsCols_in, OrdinalType startRow )
512 : CompObject(), numRowsCols_(numRowsCols_in),
513 valuesCopied_(false), values_(Source.values_) {
514 if(CV == Copy)
515 {
516 const OrdinalType numvals = (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
517 values_ = new ScalarType[numvals];
518 copyMat(Source, startRow);
519 valuesCopied_ = true;
520 }
521 else // CV == View
522 {
523 // \todo ???
524 // values_ = values_ + (stride_ * startCol) + startRow;
525 }
526}
527
528template<typename OrdinalType, typename ScalarType>
530{
531 deleteArrays();
532}
533
534//----------------------------------------------------------------------------------------------------
535// Shape methods
536//----------------------------------------------------------------------------------------------------
537
538template<typename OrdinalType, typename ScalarType>
540{
541 deleteArrays(); // Get rid of anything that might be already allocated
542 numRowsCols_ = numRowsCols_in;
543 const OrdinalType numvals = ( numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
544 values_ = new ScalarType[numvals];
545
546 putScalar();
547 valuesCopied_ = true;
548 return(0);
549}
550
551template<typename OrdinalType, typename ScalarType>
553{
554 deleteArrays(); // Get rid of anything that might be already allocated
555 numRowsCols_ = numRowsCols_in;
556 const OrdinalType numvals = ( numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
557 values_ = new ScalarType[numvals];
558 valuesCopied_ = true;
559 return(0);
560}
561
562template<typename OrdinalType, typename ScalarType>
564{
565 if(numRowsCols_in <1 ) {
566 deleteArrays();
567 return 0;
568 }
569 // Allocate space for new matrix
570 const OrdinalType numvals = ( numRowsCols_in == 1) ? 1 : 4*(numRowsCols_in - 1);
571 ScalarType* values_tmp = new ScalarType[numvals];
572 ScalarType zero = ScalarTraits<ScalarType>::zero();
573 for(OrdinalType i= 0; i<numvals ; ++i)
574 values_tmp[i] = zero;
575
576 OrdinalType min = TEUCHOS_MIN(numRowsCols_, numRowsCols_in);
577 ScalarType* dl = values_tmp;
578 ScalarType* d = values_tmp + (numRowsCols_in-1);
579 ScalarType* du = d+(numRowsCols_in);
580 ScalarType* du2 = du+(numRowsCols_in - 1);
581
582 if(values_ != 0) {
583 for(OrdinalType i = 0 ; i< min ; ++i) {
584 dl[i] = DL_[i];
585 d[i] = D_[i];
586 du[i] = DU_[i];
587 du2[i] = DU2_[i];
588 }
589 }
590
591 deleteArrays(); // Get rid of anything that might be already allocated
592 numRowsCols_ = numRowsCols_in;
593
594 values_ = values_tmp; // Set pointer to new A
595 DL_ = dl;
596 D_ = d;
597 DU_ = du;
598 DU2_ = du2;
599
600 valuesCopied_ = true;
601 return(0);
602}
603
604//----------------------------------------------------------------------------------------------------
605// Set methods
606//----------------------------------------------------------------------------------------------------
607
608template<typename OrdinalType, typename ScalarType>
610 // Set each value of the TriDi matrix to "value".
611 const OrdinalType numvals = (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
612
613 for(OrdinalType i = 0; i<numvals ; ++i)
614 {
615 values_[i] = value_in;
616 }
617 return 0;
618}
619
620// template<typename OrdinalType, typename ScalarType>
621// int SerialTriDiMatrix<OrdinalType, ScalarType>::random()
622// {
623// // Set each value of the TriDi matrix to a random value.
624// for(OrdinalType j = 0; j < numCols_; j++)
625// {
626// for(OrdinalType i = 0; i < numRowsCols_; i++)
627// {
628// values_[i + j*stride_] = ScalarTraits<ScalarType>::random();
629// }
630// }
631// return 0;
632// }
633
634template<typename OrdinalType, typename ScalarType>
637{
638 if(this == &Source)
639 return(*this); // Special case of source same as target
640 if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
641 return(*this); // Special case of both are views to same data.
642
643 // If the source is a view then we will return a view, else we will return a copy.
644 if (!Source.valuesCopied_) {
645 if(valuesCopied_) {
646 // Clean up stored data if this was previously a copy.
647 deleteArrays();
648 }
649 numRowsCols_ = Source.numRowsCols_;
650 values_ = Source.values_;
651 }
652 else {
653 // If we were a view, we will now be a copy.
654 if(!valuesCopied_) {
655 numRowsCols_ = Source.numRowsCols_;
656 const OrdinalType numvals = ( Source.numRowsCols_ == 1) ? 1 : 4*(Source.numRowsCols_ - 1);
657 if(numvals > 0) {
658 values_ = new ScalarType[numvals];
659 valuesCopied_ = true;
660 }
661 else {
662 values_ = 0;
663 }
664 }
665 // If we were a copy, we will stay a copy.
666 else {
667 // we need to allocate more space (or less space)
668 deleteArrays();
669 numRowsCols_ = Source.numRowsCols_;
670 const OrdinalType numvals = ( Source.numRowsCols_ == 1) ? 1 : 4*(Source.numRowsCols_ - 1);
671 if(numvals > 0) {
672 values_ = new ScalarType[numvals];
673 valuesCopied_ = true;
674 }
675 }
676 DL_ = values_;
677 if(values_ != 0) {
678 D_ = DL_ + (numRowsCols_-1);
679 DU_ = D_ + numRowsCols_;
680 DU2_ = DU_ + (numRowsCols_-1);
681
682 OrdinalType min = TEUCHOS_MIN(numRowsCols_, Source.numRowsCols_);
683 for(OrdinalType i = 0 ; i < min ; ++i ) {
684 D_[i] = Source.D()[i];
685 if(i< (min-1 ) )
686 {
687 DL_[i] = Source.DL()[i];
688 DU_[i] = Source.DU()[i];
689 }
690 if(i< (min-2))
691 DU2_[i] = Source.DU2()[i];
692 }
693
694 }
695 else {
696 D_ = DU_ = DU2_ = 0;
697 }
698 }
699 return(*this);
700}
701
702template<typename OrdinalType, typename ScalarType>
704{
705 // Check for compatible dimensions
706 if ((numRowsCols_ != Source.numRowsCols_) )
707 {
708 TEUCHOS_CHK_REF(*this); // Return *this without altering it.
709 }
710 copyMat(Source, 0, ScalarTraits<ScalarType>::one());
711 return(*this);
712}
713
714template<typename OrdinalType, typename ScalarType>
716{
717 // Check for compatible dimensions
718 if ((numRowsCols_ != Source.numRowsCols_) )
719 {
720 TEUCHOS_CHK_REF(*this); // Return *this without altering it.
721 }
722 copyMat(Source, 0, -ScalarTraits<ScalarType>::one());
723 return(*this);
724}
725
726template<typename OrdinalType,typename ScalarType>
728{
729 if(this == &Source)
730 return(*this); // Special case of source same as target
731 if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
732 return(*this); // Special case of both are views to same data.
733
734 // Check for compatible dimensions
735 if ((numRowsCols_ != Source.numRowsCols_) )
736 {
737 TEUCHOS_CHK_REF(*this); // Return *this without altering it.
738 }
739 copyMat(Source, 0, 0);
740 return(*this);
741}
742
743//----------------------------------------------------------------------------------------------------
744// Accessor methods
745//----------------------------------------------------------------------------------------------------
746
747template<typename OrdinalType,typename ScalarType>
748inline const ScalarType& SerialTriDiMatrix<OrdinalType,ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex) const
749{
750 OrdinalType diff = colIndex - rowIndex;
751
752 //#ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
753 checkIndex( rowIndex, colIndex );
754 //#endif
755 switch (diff) {
756 case -1:
757 return DL_[colIndex];
758 case 0:
759 return D_[colIndex];
760 case 1:
761 return DU_[rowIndex];
762 case 2:
763 return DU2_[rowIndex];
764 default:
765 TEUCHOS_TEST_FOR_EXCEPTION(true, std::out_of_range,
766 "SerialTriDiMatrix<T>::operator (row,col) "
767 "Index (" << rowIndex <<","<<colIndex<<") out of range ");
768 }
769}
770
771template<typename OrdinalType,typename ScalarType>
772inline ScalarType& Teuchos::SerialTriDiMatrix<OrdinalType,ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex)
773{
774 OrdinalType diff = colIndex - rowIndex;
775 //#ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
776 checkIndex( rowIndex, colIndex );
777 //#endif
778 switch (diff) {
779 case -1:
780 return DL_[colIndex];
781 case 0:
782 return D_[colIndex];
783 case 1:
784 return DU_[rowIndex];
785 case 2:
786 return DU2_[rowIndex];
787 default:
788 TEUCHOS_TEST_FOR_EXCEPTION(true, std::out_of_range,
789 "SerialTriDiMatrix<T>::operator (row,col) "
790 "Index (" << rowIndex <<","<<colIndex<<") out of range ");
791 }
792}
793
794//----------------------------------------------------------------------------------------------------
795// Norm methods
796//----------------------------------------------------------------------------------------------------
797
798template<typename OrdinalType,typename ScalarType>
800{
801 OrdinalType i, j;
804
805 // Fix this for Tri DI
806
807 for(j = 0; j < numRowsCols_; j++)
808 {
809 ScalarType sum = 0;
810 if(j-1>=0) sum += ScalarTraits<ScalarType>::magnitude((*this)(j-1,j));
811 sum+= ScalarTraits<ScalarType>::magnitude((*this)(j,j));
812 if(j+1<numRowsCols_) sum+= ScalarTraits<ScalarType>::magnitude((*this)(j+1,j));
814 if(absSum > anorm)
815 {
816 anorm = absSum;
817 }
818 }
819 updateFlops(numRowsCols_ * numRowsCols_);
820 return(anorm);
821}
822
823template<typename OrdinalType, typename ScalarType>
825{
826 OrdinalType i,j;
828
829 for (i = 0; i < numRowsCols_; i++) {
831 for (j=i-1; j<= i+1; j++) {
832 if(j >= 0 && j < numRowsCols_) sum += ScalarTraits<ScalarType>::magnitude((*this)(i,j));
833 }
834 anorm = TEUCHOS_MAX( anorm, sum );
835 }
836 updateFlops(numRowsCols_ * numRowsCols_);
837 return(anorm);
838}
839
840template<typename OrdinalType, typename ScalarType>
842{
843 // \todo fix this
844 OrdinalType i, j;
846 for (j = 0; j < numRowsCols_; j++) {
847 for (i = j-1; i <= j+1; i++) {
848 if(i >= 0 && i < numRowsCols_ ) anorm += ScalarTraits<ScalarType>::magnitude((*this)(i,j));
849 }
850 }
852 updateFlops( (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1) );
853 return(anorm);
854}
855
856//----------------------------------------------------------------------------------------------------
857// Comparison methods
858//----------------------------------------------------------------------------------------------------
859
860template<typename OrdinalType, typename ScalarType>
862{
863 bool allmatch = true;
864 // bool result = 1; // unused
865 if((numRowsCols_ != Operand.numRowsCols_) )
866 {
867 // result = 0; // unused
868 }
869 else
870 {
871 OrdinalType numvals = (numRowsCols_ == 1)? 1 : 4*(numRowsCols_ -1 );
872
873 for(OrdinalType i = 0; i< numvals; ++i)
874 allmatch &= (Operand.values_[i] == values_[i]);
875 }
876 return allmatch;
877}
878
879template<typename OrdinalType, typename ScalarType>
881 return !((*this) == Operand);
882}
883
884//----------------------------------------------------------------------------------------------------
885// Multiplication method
886//----------------------------------------------------------------------------------------------------
887
888template<typename OrdinalType, typename ScalarType>
890{
891 this->scale( alpha );
892 return(*this);
893}
894
895template<typename OrdinalType, typename ScalarType>
897{
898 OrdinalType i;
899 OrdinalType numvals = (numRowsCols_ == 1)? 1 : 4*(numRowsCols_ -1 );
900 for (i=0; i < numvals ; ++i ) {
901 values_[i] *= alpha;
902 }
903 return(0);
904}
905
906template<typename OrdinalType, typename ScalarType>
908{
909 OrdinalType j;
910
911 // Check for compatible dimensions
912 if ((numRowsCols_ != A.numRowsCols_) )
913 {
914 TEUCHOS_CHK_ERR(-1); // Return error
915 }
916 OrdinalType numvals = (numRowsCols_ == 1)? 1 : 4*(numRowsCols_ -1 );
917 for (j=0; j<numvals; j++) {
918 values_[j] = A.values_ * values_[j];
919 }
920 updateFlops( numvals );
921 return(0);
922}
923
924template<typename OrdinalType, typename ScalarType>
926{
927 os << std::endl;
928 if(valuesCopied_)
929 os << "A_Copied: yes" << std::endl;
930 else
931 os << "A_Copied: no" << std::endl;
932 os << "Rows and Columns: " << numRowsCols_ << std::endl;
933 if(numRowsCols_ == 0) {
934 os << "(matrix is empty, no values to display)" << std::endl;
935 return;
936 }
937 else
938 {
939 os << "DL: "<<std::endl;
940 for(int i=0;i<numRowsCols_-1;++i)
941 os << DL_[i]<<" ";
942 os << std::endl;
943 os << "D: "<<std::endl;
944 for(int i=0;i<numRowsCols_;++i)
945 os << D_[i]<<" ";
946 os << std::endl;
947 os << "DU: "<<std::endl;
948 for(int i=0;i<numRowsCols_-1;++i)
949 os << DU_[i]<<" ";
950 os << std::endl;
951 os << "DU2: "<<std::endl;
952 for(int i=0;i<numRowsCols_-2;++i)
953 os << DU2_[i]<<" ";
954 os << std::endl;
955 }
956
957 os <<" square format:"<<std::endl;
958 for(int i=0 ; i < numRowsCols_ ; ++i ) {
959 for(int j=0;j<numRowsCols_;++j) {
960 if ( j >= i-1 && j <= i+1) {
961 os << (*this)(i,j)<<" ";
962 }
963 else
964 os <<". ";
965 }
966 os << std::endl;
967 }
968
969}
970
971//----------------------------------------------------------------------------------------------------
972// Protected methods
973//----------------------------------------------------------------------------------------------------
974
975template<typename OrdinalType, typename ScalarType>
976inline void SerialTriDiMatrix<OrdinalType, ScalarType>::checkIndex( OrdinalType rowIndex, OrdinalType colIndex ) const
977{
978 OrdinalType diff = colIndex - rowIndex;
979 TEUCHOS_TEST_FOR_EXCEPTION(rowIndex < 0 || rowIndex >= numRowsCols_, std::out_of_range,
980 "SerialTriDiMatrix<T>::checkIndex: "
981 "Row index " << rowIndex << " out of range [0, "<< numRowsCols_ << "]");
982 TEUCHOS_TEST_FOR_EXCEPTION(colIndex < 0 || colIndex >= numRowsCols_,
983 std::out_of_range,
984 "SerialTriDiMatrix<T>::checkIndex: "
985 "Col index " << colIndex << " out of range [0, "<< numRowsCols_ << "]");
986 TEUCHOS_TEST_FOR_EXCEPTION(diff > 2 || diff < -1 , std::out_of_range,
987 "SerialTriDiMatrix<T>::checkIndex: "
988 "index difference " << diff << " out of range [-1, 2]");
989}
990
991template<typename OrdinalType, typename ScalarType>
993{
994 if (valuesCopied_)
995 {
996 delete [] values_;
997 values_ = 0;
998 valuesCopied_ = false;
999 }
1000}
1001
1002template<typename OrdinalType, typename ScalarType>
1004 OrdinalType startRowCol,
1005 ScalarType alpha )
1006{
1007 OrdinalType i;
1008 OrdinalType max = inputMatrix.numRowsCols_;
1009 if(max > numRowsCols_ ) max = numRowsCols_;
1010 if(startRowCol > max ) return; //
1011
1012 for(i = startRowCol ; i < max ; ++i) {
1013
1015 // diagonal
1016 D()[i] += inputMatrix.D()[i];
1017 if(i<(max-1) && (i-1) >= startRowCol) {
1018 DL()[i] += inputMatrix.DL()[i];
1019 DU()[i] += inputMatrix.DU()[i];
1020 }
1021 if(i<(max-2) && (i-2) >= startRowCol) {
1022 DU2()[i] += inputMatrix.DU2()[i];
1023 }
1024 }
1025 else {
1026 // diagonal
1027 D()[i] = inputMatrix.D()[i];
1028 if(i<(max-1) && (i-1) >= startRowCol) {
1029 DL()[i] = inputMatrix.DL()[i];
1030 DU()[i] = inputMatrix.DU()[i];
1031 }
1032 if(i<(max-2) && (i-2) >= startRowCol) {
1033 DU2()[i] = inputMatrix.DU2()[i];
1034 }
1035 }
1036 }
1037}
1038
1039}
1040
1041
1042#endif /* _TEUCHOS_SERIALTRIDIMATRIX_HPP_ */
Templated interface class to BLAS routines.
Object for storing data and providing functionality that is common to all computational classes.
Teuchos header file which uses auto-configuration information to include necessary C++ headers.
#define TEUCHOS_CHK_REF(a)
#define TEUCHOS_MIN(x, y)
#define TEUCHOS_MAX(x, y)
#define TEUCHOS_CHK_ERR(a)
Teuchos::DataAccess Mode enumerable type.
Defines basic traits for the scalar field type.
Templated BLAS wrapper.
Functionality and data that is common to all computational classes.
The base Teuchos class.
This class creates and provides basic support for TriDi matrix of templated type.
void checkIndex(OrdinalType rowIndex, OrdinalType colIndex=0) const
int reshape(OrdinalType numRowsCols)
Reshaping method for changing the size of a SerialTriDiMatrix, keeping the entries.
virtual void print(std::ostream &os) const
Print method. Defines the behavior of the std::ostream << operator inherited from the Object class.
ScalarType scalarType
Typedef for scalar type.
bool operator==(const SerialTriDiMatrix< OrdinalType, ScalarType > &Operand) const
Equality of two matrices.
void copyMat(SerialTriDiMatrix< OrdinalType, ScalarType > matrix, OrdinalType startCol, ScalarType alpha=ScalarTraits< ScalarType >::zero())
SerialTriDiMatrix< OrdinalType, ScalarType > & assign(const SerialTriDiMatrix< OrdinalType, ScalarType > &Source)
Copies values from one matrix to another.
ScalarTraits< ScalarType >::magnitudeType normOne() const
Returns the 1-norm of the matrix.
int shape(OrdinalType numRows)
Shape method for changing the size of a SerialTriDiMatrix, initializing entries to zero.
int scale(const ScalarType alpha)
Scale this matrix by alpha; *this = alpha**this.
SerialTriDiMatrix< OrdinalType, ScalarType > & operator+=(const SerialTriDiMatrix< OrdinalType, ScalarType > &Source)
Add another matrix to this matrix.
SerialTriDiMatrix< OrdinalType, ScalarType > & operator-=(const SerialTriDiMatrix< OrdinalType, ScalarType > &Source)
Subtract another matrix from this matrix.
bool empty() const
Returns the column dimension of this matrix.
OrdinalType numRowsCols() const
Returns the row dimension of this matrix.
ScalarTraits< ScalarType >::magnitudeType normInf() const
Returns the Infinity-norm of the matrix.
int shapeUninitialized(OrdinalType numRows)
Same as shape() except leaves uninitialized.
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero())
Set all values in the matrix to a constant value.
SerialTriDiMatrix< OrdinalType, ScalarType > & operator*=(const ScalarType alpha)
Scale this matrix by alpha; *this = alpha**this.
SerialTriDiMatrix< OrdinalType, ScalarType > & operator=(const SerialTriDiMatrix< OrdinalType, ScalarType > &Source)
Copies values from one matrix to another.
ScalarType & operator()(OrdinalType rowIndex, OrdinalType colIndex)
Element access method (non-const).
OrdinalType ordinalType
Typedef for ordinal type.
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
Returns the Frobenius-norm of the matrix.
bool operator!=(const SerialTriDiMatrix< OrdinalType, ScalarType > &Operand) const
Inequality of two matrices.
ScalarType * values() const
Column access method (non-const).
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Macro for throwing an exception with breakpointing to ease debugging.
Definition: PackageA.cpp:3
This structure defines some basic traits for a scalar field type.
static magnitudeType magnitude(T a)
Returns the magnitudeType of the scalar type a.
T magnitudeType
Mandatory typedef for result of magnitude.
static T zero()
Returns representation of zero for this scalar type.
static T conjugate(T a)
Returns the conjugate of the scalar type a.