Teuchos - Trilinos Tools Package  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Teuchos_SerialSymDenseMatrix.hpp
Go to the documentation of this file.
1 
2 // @HEADER
3 // ***********************************************************************
4 //
5 // Teuchos: Common Tools Package
6 // Copyright (2004) Sandia Corporation
7 //
8 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
9 // license for use of this work by or on behalf of the U.S. Government.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39 //
40 // ***********************************************************************
41 // @HEADER
42 
43 #ifndef _TEUCHOS_SERIALSYMDENSEMATRIX_HPP_
44 #define _TEUCHOS_SERIALSYMDENSEMATRIX_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 #include <utility>
56 #include <vector>
57 
120 namespace Teuchos {
121 
122 template<typename OrdinalType, typename ScalarType>
123 class SerialSymDenseMatrix : public CompObject, public BLAS<OrdinalType,ScalarType>
124 {
125  public:
126 
128  typedef OrdinalType ordinalType;
130  typedef ScalarType scalarType;
131 
133 
134 
144 
146 
156  SerialSymDenseMatrix(OrdinalType numRowsCols, bool zeroOut = true);
157 
159 
170  SerialSymDenseMatrix(DataAccess CV, bool upper, ScalarType* values, OrdinalType stride, OrdinalType numRowsCols);
171 
174 
176 
185  SerialSymDenseMatrix(DataAccess CV, const SerialSymDenseMatrix<OrdinalType, ScalarType> &Source, OrdinalType numRowCols, OrdinalType startRowCol=0);
186 
188  virtual ~SerialSymDenseMatrix ();
190 
192 
193 
195 
204  int shape(OrdinalType numRowsCols);
205 
207 
216  int shapeUninitialized(OrdinalType numRowsCols);
217 
219 
229  int reshape(OrdinalType numRowsCols);
230 
232 
234  void setLower();
235 
237 
239  void setUpper();
240 
242 
244 
245 
247 
254 
256 
261 
263 
266  SerialSymDenseMatrix<OrdinalType, ScalarType>& operator= (const ScalarType value) { putScalar(value); return(*this); }
267 
269 
274  int putScalar( const ScalarType value = Teuchos::ScalarTraits<ScalarType>::zero(), bool fullMatrix = false );
275 
277 
281 
283 
285  int random( const ScalarType bias = 0.1*Teuchos::ScalarTraits<ScalarType>::one() );
286 
288 
290 
291 
293 
300  ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex);
301 
303 
310  const ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex) const;
311 
313 
315  ScalarType* values() const { return(values_); }
316 
318 
320 
321 
323  bool upper() const {return(upper_);};
324 
326  char UPLO() const {return(UPLO_);};
328 
330 
331 
333 
340 
342 
346 
348 
352 
354 
356 
357 
359 
363 
365 
369 
371 
373 
374 
376  OrdinalType numRows() const { return(numRowCols_); }
377 
379  OrdinalType numCols() const { return(numRowCols_); }
380 
382  OrdinalType stride() const { return(stride_); }
383 
385  bool empty() const { return(numRowCols_ == 0); }
386 
388 
390 
391 
394 
397 
401 
403 
404 #ifndef TEUCHOS_HIDE_DEPRECATED_CODE
406  virtual void print(std::ostream& os) const;
407 #else
408  virtual std::ostream& print(std::ostream& os) const;
409 #endif
410 
412 
413  protected:
414 
415  // In-place scaling of matrix.
416  void scale( const ScalarType alpha );
417 
418  // Copy the values from one matrix to the other.
419  void copyMat(bool inputUpper, ScalarType* inputMatrix, OrdinalType inputStride,
420  OrdinalType numRowCols, bool outputUpper, ScalarType* outputMatrix,
421  OrdinalType outputStride, OrdinalType startRowCol,
422  ScalarType alpha = ScalarTraits<ScalarType>::zero() );
423 
424  // Copy the values from the active triangle of the matrix to the other to make the matrix full symmetric.
425  void copyUPLOMat(bool inputUpper, ScalarType* inputMatrix,
426  OrdinalType inputStride, OrdinalType inputRows);
427 
428  void deleteArrays();
429  void checkIndex( OrdinalType rowIndex, OrdinalType colIndex = 0 ) const;
430  OrdinalType numRowCols_;
431  OrdinalType stride_;
432  bool valuesCopied_;
433  ScalarType* values_;
434  bool upper_;
435  char UPLO_;
436 
437 
438 };
439 
440 //----------------------------------------------------------------------------------------------------
441 // Constructors and Destructor
442 //----------------------------------------------------------------------------------------------------
443 
444 template<typename OrdinalType, typename ScalarType>
446  : CompObject(), BLAS<OrdinalType,ScalarType>(),
447  numRowCols_(0), stride_(0), valuesCopied_(false), values_(0), upper_(false), UPLO_('L')
448 {}
449 
450 template<typename OrdinalType, typename ScalarType>
452  : CompObject(), BLAS<OrdinalType,ScalarType>(),
453  numRowCols_(numRowCols_in), stride_(numRowCols_in), valuesCopied_(false), upper_(false), UPLO_('L')
454 {
455  values_ = new ScalarType[stride_*numRowCols_];
456  valuesCopied_ = true;
457  if (zeroOut == true)
459 }
460 
461 template<typename OrdinalType, typename ScalarType>
463  DataAccess CV, bool upper_in, ScalarType* values_in, OrdinalType stride_in, OrdinalType numRowCols_in
464  )
465  : CompObject(), BLAS<OrdinalType,ScalarType>(),
466  numRowCols_(numRowCols_in), stride_(stride_in), valuesCopied_(false),
467  values_(values_in), upper_(upper_in)
468 {
469  if (upper_)
470  UPLO_ = 'U';
471  else
472  UPLO_ = 'L';
473 
474  if(CV == Copy)
475  {
476  stride_ = numRowCols_;
477  values_ = new ScalarType[stride_*numRowCols_];
478  copyMat(upper_in, values_in, stride_in, numRowCols_, upper_, values_, stride_, 0);
479  valuesCopied_ = true;
480  }
481 }
482 
483 template<typename OrdinalType, typename ScalarType>
485  : CompObject(), BLAS<OrdinalType,ScalarType>(),
486  numRowCols_(Source.numRowCols_), stride_(0), valuesCopied_(true),
487  values_(0), upper_(Source.upper_), UPLO_(Source.UPLO_)
488 {
489  if (!Source.valuesCopied_)
490  {
491  stride_ = Source.stride_;
492  values_ = Source.values_;
493  valuesCopied_ = false;
494  }
495  else
496  {
497  stride_ = numRowCols_;
498  const OrdinalType newsize = stride_ * numRowCols_;
499  if(newsize > 0) {
500  values_ = new ScalarType[newsize];
501  copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0);
502  }
503  else {
504  numRowCols_ = 0; stride_ = 0;
505  valuesCopied_ = false;
506  }
507  }
508 }
509 
510 template<typename OrdinalType, typename ScalarType>
512  DataAccess CV, const SerialSymDenseMatrix<OrdinalType,
513  ScalarType> &Source, OrdinalType numRowCols_in, OrdinalType startRowCol )
514  : CompObject(), BLAS<OrdinalType,ScalarType>(),
515  numRowCols_(numRowCols_in), stride_(Source.stride_), valuesCopied_(false), upper_(Source.upper_), UPLO_(Source.UPLO_)
516 {
517  if(CV == Copy)
518  {
519  stride_ = numRowCols_in;
520  values_ = new ScalarType[stride_ * numRowCols_in];
521  copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_in, upper_, values_, stride_, startRowCol);
522  valuesCopied_ = true;
523  }
524  else // CV == View
525  {
526  values_ = Source.values_ + (stride_ * startRowCol) + startRowCol;
527  }
528 }
529 
530 template<typename OrdinalType, typename ScalarType>
532 {
533  deleteArrays();
534 }
535 
536 //----------------------------------------------------------------------------------------------------
537 // Shape methods
538 //----------------------------------------------------------------------------------------------------
539 
540 template<typename OrdinalType, typename ScalarType>
542 {
543  deleteArrays(); // Get rid of anything that might be already allocated
544  numRowCols_ = numRowCols_in;
545  stride_ = numRowCols_;
546  values_ = new ScalarType[stride_*numRowCols_];
547  putScalar( Teuchos::ScalarTraits<ScalarType>::zero(), true );
548  valuesCopied_ = true;
549  return(0);
550 }
551 
552 template<typename OrdinalType, typename ScalarType>
554 {
555  deleteArrays(); // Get rid of anything that might be already allocated
556  numRowCols_ = numRowCols_in;
557  stride_ = numRowCols_;
558  values_ = new ScalarType[stride_*numRowCols_];
559  valuesCopied_ = true;
560  return(0);
561 }
562 
563 template<typename OrdinalType, typename ScalarType>
565 {
566  // Allocate space for new matrix
567  ScalarType* values_tmp = new ScalarType[numRowCols_in * numRowCols_in];
568  ScalarType zero = ScalarTraits<ScalarType>::zero();
569  for(OrdinalType k = 0; k < numRowCols_in * numRowCols_in; k++)
570  {
571  values_tmp[k] = zero;
572  }
573  OrdinalType numRowCols_tmp = TEUCHOS_MIN(numRowCols_, numRowCols_in);
574  if(values_ != 0)
575  {
576  copyMat(upper_, values_, stride_, numRowCols_tmp, upper_, values_tmp, numRowCols_in, 0); // Copy principal submatrix of A to new A
577  }
578  deleteArrays(); // Get rid of anything that might be already allocated
579  numRowCols_ = numRowCols_in;
580  stride_ = numRowCols_;
581  values_ = values_tmp; // Set pointer to new A
582  valuesCopied_ = true;
583  return(0);
584 }
585 
586 //----------------------------------------------------------------------------------------------------
587 // Set methods
588 //----------------------------------------------------------------------------------------------------
589 
590 template<typename OrdinalType, typename ScalarType>
592 {
593  // Do nothing if the matrix is already a lower triangular matrix
594  if (upper_ != false) {
595  copyUPLOMat( true, values_, stride_, numRowCols_ );
596  upper_ = false;
597  UPLO_ = 'L';
598  }
599 }
600 
601 template<typename OrdinalType, typename ScalarType>
603 {
604  // Do nothing if the matrix is already an upper triangular matrix
605  if (upper_ == false) {
606  copyUPLOMat( false, values_, stride_, numRowCols_ );
607  upper_ = true;
608  UPLO_ = 'U';
609  }
610 }
611 
612 template<typename OrdinalType, typename ScalarType>
613 int SerialSymDenseMatrix<OrdinalType, ScalarType>::putScalar( const ScalarType value_in, bool fullMatrix )
614 {
615  // Set each value of the dense matrix to "value".
616  if (fullMatrix == true) {
617  for(OrdinalType j = 0; j < numRowCols_; j++)
618  {
619  for(OrdinalType i = 0; i < numRowCols_; i++)
620  {
621  values_[i + j*stride_] = value_in;
622  }
623  }
624  }
625  // Set the active upper or lower triangular part of the matrix to "value"
626  else {
627  if (upper_) {
628  for(OrdinalType j = 0; j < numRowCols_; j++) {
629  for(OrdinalType i = 0; i <= j; i++) {
630  values_[i + j*stride_] = value_in;
631  }
632  }
633  }
634  else {
635  for(OrdinalType j = 0; j < numRowCols_; j++) {
636  for(OrdinalType i = j; i < numRowCols_; i++) {
637  values_[i + j*stride_] = value_in;
638  }
639  }
640  }
641  }
642  return 0;
643 }
644 
645 template<typename OrdinalType, typename ScalarType> void
648 {
649  std::swap(values_ , B.values_);
650  std::swap(numRowCols_, B.numRowCols_);
651  std::swap(stride_, B.stride_);
652  std::swap(valuesCopied_, B.valuesCopied_);
653  std::swap(upper_, B.upper_);
654  std::swap(UPLO_, B.UPLO_);
655 }
656 
657 template<typename OrdinalType, typename ScalarType>
659 {
661 
662  // Set each value of the dense matrix to a random value.
663  std::vector<MT> diagSum( numRowCols_, Teuchos::ScalarTraits<MT>::zero() );
664  if (upper_) {
665  for(OrdinalType j = 0; j < numRowCols_; j++) {
666  for(OrdinalType i = 0; i < j; i++) {
667  values_[i + j*stride_] = ScalarTraits<ScalarType>::random();
668  diagSum[i] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
669  diagSum[j] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
670  }
671  }
672  }
673  else {
674  for(OrdinalType j = 0; j < numRowCols_; j++) {
675  for(OrdinalType i = j+1; i < numRowCols_; i++) {
676  values_[i + j*stride_] = ScalarTraits<ScalarType>::random();
677  diagSum[i] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
678  diagSum[j] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
679  }
680  }
681  }
682 
683  // Set the diagonal.
684  for(OrdinalType i = 0; i < numRowCols_; i++) {
685  values_[i + i*stride_] = diagSum[i] + bias;
686  }
687  return 0;
688 }
689 
690 template<typename OrdinalType, typename ScalarType>
693 {
694  if(this == &Source)
695  return(*this); // Special case of source same as target
696  if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_)) {
697  upper_ = Source.upper_; // Might have to change the active part of the matrix.
698  return(*this); // Special case of both are views to same data.
699  }
700 
701  // If the source is a view then we will return a view, else we will return a copy.
702  if (!Source.valuesCopied_) {
703  if(valuesCopied_) {
704  // Clean up stored data if this was previously a copy.
705  deleteArrays();
706  }
707  numRowCols_ = Source.numRowCols_;
708  stride_ = Source.stride_;
709  values_ = Source.values_;
710  upper_ = Source.upper_;
711  UPLO_ = Source.UPLO_;
712  }
713  else {
714  // If we were a view, we will now be a copy.
715  if(!valuesCopied_) {
716  numRowCols_ = Source.numRowCols_;
717  stride_ = Source.numRowCols_;
718  upper_ = Source.upper_;
719  UPLO_ = Source.UPLO_;
720  const OrdinalType newsize = stride_ * numRowCols_;
721  if(newsize > 0) {
722  values_ = new ScalarType[newsize];
723  valuesCopied_ = true;
724  }
725  else {
726  values_ = 0;
727  }
728  }
729  // If we were a copy, we will stay a copy.
730  else {
731  if((Source.numRowCols_ <= stride_) && (Source.numRowCols_ == numRowCols_)) { // we don't need to reallocate
732  numRowCols_ = Source.numRowCols_;
733  upper_ = Source.upper_;
734  UPLO_ = Source.UPLO_;
735  }
736  else { // we need to allocate more space (or less space)
737  deleteArrays();
738  numRowCols_ = Source.numRowCols_;
739  stride_ = Source.numRowCols_;
740  upper_ = Source.upper_;
741  UPLO_ = Source.UPLO_;
742  const OrdinalType newsize = stride_ * numRowCols_;
743  if(newsize > 0) {
744  values_ = new ScalarType[newsize];
745  valuesCopied_ = true;
746  }
747  }
748  }
749  copyMat(Source.upper_, Source.values_, Source.stride_, Source.numRowCols_, upper_, values_, stride_, 0);
750  }
751  return(*this);
752 }
753 
754 template<typename OrdinalType, typename ScalarType>
756 {
757  this->scale(alpha);
758  return(*this);
759 }
760 
761 template<typename OrdinalType, typename ScalarType>
763 {
764  // Check for compatible dimensions
765  if ((numRowCols_ != Source.numRowCols_))
766  {
767  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
768  }
769  copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0, ScalarTraits<ScalarType>::one());
770  return(*this);
771 }
772 
773 template<typename OrdinalType, typename ScalarType>
775 {
776  // Check for compatible dimensions
777  if ((numRowCols_ != Source.numRowCols_))
778  {
779  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
780  }
781  copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0, -ScalarTraits<ScalarType>::one());
782  return(*this);
783 }
784 
785 template<typename OrdinalType, typename ScalarType>
787  if(this == &Source)
788  return(*this); // Special case of source same as target
789  if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_)) {
790  upper_ = Source.upper_; // We may have to change the active part of the matrix.
791  return(*this); // Special case of both are views to same data.
792  }
793 
794  // Check for compatible dimensions
795  if ((numRowCols_ != Source.numRowCols_))
796  {
797  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
798  }
799  copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0 );
800  return(*this);
801 }
802 
803 //----------------------------------------------------------------------------------------------------
804 // Accessor methods
805 //----------------------------------------------------------------------------------------------------
806 
807 template<typename OrdinalType, typename ScalarType>
808 inline ScalarType& SerialSymDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex)
809 {
810 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
811  checkIndex( rowIndex, colIndex );
812 #endif
813  if ( rowIndex <= colIndex ) {
814  // Accessing upper triangular part of matrix
815  if (upper_)
816  return(values_[colIndex * stride_ + rowIndex]);
817  else
818  return(values_[rowIndex * stride_ + colIndex]);
819  }
820  else {
821  // Accessing lower triangular part of matrix
822  if (upper_)
823  return(values_[rowIndex * stride_ + colIndex]);
824  else
825  return(values_[colIndex * stride_ + rowIndex]);
826  }
827 }
828 
829 template<typename OrdinalType, typename ScalarType>
830 inline const ScalarType& SerialSymDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex) const
831 {
832 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
833  checkIndex( rowIndex, colIndex );
834 #endif
835  if ( rowIndex <= colIndex ) {
836  // Accessing upper triangular part of matrix
837  if (upper_)
838  return(values_[colIndex * stride_ + rowIndex]);
839  else
840  return(values_[rowIndex * stride_ + colIndex]);
841  }
842  else {
843  // Accessing lower triangular part of matrix
844  if (upper_)
845  return(values_[rowIndex * stride_ + colIndex]);
846  else
847  return(values_[colIndex * stride_ + rowIndex]);
848  }
849 }
850 
851 //----------------------------------------------------------------------------------------------------
852 // Norm methods
853 //----------------------------------------------------------------------------------------------------
854 
855 template<typename OrdinalType, typename ScalarType>
857 {
858  return(normInf());
859 }
860 
861 template<typename OrdinalType, typename ScalarType>
863 {
864  typedef typename ScalarTraits<ScalarType>::magnitudeType MT;
865 
866  OrdinalType i, j;
867 
868  MT sum, anorm = ScalarTraits<MT>::zero();
869  ScalarType* ptr;
870 
871  if (upper_) {
872  for (j=0; j<numRowCols_; j++) {
873  sum = ScalarTraits<MT>::zero();
874  ptr = values_ + j*stride_;
875  for (i=0; i<j; i++) {
876  sum += ScalarTraits<ScalarType>::magnitude( *ptr++ );
877  }
878  ptr = values_ + j + j*stride_;
879  for (i=j; i<numRowCols_; i++) {
881  ptr += stride_;
882  }
883  anorm = TEUCHOS_MAX( anorm, sum );
884  }
885  }
886  else {
887  for (j=0; j<numRowCols_; j++) {
888  sum = ScalarTraits<MT>::zero();
889  ptr = values_ + j + j*stride_;
890  for (i=j; i<numRowCols_; i++) {
891  sum += ScalarTraits<ScalarType>::magnitude( *ptr++ );
892  }
893  ptr = values_ + j;
894  for (i=0; i<j; i++) {
896  ptr += stride_;
897  }
898  anorm = TEUCHOS_MAX( anorm, sum );
899  }
900  }
901  return(anorm);
902 }
903 
904 template<typename OrdinalType, typename ScalarType>
906 {
907  typedef typename ScalarTraits<ScalarType>::magnitudeType MT;
908 
909  OrdinalType i, j;
910 
911  MT sum = ScalarTraits<MT>::zero(), anorm = ScalarTraits<MT>::zero();
912 
913  if (upper_) {
914  for (j = 0; j < numRowCols_; j++) {
915  for (i = 0; i < j; i++) {
916  sum += ScalarTraits<ScalarType>::magnitude(2.0*values_[i+j*stride_]*values_[i+j*stride_]);
917  }
918  sum += ScalarTraits<ScalarType>::magnitude(values_[j + j*stride_]*values_[j + j*stride_]);
919  }
920  }
921  else {
922  for (j = 0; j < numRowCols_; j++) {
923  sum += ScalarTraits<ScalarType>::magnitude(values_[j + j*stride_]*values_[j + j*stride_]);
924  for (i = j+1; i < numRowCols_; i++) {
925  sum += ScalarTraits<ScalarType>::magnitude(2.0*values_[i+j*stride_]*values_[i+j*stride_]);
926  }
927  }
928  }
930  return(anorm);
931 }
932 
933 //----------------------------------------------------------------------------------------------------
934 // Comparison methods
935 //----------------------------------------------------------------------------------------------------
936 
937 template<typename OrdinalType, typename ScalarType>
939 {
940  bool result = 1;
941  if((numRowCols_ != Operand.numRowCols_)) {
942  result = 0;
943  }
944  else {
945  OrdinalType i, j;
946  for(i = 0; i < numRowCols_; i++) {
947  for(j = 0; j < numRowCols_; j++) {
948  if((*this)(i, j) != Operand(i, j)) {
949  return 0;
950  }
951  }
952  }
953  }
954  return result;
955 }
956 
957 template<typename OrdinalType, typename ScalarType>
959 {
960  return !((*this) == Operand);
961 }
962 
963 //----------------------------------------------------------------------------------------------------
964 // Multiplication method
965 //----------------------------------------------------------------------------------------------------
966 
967 template<typename OrdinalType, typename ScalarType>
968 void SerialSymDenseMatrix<OrdinalType, ScalarType>::scale( const ScalarType alpha )
969 {
970  OrdinalType i, j;
971  ScalarType* ptr;
972 
973  if (upper_) {
974  for (j=0; j<numRowCols_; j++) {
975  ptr = values_ + j*stride_;
976  for (i=0; i<= j; i++) { *ptr = alpha * (*ptr); ptr++; }
977  }
978  }
979  else {
980  for (j=0; j<numRowCols_; j++) {
981  ptr = values_ + j*stride_ + j;
982  for (i=j; i<numRowCols_; i++) { *ptr = alpha * (*ptr); ptr++; }
983  }
984  }
985 }
986 
987 /*
988 template<typename OrdinalType, typename ScalarType>
989 int SerialSymDenseMatrix<OrdinalType, ScalarType>::scale( const SerialSymDenseMatrix<OrdinalType,ScalarType>& A )
990 {
991  OrdinalType i, j;
992  ScalarType* ptr;
993 
994  // Check for compatible dimensions
995  if ((numRowCols_ != A.numRowCols_)) {
996  TEUCHOS_CHK_ERR(-1); // Return error
997  }
998 
999  if (upper_) {
1000  for (j=0; j<numRowCols_; j++) {
1001  ptr = values_ + j*stride_;
1002  for (i=0; i<= j; i++) { *ptr = A(i,j) * (*ptr); ptr++; }
1003  }
1004  }
1005  else {
1006  for (j=0; j<numRowCols_; j++) {
1007  ptr = values_ + j*stride_;
1008  for (i=j; i<numRowCols_; i++) { *ptr = A(i,j) * (*ptr); ptr++; }
1009  }
1010  }
1011 
1012  return(0);
1013 }
1014 */
1015 
1016 template<typename OrdinalType, typename ScalarType>
1017 #ifndef TEUCHOS_HIDE_DEPRECATED_CODE
1019 #else
1020 std::ostream& SerialSymDenseMatrix<OrdinalType, ScalarType>::print(std::ostream& os) const
1021 #endif
1022 {
1023  os << std::endl;
1024  if(valuesCopied_)
1025  os << "Values_copied : yes" << std::endl;
1026  else
1027  os << "Values_copied : no" << std::endl;
1028  os << "Rows / Columns : " << numRowCols_ << std::endl;
1029  os << "LDA : " << stride_ << std::endl;
1030  if (upper_)
1031  os << "Storage: Upper" << std::endl;
1032  else
1033  os << "Storage: Lower" << std::endl;
1034  if(numRowCols_ == 0) {
1035  os << "(matrix is empty, no values to display)" << std::endl;
1036  } else {
1037  for(OrdinalType i = 0; i < numRowCols_; i++) {
1038  for(OrdinalType j = 0; j < numRowCols_; j++){
1039  os << (*this)(i,j) << " ";
1040  }
1041  os << std::endl;
1042  }
1043  }
1044 #ifdef TEUCHOS_HIDE_DEPRECATED_CODE
1045  return os;
1046 #endif
1047 }
1048 
1049 //----------------------------------------------------------------------------------------------------
1050 // Protected methods
1051 //----------------------------------------------------------------------------------------------------
1052 
1053 template<typename OrdinalType, typename ScalarType>
1054 inline void SerialSymDenseMatrix<OrdinalType, ScalarType>::checkIndex( OrdinalType rowIndex, OrdinalType colIndex ) const {
1055  TEUCHOS_TEST_FOR_EXCEPTION(rowIndex < 0 || rowIndex >= numRowCols_, std::out_of_range,
1056  "SerialSymDenseMatrix<T>::checkIndex: "
1057  "Row index " << rowIndex << " out of range [0, "<< numRowCols_ << ")");
1058  TEUCHOS_TEST_FOR_EXCEPTION(colIndex < 0 || colIndex >= numRowCols_, std::out_of_range,
1059  "SerialSymDenseMatrix<T>::checkIndex: "
1060  "Col index " << colIndex << " out of range [0, "<< numRowCols_ << ")");
1061 }
1062 
1063 template<typename OrdinalType, typename ScalarType>
1064 void SerialSymDenseMatrix<OrdinalType, ScalarType>::deleteArrays(void)
1065 {
1066  if (valuesCopied_)
1067  {
1068  delete [] values_;
1069  values_ = 0;
1070  valuesCopied_ = false;
1071  }
1072 }
1073 
1074 template<typename OrdinalType, typename ScalarType>
1075 void SerialSymDenseMatrix<OrdinalType, ScalarType>::copyMat(
1076  bool inputUpper, ScalarType* inputMatrix,
1077  OrdinalType inputStride, OrdinalType numRowCols_in,
1078  bool outputUpper, ScalarType* outputMatrix,
1079  OrdinalType outputStride, OrdinalType startRowCol,
1080  ScalarType alpha
1081  )
1082 {
1083  OrdinalType i, j;
1084  ScalarType* ptr1 = 0;
1085  ScalarType* ptr2 = 0;
1086 
1087  for (j = 0; j < numRowCols_in; j++) {
1088  if (inputUpper == true) {
1089  // The input matrix is upper triangular, start at the beginning of each column.
1090  ptr2 = inputMatrix + (j + startRowCol) * inputStride + startRowCol;
1091  if (outputUpper == true) {
1092  // The output matrix matches the same pattern as the input matrix.
1093  ptr1 = outputMatrix + j*outputStride;
1094  if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
1095  for(i = 0; i <= j; i++) {
1096  *ptr1++ += alpha*(*ptr2++);
1097  }
1098  } else {
1099  for(i = 0; i <= j; i++) {
1100  *ptr1++ = *ptr2++;
1101  }
1102  }
1103  }
1104  else {
1105  // The output matrix has the opposite pattern as the input matrix.
1106  // Copy down across rows of the output matrix, but down columns of the input matrix.
1107  ptr1 = outputMatrix + j;
1108  if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
1109  for(i = 0; i <= j; i++) {
1110  *ptr1 += alpha*(*ptr2++);
1111  ptr1 += outputStride;
1112  }
1113  } else {
1114  for(i = 0; i <= j; i++) {
1115  *ptr1 = *ptr2++;
1116  ptr1 += outputStride;
1117  }
1118  }
1119  }
1120  }
1121  else {
1122  // The input matrix is lower triangular, start at the diagonal of each row.
1123  ptr2 = inputMatrix + (startRowCol+j) * inputStride + startRowCol + j;
1124  if (outputUpper == true) {
1125  // The output matrix has the opposite pattern as the input matrix.
1126  // Copy across rows of the output matrix, but down columns of the input matrix.
1127  ptr1 = outputMatrix + j*outputStride + j;
1128  if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
1129  for(i = j; i < numRowCols_in; i++) {
1130  *ptr1 += alpha*(*ptr2++);
1131  ptr1 += outputStride;
1132  }
1133  } else {
1134  for(i = j; i < numRowCols_in; i++) {
1135  *ptr1 = *ptr2++;
1136  ptr1 += outputStride;
1137  }
1138  }
1139  }
1140  else {
1141  // The output matrix matches the same pattern as the input matrix.
1142  ptr1 = outputMatrix + j*outputStride + j;
1143  if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
1144  for(i = j; i < numRowCols_in; i++) {
1145  *ptr1++ += alpha*(*ptr2++);
1146  }
1147  } else {
1148  for(i = j; i < numRowCols_in; i++) {
1149  *ptr1++ = *ptr2++;
1150  }
1151  }
1152  }
1153  }
1154  }
1155 }
1156 
1157 template<typename OrdinalType, typename ScalarType>
1158 void SerialSymDenseMatrix<OrdinalType, ScalarType>::copyUPLOMat(
1159  bool inputUpper, ScalarType* inputMatrix,
1160  OrdinalType inputStride, OrdinalType inputRows
1161  )
1162 {
1163  OrdinalType i, j;
1164  ScalarType * ptr1 = 0;
1165  ScalarType * ptr2 = 0;
1166 
1167  if (inputUpper) {
1168  for (j=1; j<inputRows; j++) {
1169  ptr1 = inputMatrix + j;
1170  ptr2 = inputMatrix + j*inputStride;
1171  for (i=0; i<j; i++) {
1172  *ptr1 = *ptr2++;
1173  ptr1+=inputStride;
1174  }
1175  }
1176  }
1177  else {
1178  for (i=1; i<inputRows; i++) {
1179  ptr1 = inputMatrix + i;
1180  ptr2 = inputMatrix + i*inputStride;
1181  for (j=0; j<i; j++) {
1182  *ptr2++ = *ptr1;
1183  ptr1+=inputStride;
1184  }
1185  }
1186  }
1187 }
1188 
1189 #ifndef TEUCHOS_HIDE_DEPRECATED_CODE
1190 template<typename OrdinalType, typename ScalarType>
1192 std::ostream& operator<< (std::ostream& os, const Teuchos::SerialSymDenseMatrix<OrdinalType, ScalarType>& obj)
1193 {
1194  obj.print (os);
1195  return os;
1196 }
1197 #endif
1198 
1200 template<typename OrdinalType, typename ScalarType>
1202 public:
1206  : obj(obj_in) {}
1207 };
1208 
1210 template<typename OrdinalType, typename ScalarType>
1211 std::ostream&
1212 operator<<(std::ostream &out,
1214 {
1215  printer.obj.print(out);
1216  return out;
1217 }
1218 
1220 template<typename OrdinalType, typename ScalarType>
1221 SerialSymDenseMatrixPrinter<OrdinalType,ScalarType>
1223 {
1225 }
1226 
1227 
1228 } // namespace Teuchos
1229 
1230 #endif /* _TEUCHOS_SERIALSYMDENSEMATRIX_HPP_ */
bool upper() const
Returns true if upper triangular part of this matrix has and will be used.
Ostream manipulator for SerialSymDenseMatrix.
T magnitudeType
Mandatory typedef for result of magnitude.
SerialSymDenseMatrix()
Default constructor; defines a zero size object.
virtual void print(std::ostream &os) const
Print method. Defines the behavior of the std::ostream &lt;&lt; operator.
Templated interface class to BLAS routines.
SerialBandDenseMatrixPrinter< OrdinalType, ScalarType > printMat(const SerialBandDenseMatrix< OrdinalType, ScalarType > &obj)
Return SerialBandDenseMatrix ostream manipulator Use as:
ScalarType * values() const
Returns the pointer to the ScalarType data array contained in the object.
ScalarType & operator()(OrdinalType rowIndex, OrdinalType colIndex)
Element access method (non-const).
void swap(SerialSymDenseMatrix< OrdinalType, ScalarType > &B)
Swap values between this matrix and incoming matrix.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Macro for throwing an exception with breakpointing to ease debugging.
int shape(OrdinalType numRowsCols)
Set dimensions of a Teuchos::SerialSymDenseMatrix object; init values to zero.
Teuchos header file which uses auto-configuration information to include necessary C++ headers...
Templated BLAS wrapper.
This class creates and provides basic support for symmetric, positive-definite dense matrices of temp...
OrdinalType numRows() const
Returns the row dimension of this matrix.
ScalarTraits< ScalarType >::magnitudeType normInf() const
Returns the Infinity-norm of the matrix.
Object for storing data and providing functionality that is common to all computational classes...
int reshape(OrdinalType numRowsCols)
Reshape a Teuchos::SerialSymDenseMatrix object.
This structure defines some basic traits for a scalar field type.
void setLower()
Specify that the lower triangle of the this matrix should be used.
void setUpper()
Specify that the upper triangle of the this matrix should be used.
Functionality and data that is common to all computational classes.
SerialSymDenseMatrix< OrdinalType, ScalarType > & operator*=(const ScalarType alpha)
Inplace scalar-matrix product A = alpha*A.
bool operator==(const SerialSymDenseMatrix< OrdinalType, ScalarType > &Operand) const
Equality of two matrices.
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero(), bool fullMatrix=false)
Set all values in the matrix to a constant value.
SerialSymDenseMatrix< OrdinalType, ScalarType > & operator+=(const SerialSymDenseMatrix< OrdinalType, ScalarType > &Source)
Add another matrix to this matrix.
ScalarTraits< ScalarType >::magnitudeType normOne() const
Returns the 1-norm of the matrix.
int random(const ScalarType bias=0.1 *Teuchos::ScalarTraits< ScalarType >::one())
Set all values in the active area (upper/lower triangle) of this matrix to be random numbers...
char UPLO() const
Returns character value of UPLO used by LAPACK routines.
virtual ~SerialSymDenseMatrix()
Teuchos::SerialSymDenseMatrix destructor.
static magnitudeType magnitude(T a)
Returns the magnitudeType of the scalar type a.
SerialSymDenseMatrix< OrdinalType, ScalarType > & operator-=(const SerialSymDenseMatrix< OrdinalType, ScalarType > &Source)
Subtract another matrix from this matrix.
OrdinalType ordinalType
Typedef for ordinal type.
SerialSymDenseMatrix< OrdinalType, ScalarType > & operator=(const SerialSymDenseMatrix< OrdinalType, ScalarType > &Source)
Copies values from one matrix to another.
int shapeUninitialized(OrdinalType numRowsCols)
Set dimensions of a Teuchos::SerialSymDenseMatrix object; don&#39;t initialize values.
Defines basic traits for the scalar field type.
static T zero()
Returns representation of zero for this scalar type.
static T random()
Returns a random number (between -one() and +one()) of this scalar type.
OrdinalType stride() const
Returns the stride between the columns of this matrix in memory.
OrdinalType numCols() const
Returns the column dimension of this matrix.
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
Returns the Frobenius-norm of the matrix.
SerialSymDenseMatrix< OrdinalType, ScalarType > & assign(const SerialSymDenseMatrix< OrdinalType, ScalarType > &Source)
Copies values from one matrix to another.
bool operator!=(const SerialSymDenseMatrix< OrdinalType, ScalarType > &Operand) const
Inequality of two matrices.
Teuchos::DataAccess Mode enumerable type.
bool empty() const
Returns whether this matrix is empty.
ScalarType scalarType
Typedef for scalar type.