QCMatrix.cpp

Go to the documentation of this file.
00001 //*****************************************************************************//
00002 //                                                                             //
00003 //   Copyright (c) 2001                                                        //
00004 //      INRIA                                                                  //
00005 //      54600 VILLERS LES NANCY                                                //
00006 //      France                                                                 //
00007 //                                                                             //
00008 //*****************************************************************************//
00009 //                                                                             //
00010 //               *** NOTICE OF PROPRIETARY INFORMATION ***                     //
00011 //                                                                             //
00012 // The information contained in this file is considered proprietary and the    //
00013 // exclusive property of  INRIA. This information may not be disclosed,        //
00014 // duplicated or used, in whole or in part, for  any purpose  whatsoever       //
00015 // without express written authorization from INRIA                            //
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   // ok, on vire explicitement les parties dynamiques.
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     //sinon l appli est out.
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   // On copie la matrice symetriques A dans une matrice pleine 
00136   // A totalement remplie.
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   // On copie les matrices symetriques A et B dans une 
00172   // matrice pleine A totalement remplie.
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     // on alloue le grand tableau
00207     elems = new QCFloat [dim*dim];
00208 
00209     // on alloue le tableau de pointeur sur les lignes ou les colonnes.
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   // on alloue le grand tableau
00229   elems = memoryInst.takeTmpMem<QCFloat>(nbElems);
00230   
00231   // on alloue le tableau de pointeur sur les lignes ou les colonnes.
00232   vectorPtr = memoryInst.takeTmpMem<QCFloat*>(dim);
00233   for (i = 0; i < dim; i++) {
00234     vectorPtr[i] = &elems[i*dim];
00235   }//for i
00236   initClean();
00237 
00238   retVal = (elems != NULL && vectorPtr != NULL);
00239 
00240   //Si retVal est vrai alors allocateByNew est faux.
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   // const_cast sinon warning a la compil a cause des 2 QCRestrict 
00257   // a la declaration de vectorPtr.
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     // sinon l appli est out.
00265     // std::cout.setf(ios::boolalpha) ;
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   // On remet allocateByNew à sa valeur par defaut.
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   //ouverture du fichier de
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   }//for i
00404   
00405   FILE* streamOfFile = NULL;
00406   string absoluteFileName = string(preFileName) + OUTPUT_SUFFIX;
00407 
00408   //ouverture du fichier de
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   //ouverture du fichier de
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     //sinon l appli est out.
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   // QCWriter&  theWriter = managerInst.getTheWriter();
00549   
00550 //   //On recupere avec des restricts pour le calcul les
00551 //   //structures qui contiennent les donnees a manipuler.
00552 //   CQCAtom* QCRestrict        atoms = managerInst.getAtoms();
00553 //   CQCAtoms& QCRestrict       system = managerInst.getAtom();
00554 //   //
00555 //   //les atomes courants.
00556 //   const CQCAtom* QCRestrict atomA;
00557 //   //
00558 //   //Le tableau de parametres.
00559 //   const typename TPManager::TModel::TParam*  QCRestrict parameters = 
00560 //     managerInst.getParameters();
00561 //   //et leurs parametres.
00562 //   const typename TPManager::TModel::TParam* QCRestrict parameterA;
00563 //   //
00564 //   int nbAtoms = system.getNbAtoms();
00565 //   int nbOAOfSystem = system.getNbOAOfSystem();
00566 //    //
00567 //   //Numero des orbitales.
00568 //   int        firstOAOfA, firstOAOfSucA;
00569 //   int i, mu, nu;
00570 //   //
00571 //   //La boucle sur les atomes.
00572 //   for (i = 0; i < nbAtoms; ++i) {
00573 //     atomA = &atoms[i];
00574 //     parameterA  = &parameters[atomA->getParamInd() ];
00575 //     firstOAOfA = atomA->getStartOrb();
00576 //     firstOAOfSucA = firstOAOfA + parameterA->getNbOA();
00577 //     theWriter.writeSubMatrixHeader(parameterA->getIdentificator(), system.getNumber(i) );
00578 //     //
00579 //     //Attention, on est a OA cte le long des colonnes de la matrice, pas le long
00580 //     // des lignes.
00581 //     for (nu = firstOAOfA; nu < firstOAOfSucA; ++nu) {
00582 //       for (mu = 0; mu < nbOAOfSystem; ++mu) {
00583 //      theWriter.writeQCFloat(vectorPtr[mu][nu]);
00584 //       }//for nu
00585 //       theWriter.writeEndl();
00586 //     }//for mu
00587 //   }//for i
00588 }

Generated on Sat Jan 28 21:07:28 2006 for QC++ by  doxygen 1.4.4