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
00038 #ifndef MAT_MatrixSymmetric
00039 #define MAT_MatrixSymmetric
00040
00041 #include <algorithm>
00042
00043 #include "MatrixBase.h"
00044 #include "Interval.h"
00045 #include "LanczosLargestMagnitudeEig.h"
00046 #include "mat_utils.h"
00047 #include "truncation.h"
00048
00049 namespace mat {
00067 template<typename Treal, typename Tmatrix>
00068 class MatrixSymmetric : public MatrixBase<Treal, Tmatrix> {
00069 public:
00070 typedef VectorGeneral<Treal, typename Tmatrix::VectorType> VectorType;
00071 typedef Treal real;
00072
00073 MatrixSymmetric()
00074 :MatrixBase<Treal, Tmatrix>() {}
00075 explicit MatrixSymmetric(const MatrixSymmetric<Treal, Tmatrix>& symm)
00076 :MatrixBase<Treal, Tmatrix>(symm) {}
00077 explicit MatrixSymmetric(const XY<Treal, MatrixSymmetric<Treal, Tmatrix> >& sm)
00078 :MatrixBase<Treal, Tmatrix>() { *this = sm.A * sm.B; }
00079 explicit MatrixSymmetric(const MatrixGeneral<Treal, Tmatrix>& matr)
00080 :MatrixBase<Treal, Tmatrix>(matr) {
00081 this->matrixPtr->nosymToSym();
00082 }
00085 #if 0
00086 template<typename Tfull>
00087 inline void assign_from_full
00088 (Tfull const* const fullmatrix,
00089 int const totnrows, int const totncols) {
00090 assert(totnrows == totncols);
00091 this->matrixPtr->assign_from_full(fullmatrix, totnrows, totncols);
00092 this->matrixPtr->nosym_to_sym();
00093 }
00094 inline void assign_from_full
00095 (Treal const* const fullmatrix,
00096 int const totnrows, int const totncols) {
00097 assert(totnrows == totncols);
00098 this->matrixPtr->assign_from_full(fullmatrix, totnrows, totncols);
00099 this->matrixPtr->nosym_to_sym();
00100 }
00101 #endif
00102
00103 inline void assignFromFull
00104 (std::vector<Treal> const & fullMat) {
00105 assert((int)fullMat.size() == this->get_nrows() * this->get_ncols());
00106 this->matrixPtr->assignFromFull(fullMat);
00107 this->matrixPtr->nosymToSym();
00108 }
00109
00110 inline void assignFromFull
00111 (std::vector<Treal> const & fullMat,
00112 std::vector<int> const & rowPermutation,
00113 std::vector<int> const & colPermutation) {
00114 assert((int)fullMat.size() == this->get_nrows() * this->get_ncols());
00115 std::vector<int> rowind(fullMat.size());
00116 std::vector<int> colind(fullMat.size());
00117 int ind = 0;
00118 for (int col = 0; col < this->get_ncols(); ++col)
00119 for (int row = 0; row < this->get_nrows(); ++row) {
00120 rowind[ind] = row;
00121 colind[ind] = col;
00122 ++ind;
00123 }
00124 this->assign_from_sparse(rowind,
00125 colind,
00126 fullMat,
00127 rowPermutation,
00128 colPermutation);
00129 this->matrixPtr->nosymToSym();
00130 }
00131
00132 inline void fullMatrix(std::vector<Treal> & fullMat) const {
00133 this->matrixPtr->syFullMatrix(fullMat);
00134 }
00141 inline void fullMatrix
00142 (std::vector<Treal> & fullMat,
00143 std::vector<int> const & rowInversePermutation,
00144 std::vector<int> const & colInversePermutation) const {
00145 std::vector<int> rowind;
00146 std::vector<int> colind;
00147 std::vector<Treal> values;
00148 get_all_values(rowind, colind, values,
00149 rowInversePermutation,
00150 colInversePermutation);
00151 fullMat.assign(this->get_nrows()*this->get_ncols(),0);
00152 assert(rowind.size() == values.size());
00153 assert(rowind.size() == colind.size());
00154 for (unsigned int ind = 0; ind < values.size(); ++ind) {
00155 assert(rowind[ind] + this->get_nrows() * colind[ind] <
00156 this->get_nrows()*this->get_ncols());
00157 fullMat[rowind[ind] + this->get_nrows() * colind[ind]] =
00158 values[ind];
00159 fullMat[colind[ind] + this->get_nrows() * rowind[ind]] =
00160 values[ind];
00161 }
00162 }
00169 inline void assign_from_sparse
00170 (std::vector<int> const & rowind,
00171 std::vector<int> const & colind,
00172 std::vector<Treal> const & values) {
00173 this->matrixPtr->syAssignFromSparse(rowind, colind, values);
00174 }
00186 inline void assign_from_sparse
00187 (std::vector<int> const & rowind,
00188 std::vector<int> const & colind,
00189 std::vector<Treal> const & values,
00190 SizesAndBlocks const & newRows,
00191 SizesAndBlocks const & newCols) {
00192 this->resetSizesAndBlocks(newRows, newCols);
00193 this->matrixPtr->syAssignFromSparse(rowind, colind, values);
00194 }
00204 inline void assign_from_sparse
00205 (std::vector<int> const & rowind,
00206 std::vector<int> const & colind,
00207 std::vector<Treal> const & values,
00208 std::vector<int> const & rowPermutation,
00209 std::vector<int> const & colPermutation) {
00210 std::vector<int> newRowind;
00211 std::vector<int> newColind;
00212 this->getPermutedAndSymmetrized(rowind, rowPermutation, newRowind,
00213 colind, colPermutation, newColind);
00214
00215 this->matrixPtr->syAssignFromSparse(newRowind, newColind, values);
00216 }
00222 inline void assign_from_sparse
00223 (std::vector<int> const & rowind,
00224 std::vector<int> const & colind,
00225 std::vector<Treal> const & values,
00226 SizesAndBlocks const & newRows,
00227 SizesAndBlocks const & newCols,
00228 std::vector<int> const & rowPermutation,
00229 std::vector<int> const & colPermutation) {
00230 this->resetSizesAndBlocks(newRows, newCols);
00231 assign_from_sparse(rowind, colind, values,
00232 rowPermutation, colPermutation);
00233 }
00241 inline void add_values
00242 (std::vector<int> const & rowind,
00243 std::vector<int> const & colind,
00244 std::vector<Treal> const & values) {
00245 this->matrixPtr->syAddValues(rowind, colind, values);
00246 }
00247
00251 inline void add_values
00252 (std::vector<int> const & rowind,
00253 std::vector<int> const & colind,
00254 std::vector<Treal> const & values,
00255 std::vector<int> const & rowPermutation,
00256 std::vector<int> const & colPermutation) {
00257 std::vector<int> newRowind;
00258 std::vector<int> newColind;
00259 this->getPermutedAndSymmetrized(rowind, rowPermutation, newRowind,
00260 colind, colPermutation, newColind);
00261 this->matrixPtr->syAddValues(newRowind, newColind, values);
00262 }
00263
00264
00265
00266 inline void get_values
00267 (std::vector<int> const & rowind,
00268 std::vector<int> const & colind,
00269 std::vector<Treal> & values) const {
00270 this->matrixPtr->syGetValues(rowind, colind, values);
00271 }
00279 inline void get_values
00280 (std::vector<int> const & rowind,
00281 std::vector<int> const & colind,
00282 std::vector<Treal> & values,
00283 std::vector<int> const & rowPermutation,
00284 std::vector<int> const & colPermutation) const {
00285 std::vector<int> newRowind;
00286 std::vector<int> newColind;
00287 this->getPermutedAndSymmetrized(rowind, rowPermutation, newRowind,
00288 colind, colPermutation, newColind);
00289 this->matrixPtr->syGetValues(newRowind, newColind, values);
00290 }
00296 inline void get_all_values
00297 (std::vector<int> & rowind,
00298 std::vector<int> & colind,
00299 std::vector<Treal> & values) const {
00300 rowind.resize(0);
00301 colind.resize(0);
00302 values.resize(0);
00303 rowind.reserve(nnz());
00304 colind.reserve(nnz());
00305 values.reserve(nnz());
00306 this->matrixPtr->syGetAllValues(rowind, colind, values);
00307 }
00313 inline void get_all_values
00314 (std::vector<int> & rowind,
00315 std::vector<int> & colind,
00316 std::vector<Treal> & values,
00317 std::vector<int> const & rowInversePermutation,
00318 std::vector<int> const & colInversePermutation) const {
00319 std::vector<int> tmpRowind;
00320 std::vector<int> tmpColind;
00321 tmpRowind.reserve(rowind.capacity());
00322 tmpColind.reserve(colind.capacity());
00323 values.resize(0);
00324 this->matrixPtr->syGetAllValues(tmpRowind, tmpColind, values);
00325 this->getPermutedAndSymmetrized(tmpRowind, rowInversePermutation, rowind,
00326 tmpColind, colInversePermutation, colind);
00327 }
00338 MatrixSymmetric<Treal, Tmatrix>&
00339 operator=(const MatrixSymmetric<Treal, Tmatrix>& symm) {
00340 MatrixBase<Treal, Tmatrix>::operator=(symm);
00341 return *this;
00342 }
00343 MatrixSymmetric<Treal, Tmatrix>&
00344 operator=(const MatrixGeneral<Treal, Tmatrix>& matr) {
00345 MatrixBase<Treal, Tmatrix>::operator=(matr);
00346 this->matrixPtr->nosymToSym();
00347 return *this;
00348 }
00349 inline MatrixSymmetric<Treal, Tmatrix>& operator=(int const k) {
00350 *this->matrixPtr = k;
00351 return *this;
00352 }
00353
00354 inline Treal frob() const {
00355 return this->matrixPtr->syFrob();
00356 }
00357 Treal mixed_norm(Treal const requestedAccuracy,
00358 int maxIter = -1) const;
00359 Treal eucl(Treal const requestedAccuracy,
00360 int maxIter = -1) const;
00361
00362 void quickEuclBounds(Treal & euclLowerBound,
00363 Treal & euclUpperBound) const {
00364 Treal frobTmp = frob();
00365 euclLowerBound = frobTmp / template_blas_sqrt( (Treal)this->get_nrows() );
00366 euclUpperBound = frobTmp;
00367 }
00368
00369
00375 static Interval<Treal>
00376 diff(const MatrixSymmetric<Treal, Tmatrix>& A,
00377 const MatrixSymmetric<Treal, Tmatrix>& B,
00378 normType const norm,
00379 Treal const requestedAccuracy);
00387 static Interval<Treal>
00388 diffIfSmall(const MatrixSymmetric<Treal, Tmatrix>& A,
00389 const MatrixSymmetric<Treal, Tmatrix>& B,
00390 normType const norm,
00391 Treal const requestedAccuracy,
00392 Treal const maxAbsVal);
00396 static inline Treal frob_diff
00397 (const MatrixSymmetric<Treal, Tmatrix>& A,
00398 const MatrixSymmetric<Treal, Tmatrix>& B) {
00399 return Tmatrix::syFrobDiff(*A.matrixPtr, *B.matrixPtr);
00400 }
00401
00405 static Treal eucl_diff
00406 (const MatrixSymmetric<Treal, Tmatrix>& A,
00407 const MatrixSymmetric<Treal, Tmatrix>& B,
00408 Treal const requestedAccuracy,
00409 int maxIter = -1);
00410
00414 static Treal mixed_diff
00415 (const MatrixSymmetric<Treal, Tmatrix>& A,
00416 const MatrixSymmetric<Treal, Tmatrix>& B,
00417 Treal const requestedAccuracy);
00418
00425 static Interval<Treal> euclDiffIfSmall
00426 (const MatrixSymmetric<Treal, Tmatrix>& A,
00427 const MatrixSymmetric<Treal, Tmatrix>& B,
00428 Treal const requestedAccuracy,
00429 Treal const maxAbsVal,
00430 VectorType * const eVecPtr = 0);
00431
00432
00443 Treal thresh(Treal const threshold,
00444 normType const norm);
00445
00452 inline Treal frob_thresh(Treal const threshold) {
00453 return 2.0 * this->matrixPtr->frob_thresh(threshold / 2);
00454 }
00455
00456 Treal eucl_thresh(Treal const threshold,
00457 MatrixTriangular<Treal, Tmatrix> const * const Zptr = NULL);
00458
00459 Treal eucl_element_level_thresh(Treal const threshold);
00460
00461 void getSizesAndBlocksForFrobNormMat
00462 ( SizesAndBlocks & rows_new, SizesAndBlocks & cols_new ) const;
00463
00464 Treal mixed_norm_thresh(Treal const threshold);
00465
00466 void simple_blockwise_frob_thresh(Treal const threshold) {
00467 this->matrixPtr->frobThreshLowestLevel(threshold*threshold, 0);
00468 }
00469
00470 inline void gershgorin(Treal& lmin, Treal& lmax) const {
00471 this->matrixPtr->sy_gershgorin(lmin, lmax);
00472 }
00473 static inline Treal trace_ab
00474 (const MatrixSymmetric<Treal, Tmatrix>& A,
00475 const MatrixSymmetric<Treal, Tmatrix>& B) {
00476 return Tmatrix::sy_trace_ab(*A.matrixPtr, *B.matrixPtr);
00477 }
00478 inline size_t nnz() const {
00479 return this->matrixPtr->sy_nnz();
00480 }
00481 inline size_t nvalues() const {
00482 return this->matrixPtr->sy_nvalues();
00483 }
00484 inline void write_to_buffer(void* buffer, const int n_bytes) const {
00485 this->write_to_buffer_base(buffer, n_bytes, matrix_symm);
00486 }
00487 inline void read_from_buffer(void* buffer, const int n_bytes) {
00488 this->read_from_buffer_base(buffer, n_bytes, matrix_symm);
00489 }
00490
00491
00493 MatrixSymmetric<Treal, Tmatrix>& operator=
00494 (XY<Treal, MatrixSymmetric<Treal, Tmatrix> > const & sm);
00496 MatrixSymmetric<Treal, Tmatrix>& operator=
00497 (const XYZpUV<Treal,
00498 MatrixSymmetric<Treal, Tmatrix>,
00499 MatrixSymmetric<Treal, Tmatrix>,
00500 Treal,
00501 MatrixSymmetric<Treal, Tmatrix> >& sm2psm);
00503 MatrixSymmetric<Treal, Tmatrix>& operator=
00504 (const XYZ<Treal,
00505 MatrixSymmetric<Treal, Tmatrix>,
00506 MatrixSymmetric<Treal, Tmatrix> >& sm2);
00508 MatrixSymmetric<Treal, Tmatrix>& operator+=
00509 (const XYZ<Treal,
00510 MatrixSymmetric<Treal, Tmatrix>,
00511 MatrixSymmetric<Treal, Tmatrix> >& sm2);
00513 MatrixSymmetric<Treal, Tmatrix>& operator=
00514 (const XYZpUV<Treal,
00515 MatrixGeneral<Treal, Tmatrix>,
00516 MatrixGeneral<Treal, Tmatrix>,
00517 Treal,
00518 MatrixSymmetric<Treal, Tmatrix> >& smmpsm);
00520 MatrixSymmetric<Treal, Tmatrix>& operator=
00521 (const XYZ<Treal,
00522 MatrixGeneral<Treal, Tmatrix>,
00523 MatrixGeneral<Treal, Tmatrix> >& smm);
00525 MatrixSymmetric<Treal, Tmatrix>& operator+=
00526 (const XYZ<Treal,
00527 MatrixGeneral<Treal, Tmatrix>,
00528 MatrixGeneral<Treal, Tmatrix> >& smm);
00529
00533 MatrixSymmetric<Treal, Tmatrix>& operator=
00534 (const XYZ<MatrixTriangular<Treal, Tmatrix>,
00535 MatrixSymmetric<Treal, Tmatrix>,
00536 MatrixTriangular<Treal, Tmatrix> >& zaz);
00537
00542 static void ssmmUpperTriangleOnly(const Treal alpha,
00543 const MatrixSymmetric<Treal, Tmatrix>& A,
00544 const MatrixSymmetric<Treal, Tmatrix>& B,
00545 const Treal beta,
00546 MatrixSymmetric<Treal, Tmatrix>& C);
00547
00548
00549
00551 MatrixSymmetric<Treal, Tmatrix>& operator=
00552 (XpY<MatrixSymmetric<Treal, Tmatrix>,
00553 MatrixSymmetric<Treal, Tmatrix> > const & mpm);
00555 MatrixSymmetric<Treal, Tmatrix>& operator=
00556 (XmY<MatrixSymmetric<Treal, Tmatrix>,
00557 MatrixSymmetric<Treal, Tmatrix> > const & mm);
00559 MatrixSymmetric<Treal, Tmatrix>& operator+=
00560 (MatrixSymmetric<Treal, Tmatrix> const & A);
00562 MatrixSymmetric<Treal, Tmatrix>& operator-=
00563 (MatrixSymmetric<Treal, Tmatrix> const & A);
00564
00566 MatrixSymmetric<Treal, Tmatrix>& operator+=
00567 (XY<Treal, MatrixSymmetric<Treal, Tmatrix> > const & sm);
00568
00570 MatrixSymmetric<Treal, Tmatrix>& operator-=
00571 (XY<Treal, MatrixSymmetric<Treal, Tmatrix> > const & sm);
00572
00573 template<typename Top>
00574 Treal accumulateWith(Top & op) {
00575 return this->matrixPtr->syAccumulateWith(op);
00576 }
00577
00578 void random() {
00579 this->matrixPtr->syRandom();
00580 }
00581
00582 void randomZeroStructure(Treal probabilityBeingZero) {
00583 this->matrixPtr->syRandomZeroStructure(probabilityBeingZero);
00584 }
00585
00590 template<typename TRule>
00591 void setElementsByRule(TRule & rule) {
00592 this->matrixPtr->sySetElementsByRule(rule);
00593 return;
00594 }
00595
00598 void transfer( MatrixSymmetric<Treal, Tmatrix> & dest ) {
00599 ValidPtr<Tmatrix>::swap( this->matrixPtr, dest.matrixPtr );
00600
00601 this->clear();
00602 }
00603
00604 template<typename Tvector>
00605 void matVecProd(Tvector & y, Tvector const & x) const {
00606 Treal const ONE = 1.0;
00607 y = (ONE * (*this) * x);
00608 }
00609
00610
00611 std::string obj_type_id() const {return "MatrixSymmetric";}
00612 protected:
00613 inline void writeToFileProt(std::ofstream & file) const {
00614 this->writeToFileBase(file, matrix_symm);
00615 }
00616 inline void readFromFileProt(std::ifstream & file) {
00617 this->readFromFileBase(file, matrix_symm);
00618 }
00619
00625 static void getPermutedAndSymmetrized
00626 (std::vector<int> const & rowind,
00627 std::vector<int> const & rowPermutation,
00628 std::vector<int> & newRowind,
00629 std::vector<int> const & colind,
00630 std::vector<int> const & colPermutation,
00631 std::vector<int> & newColind) {
00632 MatrixBase<Treal, Tmatrix>::
00633 getPermutedIndexes(rowind, rowPermutation, newRowind);
00634 MatrixBase<Treal, Tmatrix>::
00635 getPermutedIndexes(colind, colPermutation, newColind);
00636 int tmp;
00637 for (unsigned int i = 0; i < newRowind.size(); ++i) {
00638 if (newRowind[i] > newColind[i]) {
00639 tmp = newRowind[i];
00640 newRowind[i] = newColind[i];
00641 newColind[i] = tmp;
00642 }
00643 }
00644 }
00645 private:
00646 };
00647
00648 template<typename Treal, typename Tmatrix>
00649 Treal MatrixSymmetric<Treal, Tmatrix>::
00650 mixed_norm(Treal const requestedAccuracy,
00651 int maxIter) const {
00652
00653 SizesAndBlocks rows_new;
00654 SizesAndBlocks cols_new;
00655 this->getSizesAndBlocksForFrobNormMat( rows_new, cols_new );
00656
00657
00658 MatrixSymmetric<Treal, typename Tmatrix::ElementType> frobNormMat;
00659 frobNormMat.resetSizesAndBlocks(rows_new, cols_new);
00660 frobNormMat.getMatrix().syAssignFrobNormsLowestLevel(this->getMatrix());
00661 return frobNormMat.eucl(requestedAccuracy, maxIter);
00662 }
00663
00664
00665 template<typename Treal, typename Tmatrix>
00666 Treal MatrixSymmetric<Treal, Tmatrix>::
00667 eucl(Treal const requestedAccuracy,
00668 int maxIter) const {
00669 assert(requestedAccuracy >= 0);
00670
00671 Treal frobTmp = this->frob();
00672 if (frobTmp < requestedAccuracy)
00673 return (Treal)0.0;
00674 if (maxIter < 0)
00675 maxIter = this->get_nrows() * 100;
00676 VectorType guess;
00677 SizesAndBlocks cols;
00678 this->getCols(cols);
00679 guess.resetSizesAndBlocks(cols);
00680 guess.rand();
00681
00682 #if 0 // "new code"
00683 MatrixSymmetric<Treal, Tmatrix> Copy(*this);
00684 Copy.frob_thresh(requestedAccuracy / 2.0);
00685 arn::LanczosLargestMagnitudeEig
00686 <Treal, MatrixSymmetric<Treal, Tmatrix>, VectorType>
00687 lan(Copy, guess, maxIter);
00688 lan.setAbsTol( requestedAccuracy / 2.0 );
00689 #else // "old code"
00690 arn::LanczosLargestMagnitudeEig
00691 <Treal, MatrixSymmetric<Treal, Tmatrix>, VectorType>
00692 lan(*this, guess, maxIter);
00693 lan.setAbsTol( requestedAccuracy );
00694 #endif
00695 lan.run();
00696 Treal eVal = 0;
00697 Treal acc = 0;
00698 lan.getLargestMagnitudeEig(eVal, acc);
00699 return template_blas_fabs(eVal);
00700 }
00701
00702 template<typename Treal, typename Tmatrix>
00703 Interval<Treal> MatrixSymmetric<Treal, Tmatrix>::
00704 diff(const MatrixSymmetric<Treal, Tmatrix>& A,
00705 const MatrixSymmetric<Treal, Tmatrix>& B,
00706 normType const norm, Treal const requestedAccuracy) {
00707 Treal diff;
00708 Treal eNMin;
00709 switch (norm) {
00710 case frobNorm:
00711 diff = frob_diff(A, B);
00712 return Interval<Treal>(diff / template_blas_sqrt((Treal)A.get_nrows()), diff);
00713 break;
00714 case euclNorm:
00715 diff = eucl_diff(A, B, requestedAccuracy);
00716 eNMin = diff - requestedAccuracy;
00717 eNMin = eNMin >= 0 ? eNMin : 0;
00718 return Interval<Treal>(eNMin, diff + requestedAccuracy);
00719 break;
00720 default:
00721 throw Failure("MatrixSymmetric<Treal, Tmatrix>::"
00722 "diff(const MatrixSymmetric<Treal, Tmatrix>&, "
00723 "const MatrixSymmetric<Treal, Tmatrix>&, "
00724 "normType const, Treal): "
00725 "Diff not implemented for selected norm");
00726 }
00727 }
00728
00729 #if 1
00730 template<typename Treal, typename Tmatrix>
00731 Interval<Treal> MatrixSymmetric<Treal, Tmatrix>::
00732 diffIfSmall(const MatrixSymmetric<Treal, Tmatrix>& A,
00733 const MatrixSymmetric<Treal, Tmatrix>& B,
00734 normType const norm,
00735 Treal const requestedAccuracy,
00736 Treal const maxAbsVal) {
00737 Treal diff;
00738 switch (norm) {
00739 case frobNorm:
00740 {
00741 diff = frob_diff(A, B);
00742 return Interval<Treal>(diff / template_blas_sqrt((Treal)A.get_nrows()), diff);
00743 }
00744 break;
00745 case euclNorm:
00746 return euclDiffIfSmall(A, B, requestedAccuracy, maxAbsVal);
00747 break;
00748 default:
00749 throw Failure("MatrixSymmetric<Treal, Tmatrix>::"
00750 "diffIfSmall"
00751 "(const MatrixSymmetric<Treal, Tmatrix>&, "
00752 "const MatrixSymmetric<Treal, Tmatrix>&, "
00753 "normType const, Treal const, Treal const): "
00754 "Diff not implemented for selected norm");
00755 }
00756 }
00757
00758 #endif
00759
00760
00761 template<typename Treal, typename Tmatrix>
00762 Treal MatrixSymmetric<Treal, Tmatrix>::eucl_diff
00763 (const MatrixSymmetric<Treal, Tmatrix>& A,
00764 const MatrixSymmetric<Treal, Tmatrix>& B,
00765 Treal const requestedAccuracy,
00766 int maxIter) {
00767
00768 mat::DiffMatrix< MatrixSymmetric<Treal, Tmatrix>, Treal> Diff(A, B);
00769 Treal maxAbsVal = 2 * frob_diff(A,B);
00770
00771
00772 Treal relTol = template_blas_sqrt(template_blas_sqrt(mat::getMachineEpsilon<Treal>()));
00773 VectorType * const eVecPtrNotUsed = 0;
00774 Interval<Treal> euclInt =
00775 mat::euclIfSmall(Diff, requestedAccuracy, relTol, maxAbsVal, eVecPtrNotUsed, maxIter);
00776 return euclInt.midPoint();
00777 }
00778
00779 template<typename Treal, typename Tmatrix>
00780 Treal MatrixSymmetric<Treal, Tmatrix>::mixed_diff
00781 (const MatrixSymmetric<Treal, Tmatrix>& A,
00782 const MatrixSymmetric<Treal, Tmatrix>& B,
00783 Treal const requestedAccuracy) {
00784 MatrixSymmetric<Treal, typename Tmatrix::ElementType> frobNormMat;
00785 {
00786 SizesAndBlocks rows_new;
00787 SizesAndBlocks cols_new;
00788 A.getSizesAndBlocksForFrobNormMat( rows_new, cols_new );
00789 frobNormMat.resetSizesAndBlocks(rows_new, cols_new);
00790 frobNormMat.getMatrix().syAssignDiffFrobNormsLowestLevel(A.getMatrix(),B.getMatrix());
00791 }
00792 return frobNormMat.eucl(requestedAccuracy);
00793 }
00794
00795 #if 1
00796 template<typename Treal, typename Tmatrix>
00797 Interval<Treal> MatrixSymmetric<Treal, Tmatrix>::euclDiffIfSmall
00798 (const MatrixSymmetric<Treal, Tmatrix>& A,
00799 const MatrixSymmetric<Treal, Tmatrix>& B,
00800 Treal const requestedAccuracy,
00801 Treal const maxAbsVal,
00802 VectorType * const eVecPtr) {
00803
00804 mat::DiffMatrix< MatrixSymmetric<Treal, Tmatrix>, Treal> Diff(A, B);
00805
00806 Treal relTol = template_blas_sqrt(template_blas_sqrt(mat::getMachineEpsilon<Treal>()));
00807 Interval<Treal> tmpInterval = mat::euclIfSmall(Diff, requestedAccuracy, relTol, maxAbsVal, eVecPtr);
00808
00809
00810
00811
00812
00813
00814
00815
00816
00817 if ( tmpInterval.length() < 2*requestedAccuracy )
00818 return Interval<Treal>( tmpInterval.midPoint()-requestedAccuracy, tmpInterval.midPoint()+requestedAccuracy );
00819 return tmpInterval;
00820 }
00821
00822 #endif
00823
00824
00825 template<typename Treal, typename Tmatrix>
00826 Treal MatrixSymmetric<Treal, Tmatrix>::
00827 thresh(Treal const threshold,
00828 normType const norm) {
00829 switch (norm) {
00830 case frobNorm:
00831 return this->frob_thresh(threshold);
00832 break;
00833 case euclNorm:
00834 return this->eucl_thresh(threshold);
00835 break;
00836 case mixedNorm:
00837 return this->mixed_norm_thresh(threshold);
00838 break;
00839 default:
00840 throw Failure("MatrixSymmetric<Treal, Tmatrix>::"
00841 "thresh(Treal const, "
00842 "normType const): "
00843 "Thresholding not implemented for selected norm");
00844 }
00845 }
00846
00847 #if 1
00848
00849 template<typename Treal, typename Tmatrix>
00850 Treal MatrixSymmetric<Treal, Tmatrix>::
00851 eucl_thresh(Treal const threshold,
00852 MatrixTriangular<Treal, Tmatrix> const * const Zptr) {
00853 if ( Zptr == NULL ) {
00854 EuclTruncationSymm<MatrixSymmetric<Treal, Tmatrix>, Treal> TruncObj(*this);
00855 return TruncObj.run( threshold );
00856 }
00857 EuclTruncationSymmWithZ<MatrixSymmetric<Treal, Tmatrix>, MatrixTriangular<Treal, Tmatrix>, Treal> TruncObj(*this, *Zptr);
00858 return TruncObj.run( threshold );
00859 }
00860
00861 #endif
00862
00863
00864 template<typename Treal, typename Tmatrix>
00865 void MatrixSymmetric<Treal, Tmatrix>::getSizesAndBlocksForFrobNormMat
00866 ( SizesAndBlocks & rows_new, SizesAndBlocks & cols_new ) const {
00867 std::vector<int> rows_block_sizes;
00868 std::vector<int> cols_block_sizes;
00869 int n_rows;
00870 int n_cols;
00871 {
00872 SizesAndBlocks rows;
00873 SizesAndBlocks cols;
00874 this->getRows(rows);
00875 this->getCols(cols);
00876 rows.getBlockSizeVector( rows_block_sizes );
00877 cols.getBlockSizeVector( cols_block_sizes );
00878 rows_block_sizes.pop_back();
00879 cols_block_sizes.pop_back();
00880 n_rows = rows.getNTotalScalars();
00881 n_cols = cols.getNTotalScalars();
00882 int factor_rows = rows_block_sizes[rows_block_sizes.size()-1];
00883 int factor_cols = cols_block_sizes[cols_block_sizes.size()-1];
00884 for (unsigned int ind = 0; ind < rows_block_sizes.size(); ++ind)
00885 rows_block_sizes[ind] = rows_block_sizes[ind] / factor_rows;
00886 for (unsigned int ind = 0; ind < cols_block_sizes.size(); ++ind)
00887 cols_block_sizes[ind] = cols_block_sizes[ind] / factor_cols;
00888
00889
00890
00891 if (n_rows % factor_rows)
00892 n_rows = n_rows / factor_rows + 1;
00893 else
00894 n_rows = n_rows / factor_rows;
00895 if (n_cols % factor_cols)
00896 n_cols = n_cols / factor_cols + 1;
00897 else
00898 n_cols = n_cols / factor_cols;
00899 }
00900 rows_new = SizesAndBlocks( rows_block_sizes, n_rows );
00901 cols_new = SizesAndBlocks( cols_block_sizes, n_cols );
00902 }
00903
00904
00905 template<typename Treal, typename Tmatrix>
00906 Treal MatrixSymmetric<Treal, Tmatrix>::
00907 mixed_norm_thresh(Treal const threshold) {
00908 assert(threshold >= (Treal)0.0);
00909 if (threshold == (Treal)0.0)
00910 return (Treal)0;
00911
00912 SizesAndBlocks rows_new;
00913 SizesAndBlocks cols_new;
00914 this->getSizesAndBlocksForFrobNormMat( rows_new, cols_new );
00915
00916
00917 MatrixSymmetric<Treal, typename Tmatrix::ElementType> frobNormMat;
00918 frobNormMat.resetSizesAndBlocks(rows_new, cols_new);
00919
00920
00921
00922 frobNormMat.getMatrix().syAssignFrobNormsLowestLevel(this->getMatrix());
00923 EuclTruncationSymmElementLevel<MatrixSymmetric<Treal, typename Tmatrix::ElementType>, Treal> TruncObj( frobNormMat );
00924 Treal mixed_norm_result = TruncObj.run( threshold );
00925 frobNormMat.getMatrix().truncateAccordingToSparsityPattern(this->getMatrix());
00926 return mixed_norm_result;
00927 }
00928
00929
00930
00931 template<typename Treal, typename Tmatrix>
00932 inline MatrixSymmetric<Treal, Tmatrix>&
00933 MatrixSymmetric<Treal, Tmatrix>::operator=
00934 (XY<Treal, MatrixSymmetric<Treal, Tmatrix> > const & sm) {
00935
00936 if(this == &sm.B)
00937 {
00938 *this *= sm.A;
00939 return *this;
00940 }
00941 assert(!sm.tB);
00942
00943 this->matrixPtr.haveDataStructureSet(true);
00944 this->matrixPtr->assign(sm.A, *sm.B.matrixPtr);
00945 return *this;
00946 }
00947
00948 template<typename Treal, typename Tmatrix>
00949 MatrixSymmetric<Treal, Tmatrix>&
00950 MatrixSymmetric<Treal, Tmatrix>::operator=
00951 (const XYZpUV<Treal,
00952 MatrixSymmetric<Treal, Tmatrix>,
00953 MatrixSymmetric<Treal, Tmatrix>,
00954 Treal,
00955 MatrixSymmetric<Treal, Tmatrix> >& sm2psm) {
00956 assert(this != &sm2psm.B);
00957 if (this == &sm2psm.E && &sm2psm.B == &sm2psm.C) {
00958
00959 Tmatrix::sysq('U',
00960 sm2psm.A, *sm2psm.B.matrixPtr,
00961 sm2psm.D, *this->matrixPtr);
00962 return *this;
00963 }
00964 else
00965 throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
00966 "(const XYZpUV<Treal, MatrixSymmetric"
00967 "<Treal, Tmatrix> >& sm2psm) : "
00968 "D = alpha * A * B + beta * C not supported for C != D"
00969 " and for symmetric matrices not for A != B since this "
00970 "generally will result in a nonsymmetric matrix");
00971 }
00972
00973
00974 template<typename Treal, typename Tmatrix>
00975 MatrixSymmetric<Treal, Tmatrix>&
00976 MatrixSymmetric<Treal, Tmatrix>::operator=
00977 (const XYZ<Treal,
00978 MatrixSymmetric<Treal, Tmatrix>,
00979 MatrixSymmetric<Treal, Tmatrix> >& sm2) {
00980 assert(this != &sm2.B);
00981 if (&sm2.B == &sm2.C) {
00982 this->matrixPtr.haveDataStructureSet(true);
00983 Tmatrix::sysq('U', sm2.A, *sm2.B.matrixPtr, 0, *this->matrixPtr);
00984 return *this;
00985 }
00986 else
00987 throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
00988 "(const XYZ<Treal, MatrixSymmetric<Treal, Tmatrix>,"
00989 " MatrixSymmetric<Treal, Tmatrix> >& sm2) : "
00990 "Operation C = alpha * A * B with only symmetric "
00991 "matrices not supported for A != B");
00992 }
00993
00994
00995 template<typename Treal, typename Tmatrix>
00996 MatrixSymmetric<Treal, Tmatrix>&
00997 MatrixSymmetric<Treal, Tmatrix>::operator+=
00998 (const XYZ<Treal,
00999 MatrixSymmetric<Treal, Tmatrix>,
01000 MatrixSymmetric<Treal, Tmatrix> >& sm2) {
01001 assert(this != &sm2.B);
01002 if (&sm2.B == &sm2.C) {
01003 Tmatrix::sysq('U', sm2.A, *sm2.B.matrixPtr, 1, *this->matrixPtr);
01004 return *this;
01005 }
01006 else
01007 throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator+="
01008 "(const XYZ<Treal, MatrixSymmetric<Treal, Tmatrix>,"
01009 " MatrixSymmetric<Treal, Tmatrix> >& sm2) : "
01010 "Operation C += alpha * A * B with only symmetric "
01011 "matrices not supported for A != B");
01012 }
01013
01014
01015 template<typename Treal, typename Tmatrix>
01016 MatrixSymmetric<Treal, Tmatrix>&
01017 MatrixSymmetric<Treal, Tmatrix>::operator=
01018 (const XYZpUV<Treal,
01019 MatrixGeneral<Treal, Tmatrix>,
01020 MatrixGeneral<Treal, Tmatrix>,
01021 Treal,
01022 MatrixSymmetric<Treal, Tmatrix> >& smmpsm) {
01023 if (this == &smmpsm.E)
01024 if (&smmpsm.B == &smmpsm.C)
01025 if (!smmpsm.tB && smmpsm.tC) {
01026 Tmatrix::syrk('U', false,
01027 smmpsm.A, *smmpsm.B.matrixPtr,
01028 smmpsm.D, *this->matrixPtr);
01029 }
01030 else
01031 throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
01032 "(const XYZpUV<Treal, MatrixGeneral"
01033 "<Treal, Tmatrix>, "
01034 "MatrixGeneral<Treal, Tmatrix>, Treal, "
01035 "MatrixSymmetric<Treal, Tmatrix> >&) : "
01036 "C = alpha * A' * A + beta * C, not implemented"
01037 " only C = alpha * A * A' + beta * C");
01038 else
01039 throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
01040 "(const XYZpUV<"
01041 "Treal, MatrixGeneral<Treal, Tmatrix>, "
01042 "MatrixGeneral<Treal, Tmatrix>, Treal, "
01043 "MatrixSymmetric<Treal, Tmatrix> >&) : "
01044 "You are trying to call C = alpha * A * A' + beta * C"
01045 " with A and A' being different objects");
01046 else
01047 throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
01048 "(const XYZpUV<"
01049 "Treal, MatrixGeneral<Treal, Tmatrix>, "
01050 "MatrixGeneral<Treal, Tmatrix>, Treal, "
01051 "MatrixSymmetric<Treal, Tmatrix> >&) : "
01052 "D = alpha * A * A' + beta * C not supported for C != D");
01053 return *this;
01054 }
01055
01056
01057 template<typename Treal, typename Tmatrix>
01058 MatrixSymmetric<Treal, Tmatrix>&
01059 MatrixSymmetric<Treal, Tmatrix>::operator=
01060 (const XYZ<Treal,
01061 MatrixGeneral<Treal, Tmatrix>,
01062 MatrixGeneral<Treal, Tmatrix> >& smm) {
01063 if (&smm.B == &smm.C)
01064 if (!smm.tB && smm.tC) {
01065 Tmatrix::syrk('U', false,
01066 smm.A, *smm.B.matrixPtr,
01067 0, *this->matrixPtr);
01068 }
01069 else
01070 throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
01071 "(const XYZ<"
01072 "Treal, MatrixGeneral<Treal, Tmatrix>, "
01073 "MatrixGeneral<Treal, Tmatrix> >&) : "
01074 "C = alpha * A' * A, not implemented "
01075 "only C = alpha * A * A'");
01076 else
01077 throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
01078 "(const XYZ<"
01079 "Treal, MatrixGeneral<Treal, Tmatrix>, "
01080 "MatrixGeneral<Treal, Tmatrix> >&) : "
01081 "You are trying to call C = alpha * A * A' "
01082 "with A and A' being different objects");
01083 return *this;
01084 }
01085
01086
01087 template<typename Treal, typename Tmatrix>
01088 MatrixSymmetric<Treal, Tmatrix>&
01089 MatrixSymmetric<Treal, Tmatrix>::operator+=
01090 (const XYZ<Treal,
01091 MatrixGeneral<Treal, Tmatrix>,
01092 MatrixGeneral<Treal, Tmatrix> >& smm) {
01093 if (&smm.B == &smm.C)
01094 if (!smm.tB && smm.tC) {
01095 Tmatrix::syrk('U', false,
01096 smm.A, *smm.B.matrixPtr,
01097 1, *this->matrixPtr);
01098 }
01099 else
01100 throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator+="
01101 "(const XYZ<"
01102 "Treal, MatrixGeneral<Treal, Tmatrix>, "
01103 "MatrixGeneral<Treal, Tmatrix> >&) : "
01104 "C += alpha * A' * A, not implemented "
01105 "only C += alpha * A * A'");
01106 else
01107 throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator+="
01108 "(const XYZ<"
01109 "Treal, MatrixGeneral<Treal, Tmatrix>, "
01110 "MatrixGeneral<Treal, Tmatrix> >&) : "
01111 "You are trying to call C += alpha * A * A' "
01112 "with A and A' being different objects");
01113 return *this;
01114 }
01115
01116 #if 1
01117
01118
01119 template<typename Treal, typename Tmatrix>
01120 MatrixSymmetric<Treal, Tmatrix>&
01121 MatrixSymmetric<Treal, Tmatrix>::operator=
01122 (const XYZ<MatrixTriangular<Treal, Tmatrix>,
01123 MatrixSymmetric<Treal, Tmatrix>,
01124 MatrixTriangular<Treal, Tmatrix> >& zaz) {
01125 if (this == &zaz.B) {
01126 if (&zaz.A == &zaz.C) {
01127 if (zaz.tA && !zaz.tC) {
01128 Tmatrix::trsytriplemm('R', *zaz.A.matrixPtr, *this->matrixPtr);
01129 }
01130 else if (!zaz.tA && zaz.tC) {
01131 Tmatrix::trsytriplemm('L', *zaz.A.matrixPtr, *this->matrixPtr);
01132 }
01133 else
01134 throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
01135 "(const XYZ<MatrixTriangular<Treal, Tmatrix>,"
01136 "MatrixSymmetric<Treal, Tmatrix>,"
01137 "MatrixTriangular<Treal, Tmatrix> >&) : "
01138 "A = op1(Z) * A * op2(Z) : Either op1 xor op2 must be "
01139 "the transpose operation.");
01140 }
01141 else
01142 throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
01143 "(const XYZ<MatrixTriangular<Treal, Tmatrix>,"
01144 "MatrixSymmetric<Treal, Tmatrix>,"
01145 "MatrixTriangular<Treal, Tmatrix> >&) : "
01146 "A = op1(Z1) * A * op2(Z2) : Z1 and Z2 must be the same "
01147 "object");
01148 }
01149 else
01150 throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
01151 "(const XYZ<MatrixTriangular<Treal, Tmatrix>,"
01152 "MatrixSymmetric<Treal, Tmatrix>,"
01153 "MatrixTriangular<Treal, Tmatrix> >&) : "
01154 "C = op1(Z) * A * op2(Z) : A and C must be the same "
01155 "object");
01156 return *this;
01157 }
01158
01159 #endif
01160
01161
01166 template<typename Treal, typename Tmatrix>
01167 void MatrixSymmetric<Treal, Tmatrix>::
01168 ssmmUpperTriangleOnly(const Treal alpha,
01169 const MatrixSymmetric<Treal, Tmatrix>& A,
01170 const MatrixSymmetric<Treal, Tmatrix>& B,
01171 const Treal beta,
01172 MatrixSymmetric<Treal, Tmatrix>& C) {
01173 Tmatrix::ssmm_upper_tr_only(alpha, *A.matrixPtr, *B.matrixPtr,
01174 beta, *C.matrixPtr);
01175 }
01176
01177
01178
01179
01180
01181
01182
01183 template<typename Treal, typename Tmatrix>
01184 inline MatrixSymmetric<Treal, Tmatrix>&
01185 MatrixSymmetric<Treal, Tmatrix>::operator=
01186 (const XpY<MatrixSymmetric<Treal, Tmatrix>,
01187 MatrixSymmetric<Treal, Tmatrix> >& mpm) {
01188 assert(this != &mpm.A);
01189 (*this) = mpm.B;
01190 Tmatrix::add(1.0, *mpm.A.matrixPtr, *this->matrixPtr);
01191 return *this;
01192 }
01193
01194 template<typename Treal, typename Tmatrix>
01195 inline MatrixSymmetric<Treal, Tmatrix>&
01196 MatrixSymmetric<Treal, Tmatrix>::operator=
01197 (const XmY<MatrixSymmetric<Treal, Tmatrix>,
01198 MatrixSymmetric<Treal, Tmatrix> >& mmm) {
01199 assert(this != &mmm.B);
01200 (*this) = mmm.A;
01201 Tmatrix::add(-1.0, *mmm.B.matrixPtr, *this->matrixPtr);
01202 return *this;
01203 }
01204
01205 template<typename Treal, typename Tmatrix>
01206 inline MatrixSymmetric<Treal, Tmatrix>&
01207 MatrixSymmetric<Treal, Tmatrix>::operator+=
01208 (MatrixSymmetric<Treal, Tmatrix> const & A) {
01209 Tmatrix::add(1.0, *A.matrixPtr, *this->matrixPtr);
01210 return *this;
01211 }
01212
01213 template<typename Treal, typename Tmatrix>
01214 inline MatrixSymmetric<Treal, Tmatrix>&
01215 MatrixSymmetric<Treal, Tmatrix>::operator-=
01216 (MatrixSymmetric<Treal, Tmatrix> const & A) {
01217 Tmatrix::add(-1.0, *A.matrixPtr, *this->matrixPtr);
01218 return *this;
01219 }
01220
01221 template<typename Treal, typename Tmatrix>
01222 inline MatrixSymmetric<Treal, Tmatrix>&
01223 MatrixSymmetric<Treal, Tmatrix>::operator+=
01224 (XY<Treal, MatrixSymmetric<Treal, Tmatrix> > const & sm) {
01225 assert(!sm.tB);
01226 Tmatrix::add(sm.A, *sm.B.matrixPtr, *this->matrixPtr);
01227 return *this;
01228 }
01229
01230
01231 template<typename Treal, typename Tmatrix>
01232 inline MatrixSymmetric<Treal, Tmatrix>&
01233 MatrixSymmetric<Treal, Tmatrix>::operator-=
01234 (XY<Treal, MatrixSymmetric<Treal, Tmatrix> > const & sm) {
01235 assert(!sm.tB);
01236 Tmatrix::add(-sm.A, *sm.B.matrixPtr, *this->matrixPtr);
01237 return *this;
01238 }
01239
01244 template<typename Treal, typename Tmatrix, typename Top>
01245 Treal accumulate(MatrixSymmetric<Treal, Tmatrix> & A,
01246 Top & op) {
01247 return A.accumulateWith(op);
01248 }
01249
01250
01251
01252 }
01253 #endif
01254