Teuchos - Trilinos Tools Package  Version of the Day
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 
102  SerialTriDiMatrix(DataAccess CV, ScalarType* values, OrdinalType numRowsCols);
103 
105 
111 
113 
124  SerialTriDiMatrix(DataAccess CV, const SerialTriDiMatrix<OrdinalType, ScalarType> &Source, OrdinalType numRowsCols, OrdinalType startRowCols=0);
125 
127  virtual ~SerialTriDiMatrix();
129 
131 
132 
142  int shape(OrdinalType numRows);
143 
145  int shapeUninitialized(OrdinalType numRows);
146 
148 
157  int reshape(OrdinalType numRowsCols);
158 
160 
162 
163 
165 
172 
174 
180 
182 
185  SerialTriDiMatrix<OrdinalType, ScalarType>& operator= (const ScalarType value) { putScalar(value); return(*this); }
186 
188 
192  int putScalar( const ScalarType value = Teuchos::ScalarTraits<ScalarType>::zero() );
193 
195  // int random();
196 
198 
200 
201 
203 
210  ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex);
211 
213 
220  const ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex) const;
221 
223 
230  // ScalarType* operator [] (OrdinalType colIndex);
231 
233 
240  // const ScalarType* operator [] (OrdinalType colIndex) const;
241 
243 
244  ScalarType* values() const { return(values_); }
245 
246  ScalarType* D() const { return D_;}
247  ScalarType* DL() const { return DL_;}
248  ScalarType* DU() const { return DU_;}
249  ScalarType* DU2() const { return DU2_;}
250 
252 
254 
255 
257 
261 
263 
267 
269 
273 
275 
279  int scale ( const ScalarType alpha );
280 
282 
289 
291 
305  //int multiply (ETransp transa, ETransp transb, ScalarType alpha, const SerialTriDiMatrix<OrdinalType, ScalarType> &A, const SerialTriDiMatrix<OrdinalType, ScalarType> &B, ScalarType beta);
306 
308 
319  //int multiply (ESide sideA, ScalarType alpha, const SerialSymTriDiMatrix<OrdinalType, ScalarType> &A, const SerialTriDiMatrix<OrdinalType, ScalarType> &B, ScalarType beta);
320 
322 
324 
325 
327 
330  bool operator== (const SerialTriDiMatrix<OrdinalType, ScalarType> &Operand) const;
331 
333 
336  bool operator!= (const SerialTriDiMatrix<OrdinalType, ScalarType> &Operand) const;
337 
339 
341 
342 
344  OrdinalType numRowsCols() const { return(numRowsCols_); }
345 
347  // OrdinalType numCols() const { return(numRowsCols_); }
348 
350  // OrdinalType stride() const { return(stride_); }
351 
353  bool empty() const { return(numRowsCols_ == 0); }
355 
357 
358 
361 
364 
368 
370 
371  virtual void print(std::ostream& os) const;
373 
375 protected:
376  void copyMat(SerialTriDiMatrix<OrdinalType, ScalarType> matrix,
377  OrdinalType startCol,
378  ScalarType alpha = ScalarTraits<ScalarType>::zero() );
379  void deleteArrays();
380  void checkIndex( OrdinalType rowIndex, OrdinalType colIndex = 0 ) const;
381  OrdinalType numRowsCols_;
382 
383  bool valuesCopied_;
384  ScalarType* values_;
385  ScalarType* DL_;
386  ScalarType* D_;
387  ScalarType* DU_;
388  ScalarType* DU2_;
389 
390 }; // class Teuchos_SerialTriDiMatrix
391 
392 //----------------------------------------------------------------------------------------------------
393 // Constructors and Destructor
394 //----------------------------------------------------------------------------------------------------
395 
396 template<typename OrdinalType, typename ScalarType>
398  :
399  CompObject(),
400  numRowsCols_(0),
401  valuesCopied_(false),
402  values_(0),
403  DL_(NULL),
404  D_(NULL),
405  DU_(NULL),
406  DU2_(NULL)
407 {}
408 
409 template<typename OrdinalType, typename ScalarType>
410 SerialTriDiMatrix<OrdinalType, ScalarType>::SerialTriDiMatrix( OrdinalType numRowsCols_in, OrdinalType numCols_in, bool zeroOut)
411  : CompObject(), numRowsCols_(numRowsCols_in) {
412 
413  OrdinalType numvals = (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
414  values_ = new ScalarType [numvals];
415  DL_ = values_;
416  D_ = DL_ + (numRowsCols_-1);
417  DU_ = D_ + numRowsCols_;
418  DU2_ = DU_ + (numRowsCols_-1);
419 
420  valuesCopied_ = true;
421  if (zeroOut == true)
422  putScalar();
423 }
424 
425 template<typename OrdinalType, typename ScalarType>
426 SerialTriDiMatrix<OrdinalType, ScalarType>::SerialTriDiMatrix(DataAccess CV, ScalarType* values_in, OrdinalType numRowsCols_in )
427  : CompObject(), numRowsCols_(numRowsCols_in),
428  valuesCopied_(false), values_(values_in)
429 {
430  if(CV == Copy) {
431  const OrdinalType numvals = (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
432  values_ = new ScalarType[numvals];
433  DL_ = values_;
434  D_ = DL_ + (numRowsCols_-1);
435  DU_ = D_ + numRowsCols_;
436  DU2_ = DU_ + (numRowsCols_-1);
437  valuesCopied_ = true;
438 
439  for(OrdinalType i = 0 ; i < numRowsCols_ ; ++i )
440  values_[i] = values_in[i];
441  }
442 }
443 
444 template<typename OrdinalType, typename ScalarType>
445 SerialTriDiMatrix<OrdinalType, ScalarType>::SerialTriDiMatrix(const SerialTriDiMatrix<OrdinalType, ScalarType> &Source, ETransp trans) : CompObject(), BLAS<OrdinalType,ScalarType>(), numRowsCols_(0), valuesCopied_(true), values_(0)
446 {
447  if ( trans == Teuchos::NO_TRANS ) {
448  numRowsCols_ = Source.numRowsCols_;
449 
450  const OrdinalType numvals = (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
451  values_ = new ScalarType[numvals];
452  DL_ = values_;
453  D_ = DL_+ (numRowsCols_-1);
454  DU_ = D_ + numRowsCols_;
455  DU2_ = DU_ + (numRowsCols_-1);
456 
457  copyMat(Source, 0, 0);
458  }
460  {
461  numRowsCols_ = Source.numRowsCols_;
462  const OrdinalType numvals = (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
463  values_ = new ScalarType[numvals];
464  DL_ = values_;
465  D_ = DL_+(numRowsCols_-1);
466  DU_ = D_ + numRowsCols_;
467  DU2_ = DU_ + (numRowsCols_-1);
468 
469  OrdinalType min = numRowsCols_;
470  if(min > Source.numRowsCols_) min = Source.numRowsCols_;
471 
472  for(OrdinalType i = 0 ; i< min ; ++i) {
474  if(i < (min-1)) {
477  }
478  if(i < (min-2)) {
479  DU2_ = Teuchos::ScalarTraits<ScalarType>::conjugate(Source.DU2_[i]);
480  }
481  }
482  }
483  else
484  {
485  numRowsCols_ = Source.numCols_;
486  const OrdinalType numvals = (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
487  values_ = new ScalarType[numvals];
488  OrdinalType min = numRowsCols_;
489  if(min > Source.numRowsCols_) min = Source.numRowsCols_;
490  for(OrdinalType i = 0 ; i< min ; ++i) {
491  D_ = Source.D_[i];
492  if(i < (min-1)) {
493  DL_ = Source.DL_[i];
494  DU_ = Source.DU_[i];
495  }
496  if(i < (min-2)) {
497  DU2_ = Source.DU2_[i];
498  }
499  }
500  }
501 }
502 
503 template<typename OrdinalType, typename ScalarType>
506  OrdinalType numRowsCols_in, OrdinalType startRow )
507  : CompObject(), numRowsCols_(numRowsCols_in),
508  valuesCopied_(false), values_(Source.values_) {
509  if(CV == Copy)
510  {
511  const OrdinalType numvals = (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
512  values_ = new ScalarType[numvals];
513  copyMat(Source, startRow);
514  valuesCopied_ = true;
515  }
516  else // CV == View
517  {
518  // \todo ???
519  // values_ = values_ + (stride_ * startCol) + startRow;
520  }
521 }
522 
523 template<typename OrdinalType, typename ScalarType>
525 {
526  deleteArrays();
527 }
528 
529 //----------------------------------------------------------------------------------------------------
530 // Shape methods
531 //----------------------------------------------------------------------------------------------------
532 
533 template<typename OrdinalType, typename ScalarType>
534 int SerialTriDiMatrix<OrdinalType, ScalarType>::shape( OrdinalType numRowsCols_in )
535 {
536  deleteArrays(); // Get rid of anything that might be already allocated
537  numRowsCols_ = numRowsCols_in;
538  const OrdinalType numvals = ( numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
539  values_ = new ScalarType[numvals];
540 
541  putScalar();
542  valuesCopied_ = true;
543  return(0);
544 }
545 
546 template<typename OrdinalType, typename ScalarType>
548 {
549  deleteArrays(); // Get rid of anything that might be already allocated
550  numRowsCols_ = numRowsCols_in;
551  const OrdinalType numvals = ( numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
552  values_ = new ScalarType[numvals];
553  valuesCopied_ = true;
554  return(0);
555 }
556 
557 template<typename OrdinalType, typename ScalarType>
559 {
560  if(numRowsCols_in <1 ) {
561  deleteArrays();
562  return 0;
563  }
564  // Allocate space for new matrix
565  const OrdinalType numvals = ( numRowsCols_in == 1) ? 1 : 4*(numRowsCols_in - 1);
566  ScalarType* values_tmp = new ScalarType[numvals];
567  ScalarType zero = ScalarTraits<ScalarType>::zero();
568  for(OrdinalType i= 0; i<numvals ; ++i)
569  values_tmp[i] = zero;
570 
571  OrdinalType min = TEUCHOS_MIN(numRowsCols_, numRowsCols_in);
572  ScalarType* dl = values_tmp;
573  ScalarType* d = values_tmp + (numRowsCols_in-1);
574  ScalarType* du = d+(numRowsCols_in);
575  ScalarType* du2 = du+(numRowsCols_in - 1);
576 
577  if(values_ != 0) {
578  for(OrdinalType i = 0 ; i< min ; ++i) {
579  dl[i] = DL_[i];
580  d[i] = D_[i];
581  du[i] = DU_[i];
582  du2[i] = DU2_[i];
583  }
584  }
585 
586  deleteArrays(); // Get rid of anything that might be already allocated
587  numRowsCols_ = numRowsCols_in;
588 
589  values_ = values_tmp; // Set pointer to new A
590  DL_ = dl;
591  D_ = d;
592  DU_ = du;
593  DU2_ = du2;
594 
595  valuesCopied_ = true;
596  return(0);
597 }
598 
599 //----------------------------------------------------------------------------------------------------
600 // Set methods
601 //----------------------------------------------------------------------------------------------------
602 
603 template<typename OrdinalType, typename ScalarType>
605  // Set each value of the TriDi matrix to "value".
606  const OrdinalType numvals = (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1);
607 
608  for(OrdinalType i = 0; i<numvals ; ++i)
609  {
610  values_[i] = value_in;
611  }
612  return 0;
613 }
614 
615 // template<typename OrdinalType, typename ScalarType>
616 // int SerialTriDiMatrix<OrdinalType, ScalarType>::random()
617 // {
618 // // Set each value of the TriDi matrix to a random value.
619 // for(OrdinalType j = 0; j < numCols_; j++)
620 // {
621 // for(OrdinalType i = 0; i < numRowsCols_; i++)
622 // {
623 // values_[i + j*stride_] = ScalarTraits<ScalarType>::random();
624 // }
625 // }
626 // return 0;
627 // }
628 
629 template<typename OrdinalType, typename ScalarType>
632 {
633  if(this == &Source)
634  return(*this); // Special case of source same as target
635  if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
636  return(*this); // Special case of both are views to same data.
637 
638  // If the source is a view then we will return a view, else we will return a copy.
639  if (!Source.valuesCopied_) {
640  if(valuesCopied_) {
641  // Clean up stored data if this was previously a copy.
642  deleteArrays();
643  }
644  numRowsCols_ = Source.numRowsCols_;
645  values_ = Source.values_;
646  }
647  else {
648  // If we were a view, we will now be a copy.
649  if(!valuesCopied_) {
650  numRowsCols_ = Source.numRowsCols_;
651  const OrdinalType numvals = ( Source.numRowsCols_ == 1) ? 1 : 4*(Source.numRowsCols_ - 1);
652  if(numvals > 0) {
653  values_ = new ScalarType[numvals];
654  valuesCopied_ = true;
655  }
656  else {
657  values_ = 0;
658  }
659  }
660  // If we were a copy, we will stay a copy.
661  else {
662  // we need to allocate more space (or less space)
663  deleteArrays();
664  numRowsCols_ = Source.numRowsCols_;
665  const OrdinalType numvals = ( Source.numRowsCols_ == 1) ? 1 : 4*(Source.numRowsCols_ - 1);
666  if(numvals > 0) {
667  values_ = new ScalarType[numvals];
668  valuesCopied_ = true;
669  }
670  }
671  DL_ = values_;
672  if(values_ != 0) {
673  D_ = DL_ + (numRowsCols_-1);
674  DU_ = D_ + numRowsCols_;
675  DU2_ = DU_ + (numRowsCols_-1);
676 
677  OrdinalType min = TEUCHOS_MIN(numRowsCols_, Source.numRowsCols_);
678  for(OrdinalType i = 0 ; i < min ; ++i ) {
679  D_[i] = Source.D()[i];
680  if(i< (min-1 ) )
681  {
682  DL_[i] = Source.DL()[i];
683  DU_[i] = Source.DU()[i];
684  }
685  if(i< (min-2))
686  DU2_[i] = Source.DU2()[i];
687  }
688 
689  }
690  else {
691  D_ = DU_ = DU2_ = 0;
692  }
693  }
694  return(*this);
695 }
696 
697 template<typename OrdinalType, typename ScalarType>
699 {
700  // Check for compatible dimensions
701  if ((numRowsCols_ != Source.numRowsCols_) )
702  {
703  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
704  }
705  copyMat(Source, 0, ScalarTraits<ScalarType>::one());
706  return(*this);
707 }
708 
709 template<typename OrdinalType, typename ScalarType>
711 {
712  // Check for compatible dimensions
713  if ((numRowsCols_ != Source.numRowsCols_) )
714  {
715  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
716  }
717  copyMat(Source, 0, -ScalarTraits<ScalarType>::one());
718  return(*this);
719 }
720 
721 template<typename OrdinalType,typename ScalarType>
723 {
724  if(this == &Source)
725  return(*this); // Special case of source same as target
726  if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
727  return(*this); // Special case of both are views to same data.
728 
729  // Check for compatible dimensions
730  if ((numRowsCols_ != Source.numRowsCols_) )
731  {
732  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
733  }
734  copyMat(Source, 0, 0);
735  return(*this);
736 }
737 
738 //----------------------------------------------------------------------------------------------------
739 // Accessor methods
740 //----------------------------------------------------------------------------------------------------
741 
742 template<typename OrdinalType,typename ScalarType>
743 inline const ScalarType& SerialTriDiMatrix<OrdinalType,ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex) const
744 {
745  OrdinalType diff = colIndex - rowIndex;
746 
747  //#ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
748  checkIndex( rowIndex, colIndex );
749  //#endif
750  switch (diff) {
751  case -1:
752  return DL_[colIndex];
753  break;
754  case 0:
755  return D_[colIndex];
756  break;
757  case 1:
758  return DU_[rowIndex];
759  break;
760  case 2:
761  return DU2_[rowIndex];
762  break;
763  default:
764  TEUCHOS_TEST_FOR_EXCEPTION(true, std::out_of_range,
765  "SerialTriDiMatrix<T>::operator (row,col) "
766  "Index (" << rowIndex <<","<<colIndex<<") out of range ");
767  }
768  TEUCHOS_TEST_FOR_EXCEPTION(true, std::out_of_range,
769  "SerialTriDiMatrix<T>::operator (row,col) "
770  "Index (" << rowIndex <<","<<colIndex<<") out of range ");
771  return D_[0];
772 }
773 
774 template<typename OrdinalType,typename ScalarType>
775 inline ScalarType& Teuchos::SerialTriDiMatrix<OrdinalType,ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex)
776 {
777  OrdinalType diff = colIndex - rowIndex;
778  //#ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
779  checkIndex( rowIndex, colIndex );
780  //#endif
781  switch (diff) {
782  case -1:
783  return DL_[colIndex];
784  break;
785  case 0:
786  return D_[colIndex];
787  break;
788  case 1:
789  return DU_[rowIndex];
790  break;
791  case 2:
792  return DU2_[rowIndex];
793  break;
794  default:
795  TEUCHOS_TEST_FOR_EXCEPTION(true, std::out_of_range,
796  "SerialTriDiMatrix<T>::operator (row,col) "
797  "Index (" << rowIndex <<","<<colIndex<<") out of range ");
798  }
799  TEUCHOS_TEST_FOR_EXCEPTION(true, std::out_of_range,
800  "SerialTriDiMatrix<T>::operator (row,col) "
801  "Index (" << rowIndex <<","<<colIndex<<") out of range ");
802  // TEUCHOS_CHK_ERR(-1);
803  return D_[0];
804 }
805 
806 //----------------------------------------------------------------------------------------------------
807 // Norm methods
808 //----------------------------------------------------------------------------------------------------
809 
810 template<typename OrdinalType,typename ScalarType>
812 {
813  OrdinalType i, j;
816  ScalarType* ptr;
817 
818  // Fix this for Tri DI
819 
820  for(j = 0; j < numRowsCols_; j++)
821  {
822  ScalarType sum = 0;
823  if(j-1>=0) sum += ScalarTraits<ScalarType>::magnitude((*this)(j-1,j));
824  sum+= ScalarTraits<ScalarType>::magnitude((*this)(j,j));
825  if(j+1<numRowsCols_) sum+= ScalarTraits<ScalarType>::magnitude((*this)(j+1,j));
827  if(absSum > anorm)
828  {
829  anorm = absSum;
830  }
831  }
832  updateFlops(numRowsCols_ * numRowsCols_);
833  return(anorm);
834 }
835 
836 template<typename OrdinalType, typename ScalarType>
838 {
839  OrdinalType i,j;
841 
842  for (i = 0; i < numRowsCols_; i++) {
844  for (j=i-1; j<= i+1; j++) {
845  if(j >= 0 && j < numRowsCols_) sum += ScalarTraits<ScalarType>::magnitude((*this)(i,j));
846  }
847  anorm = TEUCHOS_MAX( anorm, sum );
848  }
849  updateFlops(numRowsCols_ * numRowsCols_);
850  return(anorm);
851 }
852 
853 template<typename OrdinalType, typename ScalarType>
855 {
856  // \todo fix this
857  OrdinalType i, j;
859  for (j = 0; j < numRowsCols_; j++) {
860  for (i = j-1; i <= j+1; i++) {
861  if(i >= 0 && i < numRowsCols_ ) anorm += ScalarTraits<ScalarType>::magnitude((*this)(i,j));
862  }
863  }
865  updateFlops( (numRowsCols_ == 1) ? 1 : 4*(numRowsCols_-1) );
866  return(anorm);
867 }
868 
869 //----------------------------------------------------------------------------------------------------
870 // Comparison methods
871 //----------------------------------------------------------------------------------------------------
872 
873 template<typename OrdinalType, typename ScalarType>
875 {
876  bool allmatch = true;
877  bool result = 1;
878  if((numRowsCols_ != Operand.numRowsCols_) )
879  {
880  result = 0;
881  }
882  else
883  {
884  OrdinalType numvals = (numRowsCols_ == 1)? 1 : 4*(numRowsCols_ -1 );
885 
886  for(OrdinalType i = 0; i< numvals; ++i)
887  allmatch &= (Operand.values_[i] == values_[i]);
888  }
889  return allmatch;
890 }
891 
892 template<typename OrdinalType, typename ScalarType>
894  return !((*this) == Operand);
895 }
896 
897 //----------------------------------------------------------------------------------------------------
898 // Multiplication method
899 //----------------------------------------------------------------------------------------------------
900 
901 template<typename OrdinalType, typename ScalarType>
903 {
904  this->scale( alpha );
905  return(*this);
906 }
907 
908 template<typename OrdinalType, typename ScalarType>
910 {
911  OrdinalType i;
912  OrdinalType numvals = (numRowsCols_ == 1)? 1 : 4*(numRowsCols_ -1 );
913  for (i=0; i < numvals ; ++i ) {
914  values_[i] *= alpha;
915  }
916  return(0);
917 }
918 
919 template<typename OrdinalType, typename ScalarType>
921 {
922  OrdinalType j;
923 
924  // Check for compatible dimensions
925  if ((numRowsCols_ != A.numRowsCols_) )
926  {
927  TEUCHOS_CHK_ERR(-1); // Return error
928  }
929  OrdinalType numvals = (numRowsCols_ == 1)? 1 : 4*(numRowsCols_ -1 );
930  for (j=0; j<numvals; j++) {
931  values_[j] = A.values_ * values_[j];
932  }
933  updateFlops( numvals );
934  return(0);
935 }
936 
937 template<typename OrdinalType, typename ScalarType>
939 {
940  os << std::endl;
941  if(valuesCopied_)
942  os << "A_Copied: yes" << std::endl;
943  else
944  os << "A_Copied: no" << std::endl;
945  os << "Rows and Columns: " << numRowsCols_ << std::endl;
946  if(numRowsCols_ == 0) {
947  os << "(matrix is empty, no values to display)" << std::endl;
948  return;
949  }
950  else
951  {
952  os << "DL: "<<std::endl;
953  for(int i=0;i<numRowsCols_-1;++i)
954  os << DL_[i]<<" ";
955  os << std::endl;
956  os << "D: "<<std::endl;
957  for(int i=0;i<numRowsCols_;++i)
958  os << D_[i]<<" ";
959  os << std::endl;
960  os << "DU: "<<std::endl;
961  for(int i=0;i<numRowsCols_-1;++i)
962  os << DU_[i]<<" ";
963  os << std::endl;
964  os << "DU2: "<<std::endl;
965  for(int i=0;i<numRowsCols_-2;++i)
966  os << DU2_[i]<<" ";
967  os << std::endl;
968  }
969 
970  os <<" square format:"<<std::endl;
971  for(int i=0 ; i < numRowsCols_ ; ++i ) {
972  for(int j=0;j<numRowsCols_;++j) {
973  if ( j >= i-1 && j <= i+1) {
974  os << (*this)(i,j)<<" ";
975  }
976  else
977  os <<". ";
978  }
979  os << std::endl;
980  }
981 
982 }
983 
984 //----------------------------------------------------------------------------------------------------
985 // Protected methods
986 //----------------------------------------------------------------------------------------------------
987 
988 template<typename OrdinalType, typename ScalarType>
989 inline void SerialTriDiMatrix<OrdinalType, ScalarType>::checkIndex( OrdinalType rowIndex, OrdinalType colIndex ) const
990 {
991  OrdinalType diff = colIndex - rowIndex;
992  TEUCHOS_TEST_FOR_EXCEPTION(rowIndex < 0 || rowIndex >= numRowsCols_, std::out_of_range,
993  "SerialTriDiMatrix<T>::checkIndex: "
994  "Row index " << rowIndex << " out of range [0, "<< numRowsCols_ << "]");
995  TEUCHOS_TEST_FOR_EXCEPTION(colIndex < 0 || colIndex >= numRowsCols_,
996  std::out_of_range,
997  "SerialTriDiMatrix<T>::checkIndex: "
998  "Col index " << colIndex << " out of range [0, "<< numRowsCols_ << "]");
999  TEUCHOS_TEST_FOR_EXCEPTION(diff > 2 || diff < -1 , std::out_of_range,
1000  "SerialTriDiMatrix<T>::checkIndex: "
1001  "index difference " << diff << " out of range [-1, 2]");
1002 }
1003 
1004 template<typename OrdinalType, typename ScalarType>
1006 {
1007  if (valuesCopied_)
1008  {
1009  delete [] values_;
1010  values_ = 0;
1011  valuesCopied_ = false;
1012  }
1013 }
1014 
1015 template<typename OrdinalType, typename ScalarType>
1017  OrdinalType startRowCol,
1018  ScalarType alpha )
1019 {
1020  OrdinalType i;
1021  OrdinalType max = inputMatrix.numRowsCols_;
1022  if(max > numRowsCols_ ) max = numRowsCols_;
1023  if(startRowCol > max ) return; //
1024 
1025  for(i = startRowCol ; i < max ; ++i) {
1026 
1027  if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
1028  // diagonal
1029  D()[i] += inputMatrix.D()[i];
1030  if(i<(max-1) && (i-1) >= startRowCol) {
1031  DL()[i] += inputMatrix.DL()[i];
1032  DU()[i] += inputMatrix.DU()[i];
1033  }
1034  if(i<(max-2) && (i-2) >= startRowCol) {
1035  DU2()[i] += inputMatrix.DU2()[i];
1036  }
1037  }
1038  else {
1039  // diagonal
1040  D()[i] = inputMatrix.D()[i];
1041  if(i<(max-1) && (i-1) >= startRowCol) {
1042  DL()[i] = inputMatrix.DL()[i];
1043  DU()[i] = inputMatrix.DU()[i];
1044  }
1045  if(i<(max-2) && (i-2) >= startRowCol) {
1046  DU2()[i] = inputMatrix.DU2()[i];
1047  }
1048  }
1049  }
1050 }
1051 
1052 }
1053 
1054 
1055 #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 << 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.
The Teuchos namespace contains all of the classes, structs and enums used by Teuchos, as well as a number of utility routines.
void updateFlops(int addflops) const
Increment Flop count for this object.
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.