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  virtual void print(std::ostream& os) const;
406 
408 
409  protected:
410 
411  // In-place scaling of matrix.
412  void scale( const ScalarType alpha );
413 
414  // Copy the values from one matrix to the other.
415  void copyMat(bool inputUpper, ScalarType* inputMatrix, OrdinalType inputStride,
416  OrdinalType numRowCols, bool outputUpper, ScalarType* outputMatrix,
417  OrdinalType outputStride, OrdinalType startRowCol,
418  ScalarType alpha = ScalarTraits<ScalarType>::zero() );
419 
420  // Copy the values from the active triangle of the matrix to the other to make the matrix full symmetric.
421  void copyUPLOMat(bool inputUpper, ScalarType* inputMatrix,
422  OrdinalType inputStride, OrdinalType inputRows);
423 
424  void deleteArrays();
425  void checkIndex( OrdinalType rowIndex, OrdinalType colIndex = 0 ) const;
426  OrdinalType numRowCols_;
427  OrdinalType stride_;
428  bool valuesCopied_;
429  ScalarType* values_;
430  bool upper_;
431  char UPLO_;
432 
433 
434 };
435 
436 //----------------------------------------------------------------------------------------------------
437 // Constructors and Destructor
438 //----------------------------------------------------------------------------------------------------
439 
440 template<typename OrdinalType, typename ScalarType>
442  : CompObject(), BLAS<OrdinalType,ScalarType>(),
443  numRowCols_(0), stride_(0), valuesCopied_(false), values_(0), upper_(false), UPLO_('L')
444 {}
445 
446 template<typename OrdinalType, typename ScalarType>
448  : CompObject(), BLAS<OrdinalType,ScalarType>(),
449  numRowCols_(numRowCols_in), stride_(numRowCols_in), valuesCopied_(false), upper_(false), UPLO_('L')
450 {
451  values_ = new ScalarType[stride_*numRowCols_];
452  valuesCopied_ = true;
453  if (zeroOut == true)
455 }
456 
457 template<typename OrdinalType, typename ScalarType>
459  DataAccess CV, bool upper_in, ScalarType* values_in, OrdinalType stride_in, OrdinalType numRowCols_in
460  )
461  : CompObject(), BLAS<OrdinalType,ScalarType>(),
462  numRowCols_(numRowCols_in), stride_(stride_in), valuesCopied_(false),
463  values_(values_in), upper_(upper_in)
464 {
465  if (upper_)
466  UPLO_ = 'U';
467  else
468  UPLO_ = 'L';
469 
470  if(CV == Copy)
471  {
472  stride_ = numRowCols_;
473  values_ = new ScalarType[stride_*numRowCols_];
474  copyMat(upper_in, values_in, stride_in, numRowCols_, upper_, values_, stride_, 0);
475  valuesCopied_ = true;
476  }
477 }
478 
479 template<typename OrdinalType, typename ScalarType>
481  : CompObject(), BLAS<OrdinalType,ScalarType>(),
482  numRowCols_(Source.numRowCols_), stride_(0), valuesCopied_(true),
483  values_(0), upper_(Source.upper_), UPLO_(Source.UPLO_)
484 {
485  if (!Source.valuesCopied_)
486  {
487  stride_ = Source.stride_;
488  values_ = Source.values_;
489  valuesCopied_ = false;
490  }
491  else
492  {
493  stride_ = numRowCols_;
494  const OrdinalType newsize = stride_ * numRowCols_;
495  if(newsize > 0) {
496  values_ = new ScalarType[newsize];
497  copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0);
498  }
499  else {
500  numRowCols_ = 0; stride_ = 0;
501  valuesCopied_ = false;
502  }
503  }
504 }
505 
506 template<typename OrdinalType, typename ScalarType>
508  DataAccess CV, const SerialSymDenseMatrix<OrdinalType,
509  ScalarType> &Source, OrdinalType numRowCols_in, OrdinalType startRowCol )
510  : CompObject(), BLAS<OrdinalType,ScalarType>(),
511  numRowCols_(numRowCols_in), stride_(Source.stride_), valuesCopied_(false), upper_(Source.upper_), UPLO_(Source.UPLO_)
512 {
513  if(CV == Copy)
514  {
515  stride_ = numRowCols_in;
516  values_ = new ScalarType[stride_ * numRowCols_in];
517  copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_in, upper_, values_, stride_, startRowCol);
518  valuesCopied_ = true;
519  }
520  else // CV == View
521  {
522  values_ = Source.values_ + (stride_ * startRowCol) + startRowCol;
523  }
524 }
525 
526 template<typename OrdinalType, typename ScalarType>
528 {
529  deleteArrays();
530 }
531 
532 //----------------------------------------------------------------------------------------------------
533 // Shape methods
534 //----------------------------------------------------------------------------------------------------
535 
536 template<typename OrdinalType, typename ScalarType>
538 {
539  deleteArrays(); // Get rid of anything that might be already allocated
540  numRowCols_ = numRowCols_in;
541  stride_ = numRowCols_;
542  values_ = new ScalarType[stride_*numRowCols_];
543  putScalar( Teuchos::ScalarTraits<ScalarType>::zero(), true );
544  valuesCopied_ = true;
545  return(0);
546 }
547 
548 template<typename OrdinalType, typename ScalarType>
550 {
551  deleteArrays(); // Get rid of anything that might be already allocated
552  numRowCols_ = numRowCols_in;
553  stride_ = numRowCols_;
554  values_ = new ScalarType[stride_*numRowCols_];
555  valuesCopied_ = true;
556  return(0);
557 }
558 
559 template<typename OrdinalType, typename ScalarType>
561 {
562  // Allocate space for new matrix
563  ScalarType* values_tmp = new ScalarType[numRowCols_in * numRowCols_in];
564  ScalarType zero = ScalarTraits<ScalarType>::zero();
565  for(OrdinalType k = 0; k < numRowCols_in * numRowCols_in; k++)
566  {
567  values_tmp[k] = zero;
568  }
569  OrdinalType numRowCols_tmp = TEUCHOS_MIN(numRowCols_, numRowCols_in);
570  if(values_ != 0)
571  {
572  copyMat(upper_, values_, stride_, numRowCols_tmp, upper_, values_tmp, numRowCols_in, 0); // Copy principal submatrix of A to new A
573  }
574  deleteArrays(); // Get rid of anything that might be already allocated
575  numRowCols_ = numRowCols_in;
576  stride_ = numRowCols_;
577  values_ = values_tmp; // Set pointer to new A
578  valuesCopied_ = true;
579  return(0);
580 }
581 
582 //----------------------------------------------------------------------------------------------------
583 // Set methods
584 //----------------------------------------------------------------------------------------------------
585 
586 template<typename OrdinalType, typename ScalarType>
588 {
589  // Do nothing if the matrix is already a lower triangular matrix
590  if (upper_ != false) {
591  copyUPLOMat( true, values_, stride_, numRowCols_ );
592  upper_ = false;
593  UPLO_ = 'L';
594  }
595 }
596 
597 template<typename OrdinalType, typename ScalarType>
599 {
600  // Do nothing if the matrix is already an upper triangular matrix
601  if (upper_ == false) {
602  copyUPLOMat( false, values_, stride_, numRowCols_ );
603  upper_ = true;
604  UPLO_ = 'U';
605  }
606 }
607 
608 template<typename OrdinalType, typename ScalarType>
609 int SerialSymDenseMatrix<OrdinalType, ScalarType>::putScalar( const ScalarType value_in, bool fullMatrix )
610 {
611  // Set each value of the dense matrix to "value".
612  if (fullMatrix == true) {
613  for(OrdinalType j = 0; j < numRowCols_; j++)
614  {
615  for(OrdinalType i = 0; i < numRowCols_; i++)
616  {
617  values_[i + j*stride_] = value_in;
618  }
619  }
620  }
621  // Set the active upper or lower triangular part of the matrix to "value"
622  else {
623  if (upper_) {
624  for(OrdinalType j = 0; j < numRowCols_; j++) {
625  for(OrdinalType i = 0; i <= j; i++) {
626  values_[i + j*stride_] = value_in;
627  }
628  }
629  }
630  else {
631  for(OrdinalType j = 0; j < numRowCols_; j++) {
632  for(OrdinalType i = j; i < numRowCols_; i++) {
633  values_[i + j*stride_] = value_in;
634  }
635  }
636  }
637  }
638  return 0;
639 }
640 
641 template<typename OrdinalType, typename ScalarType> void
644 {
645  std::swap(values_ , B.values_);
646  std::swap(numRowCols_, B.numRowCols_);
647  std::swap(stride_, B.stride_);
648  std::swap(valuesCopied_, B.valuesCopied_);
649  std::swap(upper_, B.upper_);
650  std::swap(UPLO_, B.UPLO_);
651 }
652 
653 template<typename OrdinalType, typename ScalarType>
655 {
657 
658  // Set each value of the dense matrix to a random value.
659  std::vector<MT> diagSum( numRowCols_, Teuchos::ScalarTraits<MT>::zero() );
660  if (upper_) {
661  for(OrdinalType j = 0; j < numRowCols_; j++) {
662  for(OrdinalType i = 0; i < j; i++) {
663  values_[i + j*stride_] = ScalarTraits<ScalarType>::random();
664  diagSum[i] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
665  diagSum[j] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
666  }
667  }
668  }
669  else {
670  for(OrdinalType j = 0; j < numRowCols_; j++) {
671  for(OrdinalType i = j+1; i < numRowCols_; i++) {
672  values_[i + j*stride_] = ScalarTraits<ScalarType>::random();
673  diagSum[i] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
674  diagSum[j] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
675  }
676  }
677  }
678 
679  // Set the diagonal.
680  for(OrdinalType i = 0; i < numRowCols_; i++) {
681  values_[i + i*stride_] = diagSum[i] + bias;
682  }
683  return 0;
684 }
685 
686 template<typename OrdinalType, typename ScalarType>
689 {
690  if(this == &Source)
691  return(*this); // Special case of source same as target
692  if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_)) {
693  upper_ = Source.upper_; // Might have to change the active part of the matrix.
694  return(*this); // Special case of both are views to same data.
695  }
696 
697  // If the source is a view then we will return a view, else we will return a copy.
698  if (!Source.valuesCopied_) {
699  if(valuesCopied_) {
700  // Clean up stored data if this was previously a copy.
701  deleteArrays();
702  }
703  numRowCols_ = Source.numRowCols_;
704  stride_ = Source.stride_;
705  values_ = Source.values_;
706  upper_ = Source.upper_;
707  UPLO_ = Source.UPLO_;
708  }
709  else {
710  // If we were a view, we will now be a copy.
711  if(!valuesCopied_) {
712  numRowCols_ = Source.numRowCols_;
713  stride_ = Source.numRowCols_;
714  upper_ = Source.upper_;
715  UPLO_ = Source.UPLO_;
716  const OrdinalType newsize = stride_ * numRowCols_;
717  if(newsize > 0) {
718  values_ = new ScalarType[newsize];
719  valuesCopied_ = true;
720  }
721  else {
722  values_ = 0;
723  }
724  }
725  // If we were a copy, we will stay a copy.
726  else {
727  if((Source.numRowCols_ <= stride_) && (Source.numRowCols_ == numRowCols_)) { // we don't need to reallocate
728  numRowCols_ = Source.numRowCols_;
729  upper_ = Source.upper_;
730  UPLO_ = Source.UPLO_;
731  }
732  else { // we need to allocate more space (or less space)
733  deleteArrays();
734  numRowCols_ = Source.numRowCols_;
735  stride_ = Source.numRowCols_;
736  upper_ = Source.upper_;
737  UPLO_ = Source.UPLO_;
738  const OrdinalType newsize = stride_ * numRowCols_;
739  if(newsize > 0) {
740  values_ = new ScalarType[newsize];
741  valuesCopied_ = true;
742  }
743  }
744  }
745  copyMat(Source.upper_, Source.values_, Source.stride_, Source.numRowCols_, upper_, values_, stride_, 0);
746  }
747  return(*this);
748 }
749 
750 template<typename OrdinalType, typename ScalarType>
752 {
753  this->scale(alpha);
754  return(*this);
755 }
756 
757 template<typename OrdinalType, typename ScalarType>
759 {
760  // Check for compatible dimensions
761  if ((numRowCols_ != Source.numRowCols_))
762  {
763  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
764  }
765  copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0, ScalarTraits<ScalarType>::one());
766  return(*this);
767 }
768 
769 template<typename OrdinalType, typename ScalarType>
771 {
772  // Check for compatible dimensions
773  if ((numRowCols_ != Source.numRowCols_))
774  {
775  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
776  }
777  copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0, -ScalarTraits<ScalarType>::one());
778  return(*this);
779 }
780 
781 template<typename OrdinalType, typename ScalarType>
783  if(this == &Source)
784  return(*this); // Special case of source same as target
785  if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_)) {
786  upper_ = Source.upper_; // We may have to change the active part of the matrix.
787  return(*this); // Special case of both are views to same data.
788  }
789 
790  // Check for compatible dimensions
791  if ((numRowCols_ != Source.numRowCols_))
792  {
793  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
794  }
795  copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0 );
796  return(*this);
797 }
798 
799 //----------------------------------------------------------------------------------------------------
800 // Accessor methods
801 //----------------------------------------------------------------------------------------------------
802 
803 template<typename OrdinalType, typename ScalarType>
804 inline ScalarType& SerialSymDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex)
805 {
806 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
807  checkIndex( rowIndex, colIndex );
808 #endif
809  if ( rowIndex <= colIndex ) {
810  // Accessing upper triangular part of matrix
811  if (upper_)
812  return(values_[colIndex * stride_ + rowIndex]);
813  else
814  return(values_[rowIndex * stride_ + colIndex]);
815  }
816  else {
817  // Accessing lower triangular part of matrix
818  if (upper_)
819  return(values_[rowIndex * stride_ + colIndex]);
820  else
821  return(values_[colIndex * stride_ + rowIndex]);
822  }
823 }
824 
825 template<typename OrdinalType, typename ScalarType>
826 inline const ScalarType& SerialSymDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex) const
827 {
828 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
829  checkIndex( rowIndex, colIndex );
830 #endif
831  if ( rowIndex <= colIndex ) {
832  // Accessing upper triangular part of matrix
833  if (upper_)
834  return(values_[colIndex * stride_ + rowIndex]);
835  else
836  return(values_[rowIndex * stride_ + colIndex]);
837  }
838  else {
839  // Accessing lower triangular part of matrix
840  if (upper_)
841  return(values_[rowIndex * stride_ + colIndex]);
842  else
843  return(values_[colIndex * stride_ + rowIndex]);
844  }
845 }
846 
847 //----------------------------------------------------------------------------------------------------
848 // Norm methods
849 //----------------------------------------------------------------------------------------------------
850 
851 template<typename OrdinalType, typename ScalarType>
853 {
854  return(normInf());
855 }
856 
857 template<typename OrdinalType, typename ScalarType>
859 {
860  typedef typename ScalarTraits<ScalarType>::magnitudeType MT;
861 
862  OrdinalType i, j;
863 
864  MT sum, anorm = ScalarTraits<MT>::zero();
865  ScalarType* ptr;
866 
867  if (upper_) {
868  for (j=0; j<numRowCols_; j++) {
869  sum = ScalarTraits<MT>::zero();
870  ptr = values_ + j*stride_;
871  for (i=0; i<j; i++) {
872  sum += ScalarTraits<ScalarType>::magnitude( *ptr++ );
873  }
874  ptr = values_ + j + j*stride_;
875  for (i=j; i<numRowCols_; i++) {
877  ptr += stride_;
878  }
879  anorm = TEUCHOS_MAX( anorm, sum );
880  }
881  }
882  else {
883  for (j=0; j<numRowCols_; j++) {
884  sum = ScalarTraits<MT>::zero();
885  ptr = values_ + j + j*stride_;
886  for (i=j; i<numRowCols_; i++) {
887  sum += ScalarTraits<ScalarType>::magnitude( *ptr++ );
888  }
889  ptr = values_ + j;
890  for (i=0; i<j; i++) {
892  ptr += stride_;
893  }
894  anorm = TEUCHOS_MAX( anorm, sum );
895  }
896  }
897  return(anorm);
898 }
899 
900 template<typename OrdinalType, typename ScalarType>
902 {
903  typedef typename ScalarTraits<ScalarType>::magnitudeType MT;
904 
905  OrdinalType i, j;
906 
907  MT sum = ScalarTraits<MT>::zero(), anorm = ScalarTraits<MT>::zero();
908 
909  if (upper_) {
910  for (j = 0; j < numRowCols_; j++) {
911  for (i = 0; i < j; i++) {
912  sum += ScalarTraits<ScalarType>::magnitude(2.0*values_[i+j*stride_]*values_[i+j*stride_]);
913  }
914  sum += ScalarTraits<ScalarType>::magnitude(values_[j + j*stride_]*values_[j + j*stride_]);
915  }
916  }
917  else {
918  for (j = 0; j < numRowCols_; j++) {
919  sum += ScalarTraits<ScalarType>::magnitude(values_[j + j*stride_]*values_[j + j*stride_]);
920  for (i = j+1; i < numRowCols_; i++) {
921  sum += ScalarTraits<ScalarType>::magnitude(2.0*values_[i+j*stride_]*values_[i+j*stride_]);
922  }
923  }
924  }
926  return(anorm);
927 }
928 
929 //----------------------------------------------------------------------------------------------------
930 // Comparison methods
931 //----------------------------------------------------------------------------------------------------
932 
933 template<typename OrdinalType, typename ScalarType>
935 {
936  bool result = 1;
937  if((numRowCols_ != Operand.numRowCols_)) {
938  result = 0;
939  }
940  else {
941  OrdinalType i, j;
942  for(i = 0; i < numRowCols_; i++) {
943  for(j = 0; j < numRowCols_; j++) {
944  if((*this)(i, j) != Operand(i, j)) {
945  return 0;
946  }
947  }
948  }
949  }
950  return result;
951 }
952 
953 template<typename OrdinalType, typename ScalarType>
955 {
956  return !((*this) == Operand);
957 }
958 
959 //----------------------------------------------------------------------------------------------------
960 // Multiplication method
961 //----------------------------------------------------------------------------------------------------
962 
963 template<typename OrdinalType, typename ScalarType>
964 void SerialSymDenseMatrix<OrdinalType, ScalarType>::scale( const ScalarType alpha )
965 {
966  OrdinalType i, j;
967  ScalarType* ptr;
968 
969  if (upper_) {
970  for (j=0; j<numRowCols_; j++) {
971  ptr = values_ + j*stride_;
972  for (i=0; i<= j; i++) { *ptr = alpha * (*ptr); ptr++; }
973  }
974  }
975  else {
976  for (j=0; j<numRowCols_; j++) {
977  ptr = values_ + j*stride_ + j;
978  for (i=j; i<numRowCols_; i++) { *ptr = alpha * (*ptr); ptr++; }
979  }
980  }
981 }
982 
983 /*
984 template<typename OrdinalType, typename ScalarType>
985 int SerialSymDenseMatrix<OrdinalType, ScalarType>::scale( const SerialSymDenseMatrix<OrdinalType,ScalarType>& A )
986 {
987  OrdinalType i, j;
988  ScalarType* ptr;
989 
990  // Check for compatible dimensions
991  if ((numRowCols_ != A.numRowCols_)) {
992  TEUCHOS_CHK_ERR(-1); // Return error
993  }
994 
995  if (upper_) {
996  for (j=0; j<numRowCols_; j++) {
997  ptr = values_ + j*stride_;
998  for (i=0; i<= j; i++) { *ptr = A(i,j) * (*ptr); ptr++; }
999  }
1000  }
1001  else {
1002  for (j=0; j<numRowCols_; j++) {
1003  ptr = values_ + j*stride_;
1004  for (i=j; i<numRowCols_; i++) { *ptr = A(i,j) * (*ptr); ptr++; }
1005  }
1006  }
1007 
1008  return(0);
1009 }
1010 */
1011 
1012 template<typename OrdinalType, typename ScalarType>
1014 {
1015  os << std::endl;
1016  if(valuesCopied_)
1017  os << "Values_copied : yes" << std::endl;
1018  else
1019  os << "Values_copied : no" << std::endl;
1020  os << "Rows / Columns : " << numRowCols_ << std::endl;
1021  os << "LDA : " << stride_ << std::endl;
1022  if (upper_)
1023  os << "Storage: Upper" << std::endl;
1024  else
1025  os << "Storage: Lower" << std::endl;
1026  if(numRowCols_ == 0) {
1027  os << "(matrix is empty, no values to display)" << std::endl;
1028  } else {
1029  for(OrdinalType i = 0; i < numRowCols_; i++) {
1030  for(OrdinalType j = 0; j < numRowCols_; j++){
1031  os << (*this)(i,j) << " ";
1032  }
1033  os << std::endl;
1034  }
1035  }
1036 }
1037 
1038 //----------------------------------------------------------------------------------------------------
1039 // Protected methods
1040 //----------------------------------------------------------------------------------------------------
1041 
1042 template<typename OrdinalType, typename ScalarType>
1043 inline void SerialSymDenseMatrix<OrdinalType, ScalarType>::checkIndex( OrdinalType rowIndex, OrdinalType colIndex ) const {
1044  TEUCHOS_TEST_FOR_EXCEPTION(rowIndex < 0 || rowIndex >= numRowCols_, std::out_of_range,
1045  "SerialSymDenseMatrix<T>::checkIndex: "
1046  "Row index " << rowIndex << " out of range [0, "<< numRowCols_ << ")");
1047  TEUCHOS_TEST_FOR_EXCEPTION(colIndex < 0 || colIndex >= numRowCols_, std::out_of_range,
1048  "SerialSymDenseMatrix<T>::checkIndex: "
1049  "Col index " << colIndex << " out of range [0, "<< numRowCols_ << ")");
1050 }
1051 
1052 template<typename OrdinalType, typename ScalarType>
1053 void SerialSymDenseMatrix<OrdinalType, ScalarType>::deleteArrays(void)
1054 {
1055  if (valuesCopied_)
1056  {
1057  delete [] values_;
1058  values_ = 0;
1059  valuesCopied_ = false;
1060  }
1061 }
1062 
1063 template<typename OrdinalType, typename ScalarType>
1064 void SerialSymDenseMatrix<OrdinalType, ScalarType>::copyMat(
1065  bool inputUpper, ScalarType* inputMatrix,
1066  OrdinalType inputStride, OrdinalType numRowCols_in,
1067  bool outputUpper, ScalarType* outputMatrix,
1068  OrdinalType outputStride, OrdinalType startRowCol,
1069  ScalarType alpha
1070  )
1071 {
1072  OrdinalType i, j;
1073  ScalarType* ptr1 = 0;
1074  ScalarType* ptr2 = 0;
1075 
1076  for (j = 0; j < numRowCols_in; j++) {
1077  if (inputUpper == true) {
1078  // The input matrix is upper triangular, start at the beginning of each column.
1079  ptr2 = inputMatrix + (j + startRowCol) * inputStride + startRowCol;
1080  if (outputUpper == true) {
1081  // The output matrix matches the same pattern as the input matrix.
1082  ptr1 = outputMatrix + j*outputStride;
1083  if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
1084  for(i = 0; i <= j; i++) {
1085  *ptr1++ += alpha*(*ptr2++);
1086  }
1087  } else {
1088  for(i = 0; i <= j; i++) {
1089  *ptr1++ = *ptr2++;
1090  }
1091  }
1092  }
1093  else {
1094  // The output matrix has the opposite pattern as the input matrix.
1095  // Copy down across rows of the output matrix, but down columns of the input matrix.
1096  ptr1 = outputMatrix + j;
1097  if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
1098  for(i = 0; i <= j; i++) {
1099  *ptr1 += alpha*(*ptr2++);
1100  ptr1 += outputStride;
1101  }
1102  } else {
1103  for(i = 0; i <= j; i++) {
1104  *ptr1 = *ptr2++;
1105  ptr1 += outputStride;
1106  }
1107  }
1108  }
1109  }
1110  else {
1111  // The input matrix is lower triangular, start at the diagonal of each row.
1112  ptr2 = inputMatrix + (startRowCol+j) * inputStride + startRowCol + j;
1113  if (outputUpper == true) {
1114  // The output matrix has the opposite pattern as the input matrix.
1115  // Copy across rows of the output matrix, but down columns of the input matrix.
1116  ptr1 = outputMatrix + j*outputStride + j;
1117  if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
1118  for(i = j; i < numRowCols_in; i++) {
1119  *ptr1 += alpha*(*ptr2++);
1120  ptr1 += outputStride;
1121  }
1122  } else {
1123  for(i = j; i < numRowCols_in; i++) {
1124  *ptr1 = *ptr2++;
1125  ptr1 += outputStride;
1126  }
1127  }
1128  }
1129  else {
1130  // The output matrix matches the same pattern as the input matrix.
1131  ptr1 = outputMatrix + j*outputStride + j;
1132  if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
1133  for(i = j; i < numRowCols_in; i++) {
1134  *ptr1++ += alpha*(*ptr2++);
1135  }
1136  } else {
1137  for(i = j; i < numRowCols_in; i++) {
1138  *ptr1++ = *ptr2++;
1139  }
1140  }
1141  }
1142  }
1143  }
1144 }
1145 
1146 template<typename OrdinalType, typename ScalarType>
1147 void SerialSymDenseMatrix<OrdinalType, ScalarType>::copyUPLOMat(
1148  bool inputUpper, ScalarType* inputMatrix,
1149  OrdinalType inputStride, OrdinalType inputRows
1150  )
1151 {
1152  OrdinalType i, j;
1153  ScalarType * ptr1 = 0;
1154  ScalarType * ptr2 = 0;
1155 
1156  if (inputUpper) {
1157  for (j=1; j<inputRows; j++) {
1158  ptr1 = inputMatrix + j;
1159  ptr2 = inputMatrix + j*inputStride;
1160  for (i=0; i<j; i++) {
1161  *ptr1 = *ptr2++;
1162  ptr1+=inputStride;
1163  }
1164  }
1165  }
1166  else {
1167  for (i=1; i<inputRows; i++) {
1168  ptr1 = inputMatrix + i;
1169  ptr2 = inputMatrix + i*inputStride;
1170  for (j=0; j<i; j++) {
1171  *ptr2++ = *ptr1;
1172  ptr1+=inputStride;
1173  }
1174  }
1175  }
1176 }
1177 
1178 #ifndef TEUCHOS_HIDE_DEPRECATED_CODE
1179 template<typename OrdinalType, typename ScalarType>
1181 std::ostream& operator<< (std::ostream& os, const Teuchos::SerialSymDenseMatrix<OrdinalType, ScalarType>& obj)
1182 {
1183  obj.print (os);
1184  return os;
1185 }
1186 #endif
1187 
1188 } // namespace Teuchos
1189 
1190 #endif /* _TEUCHOS_SERIALSYMDENSEMATRIX_HPP_ */
bool upper() const
Returns true if upper triangular part of this matrix has and will be used.
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.
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.