00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #include <ios>
00019
00020 #include "QCMatrix.hpp"
00021 #include "QCSymMatrix.hpp"
00022
00023
00024
00028 QCMatrix::QCMatrix (void) :
00029 vectorPtr(NULL),
00030 elems(NULL),
00031 dim(0),
00032 nbElems(0),
00033 allocateByNew(true) {
00034 }
00035
00036
00037
00038
00042 QCMatrix::QCMatrix (const int& dimMatrix) :
00043 vectorPtr(NULL),
00044 elems(NULL),
00045 dim(dimMatrix),
00046 nbElems(0),
00047 allocateByNew(true)
00048 {
00049 setDimAndAllocate(dim);
00050 }
00051
00052
00053
00054
00058 QCMatrix::~QCMatrix(void) {
00059
00060
00061 if (allocateByNew ) {
00062 if (elems ) {
00063 delete [] elems;
00064 elems = NULL;
00065 }
00066 if (vectorPtr ) {
00067 delete [] vectorPtr;
00068 vectorPtr = NULL;
00069 }
00070 }
00071
00072 #if defined( DEBUG_QC_MEMORY ) || defined( DEBUG_QC )
00073 if( (elems || vectorPtr) && !allocateByNew) {
00074
00075
00076 cout << "QCMatrix::~QCMatrix : " << endl
00077 << " vectorPtr : " << (QCFloat*) vectorPtr << endl
00078 << " elems : " << elems << endl
00079 << " allocateByNew : " << allocateByNew << " true : "
00080 << true << endl;
00081
00082 cout << "Free memory by new rather by givebackMemory method"
00083 << __FILE__
00084 << " at line "
00085 << __LINE__ << endl;
00086 abort();
00087 }
00088 #endif //DEBUG_QC
00089 }
00090
00091
00092
00093
00097 void
00098 QCMatrix::setDim (const int& dimMatrix) {
00099 if (dimMatrix > 0) {
00100 dim = dimMatrix;
00101 nbElems = dim * dim;
00102 }
00103 }
00104
00105
00106
00107
00111 void
00112 QCMatrix::setVectorPtr(const int indexOfTheVectorPtr,
00113 const int indexInTheElemsArray) {
00114 vectorPtr[indexOfTheVectorPtr] = elems + indexInTheElemsArray;
00115 }
00116
00117
00118
00119
00124 void
00125 QCMatrix::mult (QCMemory& QCRestrict memoryInst,
00126 const QCSymMatrix& matrixA,
00127 const QCMatrix& matrixB,
00128 const QCFloat& alpha,
00129 const QCFloat& beta)
00130 {
00131 QCMatrix workingMatrixA;
00132 workingMatrixA.setDim(matrixA.getDim() );
00133 workingMatrixA.takeTmpMem(memoryInst);
00134
00135
00136
00137 matrixA.copyInAndFill(workingMatrixA);
00138
00139 QCCBlas(dgemm)(
00140 #ifdef USE_OPTIMIZED_CBLAS
00141 CblasRowMajor,
00142 CblasNoTrans, CblasTrans,
00143 #else //USE_OPTIMIZED_CBLAS
00144 CblasTrans, CblasNoTrans,
00145 #endif //USE_OPTIMIZED_CBLAS
00146 dim, dim, dim, alpha, workingMatrixA.elems, dim, matrixB.elems,
00147 dim, beta, elems, dim);
00148
00149 workingMatrixA.giveBackTmpMem(memoryInst);
00150 }
00151
00152
00153
00154
00159 void
00160 QCMatrix::mult (QCMemory& QCRestrict memoryInst,
00161 const QCSymMatrix& matrixA, const QCSymMatrix& matrixB,
00162 const QCFloat& alpha, const QCFloat& beta) {
00163
00164 QCMatrix workingMatrixA;
00165 QCMatrix workingMatrixB;
00166 workingMatrixA.setDim(matrixA.getDim() );
00167 workingMatrixB.setDim(matrixB.getDim() );
00168 workingMatrixA.takeTmpMem(memoryInst);
00169 workingMatrixB.takeTmpMem(memoryInst);
00170
00171
00172
00173 matrixA.copyInAndFill(workingMatrixA);
00174 matrixB.copyInAndFill(workingMatrixB);
00175
00176 QCCBlas(dgemm)(
00177 #ifdef USE_OPTIMIZED_CBLAS
00178 CblasRowMajor,
00179 CblasNoTrans, CblasTrans,
00180 #else //USE_OPTIMIZED_CBLAS
00181 CblasTrans, CblasNoTrans,
00182 #endif //USE_OPTIMIZED_CBLAS
00183 dim, dim, dim, alpha,
00184 workingMatrixA.elems,
00185 dim, workingMatrixB.elems,
00186 dim, beta, elems, dim);
00187 workingMatrixB.giveBackTmpMem(memoryInst);
00188 workingMatrixA.giveBackTmpMem(memoryInst);
00189 }
00190
00191
00192
00193
00194
00198 bool
00199 QCMatrix::setDimAndAllocate (const int& dimMatrix) {
00200 int i;
00201
00202 if (dimMatrix > 0) {
00203 dim = dimMatrix;
00204 nbElems = dim * dim;
00205
00206
00207 elems = new QCFloat [dim*dim];
00208
00209
00210 vectorPtr = new QCFloat* [dim];
00211 for (i = 0; i < dim; i++) {
00212 vectorPtr[i] = &elems[i*dim];
00213 }
00214 initClean();
00215 }
00216 return (elems != NULL && vectorPtr != NULL);
00217 }
00218
00222 bool
00223 QCMatrix::takeTmpMem (QCMemory& QCRestrict memoryInst) {
00224
00225 bool retVal;
00226 int i;
00227
00228
00229 elems = memoryInst.takeTmpMem<QCFloat>(nbElems);
00230
00231
00232 vectorPtr = memoryInst.takeTmpMem<QCFloat*>(dim);
00233 for (i = 0; i < dim; i++) {
00234 vectorPtr[i] = &elems[i*dim];
00235 }
00236 initClean();
00237
00238 retVal = (elems != NULL && vectorPtr != NULL);
00239
00240
00241 allocateByNew = !retVal;
00242
00243 return retVal;
00244 }
00245
00246
00247
00248
00252 bool
00253 QCMatrix::giveBackTmpMem (QCMemory& QCRestrict memoryInst) {
00254 bool retVal = true;
00255
00256
00257
00258 if (! allocateByNew) {
00259 retVal = retVal && memoryInst.giveBackTmpMem<QCFloat* QCRestrict>(vectorPtr);
00260 retVal = retVal && memoryInst.giveBackTmpMem<QCFloat>(elems);
00261 }
00262 #if defined( DEBUG_QC_MEMORY ) || defined( DEBUG_QC )
00263 else {
00264
00265
00266 cout << "QCMatrix::giveBackTmpMem : " << endl
00267 << " vectorPtr : " << (QCFloat* )vectorPtr << endl
00268 << " elems : " << elems << endl
00269 << " allocateByNew : " << allocateByNew
00270 << " true : " << true << endl;
00271
00272 cout << "Free memory by givebackMemory method rather by ~CQCMatrix." << endl
00273 << " " << __FILE__ << " at line " << __LINE__ << endl;
00274 abort();
00275 }
00276 #endif
00277
00278
00279 allocateByNew = true;
00280
00281 return retVal;
00282 }
00286 void
00287 QCMatrix::initClean (void) {
00288 memset(elems, 0, nbElems*sizeof(QCFloat) );
00289 }
00293 void
00294 QCMatrix::print (char* matrixName) const {
00295 int i, j;
00296
00297 cout << endl << matrixName << endl << endl;
00298 for (i = 0; i < dim; i++) {
00299 for (j = 0; j < dim; j++) {
00300 cout << "elems[" << i << "]["<< j << "] = "<< vectorPtr[i][j] << endl;
00301 }
00302 }
00303 }
00304
00305 void
00306 QCMatrix::print (const std::string & matrixName, std::ofstream &out) const {
00307 int i, j;
00308
00309 out << endl << matrixName << endl << endl;
00310 for (i = 0; i < dim; i++) {
00311 for (j = 0; j < dim; j++) {
00312 out << "elems[" << i << "]["<< j << "] = "<< vectorPtr[i][j] << endl;
00313 }
00314 }
00315 }
00316
00317
00318 void
00319 QCMatrix::print (char* matrixName, const int& dim1, const int& dim2) const {
00320 int i, j;
00321
00322 cout << endl << " " << matrixName << endl;
00323 for (i = 0; i < dim1; i++) {
00324 for (j = 0; j < dim2; j++) {
00325 }
00326 }
00327 }
00328
00329
00330
00331
00336 void
00337 QCMatrix::printInFile (const char* preFileName, bool date ) {
00338 FILE* streamOfFile = NULL;
00339 int i, j;
00340 string absoluteFileName = QCCommon::outdir + PATH_SEPARATOR;
00341 if( date ){
00342 absoluteFileName += string(preFileName) + QCCommon::getTimeStr() +
00343 OUTPUT_SUFFIX;
00344
00345 } else {
00346 absoluteFileName += string(preFileName) + OUTPUT_SUFFIX;
00347 }
00348
00349
00350
00351 if ( (streamOfFile = fopen(absoluteFileName.c_str(), "a") ) ) {
00352 fprintf(streamOfFile, "dim = %5d \n ", dim);
00353 for (i = 0; i < dim; i++) {
00354 for (j = 0; j < dim; j++) {
00355 fprintf(streamOfFile, "%5d %5d %20.20e\n", i, j, vectorPtr[i][j]);
00356 }
00357 }
00358 fclose(streamOfFile);
00359 streamOfFile = NULL;
00360 }
00361 }
00366 void
00367 QCMatrix::sampleDistanceInFile (const char* preFileName) {
00368 int i;
00369 int nbInf4 = 0;
00370 int nbInf6 = 0;
00371 int nbInf8 = 0;
00372 int nbInf10 = 0;
00373 int nbInf12 = 0;
00374 int nbInf14 = 0;
00375 int nbInf16 = 0;
00376 int nbSup16 = 0;
00377
00378 for (i =0; i < nbElems; ++i) {
00379 if (QCAbs(elems[i]) > 16) {
00380 nbSup16++;
00381 }
00382 else if (QCAbs(elems[i]) > 14) {
00383 nbInf16++;
00384 }
00385 else if (QCAbs(elems[i]) > 12) {
00386 nbInf14++;
00387 }
00388 else if (QCAbs(elems[i]) > 10) {
00389 nbInf12++;
00390 }
00391 else if (QCAbs(elems[i]) > 8) {
00392 nbInf10++;
00393 }
00394 else if (QCAbs(elems[i]) > 6) {
00395 nbInf8++;
00396 }
00397 else if (QCAbs(elems[i]) > 4) {
00398 nbInf6++;
00399 }
00400 else {
00401 nbInf4++;
00402 }
00403 }
00404
00405 FILE* streamOfFile = NULL;
00406 string absoluteFileName = string(preFileName) + OUTPUT_SUFFIX;
00407
00408
00409 if ( (streamOfFile = fopen(absoluteFileName.c_str(), "a") ) ) {
00410 fprintf(streamOfFile, "\nNEW ITERATION\n");
00411 fprintf(streamOfFile, " %8d |elems| < 4\n", nbInf4);
00412 fprintf(streamOfFile, "4 < %8d |elems| < 6\n", nbInf6);
00413 fprintf(streamOfFile, "6 < %8d |elems| < 8\n", nbInf8);
00414 fprintf(streamOfFile, "8 < %8d |elems| < 10\n", nbInf10);
00415 fprintf(streamOfFile, "10 < %8d |elems| < 12\n", nbInf12);
00416 fprintf(streamOfFile, "12 < %8d |elems| < 14\n", nbInf14);
00417 fprintf(streamOfFile, "14 < %8d |elems| < 16\n", nbInf16);
00418 fprintf(streamOfFile, "16 < %8d |elems|\n", nbSup16);
00419
00420 fclose(streamOfFile);
00421 streamOfFile = NULL;
00422 }
00423 }
00428 void
00429 QCMatrix::sampleInFile (const char* preFileName) {
00430 int i;
00431 int nbSup3 = 0;
00432 int nbInf3 = 0;
00433 int nbInf4 = 0;
00434 int nbInf5 = 0;
00435 int nbInf6 = 0;
00436 int nbInf8 = 0;
00437 int nbInf10 = 0;
00438
00439 for (i =0; i < nbElems; ++i) {
00440 if (QCAbs(elems[i]) < 1e-10) {
00441 nbInf10++;
00442 }
00443 else if (QCAbs(elems[i]) < 1e-8) {
00444 nbInf8++;
00445 }
00446 else if (QCAbs(elems[i]) < 1e-6) {
00447 nbInf6++;
00448 }
00449 else if (QCAbs(elems[i]) < 1e-5) {
00450 nbInf5++;
00451 }
00452 else if (QCAbs(elems[i]) < 1e-4) {
00453 nbInf4++;
00454 }
00455 else if (QCAbs(elems[i]) < 1e-3) {
00456 nbInf3++;
00457 }
00458 else {
00459 nbSup3++;
00460 }
00461 }
00462
00463 FILE* streamOfFile = NULL;
00464 string absoluteFileName = string(preFileName) + OUTPUT_SUFFIX;
00465
00466
00467 if ( (streamOfFile = fopen(absoluteFileName.c_str(), "a") ) ) {
00468 fprintf(streamOfFile, "\nNEW ITERATION\n");
00469 fprintf(streamOfFile, "10-3 < %8d |elems| \n", nbSup3);
00470 fprintf(streamOfFile, "10-4 < %8d |elems| < 10-3\n", nbInf3);
00471 fprintf(streamOfFile, "10-5 < %8d |elems| < 10-4\n", nbInf4);
00472 fprintf(streamOfFile, "10-6 < %8d |elems| < 10-5\n", nbInf5);
00473 fprintf(streamOfFile, "10-8 < %8d |elems| < 10-6\n", nbInf6);
00474 fprintf(streamOfFile, "10-10 < %8d |elems| < 10-8\n", nbInf8);
00475 fprintf(streamOfFile, " %8d |elems| < 10-10\n", nbInf10);
00476
00477 fclose(streamOfFile);
00478 streamOfFile = NULL;
00479 }
00480 }
00485 QCFloat
00486 QCMatrix::deltaElemMax (const QCMatrix& otherMatrix) const {
00487 QCFloat deltaMax;
00488 int i;
00489 #ifdef DEBUG_QC
00490 if (dim == otherMatrix.getDim() && nbElems == otherMatrix.nbElems) {
00491 #endif //DEBUG_QC
00492 deltaMax = QC_ZERO;
00493 for (i = 0; i < nbElems; ++i) {
00494 if (QCAbs(elems[i] - otherMatrix.elems[i]) > deltaMax) {
00495 deltaMax = QCAbs(elems[i] - otherMatrix.elems[i]);
00496 }
00497 }
00498 #ifdef DEBUG_QC
00499
00500 } else {
00501
00502
00503 cout << "Pb in File " << __FILE__ << " at line " << __LINE__ << endl;
00504 abort();
00505 }
00506 #endif //DEBUG_QC
00507
00508 return deltaMax;
00509 }
00514 QCFloat
00515 QCMatrix::elemMax () const {
00516 QCFloat elemMax;
00517 int i;
00518 elemMax = QC_ZERO;
00519 for (i = 0; i < nbElems; ++i) {
00520 if (QCAbs(elems[i]) > elemMax) {
00521 elemMax = QCAbs(elems[i]);
00522 }
00523 }
00524 return elemMax;
00525 }
00529 QCFloat
00530 QCMatrix::trace () const {
00531 QCFloat trace;
00532 int i;
00533
00534 trace = QC_ZERO;
00535 for (i = 0; i < dim; i++) {
00536 trace += vectorPtr[i][i];
00537 }
00538
00539 return trace ;
00540 }
00545 template <class TPManager>
00546 void
00547 QCMatrix::writeLikeBands (TPManager& QCRestrict managerInst) const {
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 }