00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #ifndef QC_MATRIX_H_
00021 #define QC_MATRIX_H_
00022
00023 #include <string>
00024 #include <fstream>
00025
00026 #include "QCCommon.hpp"
00027 #include "QCBlas.hpp"
00028 #include "QCMemory.hpp"
00029
00030
00031
00032 class QCSymMatrix;
00033
00039 class QCMatrix
00040 {
00041
00043 public:
00044
00048 QCMatrix (void);
00049
00050
00056 QCMatrix (const int& dimMatrix);
00057
00058
00062 virtual ~QCMatrix (void);
00063
00064
00068 void setDim (const int& dimMatrix);
00069
00070
00074 bool setDimAndAllocate (const int& dimMatrix);
00075
00076
00081 void setVectorPtr (const int indexOfTheVectorPtr,
00082 const int indexInTheElemsArray);
00083
00087 bool takeTmpMem (QCMemory& QCRestrict memoryInst);
00088
00089
00093 bool giveBackTmpMem (QCMemory& QCRestrict memoryInst);
00094
00095
00099 void initClean (void);
00100
00101
00105 int getDim (void) const { return dim; }
00106 int getNbElems (void) const { return nbElems; }
00107 QCFloat* QCRestrict getElems (void) const { return elems; }
00108 QCFloat* QCRestrict * QCRestrict getVectorPtr (void) const { return vectorPtr; }
00109
00110
00114 void print (char* matrixName) const;
00115 void print (const std::string & matrixName, std::ofstream &out) const;
00116 void print (char* matrixName, const int& dim1, const int& dim2) const;
00117
00118
00123 void printInFile (const char* preFileName, bool date = true);
00124
00125
00130 void sampleInFile (const char* preFileName);
00131
00136 void sampleDistanceInFile (const char* preFileName);
00137
00138
00143 QCFloat* QCRestrict operator [] (int i) const { return vectorPtr[i]; }
00144
00145
00149 QCMatrix& QCRestrict operator = (const QCMatrix& QCRestrict aMatrix)
00150 { copy(aMatrix) ; return *this ;}
00151
00152
00157 QCFloat elemMax (void) const;
00158
00159
00164 QCFloat deltaElemMax (const QCMatrix& otherMatrix) const;
00165
00166
00171 inline
00172 void copy (const QCMatrix& otherMatrix);
00173
00174
00179 inline
00180 void mult (QCMemory& QCRestrict memoryInst,
00181 const QCMatrix& matrixA,
00182 const QCMatrix& matrixB,
00183 const QCFloat& alpha = QC_ONE,
00184 const QCFloat& beta = QC_ZERO);
00185
00186
00191 void mult (QCMemory& QCRestrict memoryInst,
00192 const QCSymMatrix& matrixA,
00193 const QCMatrix& matrixB,
00194 const QCFloat& alpha = QC_ONE,
00195 const QCFloat& beta = QC_ZERO);
00196
00197
00202 void mult (QCMemory& QCRestrict memoryInst,
00203 const QCMatrix& matrixA,
00204 const QCSymMatrix& matrixB,
00205 const QCFloat& alpha = QC_ONE,
00206 const QCFloat& beta = QC_ZERO);
00207
00208
00213 void mult (QCMemory& QCRestrict memoryInst,
00214 const QCSymMatrix& matrixA,
00215 const QCSymMatrix& matrixB,
00216 const QCFloat& alpha = QC_ONE,
00217 const QCFloat& beta = QC_ZERO);
00218
00219
00226 inline
00227 void swap (QCMatrix& otherMatrix);
00228
00229
00233 inline
00234 void extractSubMatrix (QCMatrix& QCRestrict subMatrix,
00235 const int row, const int nbRows,
00236 const int column, const int nbColumns) const;
00237
00238
00242 inline
00243 void insertSubMatrix (const QCMatrix& QCRestrict subMatrix,
00244 const int row, const int nbRows,
00245 const int column, const int nbColumns);
00246
00247
00253 inline
00254 void pack (const int& longDim);
00255
00256
00260 inline
00261 void unpack (const int& longDim);
00262
00263
00267 inline
00268 void transpose (const QCMatrix& QCRestrict matrixToTranspose);
00269
00270 inline
00271 void transpose (const QCMatrix& QCRestrict matrixToTranspose,
00272 const int nbRowInMatrixToTranspose,
00273 const int nbColInMatrixToTranspose);
00274
00275
00279 inline
00280 void calculateNumericalDerivatives (const QCMatrix& matrix2,
00281 const QCFloat& QCRestrict inv2DeltaQ);
00282
00283
00287 inline
00288 void calculateNumericalDerivatives (const QCMatrix& matrix2,
00289 const QCFloat& QCRestrict inv2DeltaQ,
00290 const int dim1,
00291 const int dim2);
00292
00293
00297 QCFloat trace (void) const;
00298
00299
00304 template <class TPManager>
00305 void writeLikeBands (TPManager& QCRestrict managerInst) const;
00306
00307
00311 inline
00312 void scale (const QCFloat & scaleFactor ) ;
00313
00314
00318 inline
00319 void scaleAndCopy (const QCFloat & scaleFactor,
00320 const QCMatrix& QCRestrict A) ;
00321
00322
00326 inline
00327 void axpy (const QCFloat & scaleFactor,
00328 const QCMatrix& QCRestrict A) ;
00329
00330
00331
00332
00333
00334
00336 protected:
00337
00338
00340 private:
00341
00342
00345 public:
00346 protected:
00347
00353 QCFloat* QCRestrict * QCRestrict vectorPtr;
00354
00355
00359 QCFloat* QCRestrict elems;
00360
00361
00365 int dim;
00366
00367
00371 int nbElems;
00372
00373
00380 bool allocateByNew;
00381
00382
00383
00384 private:
00385
00386 };
00387
00388
00389
00394 inline void
00395 QCMatrix::copy (const QCMatrix& otherMatrix) {
00396
00397 #ifdef DEBUG_QC
00398 if (dim == otherMatrix.getDim() ) {
00399 if (elems && otherMatrix.elems && elems != otherMatrix.elems ) {
00400 #endif // DEBUG_QC
00401 memcpy(elems, otherMatrix.elems, nbElems*sizeof(QCFloat) );
00402 #ifdef DEBUG_QC
00403 }
00404
00405 } else {
00406
00407
00408 cout << "Pb in File " << __FILE__ << " at line " << __LINE__ << endl;
00409 abort();
00410 }
00411 #endif //DEBUG_QC
00412 }
00413
00414
00415
00416
00424 inline void
00425 QCMatrix::mult (QCMemory& QCRestrict ,
00426 const QCMatrix& matrixA,
00427 const QCMatrix& matrixB,
00428 const QCFloat& alpha,
00429 const QCFloat& beta) {
00430 QCCBlas(dgemm)(
00431 #ifdef USE_OPTIMIZED_CBLAS
00432 CblasRowMajor,
00433 CblasNoTrans, CblasTrans,
00434 #else //USE_OPTIMIZED_CBLAS
00435 CblasTrans, CblasNoTrans,
00436 #endif //USE_OPTIMIZED_CBLAS
00437 dim, dim, dim, alpha, matrixA.elems, dim, matrixB.elems,
00438 dim, beta, elems, dim);
00439 }
00440
00441
00442
00443
00449 inline void
00450 QCMatrix::mult (QCMemory& QCRestrict memoryInst,
00451 const QCMatrix& matrixA,
00452 const QCSymMatrix& matrixB,
00453 const QCFloat& alpha,
00454 const QCFloat& beta) {
00455 mult(memoryInst, matrixB, matrixA, alpha, beta);
00456 }
00457
00458
00459
00460
00467 inline void
00468 QCMatrix::swap (QCMatrix& otherMatrix) {
00469 QCFloat* QCRestrict interElems;
00470 QCFloat* QCRestrict * QCRestrict interVectorPtr;
00471 #ifdef DEBUG_QC
00472 if (dim == otherMatrix.getDim() ) {
00473 if (elems && otherMatrix.elems && elems != otherMatrix.elems ) {
00474 #endif //DEBUG_QC
00475
00476
00477 interElems = elems;
00478 elems = otherMatrix.elems;
00479 otherMatrix.elems = interElems;
00480
00481
00482 interVectorPtr = vectorPtr;
00483 vectorPtr = otherMatrix.vectorPtr;
00484 otherMatrix.vectorPtr = interVectorPtr;
00485 #ifdef DEBUG_QC
00486 }
00487
00488 } else {
00489
00490
00491 cout << "Pb in File " << __FILE__ << " at line " << __LINE__ << endl;
00492 abort();
00493 }
00494 #endif //DEBUG_QC
00495 }
00496
00497
00498
00499
00503 inline void
00504 QCMatrix::extractSubMatrix (QCMatrix& QCRestrict subMatrix,
00505 const int row, const int nbRows,
00506 const int column, const int nbColumns) const {
00507 #ifdef OLD_JOK3
00508 int matRowInd, subMatRowInd;
00509
00510
00511
00512 for (subMatRowInd = 0, matRowInd = row;
00513 subMatRowInd < nbRows;
00514 ++subMatRowInd, ++matRowInd) {
00515 memcpy(subMatrix[subMatRowInd],
00516 &(vectorPtr[matRowInd][column]), nbColumns*sizeof(QCFloat) );
00517 }
00518 #else // OLD_JOK3
00519 int matRowInd, matColInd;
00520 int subMatRowInd, subMatColInd;
00521
00522 for (subMatRowInd = 0, matRowInd = row;
00523 subMatRowInd < nbRows;
00524 ++subMatRowInd, ++matRowInd) {
00525 for (subMatColInd = 0, matColInd = column;
00526 subMatColInd < nbColumns;
00527 ++subMatColInd, ++matColInd) {
00528 subMatrix[subMatRowInd][subMatColInd] = vectorPtr[matRowInd][matColInd];
00529 }
00530 }
00531 #endif // OLD_JOK3
00532 }
00533
00534
00535
00536
00540 inline void
00541 QCMatrix::insertSubMatrix (const QCMatrix& QCRestrict subMatrix,
00542 const int row, const int nbRows,
00543 const int column, const int nbColumns) {
00544 #ifdef OLD_JOK3
00545 int matRowInd, subMatRowInd;
00546
00547
00548
00549 for (subMatRowInd = 0, matRowInd = row;
00550 subMatRowInd < nbRows;
00551 ++subMatRowInd, ++matRowInd) {
00552 memcpy(&(vectorPtr[matRowInd][column]),
00553 subMatrix[subMatRowInd], nbColumns*sizeof(QCFloat) );
00554 }
00555 #else // OLD_JOK3
00556 int matRowInd, matColInd;
00557 int subMatRowInd, subMatColInd;
00558
00559 for (subMatRowInd = 0, matRowInd = row;
00560 subMatRowInd < nbRows;
00561 ++subMatRowInd, ++matRowInd) {
00562 for (subMatColInd = 0, matColInd = column;
00563 subMatColInd < nbColumns;
00564 ++subMatColInd, ++matColInd) {
00565 vectorPtr[matRowInd][matColInd] = subMatrix[subMatRowInd][subMatColInd];
00566 }
00567 }
00568 #endif // OLD_JOK3
00569 }
00570
00571
00572
00573
00578 inline void
00579 QCMatrix::pack (const int& longDim) {
00580 int i;
00581 for (i = 0; i < longDim; ++i) {
00582 vectorPtr[0][i] = vectorPtr[i][0];
00583 }
00584 }
00585
00586
00587 inline void
00588 QCMatrix::unpack (const int& longDim) {
00589 int i;
00590 for (i = 0; i < longDim; ++i) {
00591 vectorPtr[i][0] = vectorPtr[0][i];
00592 }
00593 }
00594
00595
00596
00597
00601 inline void
00602 QCMatrix::transpose (const QCMatrix& QCRestrict matrixToTranspose) {
00603 int i, j;
00604
00605
00606 for (i = 0; i < dim; ++i) {
00607 for (j = 0; j < dim; ++j) {
00608 vectorPtr[i][j] = matrixToTranspose[j][i];
00609 }
00610 }
00611 }
00612
00613
00614 inline void
00615 QCMatrix::transpose (const QCMatrix& QCRestrict matrixToTranspose,
00616 const int nbRowInMatrixToTranspose,
00617 const int nbColInMatrixToTranspose) {
00618 int i, j;
00619
00620
00621 for (i = 0; i < nbColInMatrixToTranspose; ++i) {
00622 for (j = 0; j < nbRowInMatrixToTranspose; ++j) {
00623 vectorPtr[i][j] = matrixToTranspose[j][i];
00624 }
00625 }
00626 }
00627
00628
00629
00630
00631
00638 inline void
00639 QCMatrix::calculateNumericalDerivatives(const QCMatrix& matrix2,
00640 const QCFloat& QCRestrict inv2DeltaQ)
00641 {
00642 int i;
00643 #ifdef DEBUG_QC
00644 if (nbElems == matrix2.nbElems) {
00645 #endif //DEBUG_QC
00646
00647
00648
00649 if (inv2DeltaQ != QC_ZERO) {
00650 for (i = 0; i < nbElems; ++i) {
00651 elems[i] = (elems[i] - matrix2.elems[i]) * inv2DeltaQ;
00652 }
00653
00654 } else {
00655 for (i = 0; i < nbElems; ++i) {
00656 elems[i] -= matrix2.elems[i];
00657 }
00658 }
00659 #ifdef DEBUG_QC
00660
00661 } else {
00662
00663
00664 cout << "Pb in File " << __FILE__ << " at line " << __LINE__ << endl;
00665 abort();
00666 }
00667 #endif //DEBUG_QC
00668 }
00669
00670
00671
00672
00679 inline void
00680 QCMatrix::calculateNumericalDerivatives (const QCMatrix& matrix2,
00681 const QCFloat& QCRestrict inv2DeltaQ,
00682 const int dim1,
00683 const int dim2)
00684 {
00685 int i,j;
00686 #ifdef DEBUG_QC
00687 if (nbElems == matrix2.nbElems) {
00688 #endif //DEBUG_QC
00689
00690
00691
00692 if (inv2DeltaQ != QC_ZERO) {
00693 for (i = 0; i < dim1; ++i) {
00694 for (j = 0; j < dim2; ++j) {
00695 vectorPtr[i][j] = ((vectorPtr[i][j] - matrix2.vectorPtr[i][j])
00696 * inv2DeltaQ);
00697 }
00698 }
00699
00700 } else {
00701 for (i = 0; i < dim1; ++i) {
00702 for (j = 0; j < dim2; ++j) {
00703 vectorPtr[i][j] -= matrix2.vectorPtr[i][j];
00704 }
00705 }
00706 }
00707 #ifdef DEBUG_QC
00708
00709 } else {
00710
00711
00712 cout << "Pb in File " << __FILE__ << " at line " << __LINE__ << endl;
00713 abort();
00714 }
00715 #endif //DEBUG_QC
00716 }
00717
00718
00719
00720
00724 inline void
00725 QCMatrix::scale (const QCFloat& scaleFactor ) {
00726 QCCBlas(dscal)(nbElems, scaleFactor, elems, 1);
00727 }
00728
00729
00730
00731
00735 inline void
00736 QCMatrix::scaleAndCopy (const QCFloat & scaleFactor,
00737 const QCMatrix& QCRestrict A) {
00738 for (int i = 0; i < nbElems ; ++i){
00739 elems[i] = scaleFactor * A.elems[i] ;
00740 }
00741 }
00742
00743
00744
00745
00749 inline void
00750 QCMatrix::axpy (const QCFloat & scaleFactor,
00751 const QCMatrix& QCRestrict A) {
00752 QCCBlas(daxpy)(nbElems, scaleFactor, const_cast<QCFloat*>(A.elems),
00753 1, const_cast<QCFloat*>(elems), 1);
00754 }
00755
00756
00757
00758
00759 #endif // QC_MATRIX_H_