QCDensityGtr.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 
00021 #include "QCCommon.hpp"
00022 #include "QCTable.hpp"
00023 #include "QCDensityGtr.hpp"
00024 #include "QCSymMatrix.hpp"
00025 #include "QCSystem.hpp"
00026 #include "QCMndo.hpp"
00027 #include "QCAm1.hpp"
00028 #include "QCPm3.hpp"
00029 #include "QCMndoParam.hpp"
00030 #include "QCAm1Param.hpp"
00031 #include "QCPm3Param.hpp"
00032 #include "QCGlobalSystem.hpp"
00033 #include "QCMDSystem.hpp"
00034 #include "QCDistMDSystem.hpp"
00035 #include "QCSCFAlgo.hpp"
00036 #include "QCDCAlgo.hpp"
00037 #include "QCManager.hpp"
00038 #include "QCSubDomain.hpp"
00039 
00040 #if defined (QC_TRACE_IN_FILE)  
00041 #include "QCTrace.hpp"
00042 #endif
00043 
00047 template <class TPMatrix>
00048 QCDensityGtr<TPMatrix>::QCDensityGtr (void) : 
00049   QCMatElemGtr<TPMatrix>() 
00050 {}
00051 
00052 
00053 
00054 
00058 template <class TPMatrix>
00059 QCDensityGtr<TPMatrix>::QCDensityGtr (int dimMatrix) : 
00060   QCMatElemGtr<TPMatrix>(dimMatrix) 
00061 {}
00062 
00063 
00064 
00065 
00069 template <class TPMatrix>
00070 QCDensityGtr<TPMatrix>::~QCDensityGtr (void) {
00071 }
00072 
00073 template <class TPMatrix>
00074 template <class TPManager>
00075 void            
00076 QCDensityGtr<TPMatrix>::initDiag (TPManager& manager,
00077                               typename TPManager::TSystem::QCIterator& iter) {
00078   QC_TRACE_INIT("BEGIN QCDensityGtr<TPMatrix>::initDiag")  ;
00079   // On recupere avec des restricts pour le calcul les
00080   // structures qui contiennent les donnees a manipuler
00081   QCSystem& QCRestrict currentSystem     = iter->getSystem();
00082   const QCSystem& QCRestrict rootSystem  = iter->getRootSystem();
00083 
00084   // Les parametres
00085   const typename TPManager::TParam *QCRestrict parameters = 
00086     manager.getParameters();
00087   
00088   // Le parametre courant
00089   const typename TPManager::TParam *QCRestrict parameterA;
00090 
00091   // Le premier index d'orbitale de A  dans la matrice complete du
00092   // systeme et l'index du premier atome qui lui succede.
00093   int firstAOofA, firstAOofSucA;
00094 
00095   QCFloat diagElem, chargeForDiagElem;
00096 
00097   int i, mu;
00098   
00099   chargeForDiagElem = ( static_cast<QCFloat>(rootSystem.getSystemCharge()) / 
00100                         rootSystem.getNbAtomicOrbitals() );
00101   this->matrix.initClean();
00102 
00103   for (i = 0; i < currentSystem.getNbAtoms(); ++i) {
00104     
00105     parameterA    = &parameters[currentSystem.getParamIndexAt(i)];
00106     firstAOofA    = currentSystem.getFirstAOAt(i);
00107     firstAOofSucA = firstAOofA + parameterA->getNbAO();
00108 
00109     // On remplit les termes diagonaux : P_mu_mu.
00110     // Chaque element de la diag porte : nb d e / nb d OA.
00111     diagElem = ( ( static_cast<QCFloat>(parameterA->getNbValenceElec()) /
00112                    parameterA->getNbAO() )
00113                  // Il faut aussi tenir compte de la charge.
00114                  - chargeForDiagElem );
00115     for (mu = firstAOofA; mu < firstAOofSucA; ++mu) {
00116       this->matrix[mu][mu] = diagElem;
00117     }
00118   }
00119 
00120   QC_TRACE_INIT("END   QCDensityGtr<TPMatrix>::initDiag")  ;
00121 }
00125 template <class TPMatrix>
00126 template <class TPMatrix2>
00127 void
00128 QCDensityGtr<TPMatrix>::computeElems (QCMemory&                 QCRestrict memory,
00129                                       const int                 nbElecOfSystem,
00130                                       QCMatElemGtr<TPMatrix2>& QCRestrict eigenVectC) {
00131 
00132   // La matrice des vecteurs propres.
00133   TPMatrix2& QCRestrict transEigenVectMatrix = eigenVectC.getMatrix();
00134   
00135   const int nbFilledOM = nbElecOfSystem / 2;
00136   const int dim        = this->matrix.getDim();
00137 
00138   int i;
00139   QCFloat* QCRestrict workingMatrix = 
00140     memory.template takeTmpMem<QCFloat>(dim*dim);
00141 
00148   // En fait, a cause du passage fortran, eigenVectMatrix contient
00149   // la matrice transposee de C et il faut mettre CblasTrans
00150   QCCBlas(dsyrk)(
00151 #ifdef USE_OPTIMIZED_CBLAS
00152                  CblasRowMajor,
00153                  CblasLower, CblasTrans,
00154 #else // USE_OPTIMIZED_CBLAS
00155                  CblasUpper, CblasNoTrans,
00156 #endif // USE_OPTIMIZED_CBLAS
00157                  dim, nbFilledOM, 
00158                  QC_TWO, transEigenVectMatrix.getElems(), dim, 
00159                  QC_ZERO, workingMatrix, dim);
00160 
00161   // Il faut maintenant recopier la workingMatrix dans la 
00162   // matrice densite.
00163   for (i = 0; i < dim; ++i) {
00164     memcpy(this->matrix[i], &(workingMatrix[i*dim]), (i+1) * sizeof(QCFloat) );
00165   }
00166   
00167   memory.template giveBackTmpMem<QCFloat>(workingMatrix);
00168 }
00172 template <class TPMatrix>
00173 template <class TPParam>
00174 void
00175 QCDensityGtr<TPMatrix>::mergeDensity (const TPParam * QCRestrict parameters,
00176                                       QCSubDomain&    QCRestrict workingDomain,
00177                                       QCSubDomain&    QCRestrict remoteDomain,
00178                                       QCSymMatrix&    QCRestrict remWeightedDensity,
00179                                       int                        remoteMapIdx,
00180                                       QCMatrix&       QCRestrict wDensityAB,
00181                                       QCMatrix&       QCRestrict interPAB,
00182                                       QCSymMatrix&    QCRestrict wDensityAA) {
00183 
00184  QC_TRACE("BEGIN QCDensityGtr<TPMatrix>::mergeDensity  ");
00185   
00186   
00187   // les parametre de A et B
00188   const TPParam * QCRestrict paramA;
00189   const TPParam * QCRestrict paramB;  
00190 
00191   // La premiere orbitale de
00192   int firstAOofA, firstAOofB;
00193   int nbAOofA, nbAOofB;
00194   int localFirstAOHigh, localFirstAOLow;
00195   int localNbAOLow, localNbAOHigh;
00196   int remoteNbAOLow, remoteNbAOHigh;
00197   int remoteLow, remoteHigh;
00198   bool atomAisHighRemote, atomAisHighLocal;
00199 
00200   // Les indices
00201   int localMu, localNu, localLda;
00202   int mu, nu, lda, a, b;
00203   //
00204   for (a = 0; a < workingDomain.getNbSharedAtoms(remoteMapIdx); ++a) {
00205     //    
00206     int localA  = workingDomain.getIndLocalIdx(remoteMapIdx, a);
00207     int remoteA = workingDomain.getIndRemoteIdx(remoteMapIdx, a);
00208     paramA      = &parameters[workingDomain.getParamIndexAt(localA)];
00209     firstAOofA  = workingDomain.getFirstAOAt(localA);
00210     nbAOofA     = paramA->getNbAO();
00211 
00212 #ifdef DEBUG_MERGE_DENSITY
00213     cout << endl << "theDiDkCellPtrIA->localNumInThis = " << localA << " pour " << table[paramA->getIdentifier()] << endl;
00214     cout << "theDiDkCellPtrIA->localNumInOther = " << remoteA << endl;
00215 #endif
00216     
00217     for (b = 0; b < a; ++b) {
00218 
00219       int localB  = workingDomain.getIndLocalIdx(remoteMapIdx, b);
00220       int remoteB = workingDomain.getIndRemoteIdx(remoteMapIdx, b);
00221       paramB      = &parameters[workingDomain.getParamIndexAt(localB)];
00222       firstAOofB  = workingDomain.getFirstAOAt(localB);
00223       nbAOofB     = paramB->getNbAO();
00224 
00225       if (firstAOofA > firstAOofB) {
00226 
00227         localFirstAOHigh = firstAOofA;
00228         localFirstAOLow  = firstAOofB;
00229         localNbAOHigh    = nbAOofA;
00230         localNbAOLow     = nbAOofB;
00231         atomAisHighLocal = true;
00232 
00233       } else {
00234         
00235         localFirstAOHigh = firstAOofB;
00236         localFirstAOLow  = firstAOofA;
00237         localNbAOHigh    = nbAOofB;
00238         localNbAOLow     = nbAOofA;
00239         atomAisHighLocal = false;
00240         
00241       }
00242 
00243       if (remoteA > remoteB) {
00244         
00245         remoteHigh        = remoteA;
00246         remoteLow         = remoteB;
00247         remoteNbAOHigh    = nbAOofA;
00248         remoteNbAOLow     = nbAOofB;
00249         atomAisHighRemote = true;
00250 
00251       } else {
00252         
00253         remoteHigh        = remoteB;
00254         remoteLow         = remoteA;
00255         remoteNbAOHigh    = nbAOofB;
00256         remoteNbAOLow     = nbAOofA;
00257         atomAisHighRemote = false;
00258         
00259       }
00260       
00261       // On extrait la sous matrice correspondante dans l'autre sous-domaine
00262       remWeightedDensity.extractSubMatrix(wDensityAB,
00263                                           remoteDomain.getFirstAOAt(remoteHigh), 
00264                                           remoteNbAOHigh,
00265                                           remoteDomain.getFirstAOAt(remoteLow),  
00266                                           remoteNbAOLow);
00267 
00268       // On  gere le cas ou les numeros ne sont pas dans le meme ordre dans
00269       // les 2 sous-systemes.
00270       if ( (atomAisHighLocal && !atomAisHighRemote) ||
00271            (!atomAisHighLocal && atomAisHighRemote) ) {
00272         
00273         interPAB.transpose(wDensityAB, 
00274                            remoteNbAOHigh, 
00275                            remoteNbAOLow);
00276         interPAB.swap(wDensityAB);
00277       }
00278       
00279       // Boucles sur les OA
00280       for (mu = 0, localMu = localFirstAOHigh; mu < localNbAOHigh; ++mu, ++localMu) {
00281         for (lda = 0, localLda = localFirstAOLow; lda < localNbAOLow; ++lda, ++localLda) {
00282 
00283 #ifdef DEBUG_MERGE_DENSITY
00284           cout << "avant : mergedDensityMatrix[ " << localMu << "][" << localLda << "] = " 
00285                << this->matrix[localMu][localLda] 
00286                << " et weightedDensityPABInOther[" << mu << "][ " << lda << "] = " 
00287                << wDensityAB[mu][lda] << endl;
00288 #endif 
00289 
00290           this->matrix[localMu][localLda] += wDensityAB[mu][lda];
00291           
00292 #ifdef DEBUG_MERGE_DENSITY
00293           cout << "apres : mergedDensityMatrix[ " << localMu << "][" << localLda << "] = " 
00294                << this->matrix[localMu][localLda] 
00295                << " et weightedDensityPABInOther[" << mu << "][ " << lda << "] = " 
00296                << wDensityAB[mu][lda] << endl;
00297 #endif 
00298 
00299         }
00300       }
00301     }
00302     // Le cas ou b == a
00303     // On extrait la sous matrice correspondante dans l'autre sous-domaine
00304     remWeightedDensity.extractSubTriangle(wDensityAA, remoteDomain.getFirstAOAt(remoteA), nbAOofA);
00305     // Boucles sur les OA
00306     for (mu = 0, localMu = firstAOofA; mu < nbAOofA; ++mu, ++localMu) {
00307       for (nu = 0, localNu = firstAOofA; nu <= mu; ++nu, ++localNu) {
00308         this->matrix[localMu][localNu] += wDensityAA[mu][nu];
00309       }
00310     }
00311   }
00312  QC_TRACE("END   QCDensityGtr<TPMatrix>::mergeDensity  ");
00313 }
00317 template <class TPMatrix>
00318 template <class TPMatrix2>
00319 void 
00320 QCDensityGtr<TPMatrix>::computeElemsDC (QCMemory&                QCRestrict memory,
00321                                         const int                nbFilledOM,
00322                                         const int                nbOccupiedOM,
00323                                         const int                nbAOCoreShell1,
00324                                         const QCFloat *          QCRestrict orbitalOccupN,
00325                                         QCMatElemGtr<TPMatrix2>& QCRestrict eigenVectCt) {
00326 
00327 
00328   // La matrice des vecteurs propres.
00329   TPMatrix2& QCRestrict transEigenVectMatrix = eigenVectCt.getMatrix();
00330   
00331   const int dim        = this->matrix.getDim();
00332   int i;
00333   QCFloat* QCRestrict workingMatrix =  memory.template takeTmpMem<QCFloat>(dim*dim);
00334 
00335 #define QC_CUT_DENSITY
00336 
00343   for (i = nbFilledOM; i < nbOccupiedOM; ++i) {
00344 
00345 #ifdef QC_CUT_DENSITY
00346     QCCBlas(dscal)(nbAOCoreShell1, sqrt(orbitalOccupN[i] / QC_TWO), 
00347                    transEigenVectMatrix[i],  1);
00348 #else // QC_CUT_DENSITY
00349     QCCBlas(dscal)(dim, sqrt(orbitalOccupN[i] / QC_TWO), 
00350                    transEigenVectMatrix[i],  1);
00351 #endif // QC_CUT_DENSITY
00352 
00353   }
00360   // En fait, a cause du passage fortran, eigenVectMatrix contient
00361   // la matrice transposee de C et il faut mettre CblasTrans
00362   QCCBlas(dsyrk)(
00363 #ifdef USE_OPTIMIZED_CBLAS
00364                  CblasRowMajor,
00365                  CblasLower, CblasTrans,
00366 #else // USE_OPTIMIZED_CBLAS
00367                  CblasUpper, CblasNoTrans,
00368 #endif // USE_OPTIMIZED_CBLAS
00369                  //On ne considere pas les parties des matrices qui concernent
00370                  //la zone de cut.
00371 #ifdef QC_CUT_DENSITY
00372                  nbAOCoreShell1, nbOccupiedOM,
00373 #else // QC_CUT_DENSITY
00374                  dim, nbOccupiedOM,
00375 #endif // QC_CUT_DENSITY
00376                  QC_TWO, transEigenVectMatrix.getElems(), dim,
00377                  QC_ZERO, workingMatrix, dim);
00378   
00379   // Il faut maintenant recopier la workingMatrix dans la
00380   // matrice densite.
00381 #ifdef QC_CUT_DENSITY
00382   for (i = 0; i < nbAOCoreShell1; ++i) {
00383 #else // QC_CUT_DENSITY
00384   for (i = 0; i < dim; ++i) {
00385 #endif // QC_CUT_DENSITY 
00386     memcpy(this->matrix[i], &(workingMatrix[i*dim]), (i+1) * sizeof(QCFloat) );
00387   }
00388   memory.template giveBackTmpMem<QCFloat>(workingMatrix);
00389 }
00390 
00394 template class QCDensityGtr<QCSymMatrix>;
00395 
00396 QCMANAGER_ITER_METH_EXPL_INST(void QCDensityGtr<QCSymMatrix>::initDiag);
00397 
00398 template void
00399 QCDensityGtr<QCSymMatrix>::computeElems<QCMatrix>(QCMemory& QCRestrict ,
00400                                                   const int ,
00401                                                   QCMatElemGtr<QCMatrix>& QCRestrict);
00402 
00403 template void 
00404 QCDensityGtr<QCSymMatrix>::computeElemsDC<QCMatrix> (QCMemory&               QCRestrict memory,
00405                                                      const int               nbFilledOM,
00406                                                      const int               nbOccupiedOM,
00407                                                      const int               nbAOCoreShell1,
00408                                                      const QCFloat *         QCRestrict orbitalOccupN,
00409                                                      QCMatElemGtr<QCMatrix>& QCRestrict eigenVectCt);
00410 
00411 
00412 QCPARAMETER_METH_EXPL_INST_PARAM(void QCDensityGtr<QCSymMatrix>::mergeDensity,
00413                            SEVEN_PARAMS(QCSubDomain& QCRestrict,
00414                                         QCSubDomain& QCRestrict,
00415                                         QCSymMatrix& QCRestrict,
00416                                         int,
00417                                         QCMatrix&    QCRestrict,
00418                                         QCMatrix&    QCRestrict,
00419                                         QCSymMatrix& QCRestrict));

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