Teuchos - Trilinos Tools Package  Version of the Day
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 
118 namespace Teuchos {
119 
120 template<typename OrdinalType, typename ScalarType>
121 class SerialSymDenseMatrix : public CompObject, public Object, public BLAS<OrdinalType,ScalarType>
122 {
123  public:
124 
126  typedef OrdinalType ordinalType;
128  typedef ScalarType scalarType;
129 
131 
132 
142 
144 
154  SerialSymDenseMatrix(OrdinalType numRowsCols, bool zeroOut = true);
155 
157 
168  SerialSymDenseMatrix(DataAccess CV, bool upper, ScalarType* values, OrdinalType stride, OrdinalType numRowsCols);
169 
172 
174 
183  SerialSymDenseMatrix(DataAccess CV, const SerialSymDenseMatrix<OrdinalType, ScalarType> &Source, OrdinalType numRowCols, OrdinalType startRowCol=0);
184 
186  virtual ~SerialSymDenseMatrix ();
188 
190 
191 
193 
202  int shape(OrdinalType numRowsCols);
203 
205 
214  int shapeUninitialized(OrdinalType numRowsCols);
215 
217 
227  int reshape(OrdinalType numRowsCols);
228 
230 
232  void setLower();
233 
235 
237  void setUpper();
238 
240 
242 
243 
245 
252 
254 
259 
261 
264  SerialSymDenseMatrix<OrdinalType, ScalarType>& operator= (const ScalarType value) { putScalar(value); return(*this); }
265 
267 
272  int putScalar( const ScalarType value = Teuchos::ScalarTraits<ScalarType>::zero(), bool fullMatrix = false );
273 
275 
277  int random( const ScalarType bias = 0.1*Teuchos::ScalarTraits<ScalarType>::one() );
278 
280 
282 
283 
285 
292  ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex);
293 
295 
302  const ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex) const;
303 
305 
307  ScalarType* values() const { return(values_); }
308 
310 
312 
313 
315  bool upper() const {return(upper_);};
316 
318  char UPLO() const {return(UPLO_);};
320 
322 
323 
325 
332 
334 
338 
340 
344 
346 
348 
349 
351 
355 
357 
361 
363 
365 
366 
368  OrdinalType numRows() const { return(numRowCols_); }
369 
371  OrdinalType numCols() const { return(numRowCols_); }
372 
374  OrdinalType stride() const { return(stride_); }
375 
377  bool empty() const { return(numRowCols_ == 0); }
378 
380 
382 
383 
386 
389 
393 
395 
396  virtual void print(std::ostream& os) const;
398 
400 
401  protected:
402 
403  // In-place scaling of matrix.
404  void scale( const ScalarType alpha );
405 
406  // Copy the values from one matrix to the other.
407  void copyMat(bool inputUpper, ScalarType* inputMatrix, OrdinalType inputStride,
408  OrdinalType numRowCols, bool outputUpper, ScalarType* outputMatrix,
409  OrdinalType outputStride, OrdinalType startRowCol,
410  ScalarType alpha = ScalarTraits<ScalarType>::zero() );
411 
412  // Copy the values from the active triangle of the matrix to the other to make the matrix full symmetric.
413  void copyUPLOMat(bool inputUpper, ScalarType* inputMatrix,
414  OrdinalType inputStride, OrdinalType inputRows);
415 
416  void deleteArrays();
417  void checkIndex( OrdinalType rowIndex, OrdinalType colIndex = 0 ) const;
418  OrdinalType numRowCols_;
419  OrdinalType stride_;
420  bool valuesCopied_;
421  ScalarType* values_;
422  bool upper_;
423  char UPLO_;
424 
425 
426 };
427 
428 //----------------------------------------------------------------------------------------------------
429 // Constructors and Destructor
430 //----------------------------------------------------------------------------------------------------
431 
432 template<typename OrdinalType, typename ScalarType>
434  : CompObject(), numRowCols_(0), stride_(0), valuesCopied_(false), values_(0), upper_(false), UPLO_('L')
435 {}
436 
437 template<typename OrdinalType, typename ScalarType>
439  : CompObject(), numRowCols_(numRowCols_in), stride_(numRowCols_in), valuesCopied_(false), upper_(false), UPLO_('L')
440 {
441  values_ = new ScalarType[stride_*numRowCols_];
442  valuesCopied_ = true;
443  if (zeroOut == true)
445 }
446 
447 template<typename OrdinalType, typename ScalarType>
449  DataAccess CV, bool upper_in, ScalarType* values_in, OrdinalType stride_in, OrdinalType numRowCols_in
450  )
451  : CompObject(), numRowCols_(numRowCols_in), stride_(stride_in), valuesCopied_(false),
452  values_(values_in), upper_(upper_in)
453 {
454  if (upper_)
455  UPLO_ = 'U';
456  else
457  UPLO_ = 'L';
458 
459  if(CV == Copy)
460  {
461  stride_ = numRowCols_;
462  values_ = new ScalarType[stride_*numRowCols_];
463  copyMat(upper_in, values_in, stride_in, numRowCols_, upper_, values_, stride_, 0);
464  valuesCopied_ = true;
465  }
466 }
467 
468 template<typename OrdinalType, typename ScalarType>
469 SerialSymDenseMatrix<OrdinalType, ScalarType>::SerialSymDenseMatrix(const SerialSymDenseMatrix<OrdinalType, ScalarType> &Source) : CompObject(), numRowCols_(Source.numRowCols_), stride_(0), valuesCopied_(true), values_(0), upper_(Source.upper_), UPLO_(Source.UPLO_)
470 {
471  if (!Source.valuesCopied_)
472  {
473  stride_ = Source.stride_;
474  values_ = Source.values_;
475  valuesCopied_ = false;
476  }
477  else
478  {
479  stride_ = numRowCols_;
480  const OrdinalType newsize = stride_ * numRowCols_;
481  if(newsize > 0) {
482  values_ = new ScalarType[newsize];
483  copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0);
484  }
485  else {
486  numRowCols_ = 0; stride_ = 0;
487  valuesCopied_ = false;
488  }
489  }
490 }
491 
492 template<typename OrdinalType, typename ScalarType>
494  DataAccess CV, const SerialSymDenseMatrix<OrdinalType,
495  ScalarType> &Source, OrdinalType numRowCols_in, OrdinalType startRowCol )
496  : CompObject(), numRowCols_(numRowCols_in), stride_(Source.stride_), valuesCopied_(false), upper_(Source.upper_), UPLO_(Source.UPLO_)
497 {
498  if(CV == Copy)
499  {
500  stride_ = numRowCols_in;
501  values_ = new ScalarType[stride_ * numRowCols_in];
502  copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_in, upper_, values_, stride_, startRowCol);
503  valuesCopied_ = true;
504  }
505  else // CV == View
506  {
507  values_ = Source.values_ + (stride_ * startRowCol) + startRowCol;
508  }
509 }
510 
511 template<typename OrdinalType, typename ScalarType>
513 {
514  deleteArrays();
515 }
516 
517 //----------------------------------------------------------------------------------------------------
518 // Shape methods
519 //----------------------------------------------------------------------------------------------------
520 
521 template<typename OrdinalType, typename ScalarType>
523 {
524  deleteArrays(); // Get rid of anything that might be already allocated
525  numRowCols_ = numRowCols_in;
526  stride_ = numRowCols_;
527  values_ = new ScalarType[stride_*numRowCols_];
529  valuesCopied_ = true;
530  return(0);
531 }
532 
533 template<typename OrdinalType, typename ScalarType>
535 {
536  deleteArrays(); // Get rid of anything that might be already allocated
537  numRowCols_ = numRowCols_in;
538  stride_ = numRowCols_;
539  values_ = new ScalarType[stride_*numRowCols_];
540  valuesCopied_ = true;
541  return(0);
542 }
543 
544 template<typename OrdinalType, typename ScalarType>
546 {
547  // Allocate space for new matrix
548  ScalarType* values_tmp = new ScalarType[numRowCols_in * numRowCols_in];
549  ScalarType zero = ScalarTraits<ScalarType>::zero();
550  for(OrdinalType k = 0; k < numRowCols_in * numRowCols_in; k++)
551  {
552  values_tmp[k] = zero;
553  }
554  OrdinalType numRowCols_tmp = TEUCHOS_MIN(numRowCols_, numRowCols_in);
555  if(values_ != 0)
556  {
557  copyMat(upper_, values_, stride_, numRowCols_tmp, upper_, values_tmp, numRowCols_in, 0); // Copy principal submatrix of A to new A
558  }
559  deleteArrays(); // Get rid of anything that might be already allocated
560  numRowCols_ = numRowCols_in;
561  stride_ = numRowCols_;
562  values_ = values_tmp; // Set pointer to new A
563  valuesCopied_ = true;
564  return(0);
565 }
566 
567 //----------------------------------------------------------------------------------------------------
568 // Set methods
569 //----------------------------------------------------------------------------------------------------
570 
571 template<typename OrdinalType, typename ScalarType>
573 {
574  // Do nothing if the matrix is already a lower triangular matrix
575  if (upper_ != false) {
576  copyUPLOMat( true, values_, stride_, numRowCols_ );
577  upper_ = false;
578  UPLO_ = 'L';
579  }
580 }
581 
582 template<typename OrdinalType, typename ScalarType>
584 {
585  // Do nothing if the matrix is already an upper triangular matrix
586  if (upper_ == false) {
587  copyUPLOMat( false, values_, stride_, numRowCols_ );
588  upper_ = true;
589  UPLO_ = 'U';
590  }
591 }
592 
593 template<typename OrdinalType, typename ScalarType>
594 int SerialSymDenseMatrix<OrdinalType, ScalarType>::putScalar( const ScalarType value_in, bool fullMatrix )
595 {
596  // Set each value of the dense matrix to "value".
597  if (fullMatrix == true) {
598  for(OrdinalType j = 0; j < numRowCols_; j++)
599  {
600  for(OrdinalType i = 0; i < numRowCols_; i++)
601  {
602  values_[i + j*stride_] = value_in;
603  }
604  }
605  }
606  // Set the active upper or lower triangular part of the matrix to "value"
607  else {
608  if (upper_) {
609  for(OrdinalType j = 0; j < numRowCols_; j++) {
610  for(OrdinalType i = 0; i <= j; i++) {
611  values_[i + j*stride_] = value_in;
612  }
613  }
614  }
615  else {
616  for(OrdinalType j = 0; j < numRowCols_; j++) {
617  for(OrdinalType i = j; i < numRowCols_; i++) {
618  values_[i + j*stride_] = value_in;
619  }
620  }
621  }
622  }
623  return 0;
624 }
625 
626 template<typename OrdinalType, typename ScalarType>
628 {
630 
631  // Set each value of the dense matrix to a random value.
632  std::vector<MT> diagSum( numRowCols_, Teuchos::ScalarTraits<MT>::zero() );
633  if (upper_) {
634  for(OrdinalType j = 0; j < numRowCols_; j++) {
635  for(OrdinalType i = 0; i < j; i++) {
636  values_[i + j*stride_] = ScalarTraits<ScalarType>::random();
637  diagSum[i] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
638  diagSum[j] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
639  }
640  }
641  }
642  else {
643  for(OrdinalType j = 0; j < numRowCols_; j++) {
644  for(OrdinalType i = j+1; i < numRowCols_; i++) {
645  values_[i + j*stride_] = ScalarTraits<ScalarType>::random();
646  diagSum[i] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
647  diagSum[j] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
648  }
649  }
650  }
651 
652  // Set the diagonal.
653  for(OrdinalType i = 0; i < numRowCols_; i++) {
654  values_[i + i*stride_] = diagSum[i] + bias;
655  }
656  return 0;
657 }
658 
659 template<typename OrdinalType, typename ScalarType>
662 {
663  if(this == &Source)
664  return(*this); // Special case of source same as target
665  if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_)) {
666  upper_ = Source.upper_; // Might have to change the active part of the matrix.
667  return(*this); // Special case of both are views to same data.
668  }
669 
670  // If the source is a view then we will return a view, else we will return a copy.
671  if (!Source.valuesCopied_) {
672  if(valuesCopied_) {
673  // Clean up stored data if this was previously a copy.
674  deleteArrays();
675  }
676  numRowCols_ = Source.numRowCols_;
677  stride_ = Source.stride_;
678  values_ = Source.values_;
679  upper_ = Source.upper_;
680  UPLO_ = Source.UPLO_;
681  }
682  else {
683  // If we were a view, we will now be a copy.
684  if(!valuesCopied_) {
685  numRowCols_ = Source.numRowCols_;
686  stride_ = Source.numRowCols_;
687  upper_ = Source.upper_;
688  UPLO_ = Source.UPLO_;
689  const OrdinalType newsize = stride_ * numRowCols_;
690  if(newsize > 0) {
691  values_ = new ScalarType[newsize];
692  valuesCopied_ = true;
693  }
694  else {
695  values_ = 0;
696  }
697  }
698  // If we were a copy, we will stay a copy.
699  else {
700  if((Source.numRowCols_ <= stride_) && (Source.numRowCols_ == numRowCols_)) { // we don't need to reallocate
701  numRowCols_ = Source.numRowCols_;
702  upper_ = Source.upper_;
703  UPLO_ = Source.UPLO_;
704  }
705  else { // we need to allocate more space (or less space)
706  deleteArrays();
707  numRowCols_ = Source.numRowCols_;
708  stride_ = Source.numRowCols_;
709  upper_ = Source.upper_;
710  UPLO_ = Source.UPLO_;
711  const OrdinalType newsize = stride_ * numRowCols_;
712  if(newsize > 0) {
713  values_ = new ScalarType[newsize];
714  valuesCopied_ = true;
715  }
716  }
717  }
718  copyMat(Source.upper_, Source.values_, Source.stride_, Source.numRowCols_, upper_, values_, stride_, 0);
719  }
720  return(*this);
721 }
722 
723 template<typename OrdinalType, typename ScalarType>
725 {
726  this->scale(alpha);
727  return(*this);
728 }
729 
730 template<typename OrdinalType, typename ScalarType>
732 {
733  // Check for compatible dimensions
734  if ((numRowCols_ != Source.numRowCols_))
735  {
736  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
737  }
738  copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0, ScalarTraits<ScalarType>::one());
739  return(*this);
740 }
741 
742 template<typename OrdinalType, typename ScalarType>
744 {
745  // Check for compatible dimensions
746  if ((numRowCols_ != Source.numRowCols_))
747  {
748  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
749  }
750  copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0, -ScalarTraits<ScalarType>::one());
751  return(*this);
752 }
753 
754 template<typename OrdinalType, typename ScalarType>
756  if(this == &Source)
757  return(*this); // Special case of source same as target
758  if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_)) {
759  upper_ = Source.upper_; // We may have to change the active part of the matrix.
760  return(*this); // Special case of both are views to same data.
761  }
762 
763  // Check for compatible dimensions
764  if ((numRowCols_ != Source.numRowCols_))
765  {
766  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
767  }
768  copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0 );
769  return(*this);
770 }
771 
772 //----------------------------------------------------------------------------------------------------
773 // Accessor methods
774 //----------------------------------------------------------------------------------------------------
775 
776 template<typename OrdinalType, typename ScalarType>
777 inline ScalarType& SerialSymDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex)
778 {
779 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
780  checkIndex( rowIndex, colIndex );
781 #endif
782  if ( rowIndex <= colIndex ) {
783  // Accessing upper triangular part of matrix
784  if (upper_)
785  return(values_[colIndex * stride_ + rowIndex]);
786  else
787  return(values_[rowIndex * stride_ + colIndex]);
788  }
789  else {
790  // Accessing lower triangular part of matrix
791  if (upper_)
792  return(values_[rowIndex * stride_ + colIndex]);
793  else
794  return(values_[colIndex * stride_ + rowIndex]);
795  }
796 }
797 
798 template<typename OrdinalType, typename ScalarType>
799 inline const ScalarType& SerialSymDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex) const
800 {
801 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
802  checkIndex( rowIndex, colIndex );
803 #endif
804  if ( rowIndex <= colIndex ) {
805  // Accessing upper triangular part of matrix
806  if (upper_)
807  return(values_[colIndex * stride_ + rowIndex]);
808  else
809  return(values_[rowIndex * stride_ + colIndex]);
810  }
811  else {
812  // Accessing lower triangular part of matrix
813  if (upper_)
814  return(values_[rowIndex * stride_ + colIndex]);
815  else
816  return(values_[colIndex * stride_ + rowIndex]);
817  }
818 }
819 
820 //----------------------------------------------------------------------------------------------------
821 // Norm methods
822 //----------------------------------------------------------------------------------------------------
823 
824 template<typename OrdinalType, typename ScalarType>
826 {
827  return(normInf());
828 }
829 
830 template<typename OrdinalType, typename ScalarType>
832 {
833  typedef typename ScalarTraits<ScalarType>::magnitudeType MT;
834 
835  OrdinalType i, j;
836 
837  MT sum, anorm = ScalarTraits<MT>::zero();
838  ScalarType* ptr;
839 
840  if (upper_) {
841  for (j=0; j<numRowCols_; j++) {
842  sum = ScalarTraits<MT>::zero();
843  ptr = values_ + j*stride_;
844  for (i=0; i<j; i++) {
845  sum += ScalarTraits<ScalarType>::magnitude( *ptr++ );
846  }
847  ptr = values_ + j + j*stride_;
848  for (i=j; i<numRowCols_; i++) {
850  ptr += stride_;
851  }
852  anorm = TEUCHOS_MAX( anorm, sum );
853  }
854  }
855  else {
856  for (j=0; j<numRowCols_; j++) {
857  sum = ScalarTraits<MT>::zero();
858  ptr = values_ + j + j*stride_;
859  for (i=j; i<numRowCols_; i++) {
860  sum += ScalarTraits<ScalarType>::magnitude( *ptr++ );
861  }
862  ptr = values_ + j;
863  for (i=0; i<j; i++) {
865  ptr += stride_;
866  }
867  anorm = TEUCHOS_MAX( anorm, sum );
868  }
869  }
870  return(anorm);
871 }
872 
873 template<typename OrdinalType, typename ScalarType>
875 {
876  typedef typename ScalarTraits<ScalarType>::magnitudeType MT;
877 
878  OrdinalType i, j;
879 
880  MT sum = ScalarTraits<MT>::zero(), anorm = ScalarTraits<MT>::zero();
881 
882  if (upper_) {
883  for (j = 0; j < numRowCols_; j++) {
884  for (i = 0; i < j; i++) {
885  sum += ScalarTraits<ScalarType>::magnitude(2.0*values_[i+j*stride_]*values_[i+j*stride_]);
886  }
887  sum += ScalarTraits<ScalarType>::magnitude(values_[j + j*stride_]*values_[j + j*stride_]);
888  }
889  }
890  else {
891  for (j = 0; j < numRowCols_; j++) {
892  sum += ScalarTraits<ScalarType>::magnitude(values_[j + j*stride_]*values_[j + j*stride_]);
893  for (i = j+1; i < numRowCols_; i++) {
894  sum += ScalarTraits<ScalarType>::magnitude(2.0*values_[i+j*stride_]*values_[i+j*stride_]);
895  }
896  }
897  }
899  return(anorm);
900 }
901 
902 //----------------------------------------------------------------------------------------------------
903 // Comparison methods
904 //----------------------------------------------------------------------------------------------------
905 
906 template<typename OrdinalType, typename ScalarType>
908 {
909  bool result = 1;
910  if((numRowCols_ != Operand.numRowCols_)) {
911  result = 0;
912  }
913  else {
914  OrdinalType i, j;
915  for(i = 0; i < numRowCols_; i++) {
916  for(j = 0; j < numRowCols_; j++) {
917  if((*this)(i, j) != Operand(i, j)) {
918  return 0;
919  }
920  }
921  }
922  }
923  return result;
924 }
925 
926 template<typename OrdinalType, typename ScalarType>
928 {
929  return !((*this) == Operand);
930 }
931 
932 //----------------------------------------------------------------------------------------------------
933 // Multiplication method
934 //----------------------------------------------------------------------------------------------------
935 
936 template<typename OrdinalType, typename ScalarType>
937 void SerialSymDenseMatrix<OrdinalType, ScalarType>::scale( const ScalarType alpha )
938 {
939  OrdinalType i, j;
940  ScalarType* ptr;
941 
942  if (upper_) {
943  for (j=0; j<numRowCols_; j++) {
944  ptr = values_ + j*stride_;
945  for (i=0; i<= j; i++) { *ptr = alpha * (*ptr); ptr++; }
946  }
947  }
948  else {
949  for (j=0; j<numRowCols_; j++) {
950  ptr = values_ + j*stride_ + j;
951  for (i=j; i<numRowCols_; i++) { *ptr = alpha * (*ptr); ptr++; }
952  }
953  }
954 }
955 
956 /*
957 template<typename OrdinalType, typename ScalarType>
958 int SerialSymDenseMatrix<OrdinalType, ScalarType>::scale( const SerialSymDenseMatrix<OrdinalType,ScalarType>& A )
959 {
960  OrdinalType i, j;
961  ScalarType* ptr;
962 
963  // Check for compatible dimensions
964  if ((numRowCols_ != A.numRowCols_)) {
965  TEUCHOS_CHK_ERR(-1); // Return error
966  }
967 
968  if (upper_) {
969  for (j=0; j<numRowCols_; j++) {
970  ptr = values_ + j*stride_;
971  for (i=0; i<= j; i++) { *ptr = A(i,j) * (*ptr); ptr++; }
972  }
973  }
974  else {
975  for (j=0; j<numRowCols_; j++) {
976  ptr = values_ + j*stride_;
977  for (i=j; i<numRowCols_; i++) { *ptr = A(i,j) * (*ptr); ptr++; }
978  }
979  }
980 
981  return(0);
982 }
983 */
984 
985 template<typename OrdinalType, typename ScalarType>
987 {
988  os << std::endl;
989  if(valuesCopied_)
990  os << "Values_copied : yes" << std::endl;
991  else
992  os << "Values_copied : no" << std::endl;
993  os << "Rows / Columns : " << numRowCols_ << std::endl;
994  os << "LDA : " << stride_ << std::endl;
995  if (upper_)
996  os << "Storage: Upper" << std::endl;
997  else
998  os << "Storage: Lower" << std::endl;
999  if(numRowCols_ == 0) {
1000  os << "(matrix is empty, no values to display)" << std::endl;
1001  } else {
1002  for(OrdinalType i = 0; i < numRowCols_; i++) {
1003  for(OrdinalType j = 0; j < numRowCols_; j++){
1004  os << (*this)(i,j) << " ";
1005  }
1006  os << std::endl;
1007  }
1008  }
1009 }
1010 
1011 //----------------------------------------------------------------------------------------------------
1012 // Protected methods
1013 //----------------------------------------------------------------------------------------------------
1014 
1015 template<typename OrdinalType, typename ScalarType>
1016 inline void SerialSymDenseMatrix<OrdinalType, ScalarType>::checkIndex( OrdinalType rowIndex, OrdinalType colIndex ) const {
1017  TEUCHOS_TEST_FOR_EXCEPTION(rowIndex < 0 || rowIndex >= numRowCols_, std::out_of_range,
1018  "SerialSymDenseMatrix<T>::checkIndex: "
1019  "Row index " << rowIndex << " out of range [0, "<< numRowCols_ << ")");
1020  TEUCHOS_TEST_FOR_EXCEPTION(colIndex < 0 || colIndex >= numRowCols_, std::out_of_range,
1021  "SerialSymDenseMatrix<T>::checkIndex: "
1022  "Col index " << colIndex << " out of range [0, "<< numRowCols_ << ")");
1023 }
1024 
1025 template<typename OrdinalType, typename ScalarType>
1027 {
1028  if (valuesCopied_)
1029  {
1030  delete [] values_;
1031  values_ = 0;
1032  valuesCopied_ = false;
1033  }
1034 }
1035 
1036 template<typename OrdinalType, typename ScalarType>
1038  bool inputUpper, ScalarType* inputMatrix,
1039  OrdinalType inputStride, OrdinalType numRowCols_in,
1040  bool outputUpper, ScalarType* outputMatrix,
1041  OrdinalType outputStride, OrdinalType startRowCol,
1042  ScalarType alpha
1043  )
1044 {
1045  OrdinalType i, j;
1046  ScalarType* ptr1 = 0;
1047  ScalarType* ptr2 = 0;
1048 
1049  for (j = 0; j < numRowCols_in; j++) {
1050  if (inputUpper == true) {
1051  // The input matrix is upper triangular, start at the beginning of each column.
1052  ptr2 = inputMatrix + (j + startRowCol) * inputStride + startRowCol;
1053  if (outputUpper == true) {
1054  // The output matrix matches the same pattern as the input matrix.
1055  ptr1 = outputMatrix + j*outputStride;
1056  if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
1057  for(i = 0; i <= j; i++) {
1058  *ptr1++ += alpha*(*ptr2++);
1059  }
1060  } else {
1061  for(i = 0; i <= j; i++) {
1062  *ptr1++ = *ptr2++;
1063  }
1064  }
1065  }
1066  else {
1067  // The output matrix has the opposite pattern as the input matrix.
1068  // Copy down across rows of the output matrix, but down columns of the input matrix.
1069  ptr1 = outputMatrix + j;
1070  if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
1071  for(i = 0; i <= j; i++) {
1072  *ptr1 += alpha*(*ptr2++);
1073  ptr1 += outputStride;
1074  }
1075  } else {
1076  for(i = 0; i <= j; i++) {
1077  *ptr1 = *ptr2++;
1078  ptr1 += outputStride;
1079  }
1080  }
1081  }
1082  }
1083  else {
1084  // The input matrix is lower triangular, start at the diagonal of each row.
1085  ptr2 = inputMatrix + (startRowCol+j) * inputStride + startRowCol + j;
1086  if (outputUpper == true) {
1087  // The output matrix has the opposite pattern as the input matrix.
1088  // Copy across rows of the output matrix, but down columns of the input matrix.
1089  ptr1 = outputMatrix + j*outputStride + j;
1090  if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
1091  for(i = j; i < numRowCols_in; i++) {
1092  *ptr1 += alpha*(*ptr2++);
1093  ptr1 += outputStride;
1094  }
1095  } else {
1096  for(i = j; i < numRowCols_in; i++) {
1097  *ptr1 = *ptr2++;
1098  ptr1 += outputStride;
1099  }
1100  }
1101  }
1102  else {
1103  // The output matrix matches the same pattern as the input matrix.
1104  ptr1 = outputMatrix + j*outputStride + j;
1105  if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
1106  for(i = j; i < numRowCols_in; i++) {
1107  *ptr1++ += alpha*(*ptr2++);
1108  }
1109  } else {
1110  for(i = j; i < numRowCols_in; i++) {
1111  *ptr1++ = *ptr2++;
1112  }
1113  }
1114  }
1115  }
1116  }
1117 }
1118 
1119 template<typename OrdinalType, typename ScalarType>
1121  bool inputUpper, ScalarType* inputMatrix,
1122  OrdinalType inputStride, OrdinalType inputRows
1123  )
1124 {
1125  OrdinalType i, j;
1126  ScalarType * ptr1 = 0;
1127  ScalarType * ptr2 = 0;
1128 
1129  if (inputUpper) {
1130  for (j=1; j<inputRows; j++) {
1131  ptr1 = inputMatrix + j;
1132  ptr2 = inputMatrix + j*inputStride;
1133  for (i=0; i<j; i++) {
1134  *ptr1 = *ptr2++;
1135  ptr1+=inputStride;
1136  }
1137  }
1138  }
1139  else {
1140  for (i=1; i<inputRows; i++) {
1141  ptr1 = inputMatrix + i;
1142  ptr2 = inputMatrix + i*inputStride;
1143  for (j=0; j<i; j++) {
1144  *ptr2++ = *ptr1;
1145  ptr1+=inputStride;
1146  }
1147  }
1148  }
1149 }
1150 
1151 } // namespace Teuchos
1152 
1153 #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 << operator inherited from the Object class...
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).
#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.
The base Teuchos class.
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.
The Teuchos namespace contains all of the classes, structs and enums used by Teuchos, as well as a number of utility routines.
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.