Teuchos Package Browser (Single Doxygen Collection)  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 //
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_
45 
49 #include "Teuchos_CompObject.hpp"
50 #include "Teuchos_BLAS.hpp"
51 #include "Teuchos_ScalarTraits.hpp"
52 #include "Teuchos_DataAccess.hpp"
53 #include "Teuchos_ConfigDefs.hpp"
54 #include "Teuchos_Assert.hpp"
55 
64 namespace Teuchos {
65 
66 template<typename OrdinalType, typename ScalarType>
67 class SerialTriDiMatrix : public CompObject, public Object, public BLAS<OrdinalType, ScalarType > {
68 public:
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 
329  bool operator== (const SerialTriDiMatrix<OrdinalType, ScalarType> &Operand) const;
330 
332 
335  bool operator!= (const SerialTriDiMatrix<OrdinalType, ScalarType> &Operand) const;
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  virtual void print(std::ostream& os) const;
372 
374 protected:
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 
395 template<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 
408 template<typename OrdinalType, typename ScalarType>
409 SerialTriDiMatrix<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 
424 template<typename OrdinalType, typename ScalarType>
425 SerialTriDiMatrix<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 
449 template<typename OrdinalType, typename ScalarType>
450 SerialTriDiMatrix<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 
508 template<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 
528 template<typename OrdinalType, typename ScalarType>
530 {
531  deleteArrays();
532 }
533 
534 //----------------------------------------------------------------------------------------------------
535 // Shape methods
536 //----------------------------------------------------------------------------------------------------
537 
538 template<typename OrdinalType, typename ScalarType>
539 int SerialTriDiMatrix<OrdinalType, ScalarType>::shape( OrdinalType numRowsCols_in )
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 
551 template<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 
562 template<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 
608 template<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 
634 template<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 
702 template<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 
714 template<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 
726 template<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 
747 template<typename OrdinalType,typename ScalarType>
748 inline 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 
771 template<typename OrdinalType,typename ScalarType>
772 inline 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 
798 template<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 
823 template<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 
840 template<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 
860 template<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 
879 template<typename OrdinalType, typename ScalarType>
881  return !((*this) == Operand);
882 }
883 
884 //----------------------------------------------------------------------------------------------------
885 // Multiplication method
886 //----------------------------------------------------------------------------------------------------
887 
888 template<typename OrdinalType, typename ScalarType>
890 {
891  this->scale( alpha );
892  return(*this);
893 }
894 
895 template<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 
906 template<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 
924 template<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 
975 template<typename OrdinalType, typename ScalarType>
976 inline 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 
991 template<typename OrdinalType, typename ScalarType>
993 {
994  if (valuesCopied_)
995  {
996  delete [] values_;
997  values_ = 0;
998  valuesCopied_ = false;
999  }
1000 }
1001 
1002 template<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 
1014  if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
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_ */
ScalarType & operator()(OrdinalType rowIndex, OrdinalType colIndex)
Element access method (non-const).
ScalarTraits< ScalarType >::magnitudeType normInf() const
Returns the Infinity-norm of the matrix.
void checkIndex(OrdinalType rowIndex, OrdinalType colIndex=0) const
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.
#define TEUCHOS_CHK_ERR(a)
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.
#define TEUCHOS_MAX(x, y)
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.
void copyMat(SerialTriDiMatrix< OrdinalType, ScalarType > matrix, OrdinalType startCol, ScalarType alpha=ScalarTraits< ScalarType >::zero())
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero())
Set all values in the matrix to a constant value.
#define TEUCHOS_CHK_REF(a)
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.
#define TEUCHOS_MIN(x, y)
int shapeUninitialized(OrdinalType numRows)
Same as shape() except leaves uninitialized.
Teuchos::DataAccess Mode enumerable type.