00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #include "QCSymMatrix.hpp"
00021 #include "QCErrorManager.hpp"
00022 #include "QCFortran.hpp"
00023 #include "QCLapack.hpp"
00024 #include "QCBlas.hpp"
00025
00026
00027
00031 QCSymMatrix::QCSymMatrix (void) :
00032 QCMatrix (),
00033 computedSystemHasChanged(true) {
00034 }
00035
00036
00037
00042 QCSymMatrix::QCSymMatrix (const int& dimMatrix) :
00049 QCMatrix (),
00050 computedSystemHasChanged(true)
00051 {
00052 setDimAndAllocate(dimMatrix);
00053 }
00054
00055
00056
00060 QCSymMatrix::~QCSymMatrix (void) {}
00061
00062
00063
00064
00068 void
00069 QCSymMatrix::setDim (const int& dimMatrix) {
00070
00071 if (dimMatrix > 0) {
00072 dim = dimMatrix;
00073 nbElems = (dim * (dim+1) )/ 2;
00074 }
00075 }
00076
00077
00078
00079
00083 bool
00084 QCSymMatrix::setDimAndAllocate (const int& dimMatrix) {
00085 int i;
00086
00087 if (dimMatrix > 0) {
00088 dim = dimMatrix;
00089 nbElems = (dim * (dim+1) )/ 2;
00090
00091
00092 elems = new QCFloat [nbElems];
00093
00094
00095 vectorPtr = new QCFloat* [dim];
00096
00097
00098
00099
00100
00101 for (i = 0; i < dim; i++) {
00102 vectorPtr[i] = &elems[(i*(i+1) )/2];
00103 }
00104 initClean();
00105 }
00106 return (elems != NULL && vectorPtr != NULL);
00107 }
00108
00109
00110
00111
00115 bool
00116 QCSymMatrix::setDimAndElems (const int& dimMatrix,
00117 QCFloat * extElems) {
00118 int i;
00119
00120 assert (dimMatrix);
00121
00122 if (dimMatrix > 0) {
00123 dim = dimMatrix;
00124 nbElems = (dim * (dim+1) )/ 2;
00125
00126
00127 elems = extElems;
00128
00129
00130 vectorPtr = new QCFloat* [dim];
00131
00132
00133
00134
00135
00136 for (i = 0; i < dim; i++) {
00137 vectorPtr[i] = &elems[(i*(i+1) )/2];
00138 }
00139
00140 }
00141 return (elems != NULL && vectorPtr != NULL);
00142 }
00143
00144
00145
00146
00150 void
00151 QCSymMatrix::unsetElems (void) {
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162 elems = NULL;
00163 dim = 0;
00164 nbElems = 0;
00165 }
00166
00167
00168
00169
00173 bool
00174 QCSymMatrix::takeTmpMem (QCMemory& QCRestrict memoryInst) {
00175 bool retVal;
00176 int i;
00177
00178
00179 elems = memoryInst.takeTmpMem<QCFloat>(nbElems);
00180
00181
00182 vectorPtr = memoryInst.takeTmpMem<QCFloat*>(dim);
00183 for (i = 0; i < dim; i++) {
00184 vectorPtr[i] = &elems[(i*(i+1) )/2];
00185 }
00186 initClean();
00187
00188 retVal = (elems != NULL && vectorPtr != NULL);
00189
00190
00191 allocateByNew = !retVal ;
00192
00193 return retVal;
00194 }
00195
00196
00197
00198
00202 void
00203 QCSymMatrix::print (char* matrixName) const {
00204 int i, j;
00205
00206 cout << endl << " " << matrixName << endl;
00207
00208 for (i = 0; i < dim; i++) {
00209 for (j = 0; j <=i; j++) {
00210 cout << "elems[" << i << "]["<< j << "] = "<< vectorPtr[i][j] << endl;
00211 }
00212 }
00213 }
00214
00215 void
00216 QCSymMatrix::print (const std::string & matrixName, std::ofstream &out) {
00217
00218 out << endl << matrixName << endl << endl;
00219 for (int i = 0; i < dim; i++) {
00220 for (int j = 0; j <= i; j++) {
00221 out << "elems[" << i << "]["<< j << "] = "<< vectorPtr[i][j] << endl;
00222 }
00223 }
00224 }
00225
00226
00227 void
00228 QCSymMatrix::print (char* matrixName, const int& dim1) const {
00229 int i, j;
00230
00231 cout << endl << " " << matrixName << endl;
00232 for (i = 0; i < dim1; i++) {
00233 for (j = 0; j <=i; j++) {
00234 cout << "elems[" << i << "]["<< j << "] = "<< vectorPtr[i][j] << endl;
00235 }
00236 }
00237 }
00238
00239
00240
00245 void
00246 QCSymMatrix::printInFile (const char* preFileName, bool date) {
00247 FILE* streamOfFile = NULL;
00248 int i, j;
00249 string absoluteFileName = QCCommon::outdir + PATH_SEPARATOR;
00250 if( date ){
00251 absoluteFileName += string(preFileName) + QCCommon::getTimeStr() +
00252 OUTPUT_SUFFIX;
00253
00254 } else {
00255 absoluteFileName += string(preFileName) + OUTPUT_SUFFIX;
00256 }
00257
00258
00259
00260 if ( (streamOfFile = fopen(absoluteFileName.c_str(), "a") ) ) {
00261 fprintf(streamOfFile, "dim = %5d \n", dim);
00262 for (i = 0; i < dim; i++) {
00263 for (j = 0; j <= i; j++) {
00264 fprintf(streamOfFile, "%5d %5d %20.20e\n", i, j, vectorPtr[i][j]);
00265 }
00266 }
00267 fclose(streamOfFile);
00268 streamOfFile = NULL;
00269 }
00270 }
00271
00272
00273
00274
00280 void
00281 QCSymMatrix::printFullInFile (const char* preFileName) {
00282 FILE* streamOfFile = NULL;
00283 int i, j;
00284
00285 string absoluteFileName = string(preFileName) + QCCommon::getTimeStr() +
00286 OUTPUT_SUFFIX;
00287
00288
00289 if ( (streamOfFile = fopen(absoluteFileName.c_str(), "a") ) ) {
00290 for (i = 0; i < dim; i++) {
00291 for (j = 0; j < dim; j++) {
00292 if (j <= i) {
00293 fprintf(streamOfFile, "%5d %5d %20.12e\n", i, j, vectorPtr[i][j]);
00294
00295 } else {
00296 fprintf(streamOfFile, "%5d %5d %20.12e\n", i, j, vectorPtr[j][i]);
00297 }
00298 }
00299 }
00300 fclose(streamOfFile);
00301 streamOfFile = NULL;
00302 }
00303 }
00304
00305
00306
00307
00311 bool
00312 QCSymMatrix::diagonalize (QCMemory& QCRestrict memoryInst,
00313 QCDiagoAlgorithm& diagoAlgorithm,
00314 const QCFloat& toleranceFactor,
00315 QCFloat* QCRestrict eigenVal,
00316 #ifdef USE_LAPACK
00317 int* QCRestrict ,
00318 #else
00319 int* QCRestrict degener,
00320 #endif
00321 QCFloat* QCRestrict eigenVectCt,
00322 const int& nbIte)
00323 {
00324
00325
00326 int i;
00327 int intRetVal = 0;
00328 bool retVal = true;
00329 QCFloat* QCRestrict workingMatrix;
00330 QCFloat* QCRestrict workingArray1;
00331 #ifdef USE_LAPACK
00332 int* QCRestrict workingArray2;
00333 int* QCRestrict workingArray3;
00334 QCFloat dummy1=0.0, dummy2=0.0;
00335 int dummy3=0, dummy4=0;
00336 int nbEigenval=0;
00345 static int nb;
00346 static int workingDim;
00347 static int workingArray1Size, workingArray2Size;
00348
00349 updateDiagonalizeLocalVar(diagoAlgorithm, dim, nb, workingDim,
00350 workingArray1Size, workingArray2Size);
00351
00352
00353 if (diagoAlgorithm == QC_CLASSIC_DIAGO) {
00354 workingMatrix = memoryInst.takeTmpMem<QCFloat>(dim*dim);
00355 workingArray1 = memoryInst.takeTmpMem<QCFloat>(workingArray1Size);
00356 workingArray2 = memoryInst.takeTmpMem<int>(5*dim);
00357 workingArray3 = memoryInst.takeTmpMem<int>(dim);
00358
00364 for (i = 0; i < dim; ++i) {
00365
00366
00367
00368 memcpy(&(workingMatrix[i*dim]), vectorPtr[i], (i+1)*sizeof(QCFloat) );
00369 }
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380 QCLapack(dsyevx)("V", "A", "U",
00381 dim, workingMatrix, dim,
00382 &dummy1, &dummy2,
00383 &dummy3, &dummy4,
00384 toleranceFactor,
00385 nbEigenval, eigenVal,
00386 eigenVectCt, dim,
00387 workingArray1, workingArray1Size,
00388 workingArray2,
00389 workingArray3,
00390 intRetVal);
00391
00392
00393 retVal = retVal && memoryInst.giveBackTmpMem<int>(workingArray3);
00394 retVal = retVal && memoryInst.giveBackTmpMem<int>(workingArray2);
00395 retVal = retVal && memoryInst.giveBackTmpMem<QCFloat>(workingArray1);
00396 retVal = retVal && memoryInst.giveBackTmpMem<QCFloat>(workingMatrix);
00397
00398 } else {
00399
00400
00401 workingArray1 = memoryInst.takeTmpMem<QCFloat>(workingArray1Size);
00402 workingArray2 = memoryInst.takeTmpMem<int>(workingArray2Size);
00403 QCLapack(dsyevd)("V", "U",
00404 dim, eigenVectCt, dim,
00405 eigenVal,
00406 workingArray1, workingArray1Size,
00407 workingArray2, workingArray2Size,
00408 intRetVal);
00409
00410 retVal = retVal && memoryInst.giveBackTmpMem<int>(workingArray2);
00411 retVal = retVal && memoryInst.giveBackTmpMem<QCFloat>(workingArray1);
00412 }
00413
00414 #endif // USE_LAPACK
00415
00416 if (intRetVal != 0) {
00417
00418
00419
00420 QCErrorManager errorManager;
00421
00422
00423 errorManager.writeMsg(QCErrorManager::ERROR,
00424 const_cast<char*>(QCErrorManager::DIAGO_FAILS),
00425 nbIte,
00426 intRetVal);
00427 }
00428
00429 retVal = retVal && (intRetVal == 0)? true : false;
00430
00431 return retVal;
00432 }
00433
00434
00435
00439 QCFloat
00440 QCSymMatrix::traceProduct (const QCSymMatrix& otherMatrix) const {
00441
00442 QCFloat trace;
00443 const QCFloat* QCRestrict otherElems = otherMatrix.getElems();
00444 int i;
00445
00446 #ifdef DEBUG_QC
00447 if (dim == otherMatrix.getDim() ||
00448 nbElems == otherMatrix.getNbElems() ) {
00449 #endif //DEBUG_QC
00450 trace = QC_ZERO;
00451 QCFloat trace1 = QC_ZERO;
00452 QCFloat trace2 = QC_ZERO;
00453 QCFloat trace3 = QC_ZERO;
00454 QCFloat trace4 = QC_ZERO;
00455 int dimMultOf4 = dim - dim%4;
00456
00457 trace1 = QCCBlas(ddot)(nbElems, elems, 1,
00458 const_cast<QCFloat*>(otherElems), 1);
00459 trace1 *= QC_TWO;
00460
00466 for (i = 0; i < dimMultOf4; i+=4) {
00467 trace1 -= vectorPtr[i][i] * otherMatrix[i][i];
00468 trace2 -= vectorPtr[i+1][i+1] * otherMatrix[i+1][i+1];
00469 trace3 -= vectorPtr[i+2][i+2] * otherMatrix[i+2][i+2];
00470 trace4 -= vectorPtr[i+3][i+3] * otherMatrix[i+3][i+3];
00471 }
00472 for (i = dimMultOf4; i < dim; ++i) {
00473 trace1 -= vectorPtr[i][i] * otherMatrix[i][i];
00474 }
00475 trace = trace1 + trace2 + trace3 + trace4;
00476 #ifdef DEBUG_QC
00477
00478 } else {
00479
00480
00481 cout << "Pb in File " << __FILE__ << " at line " << __LINE__ << endl;
00482 abort();
00483 }
00484 #endif //DEBUG_QC
00485
00486 return trace;
00487 }
00488
00489
00490
00491
00496 void
00497 QCSymMatrix::power (QCMemory& QCRestrict memoryInst,
00498 const QCSymMatrix& matrixA,
00499 const int& dimToMultiplyOfA,
00500 const QCFloat& alpha)
00501 {
00502 QCMatrix workingMatrixA;
00503 QCMatrix workingMatrixThis;
00504 workingMatrixA.setDim(matrixA.getDim() );
00505 workingMatrixThis.setDim(dim);
00506 workingMatrixA.takeTmpMem(memoryInst);
00507 workingMatrixThis.takeTmpMem(memoryInst);
00508 workingMatrixThis.initClean();
00509
00510
00511
00512 matrixA.copyInAndFill(workingMatrixA);
00513
00514 QCCBlas(dsyrk)(
00515 #ifdef USE_OPTIMIZED_CBLAS
00516 CblasRowMajor,
00517 CblasLower, CblasTrans,
00518 #else //USE_OPTIMIZED_CBLAS
00519 CblasUpper, CblasNoTrans,
00520 #endif //USE_OPTIMIZED_CBLAS
00521 dimToMultiplyOfA, dimToMultiplyOfA,
00522 alpha, workingMatrixA.getElems(), dim,
00523 QC_ZERO, workingMatrixThis.getElems(), dim);
00524
00525
00526 this->copy(workingMatrixThis);
00527
00528 workingMatrixThis.giveBackTmpMem(memoryInst);
00529 workingMatrixA.giveBackTmpMem(memoryInst);
00530 }
00531
00532
00533
00534
00538 template <class TPManager>
00539 void
00540 QCSymMatrix::writeLikeBlocks(TPManager& QCRestrict managerInst) const {
00541
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579
00580
00581
00582
00583
00584
00585
00586
00587 }
00588
00589
00590
00591
00592
00593 #ifdef USE_LAPACK
00594
00597 void
00598 QCSymMatrix::updateDiagonalizeLocalVar (QCDiagoAlgorithm& diagoAlgorithm,
00599 const int& dimOfThePartToDiago,
00600 int& nb,
00601 int& workingDim,
00602 int& workingArray1Size,
00603 int& workingArray2Size) {
00604 int ispec = 1;
00605
00606
00607
00608 computedSystemHasChanged = false;
00609
00610
00611
00612 nb = QCFortran(ilaenv)
00613 (ispec, "DSYTRD", "VAU", dimOfThePartToDiago,
00614 dimOfThePartToDiago, dimOfThePartToDiago,
00615 dimOfThePartToDiago, 6, 3);
00616
00617
00618 workingDim = ( (nb + 3) > 8) ? (nb + 3) : 8;
00619
00620 #ifndef __QC_cxx__
00621 workingArray1Size = (diagoAlgorithm == QC_CLASSIC_DIAGO) ?
00622 workingDim * dimOfThePartToDiago :
00623 1 + 6 * dimOfThePartToDiago + 2 * QCPow<2>(dimOfThePartToDiago);
00624 #else //__QC_cxx__
00625
00634 workingArray1Size = (diagoAlgorithm == QC_CLASSIC_DIAGO) ?
00635 workingDim * dimOfThePartToDiago :
00636 1 + 5 * dimOfThePartToDiago +
00637 (2 * (static_cast<int>(log(dimOfThePartToDiago)/log(QCTwo) ) + 1)
00638 * dimOfThePartToDiago)
00639 + 3 * QCPow<2>(dimOfThePartToDiago);
00640 #endif //__QC_cxx__
00641 workingArray2Size = (diagoAlgorithm == QC_CLASSIC_DIAGO) ?
00642 0 : 3 + 5 * dimOfThePartToDiago;
00643 }
00644 #endif //USE_LAPACK