00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #ifndef QC_SYMMATRIX_H_
00021 #define QC_SYMMATRIX_H_
00022
00023 #include <string>
00024 #include <fstream>
00025
00026 #include "QCMatrix.hpp"
00027
00028
00029
00035 class QCSymMatrix : public QCMatrix
00036 {
00037
00039 public:
00040
00044 QCSymMatrix (void);
00045
00046
00051 QCSymMatrix (const int& dimMatrix);
00052
00053
00057 virtual ~QCSymMatrix (void);
00058
00059
00063 void setDim (const int& dimMatrix);
00067 bool setDimAndAllocate (const int& dimMatrix);
00071 bool setDimAndElems (const int& dimMatrix, QCFloat * extElems);
00072
00076 void unsetElems (void);
00080 bool takeTmpMem (QCMemory& QCRestrict memoryInst);
00084 void print (char* matrixName) const;
00087 void print (const std::string & matrixName, std::ofstream &out) ;
00088
00089 void print (char* matrixName, const int& dim1) const;
00090
00091 void printInFile (const char* preFileName, bool date = true);
00092
00093
00098 void printFullInFile (const char* preFileName);
00099
00100
00105 QCFloat* QCRestrict operator [] (int i) const { return vectorPtr[i]; }
00106
00107
00113 bool diagonalize (QCMemory& QCRestrict memoryInst,
00114 QCDiagoAlgorithm& diagoAlgorithm,
00115 const QCFloat& toleranceFactor,
00116 QCFloat* QCRestrict eigenVal,
00117 int* QCRestrict degener,
00118 QCFloat* QCRestrict eigenVect,
00119 const int& nbIte);
00120
00121
00125 QCFloat traceProduct (const QCSymMatrix& otherMatrix) const;
00126
00127
00132 inline void power (QCMemory& QCRestrict memoryInst,
00133 const QCSymMatrix& matrixA,
00134 const QCFloat& alpha = QC_ONE);
00135
00136
00141 void power (QCMemory& QCRestrict memoryInst,
00142 const QCSymMatrix& matrixA,
00143 const int& dimToMultiplyOfA,
00144 const QCFloat& alpha = QC_ONE);
00145
00146
00156 inline void copyIn (const QCMatrix& otherMatrix) const;
00157
00158
00162 inline void copyInAndFill (const QCMatrix& otherMatrix) const;
00163
00164
00170 inline void copy (const QCMatrix& otherMatrix) const;
00171
00172
00176 inline void extractSubTriangle (QCSymMatrix& QCRestrict subMatrix,
00177 int row, int nbRows) const;
00178
00179
00184 inline void extractSubTriangleInSquare (QCMatrix& QCRestrict subSquare,
00185 int row, int nbRows) const;
00186
00187
00191 inline void insertSubTriangle (const QCSymMatrix& QCRestrict subMatrix,
00192 int row, int nbRows);
00193
00194
00198 inline void completeSubTriangle (const QCSymMatrix& QCRestrict subMatrix,
00199 int row, int nbRows);
00200
00201
00205 inline void changeComputedSystem (void) { computedSystemHasChanged = true; }
00206
00207
00212 template <class TPManager>
00213 void writeLikeBlocks(TPManager& QCRestrict managerInst) const;
00214
00215
00216
00217
00218
00219
00221 protected:
00222
00223
00225 private:
00226
00227 #ifdef USE_LAPACK
00228
00232 void updateDiagonalizeLocalVar (QCDiagoAlgorithm& diagoAlgorithm,
00233 const int& dimOfThePartToDiago,
00234 int& nb,
00235 int& workingDim,
00236 int& workingArray1Size,
00237 int& workingArray2Size);
00238 #endif //USE_LAPACK
00239
00240
00243 public:
00244 protected:
00245 private:
00246
00252 bool computedSystemHasChanged;
00253
00254 };
00255
00256
00257
00258
00264 inline void
00265 QCSymMatrix::copyIn (const QCMatrix& otherMatrix) const {
00266 int i;
00267 #ifdef DEBUG_QC
00268 if (dim == otherMatrix.getDim() ) {
00269 if (elems && otherMatrix.getElems() ) {
00270 #endif //DEBUG_QC
00271 for (i = 0; i < dim; ++i) {
00272
00273
00274
00275 memcpy(otherMatrix.getVectorPtr()[i], vectorPtr[i], (i+1)*sizeof(QCFloat) );
00276 }
00277 #ifdef DEBUG_QC
00278 }
00279 } else {
00280
00281
00282 cout << "Pb in File " << __FILE__ << " at line " << __LINE__ << endl;
00283 abort();
00284 }
00285 #endif //DEBUG_QC
00286 }
00287
00288
00289
00295 inline void
00296 QCSymMatrix::copy (const QCMatrix& otherMatrix) const {
00297 int i;
00298 #ifdef DEBUG_QC
00299 if (dim == otherMatrix.getDim() ) {
00300 if (elems && otherMatrix.getElems() ) {
00301 #endif //DEBUG_QC
00302 for (i = 0; i < dim; ++i) {
00303
00304
00305
00306 memcpy(vectorPtr[i],otherMatrix.getVectorPtr()[i], (i+1)*sizeof(QCFloat) );
00307 }
00308 #ifdef DEBUG_QC
00309 }
00310 } else {
00311
00312
00313 cout << "Pb in File " << __FILE__ << " at line " << __LINE__ << endl;
00314 abort();
00315 }
00316 #endif //DEBUG_QC
00317 }
00318
00319
00320
00327 inline void
00328 QCSymMatrix::copyInAndFill (const QCMatrix& otherMatrix) const {
00329 int i, j;
00330
00331 copyIn(otherMatrix);
00332
00333
00334 for (i = 0; i < dim; i++) {
00335 for (j = i+1; j < dim; j++) {
00336 otherMatrix[i][j] = vectorPtr[j][i];
00337 }
00338 }
00339 }
00340
00341
00342
00343
00348 inline void
00349 QCSymMatrix::extractSubTriangle (QCSymMatrix& QCRestrict subMatrix,
00350 int row, int nbRows) const {
00351 #ifdef OLD_JOK3
00352
00353
00354
00355 int matRowInd, subMatRowInd;
00356
00357 for (subMatRowInd = 0, matRowInd = row;
00358 subMatRowInd < nbRows;
00359 ++subMatRowInd, ++matRowInd) {
00360 memcpy(subMatrix[subMatRowInd],
00361 &(vectorPtr[matRowInd][row]),(subMatRowInd+1)*sizeof(QCFloat) );
00362 }
00363 #else //OLD_JOK3
00364 int matRowInd, matColInd;
00365 int subMatRowInd, subMatColInd;
00366
00367 for (subMatRowInd = 0, matRowInd = row;
00368 subMatRowInd < nbRows;
00369 ++subMatRowInd, ++matRowInd) {
00370 for (subMatColInd = 0, matColInd = row;
00371 subMatColInd <= subMatRowInd;
00372 ++subMatColInd, ++matColInd) {
00373 subMatrix[subMatRowInd][subMatColInd] = vectorPtr[matRowInd][matColInd];
00374 }
00375 }
00376 #endif //OLD_JOK3
00377 }
00378
00379
00380
00381
00386 inline void
00387 QCSymMatrix::extractSubTriangleInSquare (QCMatrix& QCRestrict subSquare,
00388 int row, int nbRows) const {
00389 int matRowInd, matColInd;
00390 int subMatRowInd, subMatColInd;
00391
00392 for (subMatRowInd = 0, matRowInd = row;
00393 subMatRowInd < nbRows;
00394 ++subMatRowInd, ++matRowInd) {
00395 for (subMatColInd = 0, matColInd = row;
00396 subMatColInd < nbRows;
00397 ++subMatColInd, ++matColInd) {
00398 if (matColInd <= matRowInd) {
00399 subSquare[subMatRowInd][subMatColInd] = vectorPtr[matRowInd][matColInd];
00400
00401 } else {
00402 subSquare[subMatRowInd][subMatColInd] = vectorPtr[matColInd][matRowInd];
00403 }
00404 }
00405 }
00406 }
00407
00408
00409
00410
00414 inline void
00415 QCSymMatrix::insertSubTriangle (const QCSymMatrix& QCRestrict subMatrix,
00416 int row, int nbRows) {
00417 #ifdef OLD_JOK3
00418
00419
00420
00421 int matRowInd, subMatRowInd;
00422 for (subMatRowInd = 0, matRowInd = row;
00423 subMatRowInd < nbRows;
00424 ++subMatRowInd, ++matRowInd) {
00425 memcpy(&(vectorPtr[matRowInd][row]),
00426
00427
00428 subMatrix[subMatRowInd], (subMatRowInd+1)*sizeof(QCFloat) );
00429 }
00430 #else //OLD_JOK3
00431 int matRowInd, matColInd;
00432 int subMatRowInd, subMatColInd;
00433
00434 for (subMatRowInd = 0, matRowInd = row;
00435 subMatRowInd < nbRows;
00436 ++subMatRowInd, ++matRowInd) {
00437 for (subMatColInd = 0, matColInd = row;
00438 subMatColInd <= subMatRowInd;
00439 ++subMatColInd, ++matColInd) {
00440 vectorPtr[matRowInd][matColInd] = subMatrix[subMatRowInd][subMatColInd];
00441 }
00442 }
00443 #endif //OLD_JOK3
00444 }
00445
00446
00447
00448
00449
00453 inline void
00454 QCSymMatrix::completeSubTriangle (const QCSymMatrix& QCRestrict subMatrix,
00455 int row, int nbRows) {
00456
00457 int matRowInd, matColInd;
00458 int subMatRowInd, subMatColInd;
00459
00460 for (subMatRowInd = 0, matRowInd = row;
00461 subMatRowInd < nbRows;
00462 ++subMatRowInd, ++matRowInd) {
00463 for (subMatColInd = 0, matColInd = row;
00464 subMatColInd <= subMatRowInd;
00465 ++subMatColInd, ++matColInd) {
00466 vectorPtr[matRowInd][matColInd] += subMatrix[subMatRowInd][subMatColInd];
00467 }
00468 }
00469 }
00470
00471
00472
00473
00477 inline void
00478 QCSymMatrix::power (QCMemory& QCRestrict memoryInst,
00479 const QCSymMatrix& matrixA,
00480 const QCFloat& alpha)
00481 {
00482 power (memoryInst, matrixA, dim, alpha);
00483 }
00484
00485
00486
00487
00488 #endif // QC_SYMMATRIX_H_