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_MATRIXGENERAL
00039 #define MAT_MATRIXGENERAL
00040 #include "MatrixBase.h"
00041 namespace mat {
00058 template<typename Treal, typename Tmatrix>
00059 class MatrixGeneral : public MatrixBase<Treal, Tmatrix> {
00060 public:
00061 typedef VectorGeneral<Treal, typename Tmatrix::VectorType> VectorType;
00062
00063 MatrixGeneral()
00064 :MatrixBase<Treal, Tmatrix>() {}
00065 explicit MatrixGeneral(const MatrixGeneral<Treal, Tmatrix>& matr)
00066 :MatrixBase<Treal, Tmatrix>(matr) {}
00067 explicit MatrixGeneral(const MatrixSymmetric<Treal, Tmatrix>& symm)
00068 :MatrixBase<Treal, Tmatrix>(symm) {
00069 this->matrixPtr->symToNosym();
00070 }
00071 explicit MatrixGeneral(const MatrixTriangular<Treal, Tmatrix>& triang)
00072 :MatrixBase<Treal, Tmatrix>(triang) {}
00075 #if 0
00076 template<typename Tfull>
00077 inline void assign_from_full
00078 (Tfull const* const fullmatrix,
00079 const int totnrows, const int totncols) {
00080 this->matrixPtr->assign_from_full(fullmatrix, totnrows, totncols);
00081 }
00082 inline void assign_from_full
00083 (Treal const* const fullmatrix,
00084 const int totnrows, const int totncols) {
00085 this->matrixPtr->assign_from_full(fullmatrix, totnrows, totncols);
00086 }
00087 #endif
00088
00089 inline void assignFromFull
00090 (std::vector<Treal> const & fullMat) {
00091 assert((int)fullMat.size() == this->get_nrows() * this->get_ncols());
00092 this->matrixPtr->assignFromFull(fullMat);
00093 }
00094
00095 inline void fullMatrix(std::vector<Treal> & fullMat) const {
00096 this->matrixPtr->fullMatrix(fullMat);
00097 }
00098
00099 inline void fullMatrix
00100 (std::vector<Treal> & fullMat,
00101 std::vector<int> const & rowInversePermutation,
00102 std::vector<int> const & colInversePermutation) const {
00103 std::vector<int> rowind;
00104 std::vector<int> colind;
00105 std::vector<Treal> values;
00106 get_all_values(rowind, colind, values,
00107 rowInversePermutation,
00108 colInversePermutation);
00109 fullMat.assign(this->get_nrows()*this->get_ncols(),0);
00110 for (unsigned int ind = 0; ind < values.size(); ++ind)
00111 fullMat[rowind[ind] + this->get_nrows() * colind[ind]] =
00112 values[ind];
00113 }
00121 inline void assign_from_sparse
00122 (std::vector<int> const & rowind,
00123 std::vector<int> const & colind,
00124 std::vector<Treal> const & values,
00125 SizesAndBlocks const & newRows,
00126 SizesAndBlocks const & newCols) {
00127 this->resetSizesAndBlocks(newRows, newCols);
00128 this->matrixPtr->assignFromSparse(rowind, colind, values);
00129 }
00137 inline void assign_from_sparse
00138 (std::vector<int> const & rowind,
00139 std::vector<int> const & colind,
00140 std::vector<Treal> const & values,
00141 std::vector<int> const & rowPermutation,
00142 std::vector<int> const & colPermutation) {
00143 std::vector<int> newRowind;
00144 std::vector<int> newColind;
00145 MatrixBase<Treal, Tmatrix>::
00146 getPermutedIndexes(rowind, rowPermutation, newRowind);
00147 MatrixBase<Treal, Tmatrix>::
00148 getPermutedIndexes(colind, colPermutation, newColind);
00149 this->matrixPtr->assignFromSparse(newRowind, newColind, values);
00150 }
00156 inline void assign_from_sparse
00157 (std::vector<int> const & rowind,
00158 std::vector<int> const & colind,
00159 std::vector<Treal> const & values,
00160 SizesAndBlocks const & newRows,
00161 SizesAndBlocks const & newCols,
00162 std::vector<int> const & rowPermutation,
00163 std::vector<int> const & colPermutation) {
00164 this->resetSizesAndBlocks(newRows, newCols);
00165 assign_from_sparse(rowind, colind, values,
00166 rowPermutation, colPermutation);
00167 }
00172 inline void get_values
00173 (std::vector<int> const & rowind,
00174 std::vector<int> const & colind,
00175 std::vector<Treal> & values) const {
00176 this->matrixPtr->getValues(rowind, colind, values);
00177 }
00185 inline void get_values
00186 (std::vector<int> const & rowind,
00187 std::vector<int> const & colind,
00188 std::vector<Treal> & values,
00189 std::vector<int> const & rowPermutation,
00190 std::vector<int> const & colPermutation) const {
00191 std::vector<int> newRowind;
00192 std::vector<int> newColind;
00193 MatrixBase<Treal, Tmatrix>::
00194 getPermutedIndexes(rowind, rowPermutation, newRowind);
00195 MatrixBase<Treal, Tmatrix>::
00196 getPermutedIndexes(colind, colPermutation, newColind);
00197 this->matrixPtr->getValues(newRowind, newColind, values);
00198 }
00203 inline void get_all_values
00204 (std::vector<int> & rowind,
00205 std::vector<int> & colind,
00206 std::vector<Treal> & values) const {
00207 rowind.resize(0);
00208 colind.resize(0);
00209 values.resize(0);
00210 this->matrixPtr->getAllValues(rowind, colind, values);
00211 }
00221 inline void get_all_values
00222 (std::vector<int> & rowind,
00223 std::vector<int> & colind,
00224 std::vector<Treal> & values,
00225 std::vector<int> const & rowInversePermutation,
00226 std::vector<int> const & colInversePermutation) const {
00227 std::vector<int> tmpRowind;
00228 std::vector<int> tmpColind;
00229 tmpRowind.reserve(rowind.capacity());
00230 tmpColind.reserve(colind.capacity());
00231 values.resize(0);
00232 this->matrixPtr->getAllValues(tmpRowind, tmpColind, values);
00233 MatrixBase<Treal, Tmatrix>::
00234 getPermutedIndexes(tmpRowind, rowInversePermutation, rowind);
00235 MatrixBase<Treal, Tmatrix>::
00236 getPermutedIndexes(tmpColind, colInversePermutation, colind);
00237
00238 }
00248 #if 0
00249 inline void fullmatrix(Treal* const full,
00250 const int totnrows, const int totncols) const {
00251 this->matrixPtr->fullmatrix(full, totnrows, totncols);
00252 }
00253 #endif
00254 MatrixGeneral<Treal, Tmatrix>&
00255 operator=(const MatrixGeneral<Treal, Tmatrix>& mat) {
00256 MatrixBase<Treal, Tmatrix>::operator=(mat);
00257 return *this;
00258 }
00259 inline MatrixGeneral<Treal, Tmatrix>&
00260 operator=(const Xtrans<MatrixGeneral<Treal, Tmatrix> >& mt) {
00261 if (mt.tA)
00262 MatrixBase<Treal, Tmatrix>::operator=(transpose(mt.A));
00263 else
00264 MatrixBase<Treal, Tmatrix>::operator=(mt.A);
00265 return *this;
00266 }
00267
00268 MatrixGeneral<Treal, Tmatrix>&
00269 operator=(const MatrixSymmetric<Treal, Tmatrix>& symm) {
00270 MatrixBase<Treal, Tmatrix>::operator=(symm);
00271 this->matrixPtr->symToNosym();
00272 return *this;
00273 }
00274 MatrixGeneral<Treal, Tmatrix>&
00275 operator=(const MatrixTriangular<Treal, Tmatrix>& triang) {
00276 MatrixBase<Treal, Tmatrix>::operator=(triang);
00277 return *this;
00278 }
00279
00280 inline MatrixGeneral<Treal, Tmatrix>& operator=(int const k) {
00281 *this->matrixPtr = k;
00282 return *this;
00283 }
00284 inline Treal frob() const {
00285 return this->matrixPtr->frob();
00286 }
00287 static inline Treal frob_diff
00288 (const MatrixGeneral<Treal, Tmatrix>& A,
00289 const MatrixGeneral<Treal, Tmatrix>& B) {
00290 return Tmatrix::frobDiff(*A.matrixPtr, *B.matrixPtr);
00291 }
00292 Treal eucl(Treal const requestedAccuracy,
00293 int maxIter = -1) const;
00294
00295
00296 void thresh(Treal const threshold,
00297 normType const norm);
00298
00299 inline void frob_thresh(Treal threshold) {
00300 this->matrixPtr->frob_thresh(threshold);
00301 }
00302 Treal eucl_thresh(Treal const threshold);
00303
00304 inline void gershgorin(Treal& lmin, Treal& lmax) {
00305 this->matrixPtr->gershgorin(lmin, lmax);
00306 }
00307 static inline Treal trace_ab
00308 (const MatrixGeneral<Treal, Tmatrix>& A,
00309 const MatrixGeneral<Treal, Tmatrix>& B) {
00310 return Tmatrix::trace_ab(*A.matrixPtr, *B.matrixPtr);
00311 }
00312 static inline Treal trace_aTb
00313 (const MatrixGeneral<Treal, Tmatrix>& A,
00314 const MatrixGeneral<Treal, Tmatrix>& B) {
00315 return Tmatrix::trace_aTb(*A.matrixPtr, *B.matrixPtr);
00316 }
00317 inline size_t nnz() const {
00318 return this->matrixPtr->nnz();
00319 }
00320 inline size_t nvalues() const {
00321 return this->matrixPtr->nvalues();
00322 }
00323
00324 inline void write_to_buffer(void* buffer, const int n_bytes) const {
00325 this->write_to_buffer_base(buffer, n_bytes, matrix_matr);
00326 }
00327 inline void read_from_buffer(void* buffer, const int n_bytes) {
00328 this->read_from_buffer_base(buffer, n_bytes, matrix_matr);
00329 }
00330
00331
00333 MatrixGeneral<Treal, Tmatrix>& operator=
00334 (const XYZ<Treal,
00335 MatrixGeneral<Treal, Tmatrix>,
00336 MatrixGeneral<Treal, Tmatrix> >& smm);
00337
00339 MatrixGeneral<Treal, Tmatrix>& operator=
00340 (const XY<MatrixGeneral<Treal, Tmatrix>,
00341 MatrixGeneral<Treal, Tmatrix> >& mm);
00342
00344 MatrixGeneral<Treal, Tmatrix>& operator+=
00345 (const XYZ<Treal,
00346 MatrixGeneral<Treal, Tmatrix>,
00347 MatrixGeneral<Treal, Tmatrix> >& smm);
00348
00350 MatrixGeneral<Treal, Tmatrix>& operator=
00351 (const XYZpUV<Treal,
00352 MatrixGeneral<Treal, Tmatrix>,
00353 MatrixGeneral<Treal, Tmatrix>,
00354 Treal,
00355 MatrixGeneral<Treal, Tmatrix> >& smmpsm);
00356
00358 MatrixGeneral<Treal, Tmatrix>& operator=
00359 (XpY<MatrixGeneral<Treal, Tmatrix>,
00360 MatrixGeneral<Treal, Tmatrix> > const & mpm);
00362 MatrixGeneral<Treal, Tmatrix>& operator+=
00363 (MatrixGeneral<Treal, Tmatrix> const & A);
00364 MatrixGeneral<Treal, Tmatrix>& operator-=
00365 (MatrixGeneral<Treal, Tmatrix> const & A);
00367 MatrixGeneral<Treal, Tmatrix>& operator+=
00368 (XY<Treal, MatrixGeneral<Treal, Tmatrix> > const & sm);
00369
00370
00371
00373 MatrixGeneral<Treal, Tmatrix>& operator=
00374 (const XYZ<Treal,
00375 MatrixSymmetric<Treal, Tmatrix>,
00376 MatrixGeneral<Treal, Tmatrix> >& smm);
00378 MatrixGeneral<Treal, Tmatrix>& operator=
00379 (const XY<MatrixSymmetric<Treal, Tmatrix>,
00380 MatrixGeneral<Treal, Tmatrix> >& mm);
00382 MatrixGeneral<Treal, Tmatrix>& operator+=
00383 (const XYZ<Treal,
00384 MatrixSymmetric<Treal, Tmatrix>,
00385 MatrixGeneral<Treal, Tmatrix> >& smm);
00387 MatrixGeneral<Treal, Tmatrix>& operator=
00388 (const XYZpUV<Treal,
00389 MatrixSymmetric<Treal, Tmatrix>,
00390 MatrixGeneral<Treal, Tmatrix>,
00391 Treal,
00392 MatrixGeneral<Treal, Tmatrix> >& smmpsm);
00394 MatrixGeneral<Treal, Tmatrix>& operator=
00395 (const XYZ<Treal,
00396 MatrixGeneral<Treal, Tmatrix>,
00397 MatrixSymmetric<Treal, Tmatrix> >& smm);
00399 MatrixGeneral<Treal, Tmatrix>& operator=
00400 (const XY<MatrixGeneral<Treal, Tmatrix>,
00401 MatrixSymmetric<Treal, Tmatrix> >& mm);
00403 MatrixGeneral<Treal, Tmatrix>& operator+=
00404 (const XYZ<Treal,
00405 MatrixGeneral<Treal, Tmatrix>,
00406 MatrixSymmetric<Treal, Tmatrix> >& smm);
00408 MatrixGeneral<Treal, Tmatrix>& operator=
00409 (const XYZpUV<Treal,
00410 MatrixGeneral<Treal, Tmatrix>,
00411 MatrixSymmetric<Treal, Tmatrix>,
00412 Treal,
00413 MatrixGeneral<Treal, Tmatrix> >& smmpsm);
00415 MatrixGeneral<Treal, Tmatrix>& operator=
00416 (const XYZ<Treal,
00417 MatrixSymmetric<Treal, Tmatrix>,
00418 MatrixSymmetric<Treal, Tmatrix> >& smm);
00420 MatrixGeneral<Treal, Tmatrix>& operator=
00421 (const XY<MatrixSymmetric<Treal, Tmatrix>,
00422 MatrixSymmetric<Treal, Tmatrix> >& mm);
00424 MatrixGeneral<Treal, Tmatrix>& operator+=
00425 (const XYZ<Treal,
00426 MatrixSymmetric<Treal, Tmatrix>,
00427 MatrixSymmetric<Treal, Tmatrix> >& smm);
00429 MatrixGeneral<Treal, Tmatrix>& operator=
00430 (const XYZpUV<Treal,
00431 MatrixSymmetric<Treal, Tmatrix>,
00432 MatrixSymmetric<Treal, Tmatrix>,
00433 Treal,
00434 MatrixGeneral<Treal, Tmatrix> >& smmpsm);
00435
00436
00438 MatrixGeneral<Treal, Tmatrix>& operator=
00439 (const XYZ<Treal,
00440 MatrixTriangular<Treal, Tmatrix>,
00441 MatrixGeneral<Treal, Tmatrix> >& smm);
00443 MatrixGeneral<Treal, Tmatrix>& operator=
00444 (const XYZ<Treal,
00445 MatrixGeneral<Treal, Tmatrix>,
00446 MatrixTriangular<Treal, Tmatrix> >& smm);
00447
00448 void random() {
00449 this->matrixPtr->random();
00450 }
00451
00452 void randomZeroStructure(Treal probabilityBeingZero) {
00453 this->matrixPtr->randomZeroStructure(probabilityBeingZero);
00454 }
00455
00456 template<typename TRule>
00457 void setElementsByRule(TRule & rule) {
00458 this->matrixPtr->setElementsByRule(rule);
00459 return;
00460 }
00461
00462 std::string obj_type_id() const {return "MatrixGeneral";}
00463 protected:
00464 inline void writeToFileProt(std::ofstream & file) const {
00465 this->writeToFileBase(file, matrix_matr);
00466 }
00467 inline void readFromFileProt(std::ifstream & file) {
00468 this->readFromFileBase(file, matrix_matr);
00469 }
00470 private:
00471
00472 };
00473
00474 template<typename Treal, typename Tmatrix>
00475 Treal MatrixGeneral<Treal, Tmatrix>::
00476 eucl(Treal const requestedAccuracy,
00477 int maxIter) const {
00478 VectorType guess;
00479 SizesAndBlocks cols;
00480 this->getCols(cols);
00481 guess.resetSizesAndBlocks(cols);
00482 guess.rand();
00483 mat::ATAMatrix<MatrixGeneral<Treal, Tmatrix>, Treal> ata(*this);
00484 if (maxIter < 0)
00485 maxIter = this->get_nrows() * 100;
00486 arn::LanczosLargestMagnitudeEig
00487 <Treal, ATAMatrix<MatrixGeneral<Treal, Tmatrix>, Treal>, VectorType>
00488 lan(ata, guess, maxIter);
00489 lan.setRelTol( 100 * mat::getMachineEpsilon<Treal>() );
00490 lan.run();
00491 Treal eVal = 0;
00492 Treal acc = 0;
00493 lan.getLargestMagnitudeEig(eVal, acc);
00494 Interval<Treal> euclInt( sqrt(eVal-acc),
00495 sqrt(eVal+acc) );
00496 if ( euclInt.low() < 0 )
00497 euclInt = Interval<Treal>( 0, sqrt(eVal+acc) );
00498 if ( euclInt.length() / 2.0 > requestedAccuracy ) {
00499 std::cout << "req: " << requestedAccuracy
00500 << " obt: " << euclInt.length() / 2.0 << std::endl;
00501 throw std::runtime_error("Desired accuracy not obtained in Lanczos.");
00502 }
00503 return euclInt.midPoint();
00504 }
00505
00506
00507 template<typename Treal, typename Tmatrix>
00508 void MatrixGeneral<Treal, Tmatrix>::
00509 thresh(Treal const threshold,
00510 normType const norm) {
00511 switch (norm) {
00512 case frobNorm:
00513 this->frob_thresh(threshold);
00514 break;
00515 default:
00516 throw Failure("MatrixGeneral<Treal, Tmatrix>::"
00517 "thresh(Treal const, "
00518 "normType const): "
00519 "Thresholding not imlpemented for selected norm");
00520 }
00521 }
00522
00523 template<typename Treal, typename Tmatrix>
00524 Treal MatrixGeneral<Treal, Tmatrix>::
00525 eucl_thresh(Treal const threshold) {
00526 EuclTruncationGeneral<MatrixGeneral<Treal, Tmatrix>, Treal> TruncObj( *this );
00527 return TruncObj.run( threshold );
00528 }
00529
00530
00531
00532
00533 template<typename Treal, typename Tmatrix>
00534 MatrixGeneral<Treal, Tmatrix>&
00535 MatrixGeneral<Treal, Tmatrix>::operator=
00536 (const XYZ<Treal,
00537 MatrixGeneral<Treal, Tmatrix>,
00538 MatrixGeneral<Treal, Tmatrix> >& smm) {
00539 assert(this != &smm.B && this != &smm.C);
00540 this->matrixPtr.haveDataStructureSet(true);
00541 Tmatrix::gemm(smm.tB, smm.tC, smm.A,
00542 *smm.B.matrixPtr, *smm.C.matrixPtr,
00543 0, *this->matrixPtr);
00544 return *this;
00545 }
00546
00547
00548 template<typename Treal, typename Tmatrix>
00549 MatrixGeneral<Treal, Tmatrix>&
00550 MatrixGeneral<Treal, Tmatrix>::operator=
00551 (const XY<MatrixGeneral<Treal, Tmatrix>,
00552 MatrixGeneral<Treal, Tmatrix> >& mm) {
00553 assert(this != &mm.A && this != &mm.B);
00554 Tmatrix::gemm(mm.tA, mm.tB, 1.0,
00555 *mm.A.matrixPtr, *mm.B.matrixPtr,
00556 0, *this->matrixPtr);
00557 return *this;
00558 }
00559
00560
00561 template<typename Treal, typename Tmatrix>
00562 MatrixGeneral<Treal, Tmatrix>&
00563 MatrixGeneral<Treal, Tmatrix>::operator+=
00564 (const XYZ<Treal,
00565 MatrixGeneral<Treal, Tmatrix>,
00566 MatrixGeneral<Treal, Tmatrix> >& smm) {
00567 assert(this != &smm.B && this != &smm.C);
00568 Tmatrix::gemm(smm.tB, smm.tC, smm.A,
00569 *smm.B.matrixPtr, *smm.C.matrixPtr,
00570 1, *this->matrixPtr);
00571 return *this;
00572 }
00573
00574
00575 template<typename Treal, typename Tmatrix>
00576 MatrixGeneral<Treal, Tmatrix>&
00577 MatrixGeneral<Treal, Tmatrix>::operator=
00578 (const XYZpUV<Treal,
00579 MatrixGeneral<Treal, Tmatrix>,
00580 MatrixGeneral<Treal, Tmatrix>,
00581 Treal,
00582 MatrixGeneral<Treal, Tmatrix> >& smmpsm) {
00583 assert(this != &smmpsm.B && this != &smmpsm.C);
00584 assert(!smmpsm.tE);
00585 if (this == &smmpsm.E)
00586 Tmatrix::gemm(smmpsm.tB, smmpsm.tC, smmpsm.A,
00587 *smmpsm.B.matrixPtr, *smmpsm.C.matrixPtr,
00588 smmpsm.D, *this->matrixPtr);
00589 else
00590 throw Failure("MatrixGeneral<Treal, Tmatrix>::operator="
00591 "(const XYZpUV<Treal, MatrixGeneral"
00592 "<Treal, Tmatrix> >&) : D = alpha "
00593 "* op(A) * op(B) + beta * C not supported for C != D");
00594 return *this;
00595 }
00596
00597
00598
00599 template<typename Treal, typename Tmatrix>
00600 inline MatrixGeneral<Treal, Tmatrix>&
00601 MatrixGeneral<Treal, Tmatrix>::operator=
00602 (const XpY<MatrixGeneral<Treal, Tmatrix>,
00603 MatrixGeneral<Treal, Tmatrix> >& mpm) {
00604 assert(this != &mpm.A);
00605 (*this) = mpm.B;
00606 Tmatrix::add(1.0, *mpm.A.matrixPtr, *this->matrixPtr);
00607 return *this;
00608 }
00609
00610 template<typename Treal, typename Tmatrix>
00611 inline MatrixGeneral<Treal, Tmatrix>&
00612 MatrixGeneral<Treal, Tmatrix>::operator+=
00613 (MatrixGeneral<Treal, Tmatrix> const & A) {
00614 Tmatrix::add(1.0, *A.matrixPtr, *this->matrixPtr);
00615 return *this;
00616 }
00617
00618 template<typename Treal, typename Tmatrix>
00619 inline MatrixGeneral<Treal, Tmatrix>&
00620 MatrixGeneral<Treal, Tmatrix>::operator-=
00621 (MatrixGeneral<Treal, Tmatrix> const & A) {
00622 Tmatrix::add(-1.0, *A.matrixPtr, *this->matrixPtr);
00623 return *this;
00624 }
00625
00626 template<typename Treal, typename Tmatrix>
00627 inline MatrixGeneral<Treal, Tmatrix>&
00628 MatrixGeneral<Treal, Tmatrix>::operator+=
00629 (XY<Treal, MatrixGeneral<Treal, Tmatrix> > const & sm) {
00630 assert(!sm.tB);
00631 Tmatrix::add(sm.A, *sm.B.matrixPtr, *this->matrixPtr);
00632 return *this;
00633 }
00634
00635
00636
00637
00638 template<typename Treal, typename Tmatrix>
00639 MatrixGeneral<Treal, Tmatrix>&
00640 MatrixGeneral<Treal, Tmatrix>::operator=
00641 (const XYZ<Treal,
00642 MatrixSymmetric<Treal, Tmatrix>,
00643 MatrixGeneral<Treal, Tmatrix> >& smm) {
00644 assert(this != &smm.C);
00645 assert(!smm.tB && !smm.tC);
00646 this->matrixPtr.haveDataStructureSet(true);
00647 Tmatrix::symm('L', 'U', smm.A,
00648 *smm.B.matrixPtr, *smm.C.matrixPtr,
00649 0, *this->matrixPtr);
00650 return *this;
00651 }
00652
00653
00654 template<typename Treal, typename Tmatrix>
00655 MatrixGeneral<Treal, Tmatrix>&
00656 MatrixGeneral<Treal, Tmatrix>::operator=
00657 (const XY<MatrixSymmetric<Treal, Tmatrix>,
00658 MatrixGeneral<Treal, Tmatrix> >& mm) {
00659 assert(this != &mm.B);
00660 assert(!mm.tB);
00661 Tmatrix::symm('L', 'U', 1.0,
00662 *mm.A.matrixPtr, *mm.B.matrixPtr,
00663 0, *this->matrixPtr);
00664 return *this;
00665 }
00666
00667
00668 template<typename Treal, typename Tmatrix>
00669 MatrixGeneral<Treal, Tmatrix>&
00670 MatrixGeneral<Treal, Tmatrix>::operator+=
00671 (const XYZ<Treal,
00672 MatrixSymmetric<Treal, Tmatrix>,
00673 MatrixGeneral<Treal, Tmatrix> >& smm) {
00674 assert(this != &smm.C);
00675 assert(!smm.tB && !smm.tC);
00676 Tmatrix::symm('L', 'U', smm.A,
00677 *smm.B.matrixPtr, *smm.C.matrixPtr,
00678 1, *this->matrixPtr);
00679 return *this;
00680 }
00681
00682
00683 template<typename Treal, typename Tmatrix>
00684 MatrixGeneral<Treal, Tmatrix>&
00685 MatrixGeneral<Treal, Tmatrix>::operator=
00686 (const XYZpUV<Treal,
00687 MatrixSymmetric<Treal, Tmatrix>,
00688 MatrixGeneral<Treal, Tmatrix>,
00689 Treal,
00690 MatrixGeneral<Treal, Tmatrix> >& smmpsm) {
00691 assert(this != &smmpsm.C);
00692 assert(!smmpsm.tB && !smmpsm.tC && !smmpsm.tE);
00693 if (this == &smmpsm.E)
00694 Tmatrix::symm('L', 'U', smmpsm.A,
00695 *smmpsm.B.matrixPtr, *smmpsm.C.matrixPtr,
00696 smmpsm.D, *this->matrixPtr);
00697 else
00698 throw Failure("MatrixGeneral<Treal, Tmatrix>::operator="
00699 "(const XYZpUV<Treal, MatrixGeneral"
00700 "<Treal, Tmatrix>, MatrixSymmetric<Treal, "
00701 "Tmatrix>, Treal, MatrixGeneral"
00702 "<Treal, Tmatrix> >&) "
00703 ": D = alpha * A * B + beta * C (with A symmetric)"
00704 " not supported for C != D");
00705 return *this;
00706 }
00707
00708
00709 template<typename Treal, typename Tmatrix>
00710 MatrixGeneral<Treal, Tmatrix>&
00711 MatrixGeneral<Treal, Tmatrix>::operator=
00712 (const XYZ<Treal,
00713 MatrixGeneral<Treal, Tmatrix>,
00714 MatrixSymmetric<Treal, Tmatrix> >& smm) {
00715 assert(this != &smm.B);
00716 assert(!smm.tB && !smm.tC);
00717 this->matrixPtr.haveDataStructureSet(true);
00718 Tmatrix::symm('R', 'U', smm.A,
00719 *smm.C.matrixPtr, *smm.B.matrixPtr,
00720 0, *this->matrixPtr);
00721 return *this;
00722 }
00723
00724
00725 template<typename Treal, typename Tmatrix>
00726 MatrixGeneral<Treal, Tmatrix>&
00727 MatrixGeneral<Treal, Tmatrix>::operator=
00728 (const XY<MatrixGeneral<Treal, Tmatrix>,
00729 MatrixSymmetric<Treal, Tmatrix> >& mm) {
00730 assert(this != &mm.A);
00731 assert(!mm.tA && !mm.tB);
00732 Tmatrix::symm('R', 'U', 1.0,
00733 *mm.B.matrixPtr, *mm.A.matrixPtr,
00734 0, *this->matrixPtr);
00735 return *this;
00736 }
00737
00738
00739 template<typename Treal, typename Tmatrix>
00740 MatrixGeneral<Treal, Tmatrix>&
00741 MatrixGeneral<Treal, Tmatrix>::operator+=
00742 (const XYZ<Treal,
00743 MatrixGeneral<Treal, Tmatrix>,
00744 MatrixSymmetric<Treal, Tmatrix> >& smm) {
00745 assert(this != &smm.B);
00746 assert(!smm.tB && !smm.tC);
00747 Tmatrix::symm('R', 'U', smm.A,
00748 *smm.C.matrixPtr, *smm.B.matrixPtr,
00749 1, *this->matrixPtr);
00750 return *this;
00751 }
00752
00753
00754 template<typename Treal, typename Tmatrix>
00755 MatrixGeneral<Treal, Tmatrix>&
00756 MatrixGeneral<Treal, Tmatrix>::operator=
00757 (const XYZpUV<Treal,
00758 MatrixGeneral<Treal, Tmatrix>,
00759 MatrixSymmetric<Treal, Tmatrix>,
00760 Treal,
00761 MatrixGeneral<Treal, Tmatrix> >& smmpsm) {
00762 assert(this != &smmpsm.B);
00763 assert(!smmpsm.tB && !smmpsm.tC && !smmpsm.tE);
00764 if (this == &smmpsm.E)
00765 Tmatrix::symm('R', 'U', smmpsm.A,
00766 *smmpsm.C.matrixPtr, *smmpsm.B.matrixPtr,
00767 smmpsm.D, *this->matrixPtr);
00768 else
00769 throw Failure("MatrixGeneral<Treal, Tmatrix>::operator="
00770 "(const XYZpUV<Treal, MatrixSymmetric"
00771 "<Treal, Tmatrix>, MatrixGeneral<Treal, "
00772 "Tmatrix>, Treal, MatrixGeneral"
00773 "<Treal, Tmatrix> >&) "
00774 ": D = alpha * B * A + beta * C (with A symmetric)"
00775 " not supported for C != D");
00776 return *this;
00777 }
00778
00779
00781 template<typename Treal, typename Tmatrix>
00782 MatrixGeneral<Treal, Tmatrix>&
00783 MatrixGeneral<Treal, Tmatrix>::operator=
00784 (const XYZ<Treal,
00785 MatrixSymmetric<Treal, Tmatrix>,
00786 MatrixSymmetric<Treal, Tmatrix> >& smm) {
00787 assert(!smm.tB && !smm.tC);
00788 this->matrixPtr.haveDataStructureSet(true);
00789 Tmatrix::ssmm(smm.A,
00790 *smm.B.matrixPtr,
00791 *smm.C.matrixPtr,
00792 0,
00793 *this->matrixPtr);
00794 return *this;
00795 }
00796
00798 template<typename Treal, typename Tmatrix>
00799 MatrixGeneral<Treal, Tmatrix>&
00800 MatrixGeneral<Treal, Tmatrix>::operator=
00801 (const XY<MatrixSymmetric<Treal, Tmatrix>,
00802 MatrixSymmetric<Treal, Tmatrix> >& mm) {
00803 assert(!mm.tB);
00804 Tmatrix::ssmm(1.0,
00805 *mm.A.matrixPtr,
00806 *mm.B.matrixPtr,
00807 0,
00808 *this->matrixPtr);
00809 return *this;
00810 }
00811
00813 template<typename Treal, typename Tmatrix>
00814 MatrixGeneral<Treal, Tmatrix>&
00815 MatrixGeneral<Treal, Tmatrix>::operator+=
00816 (const XYZ<Treal,
00817 MatrixSymmetric<Treal, Tmatrix>,
00818 MatrixSymmetric<Treal, Tmatrix> >& smm) {
00819 assert(!smm.tB && !smm.tC);
00820 Tmatrix::ssmm(smm.A,
00821 *smm.B.matrixPtr,
00822 *smm.C.matrixPtr,
00823 1,
00824 *this->matrixPtr);
00825 return *this;
00826 }
00827
00828
00830 template<typename Treal, typename Tmatrix>
00831 MatrixGeneral<Treal, Tmatrix>&
00832 MatrixGeneral<Treal, Tmatrix>::operator=
00833 (const XYZpUV<Treal,
00834 MatrixSymmetric<Treal, Tmatrix>,
00835 MatrixSymmetric<Treal, Tmatrix>,
00836 Treal,
00837 MatrixGeneral<Treal, Tmatrix> >& smmpsm) {
00838 assert(!smmpsm.tB && !smmpsm.tC && !smmpsm.tE);
00839 if (this == &smmpsm.E)
00840 Tmatrix::ssmm(smmpsm.A,
00841 *smmpsm.B.matrixPtr,
00842 *smmpsm.C.matrixPtr,
00843 smmpsm.D,
00844 *this->matrixPtr);
00845 else
00846 throw Failure("MatrixGeneral<Treal, Tmatrix>::"
00847 "operator=(const XYZpUV<"
00848 "Treal, MatrixSymmetric<Treal, Tmatrix>, "
00849 "MatrixSymmetric<Treal, Tmatrix>, Treal, "
00850 "MatrixGeneral<Treal, Tmatrix> >&) "
00851 ": D = alpha * A * B + beta * C (with A and B symmetric)"
00852 " not supported for C != D");
00853 return *this;
00854 }
00855
00856
00857
00858
00859
00860
00861 template<typename Treal, typename Tmatrix>
00862 MatrixGeneral<Treal, Tmatrix>&
00863 MatrixGeneral<Treal, Tmatrix>::operator=
00864 (const XYZ<Treal,
00865 MatrixTriangular<Treal, Tmatrix>,
00866 MatrixGeneral<Treal, Tmatrix> >& smm) {
00867 assert(!smm.tC);
00868 if (this == &smm.C)
00869 Tmatrix::trmm('L', 'U', smm.tB, smm.A,
00870 *smm.B.matrixPtr, *this->matrixPtr);
00871 else
00872 throw Failure("MatrixGeneral<Treal, Tmatrix>::operator="
00873 "(const XYZ<Treal, MatrixTriangular"
00874 "<Treal, Tmatrix>, MatrixGeneral<Treal,"
00875 " Tmatrix> >& : D = alpha * op(A) * B (with"
00876 " A upper triangular) not supported for B != D");
00877 return *this;
00878 }
00879
00880
00881
00882 template<typename Treal, typename Tmatrix>
00883 MatrixGeneral<Treal, Tmatrix>&
00884 MatrixGeneral<Treal, Tmatrix>::operator=
00885 (const XYZ<Treal,
00886 MatrixGeneral<Treal, Tmatrix>,
00887 MatrixTriangular<Treal, Tmatrix> >& smm) {
00888 assert(!smm.tB);
00889 if (this == &smm.B)
00890 Tmatrix::trmm('R', 'U', smm.tC, smm.A,
00891 *smm.C.matrixPtr, *this->matrixPtr);
00892 else
00893 throw Failure("MatrixGeneral<Treal, Tmatrix>::operator="
00894 "(const XYZ<Treal, MatrixGeneral"
00895 "<Treal, Tmatrix>, MatrixTriangular<Treal,"
00896 " Tmatrix> >& : D = alpha * A * op(B) (with"
00897 " B upper triangular) not supported for A != D");
00898 return *this;
00899 }
00900
00901
00902
00903 template<typename Treal, typename Tmatrix>
00904 Treal trace(const XYZ<Treal,
00905 MatrixGeneral<Treal, Tmatrix>,
00906 MatrixGeneral<Treal, Tmatrix> >& smm) {
00907 if ((!smm.tB && !smm.tC) || (smm.tB && smm.tC))
00908 return smm.A * MatrixGeneral<Treal, Tmatrix>::
00909 trace_ab(smm.B, smm.C);
00910 else if (smm.tB)
00911 return smm.A * MatrixGeneral<Treal, Tmatrix>::
00912 trace_aTb(smm.B, smm.C);
00913 else
00914 return smm.A * MatrixGeneral<Treal, Tmatrix>::
00915 trace_aTb(smm.C, smm.B);
00916 }
00917
00918
00919 }
00920 #endif
00921
00922