Teuchos - Trilinos Tools Package  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Teuchos_SerialTriDiMatrix.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_SERIALTRIDIMATRIX_HPP_
11 #define _TEUCHOS_SERIALTRIDIMATRIX_HPP_
12 
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_Assert.hpp"
22 
31 namespace Teuchos {
32 
33 template<typename OrdinalType, typename ScalarType>
34 class SerialTriDiMatrix : public CompObject, public Object, public BLAS<OrdinalType, ScalarType > {
35 public:
36 
38  typedef OrdinalType ordinalType;
40  typedef ScalarType scalarType;
41 
43 
44 
46 
50 
52 
60  SerialTriDiMatrix(OrdinalType numRows, OrdinalType numCols, bool zeroOut = true);
61 
63 
68  SerialTriDiMatrix(DataAccess CV, ScalarType* values, OrdinalType numRowsCols);
69 
71 
77 
79 
90  SerialTriDiMatrix(DataAccess CV, const SerialTriDiMatrix<OrdinalType, ScalarType> &Source, OrdinalType numRowsCols, OrdinalType startRowCols=0);
91 
93  virtual ~SerialTriDiMatrix();
95 
97 
98 
108  int shape(OrdinalType numRows);
109 
111  int shapeUninitialized(OrdinalType numRows);
112 
114 
123  int reshape(OrdinalType numRowsCols);
124 
126 
128 
129 
131 
138 
140 
146 
148 
151  SerialTriDiMatrix<OrdinalType, ScalarType>& operator= (const ScalarType value) { putScalar(value); return(*this); }
152 
154 
158  int putScalar( const ScalarType value = Teuchos::ScalarTraits<ScalarType>::zero() );
159 
161  // int random();
162 
164 
166 
167 
169 
176  ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex);
177 
179 
186  const ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex) const;
187 
189 
196  // ScalarType* operator [] (OrdinalType colIndex);
197 
199 
206  // const ScalarType* operator [] (OrdinalType colIndex) const;
207 
209 
210  ScalarType* values() const { return(values_); }
211 
212  ScalarType* D() const { return D_;}
213  ScalarType* DL() const { return DL_;}
214  ScalarType* DU() const { return DU_;}
215  ScalarType* DU2() const { return DU2_;}
216 
218 
220 
221 
223 
226  SerialTriDiMatrix<OrdinalType, ScalarType>& operator+= (const SerialTriDiMatrix<OrdinalType, ScalarType>& Source);
227 
229 
232  SerialTriDiMatrix<OrdinalType, ScalarType>& operator-= (const SerialTriDiMatrix<OrdinalType, ScalarType>& Source);
233 
235 
238  SerialTriDiMatrix<OrdinalType, ScalarType>& operator*= (const ScalarType alpha);
239 
241 
245  int scale ( const ScalarType alpha );
246 
248 
254  int scale ( const SerialTriDiMatrix<OrdinalType, ScalarType>& A );
255 
257 
271  //int multiply (ETransp transa, ETransp transb, ScalarType alpha, const SerialTriDiMatrix<OrdinalType, ScalarType> &A, const SerialTriDiMatrix<OrdinalType, ScalarType> &B, ScalarType beta);
272 
274 
285  //int multiply (ESide sideA, ScalarType alpha, const SerialSymTriDiMatrix<OrdinalType, ScalarType> &A, const SerialTriDiMatrix<OrdinalType, ScalarType> &B, ScalarType beta);
286 
288 
290 
291 
293 
296  bool operator== (const SerialTriDiMatrix<OrdinalType, ScalarType> &Operand) const;
297 
299 
302  bool operator!= (const SerialTriDiMatrix<OrdinalType, ScalarType> &Operand) const;
303 
305 
307 
308 
310  OrdinalType numRowsCols() const { return(numRowsCols_); }
311 
313  // OrdinalType numCols() const { return(numRowsCols_); }
314 
316  // OrdinalType stride() const { return(stride_); }
317 
319  bool empty() const { return(numRowsCols_ == 0); }
321 
323 
324 
327 
330 
334 
336 
337  virtual void print(std::ostream& os) const;
339 
341 protected:
342  void copyMat(SerialTriDiMatrix<OrdinalType, ScalarType> matrix,
343  OrdinalType startCol,
344  ScalarType alpha = ScalarTraits<ScalarType>::zero() );
345  void deleteArrays();
346  void checkIndex( OrdinalType rowIndex, OrdinalType colIndex = 0 ) const;
347  OrdinalType numRowsCols_;
348 
349  bool valuesCopied_;
350  ScalarType* values_;
351  ScalarType* DL_;
352  ScalarType* D_;
353  ScalarType* DU_;
354  ScalarType* DU2_;
355 
356 }; // class Teuchos_SerialTriDiMatrix
357 
358 //----------------------------------------------------------------------------------------------------
359 // Constructors and Destructor
360 //----------------------------------------------------------------------------------------------------
361 
362 template<typename OrdinalType, typename ScalarType>
364  :
365  CompObject(),
366  numRowsCols_(0),
367  valuesCopied_(false),
368  values_(0),
369  DL_(NULL),
370  D_(NULL),
371  DU_(NULL),
372  DU2_(NULL)
373 {}
374 
375 template<typename OrdinalType, typename ScalarType>
376 SerialTriDiMatrix<OrdinalType, ScalarType>::SerialTriDiMatrix( OrdinalType numRowsCols_in, OrdinalType /* numCols_in */, bool zeroOut)
377  : CompObject(), numRowsCols_(numRowsCols_in) {
378 
379  OrdinalType numvals = (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
380  values_ = new ScalarType [numvals];
381  DL_ = values_;
382  D_ = DL_ + (numRowsCols_-1);
383  DU_ = D_ + numRowsCols_;
384  DU2_ = DU_ + (numRowsCols_-1);
385 
386  valuesCopied_ = true;
387  if (zeroOut == true)
388  putScalar();
389 }
390 
391 template<typename OrdinalType, typename ScalarType>
392 SerialTriDiMatrix<OrdinalType, ScalarType>::SerialTriDiMatrix(DataAccess CV, ScalarType* values_in, OrdinalType numRowsCols_in )
393  : CompObject(), numRowsCols_(numRowsCols_in),
394  valuesCopied_(false), values_(values_in)
395 {
396  const OrdinalType numvals = (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
397  if(CV == Copy) {
398  values_ = new ScalarType[numvals];
399  valuesCopied_ = true;
400  }
401  else //CV == View
402  {
403  values_ = values_in;
404  valuesCopied_ = false;
405  }
406  DL_ = values_;
407  D_ = DL_ + (numRowsCols_-1);
408  DU_ = D_ + numRowsCols_;
409  DU2_ = DU_ + (numRowsCols_-1);
410  if(CV == Copy) {
411  for(OrdinalType i = 0 ; i < numRowsCols_ ; ++i )
412  values_[i] = values_in[i];
413  }
414 }
415 
416 template<typename OrdinalType, typename ScalarType>
417 SerialTriDiMatrix<OrdinalType, ScalarType>::SerialTriDiMatrix(const SerialTriDiMatrix<OrdinalType, ScalarType> &Source, ETransp trans) : CompObject(), BLAS<OrdinalType,ScalarType>(), numRowsCols_(0), valuesCopied_(true), values_(0)
418 {
419  if ( trans == Teuchos::NO_TRANS ) {
420  numRowsCols_ = Source.numRowsCols_;
421 
422  const OrdinalType numvals = (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
423  values_ = new ScalarType[numvals];
424  DL_ = values_;
425  D_ = DL_+ (numRowsCols_-1);
426  DU_ = D_ + numRowsCols_;
427  DU2_ = DU_ + (numRowsCols_-1);
428 
429  copyMat(Source, 0, 0);
430  }
432  {
433  numRowsCols_ = Source.numRowsCols_;
434  const OrdinalType numvals = (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
435  values_ = new ScalarType[numvals];
436  DL_ = values_;
437  D_ = DL_+(numRowsCols_-1);
438  DU_ = D_ + numRowsCols_;
439  DU2_ = DU_ + (numRowsCols_-1);
440 
441  OrdinalType min = numRowsCols_;
442  if(min > Source.numRowsCols_) min = Source.numRowsCols_;
443 
444  for(OrdinalType i = 0 ; i< min ; ++i) {
445  D_[i] = Teuchos::ScalarTraits<ScalarType>::conjugate(Source.D_[i]);
446  if(i < (min-1)) {
447  DL_[i] = Teuchos::ScalarTraits<ScalarType>::conjugate(Source.DL_[i]);
448  DU_[i] = Teuchos::ScalarTraits<ScalarType>::conjugate(Source.DU_[i]);
449  }
450  if(i < (min-2)) {
451  DU2_[i] = Teuchos::ScalarTraits<ScalarType>::conjugate(Source.DU2_[i]);
452  }
453  }
454  }
455  else
456  {
457  numRowsCols_ = Source.numRowsCols_;
458  const OrdinalType numvals = (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
459  values_ = new ScalarType[numvals];
460  OrdinalType min = numRowsCols_;
461  if(min > Source.numRowsCols_) min = Source.numRowsCols_;
462  for(OrdinalType i = 0 ; i< min ; ++i) {
463  D_[i] = Source.D_[i];
464  if(i < (min-1)) {
465  DL_[i] = Source.DL_[i];
466  DU_[i] = Source.DU_[i];
467  }
468  if(i < (min-2)) {
469  DU2_[i] = Source.DU2_[i];
470  }
471  }
472  }
473 }
474 
475 template<typename OrdinalType, typename ScalarType>
478  OrdinalType numRowsCols_in, OrdinalType startRow )
479  : CompObject(), numRowsCols_(numRowsCols_in),
480  valuesCopied_(false), values_(Source.values_) {
481  if(CV == Copy)
482  {
483  const OrdinalType numvals = (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
484  values_ = new ScalarType[numvals];
485  copyMat(Source, startRow);
486  valuesCopied_ = true;
487  }
488  else // CV == View
489  {
490  // \todo ???
491  // values_ = values_ + (stride_ * startCol) + startRow;
492  }
493 }
494 
495 template<typename OrdinalType, typename ScalarType>
497 {
498  deleteArrays();
499 }
500 
501 //----------------------------------------------------------------------------------------------------
502 // Shape methods
503 //----------------------------------------------------------------------------------------------------
504 
505 template<typename OrdinalType, typename ScalarType>
506 int SerialTriDiMatrix<OrdinalType, ScalarType>::shape( OrdinalType numRowsCols_in )
507 {
508  deleteArrays(); // Get rid of anything that might be already allocated
509  numRowsCols_ = numRowsCols_in;
510  const OrdinalType numvals = ( numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
511  values_ = new ScalarType[numvals];
512 
513  putScalar();
514  valuesCopied_ = true;
515  return(0);
516 }
517 
518 template<typename OrdinalType, typename ScalarType>
520 {
521  deleteArrays(); // Get rid of anything that might be already allocated
522  numRowsCols_ = numRowsCols_in;
523  const OrdinalType numvals = ( numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
524  values_ = new ScalarType[numvals];
525  valuesCopied_ = true;
526  return(0);
527 }
528 
529 template<typename OrdinalType, typename ScalarType>
531 {
532  if(numRowsCols_in <1 ) {
533  deleteArrays();
534  return 0;
535  }
536  // Allocate space for new matrix
537  const OrdinalType numvals = ( numRowsCols_in == 1) ? 1 : 4*(numRowsCols_in - 1);
538  ScalarType* values_tmp = new ScalarType[numvals];
539  ScalarType zero = ScalarTraits<ScalarType>::zero();
540  for(OrdinalType i= 0; i<numvals ; ++i)
541  values_tmp[i] = zero;
542 
543  OrdinalType min = TEUCHOS_MIN(numRowsCols_, numRowsCols_in);
544  ScalarType* dl = values_tmp;
545  ScalarType* d = values_tmp + (numRowsCols_in-1);
546  ScalarType* du = d+(numRowsCols_in);
547  ScalarType* du2 = du+(numRowsCols_in - 1);
548 
549  if(values_ != 0) {
550  for(OrdinalType i = 0 ; i< min ; ++i) {
551  dl[i] = DL_[i];
552  d[i] = D_[i];
553  du[i] = DU_[i];
554  du2[i] = DU2_[i];
555  }
556  }
557 
558  deleteArrays(); // Get rid of anything that might be already allocated
559  numRowsCols_ = numRowsCols_in;
560 
561  values_ = values_tmp; // Set pointer to new A
562  DL_ = dl;
563  D_ = d;
564  DU_ = du;
565  DU2_ = du2;
566 
567  valuesCopied_ = true;
568  return(0);
569 }
570 
571 //----------------------------------------------------------------------------------------------------
572 // Set methods
573 //----------------------------------------------------------------------------------------------------
574 
575 template<typename OrdinalType, typename ScalarType>
577  // Set each value of the TriDi matrix to "value".
578  const OrdinalType numvals = (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
579 
580  for(OrdinalType i = 0; i<numvals ; ++i)
581  {
582  values_[i] = value_in;
583  }
584  return 0;
585 }
586 
587 // template<typename OrdinalType, typename ScalarType>
588 // int SerialTriDiMatrix<OrdinalType, ScalarType>::random()
589 // {
590 // // Set each value of the TriDi matrix to a random value.
591 // for(OrdinalType j = 0; j < numCols_; j++)
592 // {
593 // for(OrdinalType i = 0; i < numRowsCols_; i++)
594 // {
595 // values_[i + j*stride_] = ScalarTraits<ScalarType>::random();
596 // }
597 // }
598 // return 0;
599 // }
600 
601 template<typename OrdinalType, typename ScalarType>
604 {
605  if(this == &Source)
606  return(*this); // Special case of source same as target
607  if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
608  return(*this); // Special case of both are views to same data.
609 
610  // If the source is a view then we will return a view, else we will return a copy.
611  if (!Source.valuesCopied_) {
612  if(valuesCopied_) {
613  // Clean up stored data if this was previously a copy.
614  deleteArrays();
615  }
616  numRowsCols_ = Source.numRowsCols_;
617  values_ = Source.values_;
618  }
619  else {
620  // If we were a view, we will now be a copy.
621  if(!valuesCopied_) {
622  numRowsCols_ = Source.numRowsCols_;
623  const OrdinalType numvals = ( Source.numRowsCols_ == 1) ? 1 : 4*(Source.numRowsCols_ - 1);
624  if(numvals > 0) {
625  values_ = new ScalarType[numvals];
626  valuesCopied_ = true;
627  }
628  else {
629  values_ = 0;
630  }
631  }
632  // If we were a copy, we will stay a copy.
633  else {
634  // we need to allocate more space (or less space)
635  deleteArrays();
636  numRowsCols_ = Source.numRowsCols_;
637  const OrdinalType numvals = ( Source.numRowsCols_ == 1) ? 1 : 4*(Source.numRowsCols_ - 1);
638  if(numvals > 0) {
639  values_ = new ScalarType[numvals];
640  valuesCopied_ = true;
641  }
642  }
643  DL_ = values_;
644  if(values_ != 0) {
645  D_ = DL_ + (numRowsCols_-1);
646  DU_ = D_ + numRowsCols_;
647  DU2_ = DU_ + (numRowsCols_-1);
648 
649  OrdinalType min = TEUCHOS_MIN(numRowsCols_, Source.numRowsCols_);
650  for(OrdinalType i = 0 ; i < min ; ++i ) {
651  D_[i] = Source.D()[i];
652  if(i< (min-1 ) )
653  {
654  DL_[i] = Source.DL()[i];
655  DU_[i] = Source.DU()[i];
656  }
657  if(i< (min-2))
658  DU2_[i] = Source.DU2()[i];
659  }
660 
661  }
662  else {
663  D_ = DU_ = DU2_ = 0;
664  }
665  }
666  return(*this);
667 }
668 
669 template<typename OrdinalType, typename ScalarType>
671 {
672  // Check for compatible dimensions
673  if ((numRowsCols_ != Source.numRowsCols_) )
674  {
675  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
676  }
677  copyMat(Source, 0, ScalarTraits<ScalarType>::one());
678  return(*this);
679 }
680 
681 template<typename OrdinalType, typename ScalarType>
683 {
684  // Check for compatible dimensions
685  if ((numRowsCols_ != Source.numRowsCols_) )
686  {
687  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
688  }
689  copyMat(Source, 0, -ScalarTraits<ScalarType>::one());
690  return(*this);
691 }
692 
693 template<typename OrdinalType,typename ScalarType>
695 {
696  if(this == &Source)
697  return(*this); // Special case of source same as target
698  if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
699  return(*this); // Special case of both are views to same data.
700 
701  // Check for compatible dimensions
702  if ((numRowsCols_ != Source.numRowsCols_) )
703  {
704  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
705  }
706  copyMat(Source, 0, 0);
707  return(*this);
708 }
709 
710 //----------------------------------------------------------------------------------------------------
711 // Accessor methods
712 //----------------------------------------------------------------------------------------------------
713 
714 template<typename OrdinalType,typename ScalarType>
715 inline const ScalarType& SerialTriDiMatrix<OrdinalType,ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex) const
716 {
717  OrdinalType diff = colIndex - rowIndex;
718 
719  //#ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
720  checkIndex( rowIndex, colIndex );
721  //#endif
722  switch (diff) {
723  case -1:
724  return DL_[colIndex];
725  case 0:
726  return D_[colIndex];
727  case 1:
728  return DU_[rowIndex];
729  case 2:
730  return DU2_[rowIndex];
731  default:
732  TEUCHOS_TEST_FOR_EXCEPTION(true, std::out_of_range,
733  "SerialTriDiMatrix<T>::operator (row,col) "
734  "Index (" << rowIndex <<","<<colIndex<<") out of range ");
735  }
736 }
737 
738 template<typename OrdinalType,typename ScalarType>
739 inline ScalarType& Teuchos::SerialTriDiMatrix<OrdinalType,ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex)
740 {
741  OrdinalType diff = colIndex - rowIndex;
742  //#ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
743  checkIndex( rowIndex, colIndex );
744  //#endif
745  switch (diff) {
746  case -1:
747  return DL_[colIndex];
748  case 0:
749  return D_[colIndex];
750  case 1:
751  return DU_[rowIndex];
752  case 2:
753  return DU2_[rowIndex];
754  default:
755  TEUCHOS_TEST_FOR_EXCEPTION(true, std::out_of_range,
756  "SerialTriDiMatrix<T>::operator (row,col) "
757  "Index (" << rowIndex <<","<<colIndex<<") out of range ");
758  }
759 }
760 
761 //----------------------------------------------------------------------------------------------------
762 // Norm methods
763 //----------------------------------------------------------------------------------------------------
764 
765 template<typename OrdinalType,typename ScalarType>
767 {
768  OrdinalType i, j;
771 
772  // Fix this for Tri DI
773 
774  for(j = 0; j < numRowsCols_; j++)
775  {
776  ScalarType sum = 0;
777  if(j-1>=0) sum += ScalarTraits<ScalarType>::magnitude((*this)(j-1,j));
778  sum+= ScalarTraits<ScalarType>::magnitude((*this)(j,j));
779  if(j+1<numRowsCols_) sum+= ScalarTraits<ScalarType>::magnitude((*this)(j+1,j));
781  if(absSum > anorm)
782  {
783  anorm = absSum;
784  }
785  }
786  updateFlops(numRowsCols_ * numRowsCols_);
787  return(anorm);
788 }
789 
790 template<typename OrdinalType, typename ScalarType>
792 {
793  OrdinalType i,j;
795 
796  for (i = 0; i < numRowsCols_; i++) {
798  for (j=i-1; j<= i+1; j++) {
799  if(j >= 0 && j < numRowsCols_) sum += ScalarTraits<ScalarType>::magnitude((*this)(i,j));
800  }
801  anorm = TEUCHOS_MAX( anorm, sum );
802  }
803  updateFlops(numRowsCols_ * numRowsCols_);
804  return(anorm);
805 }
806 
807 template<typename OrdinalType, typename ScalarType>
809 {
810  // \todo fix this
811  OrdinalType i, j;
813  for (j = 0; j < numRowsCols_; j++) {
814  for (i = j-1; i <= j+1; i++) {
815  if(i >= 0 && i < numRowsCols_ ) anorm += ScalarTraits<ScalarType>::magnitude((*this)(i,j));
816  }
817  }
819  updateFlops( (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1) );
820  return(anorm);
821 }
822 
823 //----------------------------------------------------------------------------------------------------
824 // Comparison methods
825 //----------------------------------------------------------------------------------------------------
826 
827 template<typename OrdinalType, typename ScalarType>
829 {
830  bool allmatch = true;
831  // bool result = 1; // unused
832  if((numRowsCols_ != Operand.numRowsCols_) )
833  {
834  // result = 0; // unused
835  }
836  else
837  {
838  OrdinalType numvals = (numRowsCols_ == 1)? 1 : 4*(numRowsCols_ -1 );
839 
840  for(OrdinalType i = 0; i< numvals; ++i)
841  allmatch &= (Operand.values_[i] == values_[i]);
842  }
843  return allmatch;
844 }
845 
846 template<typename OrdinalType, typename ScalarType>
848  return !((*this) == Operand);
849 }
850 
851 //----------------------------------------------------------------------------------------------------
852 // Multiplication method
853 //----------------------------------------------------------------------------------------------------
854 
855 template<typename OrdinalType, typename ScalarType>
857 {
858  this->scale( alpha );
859  return(*this);
860 }
861 
862 template<typename OrdinalType, typename ScalarType>
864 {
865  OrdinalType i;
866  OrdinalType numvals = (numRowsCols_ == 1)? 1 : 4*(numRowsCols_ -1 );
867  for (i=0; i < numvals ; ++i ) {
868  values_[i] *= alpha;
869  }
870  return(0);
871 }
872 
873 template<typename OrdinalType, typename ScalarType>
875 {
876  OrdinalType j;
877 
878  // Check for compatible dimensions
879  if ((numRowsCols_ != A.numRowsCols_) )
880  {
881  TEUCHOS_CHK_ERR(-1); // Return error
882  }
883  OrdinalType numvals = (numRowsCols_ == 1)? 1 : 4*(numRowsCols_ -1 );
884  for (j=0; j<numvals; j++) {
885  values_[j] = A.values_ * values_[j];
886  }
887  updateFlops( numvals );
888  return(0);
889 }
890 
891 template<typename OrdinalType, typename ScalarType>
893 {
894  os << std::endl;
895  if(valuesCopied_)
896  os << "A_Copied: yes" << std::endl;
897  else
898  os << "A_Copied: no" << std::endl;
899  os << "Rows and Columns: " << numRowsCols_ << std::endl;
900  if(numRowsCols_ == 0) {
901  os << "(matrix is empty, no values to display)" << std::endl;
902  return;
903  }
904  else
905  {
906  os << "DL: "<<std::endl;
907  for(int i=0;i<numRowsCols_-1;++i)
908  os << DL_[i]<<" ";
909  os << std::endl;
910  os << "D: "<<std::endl;
911  for(int i=0;i<numRowsCols_;++i)
912  os << D_[i]<<" ";
913  os << std::endl;
914  os << "DU: "<<std::endl;
915  for(int i=0;i<numRowsCols_-1;++i)
916  os << DU_[i]<<" ";
917  os << std::endl;
918  os << "DU2: "<<std::endl;
919  for(int i=0;i<numRowsCols_-2;++i)
920  os << DU2_[i]<<" ";
921  os << std::endl;
922  }
923 
924  os <<" square format:"<<std::endl;
925  for(int i=0 ; i < numRowsCols_ ; ++i ) {
926  for(int j=0;j<numRowsCols_;++j) {
927  if ( j >= i-1 && j <= i+1) {
928  os << (*this)(i,j)<<" ";
929  }
930  else
931  os <<". ";
932  }
933  os << std::endl;
934  }
935 
936 }
937 
938 //----------------------------------------------------------------------------------------------------
939 // Protected methods
940 //----------------------------------------------------------------------------------------------------
941 
942 template<typename OrdinalType, typename ScalarType>
943 inline void SerialTriDiMatrix<OrdinalType, ScalarType>::checkIndex( OrdinalType rowIndex, OrdinalType colIndex ) const
944 {
945  OrdinalType diff = colIndex - rowIndex;
946  TEUCHOS_TEST_FOR_EXCEPTION(rowIndex < 0 || rowIndex >= numRowsCols_, std::out_of_range,
947  "SerialTriDiMatrix<T>::checkIndex: "
948  "Row index " << rowIndex << " out of range [0, "<< numRowsCols_ << "]");
949  TEUCHOS_TEST_FOR_EXCEPTION(colIndex < 0 || colIndex >= numRowsCols_,
950  std::out_of_range,
951  "SerialTriDiMatrix<T>::checkIndex: "
952  "Col index " << colIndex << " out of range [0, "<< numRowsCols_ << "]");
953  TEUCHOS_TEST_FOR_EXCEPTION(diff > 2 || diff < -1 , std::out_of_range,
954  "SerialTriDiMatrix<T>::checkIndex: "
955  "index difference " << diff << " out of range [-1, 2]");
956 }
957 
958 template<typename OrdinalType, typename ScalarType>
959 void SerialTriDiMatrix<OrdinalType, ScalarType>::deleteArrays(void)
960 {
961  if (valuesCopied_)
962  {
963  delete [] values_;
964  values_ = 0;
965  valuesCopied_ = false;
966  }
967 }
968 
969 template<typename OrdinalType, typename ScalarType>
970 void SerialTriDiMatrix<OrdinalType, ScalarType>::copyMat(SerialTriDiMatrix<OrdinalType, ScalarType> inputMatrix,
971  OrdinalType startRowCol,
972  ScalarType alpha )
973 {
974  OrdinalType i;
975  OrdinalType max = inputMatrix.numRowsCols_;
976  if(max > numRowsCols_ ) max = numRowsCols_;
977  if(startRowCol > max ) return; //
978 
979  for(i = startRowCol ; i < max ; ++i) {
980 
981  if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
982  // diagonal
983  D()[i] += inputMatrix.D()[i];
984  if(i<(max-1) && (i-1) >= startRowCol) {
985  DL()[i] += inputMatrix.DL()[i];
986  DU()[i] += inputMatrix.DU()[i];
987  }
988  if(i<(max-2) && (i-2) >= startRowCol) {
989  DU2()[i] += inputMatrix.DU2()[i];
990  }
991  }
992  else {
993  // diagonal
994  D()[i] = inputMatrix.D()[i];
995  if(i<(max-1) && (i-1) >= startRowCol) {
996  DL()[i] = inputMatrix.DL()[i];
997  DU()[i] = inputMatrix.DU()[i];
998  }
999  if(i<(max-2) && (i-2) >= startRowCol) {
1000  DU2()[i] = inputMatrix.DU2()[i];
1001  }
1002  }
1003  }
1004 }
1005 
1006 }
1007 
1008 
1009 #endif /* _TEUCHOS_SERIALTRIDIMATRIX_HPP_ */
ScalarType & operator()(OrdinalType rowIndex, OrdinalType colIndex)
Element access method (non-const).
ScalarTraits< ScalarType >::magnitudeType normInf() const
Returns the Infinity-norm of the matrix.
int shape(OrdinalType numRows)
Shape method for changing the size of a SerialTriDiMatrix, initializing entries to zero...
OrdinalType numRowsCols() const
Returns the row dimension of this matrix.
ScalarType scalarType
Typedef for scalar type.
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.
virtual void print(std::ostream &os) const
Print method. Defines the behavior of the std::ostream &lt;&lt; operator inherited from the Object class...
Teuchos header file which uses auto-configuration information to include necessary C++ headers...
SerialTriDiMatrix< OrdinalType, ScalarType > & operator=(const SerialTriDiMatrix< OrdinalType, ScalarType > &Source)
Copies values from one matrix to another.
Templated BLAS wrapper.
SerialTriDiMatrix< OrdinalType, ScalarType > & assign(const SerialTriDiMatrix< OrdinalType, ScalarType > &Source)
Copies values from one matrix to another.
Object for storing data and providing functionality that is common to all computational classes...
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
Returns the Frobenius-norm of the matrix.
This structure defines some basic traits for a scalar field type.
This class creates and provides basic support for TriDi matrix of templated type. ...
ScalarType * values() const
Column access method (non-const).
Functionality and data that is common to all computational classes.
static T conjugate(T a)
Returns the conjugate of the scalar type a.
ScalarTraits< ScalarType >::magnitudeType normOne() const
Returns the 1-norm of the matrix.
The base Teuchos class.
int scale(const ScalarType alpha)
Scale this matrix by alpha; *this = alpha**this.
bool operator==(const SerialTriDiMatrix< OrdinalType, ScalarType > &Operand) const
Equality of two matrices.
SerialTriDiMatrix< OrdinalType, ScalarType > & operator*=(const ScalarType alpha)
Scale this matrix by alpha; *this = alpha**this.
bool empty() const
Returns the column dimension of this matrix.
static magnitudeType magnitude(T a)
Returns the magnitudeType of the scalar type a.
SerialTriDiMatrix()
Default Constructor.
SerialTriDiMatrix< OrdinalType, ScalarType > & operator+=(const SerialTriDiMatrix< OrdinalType, ScalarType > &Source)
Add another matrix to this matrix.
OrdinalType ordinalType
Typedef for ordinal type.
Defines basic traits for the scalar field type.
static T zero()
Returns representation of zero for this scalar type.
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero())
Set all values in the matrix to a constant value.
int reshape(OrdinalType numRowsCols)
Reshaping method for changing the size of a SerialTriDiMatrix, keeping the entries.
bool operator!=(const SerialTriDiMatrix< OrdinalType, ScalarType > &Operand) const
Inequality of two matrices.
SerialTriDiMatrix< OrdinalType, ScalarType > & operator-=(const SerialTriDiMatrix< OrdinalType, ScalarType > &Source)
Subtract another matrix from this matrix.
int shapeUninitialized(OrdinalType numRows)
Same as shape() except leaves uninitialized.
Teuchos::DataAccess Mode enumerable type.