00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048 #ifndef _TEUCHOS_SERIALDENSEMATRIX_HPP_
00049 #define _TEUCHOS_SERIALDENSEMATRIX_HPP_
00050
00054 #include "Teuchos_CompObject.hpp"
00055 #include "Teuchos_BLAS.hpp"
00056 #include "Teuchos_ScalarTraits.hpp"
00057 #include "Teuchos_DataAccess.hpp"
00058 #include "Teuchos_ConfigDefs.hpp"
00059 #include "Teuchos_TestForException.hpp"
00060 #include "Teuchos_SerialSymDenseMatrix.hpp"
00061
00070 namespace Teuchos {
00071
00072 template<typename OrdinalType, typename ScalarType>
00073 class SerialDenseMatrix : public CompObject, public Object, public BLAS<OrdinalType, ScalarType>
00074 {
00075 public:
00076
00078
00079
00081
00084 SerialDenseMatrix();
00085
00087
00095 SerialDenseMatrix(OrdinalType numRows, OrdinalType numCols, bool zeroOut = true);
00096
00098
00106 SerialDenseMatrix(DataAccess CV, ScalarType* values, OrdinalType stride, OrdinalType numRows, OrdinalType numCols);
00107
00109
00114 SerialDenseMatrix(const SerialDenseMatrix<OrdinalType, ScalarType> &Source, ETransp trans = Teuchos::NO_TRANS);
00115
00117
00129 SerialDenseMatrix(DataAccess CV, const SerialDenseMatrix<OrdinalType, ScalarType> &Source, OrdinalType numRows, OrdinalType numCols, OrdinalType startRow=0, OrdinalType startCol=0);
00130
00132 virtual ~SerialDenseMatrix();
00134
00136
00137
00138
00148 int shape(OrdinalType numRows, OrdinalType numCols);
00149
00151 int shapeUninitialized(OrdinalType numRows, OrdinalType numCols);
00152
00154
00164 int reshape(OrdinalType numRows, OrdinalType numCols);
00165
00166
00168
00170
00171
00173
00179 SerialDenseMatrix<OrdinalType, ScalarType>& operator= (const SerialDenseMatrix<OrdinalType, ScalarType>& Source);
00180
00182
00187 SerialDenseMatrix<OrdinalType, ScalarType>& assign (const SerialDenseMatrix<OrdinalType, ScalarType>& Source);
00188
00190
00194 int putScalar( const ScalarType value = Teuchos::ScalarTraits<ScalarType>::zero() );
00195
00197 int random();
00198
00200
00202
00203
00205
00212 ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex);
00213
00215
00222 const ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex) const;
00223
00225
00232 ScalarType* operator [] (OrdinalType colIndex);
00233
00235
00242 const ScalarType* operator [] (OrdinalType colIndex) const;
00243
00245
00246 ScalarType* values() const { return(values_); }
00247
00249
00251
00252
00254
00257 SerialDenseMatrix<OrdinalType, ScalarType>& operator+= (const SerialDenseMatrix<OrdinalType, ScalarType>& Source);
00258
00260
00263 SerialDenseMatrix<OrdinalType, ScalarType>& operator-= (const SerialDenseMatrix<OrdinalType, ScalarType>& Source);
00264
00266
00269 SerialDenseMatrix<OrdinalType, ScalarType>& operator*= (const ScalarType alpha);
00270
00272
00276 int scale ( const ScalarType alpha );
00277
00279
00285 int scale ( const SerialDenseMatrix<OrdinalType, ScalarType>& A );
00286
00288
00302 int multiply (ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix<OrdinalType, ScalarType> &A, const SerialDenseMatrix<OrdinalType, ScalarType> &B, ScalarType beta);
00303
00305
00316 int multiply (ESide sideA, ScalarType alpha, const SerialSymDenseMatrix<OrdinalType, ScalarType> &A, const SerialDenseMatrix<OrdinalType, ScalarType> &B, ScalarType beta);
00317
00319
00321
00322
00324
00327 bool operator== (const SerialDenseMatrix<OrdinalType, ScalarType> &Operand);
00328
00330
00333 bool operator!= (const SerialDenseMatrix<OrdinalType, ScalarType> &Operand);
00334
00336
00338
00339
00341 OrdinalType numRows() const { return(numRows_); }
00342
00344 OrdinalType numCols() const { return(numCols_); }
00345
00347 OrdinalType stride() const { return(stride_); }
00348
00350 bool empty() const { return(numRows_ == 0 || numCols_ == 0); }
00352
00354
00355
00357 typename ScalarTraits<ScalarType>::magnitudeType normOne() const;
00358
00360 typename ScalarTraits<ScalarType>::magnitudeType normInf() const;
00361
00363 typename ScalarTraits<ScalarType>::magnitudeType normFrobenius() const;
00365
00367
00368
00369 virtual void print(std::ostream& os) const;
00370
00372 protected:
00373 void copyMat(ScalarType* inputMatrix, OrdinalType strideInput,
00374 OrdinalType numRows, OrdinalType numCols, ScalarType* outputMatrix,
00375 OrdinalType strideOutput, OrdinalType startRow, OrdinalType startCol,
00376 ScalarType alpha = ScalarTraits<ScalarType>::zero() );
00377 void deleteArrays();
00378 void checkIndex( OrdinalType rowIndex, OrdinalType colIndex = 0 ) const;
00379 OrdinalType numRows_;
00380 OrdinalType numCols_;
00381 OrdinalType stride_;
00382 bool valuesCopied_;
00383 ScalarType* values_;
00384
00385 };
00386
00387
00388
00389
00390
00391 template<typename OrdinalType, typename ScalarType>
00392 SerialDenseMatrix<OrdinalType, ScalarType>::SerialDenseMatrix()
00393 : CompObject(), Object("Teuchos::SerialDenseMatrix"), numRows_(0), numCols_(0), stride_(0), valuesCopied_(false), values_(0)
00394 {}
00395
00396 template<typename OrdinalType, typename ScalarType>
00397 SerialDenseMatrix<OrdinalType, ScalarType>::SerialDenseMatrix(
00398 OrdinalType numRows_in, OrdinalType numCols_in, bool zeroOut
00399 )
00400 : CompObject(), Object("Teuchos::SerialDenseMatrix"), numRows_(numRows_in), numCols_(numCols_in), stride_(numRows_in)
00401 {
00402 values_ = new ScalarType[stride_*numCols_];
00403 valuesCopied_ = true;
00404 if (zeroOut == true)
00405 putScalar();
00406 }
00407
00408 template<typename OrdinalType, typename ScalarType>
00409 SerialDenseMatrix<OrdinalType, ScalarType>::SerialDenseMatrix(
00410 DataAccess CV, ScalarType* values_in, OrdinalType stride_in, OrdinalType numRows_in,
00411 OrdinalType numCols_in
00412 )
00413 : CompObject(), Object("Teuchos::SerialDenseMatrix"), numRows_(numRows_in), numCols_(numCols_in), stride_(stride_in),
00414 valuesCopied_(false), values_(values_in)
00415 {
00416 if(CV == Copy)
00417 {
00418 stride_ = numRows_;
00419 values_ = new ScalarType[stride_*numCols_];
00420 copyMat(values_in, stride_in, numRows_, numCols_, values_, stride_, 0, 0, false);
00421 valuesCopied_ = true;
00422 }
00423 }
00424
00425 template<typename OrdinalType, typename ScalarType>
00426 SerialDenseMatrix<OrdinalType, ScalarType>::SerialDenseMatrix(const SerialDenseMatrix<OrdinalType, ScalarType> &Source, ETransp trans) : CompObject(), Object("Teuchos::SerialDenseMatrix"), numRows_(0), numCols_(0), stride_(0), valuesCopied_(true), values_(0)
00427 {
00428 if ( trans == Teuchos::NO_TRANS )
00429 {
00430 numRows_ = Source.numRows_;
00431 numCols_ = Source.numCols_;
00432 stride_ = numRows_;
00433 values_ = new ScalarType[stride_*numCols_];
00434 copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0, false);
00435 }
00436 else if ( trans == Teuchos::CONJ_TRANS && ScalarTraits<ScalarType>::isComplex )
00437 {
00438 numRows_ = Source.numCols_;
00439 numCols_ = Source.numRows_;
00440 stride_ = numRows_;
00441 values_ = new ScalarType[stride_*numCols_];
00442 for (OrdinalType j=0; j<numCols_; j++) {
00443 for (OrdinalType i=0; i<numRows_; i++) {
00444 values_[j*stride_ + i] = Teuchos::ScalarTraits<ScalarType>::conjugate(Source.values_[i*Source.stride_ + j]);
00445 }
00446 }
00447 }
00448 else
00449 {
00450 numRows_ = Source.numCols_;
00451 numCols_ = Source.numRows_;
00452 stride_ = numRows_;
00453 values_ = new ScalarType[stride_*numCols_];
00454 for (OrdinalType j=0; j<numCols_; j++) {
00455 for (OrdinalType i=0; i<numRows_; i++) {
00456 values_[j*stride_ + i] = Source.values_[i*Source.stride_ + j];
00457 }
00458 }
00459 }
00460 }
00461
00462 template<typename OrdinalType, typename ScalarType>
00463 SerialDenseMatrix<OrdinalType, ScalarType>::SerialDenseMatrix(
00464 DataAccess CV, const SerialDenseMatrix<OrdinalType, ScalarType> &Source,
00465 OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType startRow,
00466 OrdinalType startCol
00467 )
00468 : CompObject(), Object("Teuchos::SerialDenseMatrix"), numRows_(numRows_in), numCols_(numCols_in), stride_(Source.stride_),
00469 valuesCopied_(false), values_(Source.values_)
00470 {
00471 if(CV == Copy)
00472 {
00473 stride_ = numRows_in;
00474 values_ = new ScalarType[stride_ * numCols_in];
00475 copyMat(Source.values_, Source.stride_, numRows_in, numCols_in, values_, stride_, startRow, startCol, false);
00476 valuesCopied_ = true;
00477 }
00478 else
00479 {
00480 values_ = values_ + (stride_ * startCol) + startRow;
00481 }
00482 }
00483
00484 template<typename OrdinalType, typename ScalarType>
00485 SerialDenseMatrix<OrdinalType, ScalarType>::~SerialDenseMatrix()
00486 {
00487 deleteArrays();
00488 }
00489
00490
00491
00492
00493
00494 template<typename OrdinalType, typename ScalarType>
00495 int SerialDenseMatrix<OrdinalType, ScalarType>::shape(
00496 OrdinalType numRows_in, OrdinalType numCols_in
00497 )
00498 {
00499 deleteArrays();
00500 numRows_ = numRows_in;
00501 numCols_ = numCols_in;
00502 stride_ = numRows_;
00503 values_ = new ScalarType[stride_*numCols_];
00504 putScalar();
00505 valuesCopied_ = true;
00506 return(0);
00507 }
00508
00509 template<typename OrdinalType, typename ScalarType>
00510 int SerialDenseMatrix<OrdinalType, ScalarType>::shapeUninitialized(
00511 OrdinalType numRows_in, OrdinalType numCols_in
00512 )
00513 {
00514 deleteArrays();
00515 numRows_ = numRows_in;
00516 numCols_ = numCols_in;
00517 stride_ = numRows_;
00518 values_ = new ScalarType[stride_*numCols_];
00519 valuesCopied_ = true;
00520 return(0);
00521 }
00522
00523 template<typename OrdinalType, typename ScalarType>
00524 int SerialDenseMatrix<OrdinalType, ScalarType>::reshape(
00525 OrdinalType numRows_in, OrdinalType numCols_in
00526 )
00527 {
00528
00529 ScalarType* values_tmp = new ScalarType[numRows_in * numCols_in];
00530 ScalarType zero = ScalarTraits<ScalarType>::zero();
00531 for(OrdinalType k = 0; k < numRows_in * numCols_in; k++)
00532 {
00533 values_tmp[k] = zero;
00534 }
00535 OrdinalType numRows_tmp = TEUCHOS_MIN(numRows_, numRows_in);
00536 OrdinalType numCols_tmp = TEUCHOS_MIN(numCols_, numCols_in);
00537 if(values_ != 0)
00538 {
00539 copyMat(values_, stride_, numRows_tmp, numCols_tmp, values_tmp,
00540 numRows_in, 0, 0, false);
00541 }
00542 deleteArrays();
00543 numRows_ = numRows_in;
00544 numCols_ = numCols_in;
00545 stride_ = numRows_;
00546 values_ = values_tmp;
00547 valuesCopied_ = true;
00548 return(0);
00549 }
00550
00551
00552
00553
00554
00555 template<typename OrdinalType, typename ScalarType>
00556 int SerialDenseMatrix<OrdinalType, ScalarType>::putScalar( const ScalarType value_in )
00557 {
00558
00559 for(OrdinalType j = 0; j < numCols_; j++)
00560 {
00561 for(OrdinalType i = 0; i < numRows_; i++)
00562 {
00563 values_[i + j*stride_] = value_in;
00564 }
00565 }
00566 return 0;
00567 }
00568
00569 template<typename OrdinalType, typename ScalarType>
00570 int SerialDenseMatrix<OrdinalType, ScalarType>::random()
00571 {
00572
00573 for(OrdinalType j = 0; j < numCols_; j++)
00574 {
00575 for(OrdinalType i = 0; i < numRows_; i++)
00576 {
00577 values_[i + j*stride_] = ScalarTraits<ScalarType>::random();
00578 }
00579 }
00580 return 0;
00581 }
00582
00583 template<typename OrdinalType, typename ScalarType>
00584 SerialDenseMatrix<OrdinalType,ScalarType>&
00585 SerialDenseMatrix<OrdinalType, ScalarType>::operator=(
00586 const SerialDenseMatrix<OrdinalType,ScalarType>& Source
00587 )
00588 {
00589 if(this == &Source)
00590 return(*this);
00591 if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
00592 return(*this);
00593
00594
00595 if (!Source.valuesCopied_) {
00596 if(valuesCopied_) {
00597
00598 deleteArrays();
00599 }
00600 numRows_ = Source.numRows_;
00601 numCols_ = Source.numCols_;
00602 stride_ = Source.stride_;
00603 values_ = Source.values_;
00604 }
00605 else {
00606
00607 if(!valuesCopied_) {
00608 numRows_ = Source.numRows_;
00609 numCols_ = Source.numCols_;
00610 stride_ = Source.numRows_;
00611 const OrdinalType newsize = stride_ * numCols_;
00612 if(newsize > 0) {
00613 values_ = new ScalarType[newsize];
00614 valuesCopied_ = true;
00615 }
00616 else {
00617 values_ = 0;
00618 }
00619 }
00620
00621 else {
00622 if((Source.numRows_ <= stride_) && (Source.numCols_ == numCols_)) {
00623 numRows_ = Source.numRows_;
00624 numCols_ = Source.numCols_;
00625 }
00626 else {
00627 deleteArrays();
00628 numRows_ = Source.numRows_;
00629 numCols_ = Source.numCols_;
00630 stride_ = Source.numRows_;
00631 const OrdinalType newsize = stride_ * numCols_;
00632 if(newsize > 0) {
00633 values_ = new ScalarType[newsize];
00634 valuesCopied_ = true;
00635 }
00636 }
00637 }
00638 copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0, false);
00639 }
00640 return(*this);
00641 }
00642
00643 template<typename OrdinalType, typename ScalarType>
00644 SerialDenseMatrix<OrdinalType, ScalarType>& SerialDenseMatrix<OrdinalType, ScalarType>::operator+= (const SerialDenseMatrix<OrdinalType,ScalarType>& Source )
00645 {
00646
00647 if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_))
00648 {
00649 TEUCHOS_CHK_REF(*this);
00650 }
00651 copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0, 1.0);
00652 return(*this);
00653 }
00654
00655 template<typename OrdinalType, typename ScalarType>
00656 SerialDenseMatrix<OrdinalType, ScalarType>& SerialDenseMatrix<OrdinalType, ScalarType>::operator-= (const SerialDenseMatrix<OrdinalType,ScalarType>& Source )
00657 {
00658
00659 if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_))
00660 {
00661 TEUCHOS_CHK_REF(*this);
00662 }
00663 copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0, -1.0);
00664 return(*this);
00665 }
00666
00667 template<typename OrdinalType, typename ScalarType>
00668 SerialDenseMatrix<OrdinalType,ScalarType>& SerialDenseMatrix<OrdinalType, ScalarType>::assign (const SerialDenseMatrix<OrdinalType,ScalarType>& Source) {
00669 if(this == &Source)
00670 return(*this);
00671 if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
00672 return(*this);
00673
00674
00675 if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_))
00676 {
00677 TEUCHOS_CHK_REF(*this);
00678 }
00679 copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0 );
00680 return(*this);
00681 }
00682
00683
00684
00685
00686
00687 template<typename OrdinalType, typename ScalarType>
00688 inline ScalarType& SerialDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex)
00689 {
00690 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
00691 checkIndex( rowIndex, colIndex );
00692 #endif
00693 return(values_[colIndex * stride_ + rowIndex]);
00694 }
00695
00696 template<typename OrdinalType, typename ScalarType>
00697 inline const ScalarType& SerialDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex) const
00698 {
00699 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
00700 checkIndex( rowIndex, colIndex );
00701 #endif
00702 return(values_[colIndex * stride_ + rowIndex]);
00703 }
00704
00705 template<typename OrdinalType, typename ScalarType>
00706 inline const ScalarType* SerialDenseMatrix<OrdinalType, ScalarType>::operator [] (OrdinalType colIndex) const
00707 {
00708 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
00709 checkIndex( 0, colIndex );
00710 #endif
00711 return(values_ + colIndex * stride_);
00712 }
00713
00714 template<typename OrdinalType, typename ScalarType>
00715 inline ScalarType* SerialDenseMatrix<OrdinalType, ScalarType>::operator [] (OrdinalType colIndex)
00716 {
00717 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
00718 checkIndex( 0, colIndex );
00719 #endif
00720 return(values_ + colIndex * stride_);
00721 }
00722
00723
00724
00725
00726
00727 template<typename OrdinalType, typename ScalarType>
00728 typename ScalarTraits<ScalarType>::magnitudeType SerialDenseMatrix<OrdinalType, ScalarType>::normOne() const
00729 {
00730 OrdinalType i, j;
00731 typename ScalarTraits<ScalarType>::magnitudeType anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00732 typename ScalarTraits<ScalarType>::magnitudeType absSum = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00733 ScalarType* ptr;
00734 for(j = 0; j < numCols_; j++)
00735 {
00736 ScalarType sum = 0;
00737 ptr = values_ + j * stride_;
00738 for(i = 0; i < numRows_; i++)
00739 {
00740 sum += ScalarTraits<ScalarType>::magnitude(*ptr++);
00741 }
00742 absSum = ScalarTraits<ScalarType>::magnitude(sum);
00743 if(absSum > anorm)
00744 {
00745 anorm = absSum;
00746 }
00747 }
00748 updateFlops(numRows_ * numCols_);
00749 return(anorm);
00750 }
00751
00752 template<typename OrdinalType, typename ScalarType>
00753 typename ScalarTraits<ScalarType>::magnitudeType SerialDenseMatrix<OrdinalType, ScalarType>::normInf() const
00754 {
00755 OrdinalType i, j;
00756 typename ScalarTraits<ScalarType>::magnitudeType sum, anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00757
00758 for (i = 0; i < numRows_; i++) {
00759 sum = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00760 for (j=0; j< numCols_; j++) {
00761 sum += ScalarTraits<ScalarType>::magnitude(*(values_+i+j*stride_));
00762 }
00763 anorm = TEUCHOS_MAX( anorm, sum );
00764 }
00765 updateFlops(numRows_ * numCols_);
00766 return(anorm);
00767 }
00768
00769 template<typename OrdinalType, typename ScalarType>
00770 typename ScalarTraits<ScalarType>::magnitudeType SerialDenseMatrix<OrdinalType, ScalarType>::normFrobenius() const
00771 {
00772 OrdinalType i, j;
00773 typename ScalarTraits<ScalarType>::magnitudeType anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00774 for (j = 0; j < numCols_; j++) {
00775 for (i = 0; i < numRows_; i++) {
00776 anorm += ScalarTraits<ScalarType>::magnitude(values_[i+j*stride_]*values_[i+j*stride_]);
00777 }
00778 }
00779 anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::squareroot(anorm));
00780 updateFlops(numRows_ * numCols_);
00781 return(anorm);
00782 }
00783
00784
00785
00786
00787
00788 template<typename OrdinalType, typename ScalarType>
00789 bool SerialDenseMatrix<OrdinalType, ScalarType>::operator== (const SerialDenseMatrix<OrdinalType, ScalarType> &Operand)
00790 {
00791 bool result = 1;
00792 if((numRows_ != Operand.numRows_) || (numCols_ != Operand.numCols_))
00793 {
00794 result = 0;
00795 }
00796 else
00797 {
00798 OrdinalType i, j;
00799 for(i = 0; i < numRows_; i++)
00800 {
00801 for(j = 0; j < numCols_; j++)
00802 {
00803 if((*this)(i, j) != Operand(i, j))
00804 {
00805 return 0;
00806 }
00807 }
00808 }
00809 }
00810 return result;
00811 }
00812
00813 template<typename OrdinalType, typename ScalarType>
00814 bool SerialDenseMatrix<OrdinalType, ScalarType>::operator!= (const SerialDenseMatrix<OrdinalType, ScalarType> &Operand)
00815 {
00816 return !((*this) == Operand);
00817 }
00818
00819
00820
00821
00822
00823 template<typename OrdinalType, typename ScalarType>
00824 SerialDenseMatrix<OrdinalType, ScalarType>& SerialDenseMatrix<OrdinalType, ScalarType>::operator*= (const ScalarType alpha )
00825 {
00826 this->scale( alpha );
00827 return(*this);
00828 }
00829
00830 template<typename OrdinalType, typename ScalarType>
00831 int SerialDenseMatrix<OrdinalType, ScalarType>::scale( const ScalarType alpha )
00832 {
00833 OrdinalType i, j;
00834 ScalarType* ptr;
00835
00836 for (j=0; j<numCols_; j++) {
00837 ptr = values_ + j*stride_;
00838 for (i=0; i<numRows_; i++) { *ptr = alpha * (*ptr); ptr++; }
00839 }
00840 updateFlops( numRows_*numCols_ );
00841 return(0);
00842 }
00843
00844 template<typename OrdinalType, typename ScalarType>
00845 int SerialDenseMatrix<OrdinalType, ScalarType>::scale( const SerialDenseMatrix<OrdinalType,ScalarType>& A )
00846 {
00847 OrdinalType i, j;
00848 ScalarType* ptr;
00849
00850
00851 if ((numRows_ != A.numRows_) || (numCols_ != A.numCols_))
00852 {
00853 TEUCHOS_CHK_ERR(-1);
00854 }
00855 for (j=0; j<numCols_; j++) {
00856 ptr = values_ + j*stride_;
00857 for (i=0; i<numRows_; i++) { *ptr = A(i,j) * (*ptr); ptr++; }
00858 }
00859 updateFlops( numRows_*numCols_ );
00860 return(0);
00861 }
00862
00863 template<typename OrdinalType, typename ScalarType>
00864 int SerialDenseMatrix<OrdinalType, ScalarType>::multiply(ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix<OrdinalType, ScalarType> &A, const SerialDenseMatrix<OrdinalType, ScalarType> &B, ScalarType beta)
00865 {
00866
00867 OrdinalType A_nrows = (ETranspChar[transa]!='N') ? A.numCols() : A.numRows();
00868 OrdinalType A_ncols = (ETranspChar[transa]!='N') ? A.numRows() : A.numCols();
00869 OrdinalType B_nrows = (ETranspChar[transb]!='N') ? B.numCols() : B.numRows();
00870 OrdinalType B_ncols = (ETranspChar[transb]!='N') ? B.numRows() : B.numCols();
00871 if ((numRows_ != A_nrows) || (A_ncols != B_nrows) || (numCols_ != B_ncols))
00872 {
00873 TEUCHOS_CHK_ERR(-1);
00874 }
00875
00876 this->GEMM(transa, transb, numRows_, numCols_, A_ncols, alpha, A.values(), A.stride(), B.values(), B.stride(), beta, values_, stride_);
00877 double nflops = 2 * numRows_;
00878 nflops *= numCols_;
00879 nflops *= A_ncols;
00880 updateFlops(nflops);
00881 return(0);
00882 }
00883
00884 template<typename OrdinalType, typename ScalarType>
00885 int SerialDenseMatrix<OrdinalType, ScalarType>::multiply (ESide sideA, ScalarType alpha, const SerialSymDenseMatrix<OrdinalType, ScalarType> &A, const SerialDenseMatrix<OrdinalType, ScalarType> &B, ScalarType beta)
00886 {
00887
00888 OrdinalType A_nrows = A.numRows(), A_ncols = A.numCols();
00889 OrdinalType B_nrows = B.numRows(), B_ncols = B.numCols();
00890
00891 if (ESideChar[sideA]=='L') {
00892 if ((numRows_ != A_nrows) || (A_ncols != B_nrows) || (numCols_ != B_ncols)) {
00893 TEUCHOS_CHK_ERR(-1);
00894 }
00895 } else {
00896 if ((numRows_ != B_nrows) || (B_ncols != A_nrows) || (numCols_ != A_ncols)) {
00897 TEUCHOS_CHK_ERR(-1);
00898 }
00899 }
00900
00901
00902 EUplo uplo = (A.upper() ? Teuchos::UPPER_TRI : Teuchos::LOWER_TRI);
00903 this->SYMM(sideA, uplo, numRows_, numCols_, alpha, A.values(), A.stride(), B.values(), B.stride(), beta, values_, stride_);
00904 double nflops = 2 * numRows_;
00905 nflops *= numCols_;
00906 nflops *= A_ncols;
00907 updateFlops(nflops);
00908 return(0);
00909 }
00910
00911 template<typename OrdinalType, typename ScalarType>
00912 void SerialDenseMatrix<OrdinalType, ScalarType>::print(std::ostream& os) const
00913 {
00914 os << std::endl;
00915 if(valuesCopied_)
00916 os << "Values_copied : yes" << std::endl;
00917 else
00918 os << "Values_copied : no" << std::endl;
00919 os << "Rows : " << numRows_ << std::endl;
00920 os << "Columns : " << numCols_ << std::endl;
00921 os << "LDA : " << stride_ << std::endl;
00922 if(numRows_ == 0 || numCols_ == 0) {
00923 os << "(matrix is empty, no values to display)" << std::endl;
00924 } else {
00925 for(OrdinalType i = 0; i < numRows_; i++) {
00926 for(OrdinalType j = 0; j < numCols_; j++){
00927 os << (*this)(i,j) << " ";
00928 }
00929 os << std::endl;
00930 }
00931 }
00932 }
00933
00934
00935
00936
00937
00938 template<typename OrdinalType, typename ScalarType>
00939 inline void SerialDenseMatrix<OrdinalType, ScalarType>::checkIndex( OrdinalType rowIndex, OrdinalType colIndex ) const {
00940 TEST_FOR_EXCEPTION(rowIndex < 0 || rowIndex >= numRows_, std::out_of_range,
00941 "SerialDenseMatrix<T>::checkIndex: "
00942 "Row index " << rowIndex << " out of range [0, "<< numRows_ << ")");
00943 TEST_FOR_EXCEPTION(colIndex < 0 || colIndex >= numCols_, std::out_of_range,
00944 "SerialDenseMatrix<T>::checkIndex: "
00945 "Col index " << colIndex << " out of range [0, "<< numCols_ << ")");
00946 }
00947
00948 template<typename OrdinalType, typename ScalarType>
00949 void SerialDenseMatrix<OrdinalType, ScalarType>::deleteArrays(void)
00950 {
00951 if (valuesCopied_)
00952 {
00953 delete [] values_;
00954 values_ = 0;
00955 valuesCopied_ = false;
00956 }
00957 }
00958
00959 template<typename OrdinalType, typename ScalarType>
00960 void SerialDenseMatrix<OrdinalType, ScalarType>::copyMat(
00961 ScalarType* inputMatrix, OrdinalType strideInput, OrdinalType numRows_in,
00962 OrdinalType numCols_in, ScalarType* outputMatrix, OrdinalType strideOutput,
00963 OrdinalType startRow, OrdinalType startCol, ScalarType alpha
00964 )
00965 {
00966 OrdinalType i, j;
00967 ScalarType* ptr1 = 0;
00968 ScalarType* ptr2 = 0;
00969 for(j = 0; j < numCols_in; j++) {
00970 ptr1 = outputMatrix + (j * strideOutput);
00971 ptr2 = inputMatrix + (j + startCol) * strideInput + startRow;
00972 if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
00973 for(i = 0; i < numRows_in; i++)
00974 {
00975 *ptr1++ += alpha*(*ptr2++);
00976 }
00977 } else {
00978 for(i = 0; i < numRows_in; i++)
00979 {
00980 *ptr1++ = *ptr2++;
00981 }
00982 }
00983 }
00984 }
00985
00986 }
00987
00988
00989 #endif