QCSymMatrix.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 
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     // on alloue le grand tableau.
00092     elems = new QCFloat [nbElems];
00093 
00094     //on alloue le tableau de pointeur sur les lignes ou les colonnes.
00095     vectorPtr = new QCFloat* [dim];
00096 
00097     //Attention, la derniere ligne, est la plus longue.
00098     //La premiere est la plus courte, elle contient seulement
00099     // l element diagonal.
00100     //Ensuite, on va crescendo.
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     // on positionne le grand tableau.
00127     elems = extElems;
00128 
00129     // on alloue le tableau de pointeur sur les lignes ou les colonnes.
00130     vectorPtr = new QCFloat* [dim];
00131 
00132     // Attention, la derniere ligne, est la plus longue.
00133     // La premiere est la plus courte, elle contient seulement
00134     // l element diagonal.
00135     // Ensuite, on va crescendo.
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   //   for (int i = 0; i < dim; i++) {
00155   //     vectorPtr[i] = NULL;
00156   //   }
00157 
00158   // on libere le tableau de pointeur sur les lignes ou les colonnes.
00159   //  delete [] vectorPtr;
00160 
00161   // on deconecte le grand tableau.
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   //on alloue le grand tableau
00179   elems = memoryInst.takeTmpMem<QCFloat>(nbElems);
00180   
00181   //on alloue le tableau de pointeur sur les lignes ou les colonnes.
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   //Si retVal est vrai alors allocateByNew est faux.
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   //ouverture du fichier de
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   //ouverture du fichier de
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       /*degener*/,
00318 #else //USE_LAPACK
00319                           int* QCRestrict       degener,
00320 #endif //USE_LAPACK
00321                           QCFloat* QCRestrict   eigenVectCt,
00322                           const int&            nbIte) 
00323 {
00324   //  std::cout << " Diagonalisation depuis QCSymMAtrix"<<std::endl;
00325   // On declare les autres structures dont on a besoin pour la diago.
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   //if (computedSystemHasChanged) {
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       // On va crescendo en copiant car c est le lower triangle qui
00367       // est dans la matrice symetrique.
00368       memcpy(&(workingMatrix[i*dim]), vectorPtr[i], (i+1)*sizeof(QCFloat) );
00369     }
00370  //    cout << "  QCLapack(dsyevx)(V, A, U,"<< "N="<<dim 
00371 //       << "          A "<< workingMatrix[0] << " LDA " << dim 
00372 //       << "          ABSTOL " << toleranceFactor << endl
00373 //       << "          M      " << nbEigenval << endl
00374 //       << "          eigenVal        " << eigenVal[0] << endl
00375 //       << "          eigenVectCt     " << eigenVectCt  << "  LDZ "<< dim <<endl
00376 //       << "          workingArray1   " << workingArray1[0] << "   LWORK "<<  workingArray1Size <<endl
00377 //       << "          workingArray2   " << workingArray2[0] << endl
00378 //       << "          workingArray3   " << workingArray3[0] << endl
00379 //       << "          INFO       " << intRetVal << endl;
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     //    cout << "                       .... ) end" <<std::endl;
00392     // et on rend la memoire temporaire.
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 { //diagoAlgorithm == QC_DC_DIAGO
00399 
00400     // on recupere la memoire temporaire.
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     // et on rend la memoire temporaire.
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     // on cree le gestionnaire d erreur qui va envoyer
00419     // le message d erreur.
00420     QCErrorManager errorManager;   
00421 
00422     //Attention, la liste d arguments est variable.
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       //sinon l appli est out.
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   // On copie les matrices symetriques A et B dans une matrice pleine 
00511   // A totalement remplie.
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   // On copie le resultat dans une la QCSymMatrix this.
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   // CQCWriter& theWriter = managerInst.getTheWriter();
00542 //   //
00543 //   //On recupere avec des restricts pour le calcul les
00544 //   //structures qui contiennent les donnees a manipuler.
00545 //   CQCAtom* QCRestrict        atoms = managerInst.getAtoms();
00546 //   CQCAtoms& QCRestrict       system = managerInst.getAtom();
00547 //   //
00548 //   //les atomes courants.
00549 //   const CQCAtom* QCRestrict atomA;
00550 //   const CQCAtom* QCRestrict atomB;
00551 //   //
00552 //   //Le tableau de parametres.
00553 //   const typename TPManager::TModel::TParam*  QCRestrict parameters = 
00554 //     managerInst.getParameters();
00555 //   //et leurs parametres.
00556 //   const typename TPManager::TModel::TParam* QCRestrict parameterA;
00557 //   const typename TPManager::TModel::TParam* QCRestrict parameterB;
00558 //   //
00559 //   //Numero des orbitales.
00560 //   int        firstOAOfA, firstOAOfSucA, firstOAOfB, firstOAOfSucB;
00561 //   int i, j, mu, nu;
00562 //   //
00563 //   //La double boucle sur les atomes.
00564 //   //je prends  j <= i.
00565 //   for (i = 0; i < system.getNbAtoms(); ++i) {
00566 //     atomA = &atoms[i];
00567 //     parameterA  = &parameters[atomA->getParamInd() ];
00568 //     firstOAOfA = atomA->getStartOrb();
00569 //     firstOAOfSucA = firstOAOfA + parameterA->getNbOA();
00570 //     //
00571 //     //Boucle sur les B pour remplir les H_mu_nu.
00572 //     for (j = 0; j <= i; ++j) {
00573 //       atomB = &atoms[j];
00574 //       parameterB  = &parameters[atomB->getParamInd()];
00575 //       firstOAOfB = atomB->getStartOrb();
00576 //       firstOAOfSucB = firstOAOfB + parameterB->getNbOA();
00577 //       theWriter.writeSubMatrixHeader(parameterA->getIdentificator(), system.getNumber(i), 
00578 //                                    parameterB->getIdentificator(), system.getNumber(j) );
00579 //       for (mu = firstOAOfA; mu < firstOAOfSucA; ++mu) {
00580 //      for (nu = firstOAOfB; nu < firstOAOfSucB && nu <= mu; ++nu) {
00581 //        theWriter.writeQCFloat(vectorPtr[mu][nu]);
00582 //      }//for nu
00583 //      theWriter.writeEndl();
00584 //       }//for mu
00585 //     }//for j
00586 //   }//for i
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   // On met justCreated a false pour initialiser une seule fois les 
00607   // variables statiques de diagonalise.
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

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