QCDistMDSystem.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 <algorithm>
00021 #include <iomanip>
00022 #include <float.h>
00023 
00024 #include "QCDistMDSystem.hpp"
00025 
00026 #include "QCCommon.hpp"
00027 #include "QCMacro.hpp"
00028 #include "QCTrace.hpp"
00029 #include "QCChrono.hpp"
00030 #include "QCTools.hpp"
00031 
00032 
00033 #include "QCMndo.hpp"
00034 #include "QCAm1.hpp"
00035 #include "QCPm3.hpp"
00036 #include "QCMndoParam.hpp"
00037 #include "QCAm1Param.hpp"
00038 #include "QCPm3Param.hpp"
00039 #include "QCDCAlgo.hpp"
00040 #include "QCManager.hpp"
00041 #include "QCIntgReader.hpp"
00042 
00043 
00044 #ifdef HAVE_MPI
00045 #  include <mpi.h>
00046 #endif
00047 
00051 #define  D_NB_REQUESTS   4
00052 #define  D_INFO_SIZE     3
00053 #define  D_INFO_SD       0
00054 #define  D_INFO_DIM      1
00055 
00056 //#define MPICL_TRACE
00057 
00058 #if defined(MPICL_TRACE) && defined(__QC_xlC__)
00059 
00060 #  include "pcontrol.h"
00061 
00062 //  extern "C" {
00063 
00064 //    extern void tracenode(int tracesize, int flush, int sync);
00065 //    extern void tracefiles(char *tempfile, char *permfile, int verbose);
00066 //    extern void tracelevel(int mpi, int user, int trace);
00067 
00068 //  }
00069 
00070 
00071 #endif
00072 
00073 static void 
00074 compute_min_max (QCFloatFloat * in, QCFloatFloat * inout, int  * len,
00075                  MPI_Datatype * dptr);
00076 
00080 template <class TPSolver>
00081 QCDistMDSystem<TPSolver>::QCDistMDSystem (void) : QCMDSystem<TPSolver>(),
00082                                                   myrank(-1), cartrank(-1), nproc(0),
00083                                                   maxOvrlpSize(0), maxOvrlpSDs(0),
00084                                                   maxNbAtom(0), maxNbAO(0), maxDensitySize(0),
00085                                                   pipeline(NULL), atomsBuffers(NULL),
00086                                                   orderBuffers(NULL), ovrlpBuffers(NULL),
00087                                                   zoneBuffers(NULL), blockBuffers(NULL),
00088                                                   densBuffers(NULL), sdsFromProc(NULL),
00089                                                   sdsIntoProc(NULL),
00090                                                   _nbLocalContrib(0),_localContrib(NULL),
00091                                                   remSendsIdx(0), remRecvsIdx(0), rbuffer(NULL), sbuffer(NULL),
00092                                                   rcounts(NULL),  displs(NULL),  nbOM(NULL),
00093                                                   totalAO(0),_totalAOAll(0)
00094                                                   
00095 {}
00096 template <class TPSolver>
00097 QCDistMDSystem<TPSolver>::QCDistMDSystem (int& argc, char **&  argv) : QCMDSystem<TPSolver>(),
00098                                                                        myrank(-1), cartrank(-1), nproc(0),
00099                                                                        maxOvrlpSize(0), maxOvrlpSDs(0),
00100                                                                        maxNbAtom(0), maxNbAO(0), maxDensitySize(0),
00101                                                                        pipeline(NULL), atomsBuffers(NULL),
00102                                                                        orderBuffers(NULL), ovrlpBuffers(NULL),
00103                                                                        zoneBuffers(NULL), blockBuffers(NULL),
00104                                                                        densBuffers(NULL), sdsFromProc(NULL),
00105                                                                        sdsIntoProc(NULL), 
00106                                                                        _nbLocalContrib(0),_localContrib(NULL),
00107                                                                        remSendsIdx(0),remRecvsIdx(0), rbuffer(NULL), sbuffer(NULL),
00108                                                                        rcounts(NULL),  displs(NULL), nbOM(NULL),
00109                                                                        totalAO(0),_totalAOAll(0)
00110 {
00111   QC_TRACE_INIT("BEGIN QCDistMDSystem<TPSolver>::QCDistMDSystem (....) ");
00112   this->initMPI(argc,argv);
00113   QC_TRACE_INIT("END   QCDistMDSystem<TPSolver>::QCDistMDSystem (....) ");
00114 }
00118 template <class TPSolver>
00119 QCDistMDSystem<TPSolver>::~QCDistMDSystem (void) 
00120 {
00121   QC_TRACE_END("BEGIN QCDistMDSystem<TPSolver>::~QCDistMDSystem  "<<myrank<<")");
00122   MPI_Barrier(QC_COMM_WORLD);
00123 
00124   this->freeBuffers();
00125   
00126   if (sdsFromProc) {
00127     delete [] sdsFromProc; sdsFromProc = NULL;
00128   }
00129   if (sdsIntoProc) {
00130     delete [] sdsIntoProc; sdsIntoProc = NULL;
00131   }
00132   if (rbuffer) {
00133     delete [] rbuffer; rbuffer = NULL;
00134   }
00135   if (sbuffer) {
00136     delete [] sbuffer; sbuffer = NULL;
00137   }
00138   if (rcounts) {
00139     delete [] rcounts; rcounts = NULL;
00140   }
00141   if (displs) {
00142     delete [] displs; displs = NULL;
00143   }
00144   if (nbOM) {
00145     delete [] nbOM; nbOM = NULL;
00146   }
00147   if (_localContrib) {
00148     delete [] _localContrib; _localContrib = NULL;
00149   }
00150 
00152   MPI_Op_free(&QC_MIN_MAX);
00153   MPI_Type_free(&stype);
00154   MPI_Type_free(&QC_2DOUBLE);
00155 
00157   MPI_Barrier(QC_COMM_WORLD);
00158   MPI_Finalize();
00159   QC_TRACE_END("INT   QCDistMDSystem<TPSolver>::~QCDistMDSystem  "<<myrank<<")");
00160 }
00161 
00165 // template <class TPSolver>
00166 // void
00167 // QCDistMDSystem<TPSolver>::readPartition (const string& name,
00168 //                                       const string& path) {
00169 //   QCReader reader;
00170 //   reader.readPartitionFile(name, *this, path);
00171 // }
00172 
00173 
00174 
00175 
00179 template <class TPSolver>
00180 void QCDistMDSystem<TPSolver>::allocDomains (int nb) {
00181   QC_TRACE_INIT("BEGIN QCDistMDSystem<TPSolver>::allocDomains");
00182 
00183   int q  =  nb / nproc;
00184   int r  =  nb % nproc;
00185   this->cptIndex = 0 ;
00186   QCMDSystem<TPSolver>::nbDomains      = q;
00187   if (myrank < r) {
00188     ++QCMDSystem<TPSolver>::nbDomains ;
00189     QCMDSystem<TPSolver>::offset = myrank * (q+1);
00190     
00191   } else {
00192     QCMDSystem<TPSolver>::offset = myrank * q + r;
00193   }
00194   
00195   if (QCMDSystem<TPSolver>::nbDomains > 0) {
00196 
00197     QCMDSystem<TPSolver>::allocDomains(QCMDSystem<TPSolver>::nbDomains, QCMDSystem<TPSolver>::offset);
00198     QCMDSystem<TPSolver>::totalNbDomains = nb;
00199 
00200   } else {
00201     
00202     cerr << "PROC " << myrank  << ": ** FATAL ERROR **: AT LEAST ONE DOMAIN PER PROCESS" << std::endl;
00203     exit(EXIT_FAILURE);
00204   }
00205   QC_TRACE_INIT("END   QCDistMDSystem<TPSolver>::allocDomains");
00206 }
00207 
00211 template <class TPSolver>
00212 void QCDistMDSystem<TPSolver>::allocStructures (int sbdId, int nbatoms,
00213                                                 int nbAO, int nbovlp, int *sbdIds, int *nbshared) {
00214 
00215   QC_TRACE_INIT("BEGIN QCDistMDSystem<TPSolver>::allocStructures (process "<<myrank<<")");
00216   
00217   if (isMine(sbdId)) {
00218     int localIndex           = sbdId - QCMDSystem<TPSolver>::offset ;
00219     _globalToLocalNum[sbdId] = localIndex;
00220     QC_TRACE_INIT("               do something on subdomain "<< sbdId);
00221     
00222     QCMDSystem<TPSolver>::domains[localIndex].allocStructures(nbatoms,nbAO, nbovlp, sbdIds, nbshared);
00223     QCMDSystem<TPSolver>::domains[localIndex].setGlobalIndex(sbdId);
00224     
00225   }
00226   else{
00227     QC_TRACE_INIT("               nothing to do with subdomain "<< sbdId);
00228 
00229   }
00230   QC_TRACE_INIT("END   QCDistMDSystem<TPSolver>::allocStructures ");
00231 }
00235 template <class TPSolver>
00236 void 
00237 QCDistMDSystem<TPSolver>::addAtom (const QCPoint3D& coords, int type, int globalIdx,
00238                                    int domainIdx, QCSubDomainZone  zone, int  atomIdx) {
00239   //  QC_TRACE("BEGIN   QCDistMDSystem<TPSolver>::addATom ");
00240   if (isMine(domainIdx)) {
00241     int localIdx = _globalToLocalNum[domainIdx] ;
00242     QCMDSystem<TPSolver>::addAtom(coords,type, globalIdx,localIdx, zone, atomIdx);
00243   }
00244   //  QC_TRACE("END   QCDistMDSystem<TPSolver>::addAtom ");
00245 }
00249 template <class TPSolver>
00250 void
00251 QCDistMDSystem<TPSolver>::fillMaps (int nbsbd, const int * domainIdxs, const QCSubDomainZone * zones,
00252                                     const int * atomIdxs) {
00253   if (nbsbd > 1) {
00254     for (int i=0; i < nbsbd-1; ++i) {
00255       for (int j=i+1; j < nbsbd; ++j) {
00256         if ( isMine(domainIdxs[i]) &&  zones[i] != QC_CORE &&   zones[j] != QC_SHELL2 ) {
00257           QCMDSystem<TPSolver>::domains[domainIdxs[i] - QCMDSystem<TPSolver>::offset].addOverlpIndirectionInfo(domainIdxs[j],
00258                                                                    atomIdxs[i],
00259                                                                    atomIdxs[j],
00260                                                                    zones[j]);
00261           
00262         }
00263         if ( isMine(domainIdxs[j]) &&  zones[j] != QC_CORE &&    zones[i] != QC_SHELL2) {
00264           QCMDSystem<TPSolver>::domains[domainIdxs[j] - QCMDSystem<TPSolver>::offset].addOverlpIndirectionInfo(domainIdxs[i],
00265                                                                    atomIdxs[j],
00266                                                                    atomIdxs[i],
00267                                                                    zones[i]);
00268           
00269         }
00270       }
00271     }
00272   }
00273 }
00274 
00275 
00279 template <class TPSolver>
00280 void
00281 QCDistMDSystem<TPSolver>::fillMaps ( const vector< QCAtomIn > & v ){
00282 //
00283   int nbsbd = v.size() ; 
00284   if (nbsbd > 1) {
00285     for (int i=0; i < nbsbd-1; ++i) {
00286       for (int j=i+1; j < nbsbd; ++j) {
00287         if ( isMine(v[i].numDomain) &&  v[i].typeZone != QC_CORE &&  v[j].typeZone != QC_SHELL2 ) {
00288           QCMDSystem<TPSolver>::domains[v[i].numDomain - QCMDSystem<TPSolver>::offset].addOverlpIndirectionInfo(v[j].numDomain,
00289                                                                 v[i].localNum, v[j].localNum,  v[j].typeZone);
00290           
00291         }
00292         if ( isMine(v[j].numDomain) &&  v[j].typeZone != QC_CORE &&   v[i].typeZone != QC_SHELL2) {
00293           QCMDSystem<TPSolver>::domains[v[j].numDomain - QCMDSystem<TPSolver>::offset].addOverlpIndirectionInfo(v[i].numDomain,
00294                                                                  v[j].localNum, v[i].localNum, v[i].typeZone);
00295           
00296         }
00297       }
00298     }
00299   }
00300 }
00301 
00305 template <class TPSolver>
00306 void QCDistMDSystem<TPSolver>::initMPI(int& argc, char **&  argv)
00307 {
00308   QC_TRACE_INIT("BEGIN QCDistMDSystem<TPSolver>::initMPI ");
00309   MPI_Group QC_GROUP;
00310   MPI_Init(&argc, &argv);
00311   MPI_Comm_group (MPI_COMM_WORLD, &QC_GROUP);
00312 
00314   //MPI_Group_excl(QC_GROUP, 1, 0, &QC_NEW_GROUP);
00315 
00316   MPI_Comm_create(MPI_COMM_WORLD,  QC_GROUP, &this->QC_COMM_WORLD);
00317 
00318   MPI_Comm_size(QC_COMM_WORLD, &this->nproc);
00319   MPI_Comm_rank(QC_COMM_WORLD, &this->myrank);
00320   //
00321 
00322 #if defined(MPICL_TRACE) && defined(__QC_xlC__)
00323   /* enable tracing */
00324   ostringstream ss;
00325   ss << QCCommon::outdir << PATH_SEPARATOR   << "trace.trf";
00326   
00327   MPI_Pcontrol(TRACEFILES, "", ss.str().c_str(), 0);
00328   MPI_Pcontrol(TRACELEVEL, 1, 1, 0);
00329   MPI_Pcontrol(TRACENODE, 1000000, 0, 1);
00330   MPI_Pcontrol(TRACESTATISTICS, 25, 1, 1, 1, 1);
00331 #endif
00332   QC_TRACE_INIT("END   QCDistMDSystem<TPSolver>::initMPI ");
00333 }
00334 
00335 template <class TPSolver>
00336 template <class TPManager>
00337 void QCDistMDSystem<TPSolver>::init (TPManager&    manager, const string& path)
00338 {
00339   QC_TRACE_INIT("BEGIN QCDistMDSystem<TPSolver>::init ");
00340   //
00341   //  ***************************************************
00342   //
00343   QCPartitioner<TQCDistMDSystem> * partitioner = manager.getPartitioner();
00344   
00345   if (manager.getGeneralData().isReadPartitionFromFile()) { 
00346     
00347     partitioner->readFromFile( (*this), const_cast<QCGeneralData&>(manager.getGeneralData()),
00348                                manager.getFiles(),  path);
00349   } 
00350   else { 
00351     partitioner->initData(const_cast<QCGeneralData&>(manager.getGeneralData()), manager.getFiles());
00352    
00353     if(partitioner->getPartitionType() == QC_DIXON_PART  &&   partitioner->getNbPartitions()!= this->getNbAA() ) {
00354        std::cerr << "Error For Dixon partitioning the number of partitions ("<<partitioner->getNbPartitions()<<")"
00355                  << " must be equal to the number of amino acids ("<<this->getNbAA() <<")"<<std::endl;
00356        exit(EXIT_FAILURE);
00357     }
00358     partitioner->partitioning(*this);
00359   }
00360   QCMDSystem<TPSolver>::_partitionType = partitioner->getPartitionType() ;
00361   //  std::cout << "      QCMDSystem<TPSolver>::nbDomain : " <<    QCMDSystem<TPSolver>::nbDomains <<std::endl    ;
00362   //
00363   //
00365   //
00366   int localMaxOvrlp       = 0;
00367   int sumOvrlpSD          = 0;
00368   int localMaxNbAtoms     = 0;
00369   int localMaxDensitySize = 0;
00370   int localMaxNbAO        = 0;
00371   //
00372   for (int i = 0; i < QCMDSystem<TPSolver>::nbDomains; ++i) {
00373     QCMDSystem<TPSolver>::domains[i].sortOverlapments();
00374     QCMDSystem<TPSolver>::domains[i].setSystemNumbers(manager.getParameters());
00375     //
00376     //    QCMDSystem<TPSolver>::domains[i].updateInternalStructure();
00377     //
00378     QCMDSystem<TPSolver>::domains[i].allocateMatrices(manager, 
00379                                                       QCMDSystem<TPSolver>::domains[i].getNbAtomicOrbitals());
00380     QCMDSystem<TPSolver>::domains[i].fillWeight(manager.getGeneralData());
00381     //
00382     int ovrlpsize   = QCMDSystem<TPSolver>::domains[i].getOverlapMapSize();
00383     int ovrlpsds    = QCMDSystem<TPSolver>::domains[i].getNbOverlpSD();
00384     int nbatoms     = QCMDSystem<TPSolver>::domains[i].getNbAtoms();
00385     int nbao        = QCMDSystem<TPSolver>::domains[i].getNbAtomicOrbitals();
00386     int densitysize = QCMDSystem<TPSolver>::domains[i].getDensityBlockCount();
00387     //
00388 #ifdef TRACE_MPI
00389     std::cout << "    SD  " << i <<"     ovrlpsize : " << ovrlpsize << " ovrlpsds " 
00390               << ovrlpsds <<" nbatoms " << nbatoms << " nbao " << nbao 
00391               << " densitysize " << densitysize <<std::endl ;
00392 #endif
00393     if (localMaxOvrlp < ovrlpsize) { localMaxOvrlp = ovrlpsize; }
00394     if (localMaxNbAtoms < nbatoms) { localMaxNbAtoms = nbatoms; }
00395     if (localMaxDensitySize < densitysize) { localMaxDensitySize = densitysize;  }
00396     if (localMaxNbAO < nbao) { localMaxNbAO = nbao; }
00397     //
00398     sumOvrlpSD    += ovrlpsds;
00399   }
00400   //
00401   //   La version DISTRIBUEE
00402   //
00403   const int NB = 5 ;
00404 
00405   int sendbuf [NB], recvbuf [NB];
00406   sendbuf[0] = localMaxOvrlp;
00407   sendbuf[1] = sumOvrlpSD;
00408   sendbuf[2] = localMaxNbAtoms;
00409   sendbuf[3] = localMaxDensitySize;
00410   sendbuf[4] = localMaxNbAO;
00411 #ifdef  TRACE_MPI
00412   std::cout << "    Proc  " << myrank <<"  local    localMaxOvrlp : " << localMaxOvrlp << " sumOvrlpSD " 
00413             << sumOvrlpSD <<" localMaxNbAtoms " << localMaxNbAtoms << " localMaxNbAO " << localMaxNbAO 
00414             << " localMaxDensitySize " << localMaxDensitySize <<std::endl ;
00415 #endif
00416  
00417   //
00418   MPI_Allreduce(sendbuf, recvbuf, NB, MPI_INT, MPI_MAX, QC_COMM_WORLD);
00419   //
00420   maxOvrlpSize   = recvbuf[0];
00421   maxOvrlpSDs    = recvbuf[1];
00422   maxNbAtom      = recvbuf[2];
00423   maxDensitySize = recvbuf[3];
00424   maxNbAO        = recvbuf[4];
00425 
00426 #ifdef  TRACE_MPI
00427   std::cout << "    Proc  " << myrank <<"  global    globalMaxOvrlp : " << maxOvrlpSize << " maxOvrlpSDs " 
00428             << maxOvrlpSDs <<" globalMaxNbAtoms " << maxNbAtom << " globalMaxNbAO " << maxNbAO 
00429             << " globalMaxDensitySize " << maxDensitySize <<std::endl ;
00430 #endif
00431  
00432 //   // Allocation des buffers
00433 //   this->fillSDsCommArray();
00434 
00435 // #ifdef QC_DEBUG_DENSITY_SR
00436 //   std::cout << std::endl;
00437 //   printProcArray(sdsFromProc, false);
00438 //   printProcArray(sdsIntoProc, true);
00439 //   std::cout << std::endl;
00440 // #endif
00441 //   // Echange globale des listes de reception
00442 //   MPI_Alltoall(sdsFromProc, (maxOvrlpSDs + SDS_HEAD), MPI_INT,
00443 //             sdsIntoProc, (maxOvrlpSDs + SDS_HEAD), MPI_INT, QC_COMM_WORLD);
00444 //   //
00445 //   // Pour plus de facilite d'utilisation on construit
00446 //   // deux vecteurs correspondant aux echanges de messages 
00447 //   // a effectuer
00448 //   this->fillSDsCommVect();
00449   this->fillCommunicationStruct(maxOvrlpSDs) ;
00450   //
00451   this->allocBuffers(true);
00452   //
00453   MPI_Type_contiguous(2, MPI_DOUBLE, &QC_2DOUBLE);
00454   MPI_Type_commit(&QC_2DOUBLE);
00455   //
00456   MPI_Op_create(reinterpret_cast<MPI_User_function *>(compute_min_max),false,&QC_MIN_MAX);
00460   this->initFermiSeq();
00461 
00462 #if defined (QC_TRACE_IN_FILE)  
00463      std::cout << " Constante apres initialisation de la method QCDistMDSystem<TPSolver>::init" <<std::endl
00464                << " maxOvrlpSize   : " <<maxOvrlpSize<<std::endl
00465                << " maxOvrlpSDs    : "<<maxOvrlpSDs <<std::endl
00466                << " maxNbAtom      : "<<maxNbAtom<<std::endl
00467                << " maxDensitySize : "<<maxDensitySize<<std::endl
00468                << " maxNbAO        : "<<maxNbAO<<std::endl
00469                << " totalAO        : "<<totalAO<<std::endl
00470                << " totalAOAll     : "<<_totalAOAll<<std::endl;
00471    
00472 #endif
00473   QC_TRACE_INIT("END   QCDistMDSystem<TPSolver>::init ");
00474 }
00478 template <class TPSolver>
00479 template <class TPManager>
00480 void 
00481 QCDistMDSystem<TPSolver>::completeHamiltonMatrices (TPManager& manager) {
00482 
00483   QC_TRACE("BEGIN QCDistMDSystem<TPSolver>::completeFockMatrices ");
00484   typedef typename TPManager::TModel         TModel;
00485   typedef typename TPManager::TModel::TParam TParam;
00486   
00491   TModel&         QCRestrict model       = manager.getModel();
00492   const TParam *  QCRestrict params      = manager.getParameters();
00493   
00494   // matrices de Hamilton locales a cette methode
00495   QCSymMatrix& QCRestrict hamiltonHAA = model.getSpWorkingAA();
00496 
00497 
00498 #ifdef QC_VERBOSE
00499   std::cout << "PROC " << myrank  << ": * HAMILTON CONTRIBUTIONS *" << std::endl;
00500 #endif
00501 
00502   // Remplissage des buffers
00503   for (int i=0; i < QCMDSystem<TPSolver>::nbDomains; ++i) {
00504     
00505     QCMDSystem<TPSolver>::domains[i].copyAtomOverlpInto(atomsBuffers[i+1],
00506                                                         orderBuffers[i+1],
00507                                                         ovrlpBuffers[i+1],
00508                                                         zoneBuffers [i+1]);
00509   }
00510   //
00511   // On ajuste les choses pour MPI
00512   //
00513   int sendIdx = QCMDSystem<TPSolver>::nbDomains, recvIdx = 0;
00514 
00515   const int H_NB_REQUESTS = 8 ;
00516 
00517   MPI_Request requests[H_NB_REQUESTS];
00518   MPI_Status  statuses[H_NB_REQUESTS];
00519 
00520 
00521   QCFloat  * atomSendBuf , * atomRecvBuf;
00522   int      * orderSendBuf, * orderRecvBuf, *ovrlpSendBuf, * ovrlpRecvBuf;
00523   QCSubDomainZone      * zoneSendBuf,  * zoneRecvBuf;
00524 
00525   QCChrono recvChrono;
00526 
00527   // On commence le pipelining
00528   for (int step=0; step < QCMDSystem<TPSolver>::totalNbDomains-1; ++step) {
00529 
00530     // initialisation des pointeurs les sends
00531     atomSendBuf  = atomsBuffers[sendIdx];
00532     orderSendBuf = orderBuffers[sendIdx];
00533     ovrlpSendBuf = ovrlpBuffers[sendIdx];
00534     zoneSendBuf  = zoneBuffers [sendIdx];
00535     
00536     
00537     // les recvs
00538     atomRecvBuf  = atomsBuffers[recvIdx];
00539     orderRecvBuf = orderBuffers[recvIdx];
00540     ovrlpRecvBuf = ovrlpBuffers[recvIdx];
00541     zoneRecvBuf  = zoneBuffers [recvIdx];
00542     
00543 
00544     int fromId, toId;
00545 
00546     fromId = (myrank + nproc - 1)%nproc;
00547     toId   = (myrank + 1)%nproc;
00548 
00550     MPI_Isend (atomSendBuf, maxNbAtom * QCAtoms::DIMENSION + ATOM_OFFSET,
00551                MPI_DOUBLE, toId, HCONTRIB_ATOM, QC_COMM_WORLD, &requests[0]);
00552     MPI_Isend (orderSendBuf, maxNbAtom * ORDR_CELL_SIZE,
00553                MPI_INT,   toId, HCONTRIB_ORDER, QC_COMM_WORLD, &requests[1]);
00554     MPI_Isend (ovrlpSendBuf, maxOvrlpSize,
00555                MPI_INT,   toId, HCONTRIB_OVRLP, QC_COMM_WORLD, &requests[2]);
00556     MPI_Isend (zoneSendBuf, maxNbAtom,
00557                MPI_INT,   toId, HCONTRIB_ZONE,  QC_COMM_WORLD, &requests[3]);
00558 
00559 
00561     MPI_Irecv (atomRecvBuf, maxNbAtom * QCAtoms::DIMENSION + ATOM_OFFSET,
00562                MPI_DOUBLE, fromId, HCONTRIB_ATOM, 
00563                QC_COMM_WORLD, &requests[H_NB_REQUESTS/2]);
00564     MPI_Irecv (orderRecvBuf, maxNbAtom * ORDR_CELL_SIZE,
00565                MPI_INT,   fromId, HCONTRIB_ORDER, 
00566                QC_COMM_WORLD, &requests[H_NB_REQUESTS/2 + 1]);
00567     MPI_Irecv (ovrlpRecvBuf, maxOvrlpSize,
00568                MPI_INT,   fromId, HCONTRIB_OVRLP, 
00569                QC_COMM_WORLD, &requests[H_NB_REQUESTS/2 + 2]);
00570     MPI_Irecv (zoneRecvBuf, maxNbAtom,
00571                MPI_INT,   fromId, HCONTRIB_ZONE,  
00572                QC_COMM_WORLD, &requests[H_NB_REQUESTS/2 + 3]);
00573 
00574 
00575 
00576 
00577     for (int i=1; i < QCMDSystem<TPSolver>::nbDomains; ++i) {
00578       int currIdx = (recvIdx+i)%(QCMDSystem<TPSolver>::nbDomains+1);
00579       pipeline[i].setAtomOverlpFrom(atomsBuffers[currIdx], orderBuffers[currIdx],
00580                                     ovrlpBuffers[currIdx], zoneBuffers [currIdx]);
00581       
00582 #ifdef QC_DEBUG_HAMILTON_CONTRIB_LEVEL2
00583       std::cout << "PROC " << myrank  << ": STEP " << step << " OF HAMILTON CONTRIB ( ): SD `"
00584            << QCMDSystem<TPSolver>::domains[i].getId()  << "' <= SD `" << pipeline[i].getId() 
00585            << "'" << std::endl;
00586 #endif
00587 
00588       QCMDSystem<TPSolver>::domains[i].getHamiltonH().completeElems (// Les params
00589                                                params,
00590                                                // le domaine courant
00591                                                QCMDSystem<TPSolver>::domains[i],
00592                                                // le domaine contribuant
00593                                                pipeline[i],
00594                                                // autres
00595                                                QCMDSystem<TPSolver>::domains[i].getRepInterElecIntegrals(),
00596                                                hamiltonHAA,
00597                                                manager.getGeneralData(),
00598                                                step,
00599                                                QCMDSystem<TPSolver>::totalNbDomains);
00600       
00601     }
00602 
00603     
00604     // On attend la reception des tableaux
00605     recvChrono.start();
00606 
00607     MPI_Waitall(H_NB_REQUESTS/2, requests + H_NB_REQUESTS/2, 
00608                 statuses + H_NB_REQUESTS/2);
00609     
00610     recvChrono.pause();
00611 
00612 
00613     int currIdx = (recvIdx)%(QCMDSystem<TPSolver>::nbDomains+1);
00614     pipeline[0].setAtomOverlpFrom(atomsBuffers[currIdx], orderBuffers[currIdx],
00615                                   ovrlpBuffers[currIdx], zoneBuffers [currIdx]);
00616 
00617 #ifdef QC_DEBUG_HAMILTON_CONTRIB_LEVEL2
00618     std::cout << "PROC " << myrank       << ": STEP " << step  << " OF HAMILTON CONTRIB (*): SD `"
00619          << QCMDSystem<TPSolver>::domains[0].getId() << "' <= SD `" << pipeline[0].getId() 
00620          << "' WAIT=" << setprecision(6) << recvChrono.getvalsec() << " sec" << std::endl;
00621 #endif
00622     
00623     QCMDSystem<TPSolver>::domains[0].getHamiltonH().completeElems (// Les params
00624                                              params,
00625                                              // le domaine courant
00626                                              QCMDSystem<TPSolver>::domains[0],
00627                                              // le domaine contribuant
00628                                              pipeline[0],
00629                                              // autres
00630                                              QCMDSystem<TPSolver>::domains[0].getRepInterElecIntegrals(),
00631                                              hamiltonHAA,
00632                                              manager.getGeneralData(),
00633                                              step,
00634                                              QCMDSystem<TPSolver>::totalNbDomains);
00635 
00636 
00637     // On attend les envois
00638     MPI_Waitall(H_NB_REQUESTS/2, requests, statuses);
00639     
00640 
00641     // mise a jour des indexes pour les buffers
00642     // d'envoi et de reception
00643     recvIdx = (recvIdx + QCMDSystem<TPSolver>::nbDomains) % (QCMDSystem<TPSolver>::nbDomains+1);
00644     sendIdx = (sendIdx + QCMDSystem<TPSolver>::nbDomains) % (QCMDSystem<TPSolver>::nbDomains+1);
00645 
00646 
00647   }
00648 
00649 
00650 #ifdef QC_OUTPUT_HAMILTON_CONTRIB 
00651   for (int i=0; i < QCMDSystem<TPSolver>::nbDomains; ++i) {
00652     ostringstream osstr;
00653     osstr << "hamilton_dist_mine_" << QCMDSystem<TPSolver>::domains[i].getId(); 
00654     QCMDSystem<TPSolver>::domains[i].getHamiltonH().getMatrix().printInFile(osstr.str().c_str());
00655   }
00656 #endif
00657 }
00658 
00659 
00660 
00661 
00665 template <class TPSolver>
00666 template <class TPManager>
00667 QCFloat 
00668 QCDistMDSystem<TPSolver>::completeFockMatricesAndElecEnergy (TPManager& manager, bool isFirstCall, 
00669                                                              int /*iter*/) {
00670 
00671   QC_TRACE("BEGIN QCDistMDSystem<TPSolver>::completeFockMatricesAndElecEnergy ");
00672 
00673   typedef typename TPManager::TModel         TModel;
00674   typedef typename TPManager::TModel::TParam TParam;
00675 
00676 #if defined(MPICL_TRACE) && defined(__QC_xlC__)
00677   int work,  task = iter;
00678   MPI_Pcontrol(TRACEEVENT, "entry", task, 0);
00679 #endif
00680 
00681 
00686   TModel&         QCRestrict model  = manager.getModel();
00687   const TParam *  QCRestrict params = manager.getParameters();
00688 
00689   // Les parametrage de l'application
00690   const QCGeneralData& data = manager.getGeneralData();
00691 
00692   // matrices de Hamilton locales a cette methode
00693   QCSymMatrix& QCRestrict fockFAA    = model.getSpWorkingAA();
00694 
00695   // l'energie totale
00696   QCFloat  domainElecEnergy;
00697 
00698 
00699 #ifdef QC_VERBOSE
00700   std::cout << "PROC " << myrank << ": * FOCK CONTRIBUTIONS *" << std::endl;
00701 #endif
00702 
00703 
00704   // Remplissage des buffers
00705   for (int i=0; i < QCMDSystem<TPSolver>::nbDomains; ++i) {
00706 
00707     QCMDSystem<TPSolver>::domains[i].copyAtomOverlpInto(atomsBuffers[i+1],
00708                                                         orderBuffers[i+1],
00709                                                         ovrlpBuffers[i+1],
00710                                                         zoneBuffers [i+1],
00711                                                         blockBuffers[i+1]);
00712     QCMDSystem<TPSolver>::domains[i].initWeightedDensity(params);
00713     QCMDSystem<TPSolver>::domains[i].fillDensityArray(params, densBuffers[i+1], false);
00714   }
00715   // On ajouste les choses pou MPI
00716   int sendIdx = QCMDSystem<TPSolver>::nbDomains;
00717   int recvIdx = 0;
00718 
00719 #define  F_NB_REQUESTS  12
00720 
00721   MPI_Request requests[F_NB_REQUESTS];
00722   MPI_Status  statuses[F_NB_REQUESTS];
00723 
00724 
00725   QCFloat  * atomSendBuf , * atomRecvBuf;
00726   int      * orderSendBuf, * orderRecvBuf;
00727   int      * ovrlpSendBuf, * ovrlpRecvBuf;
00728   QCSubDomainZone      * zoneSendBuf,  * zoneRecvBuf;
00729   int      * blockSendBuf, * blockRecvBuf;
00730   QCFloat  * densSendBuf,  * densRecvBuf;
00731 
00732   QCChrono recvChrono;
00733 
00734 
00735   // Try synchronizing processes
00736   //  MPI_Barrier(QC_COMM_WORLD);
00737 
00738 
00739   // On commence le pipelining
00740   for (int step=0; step < QCMDSystem<TPSolver>::totalNbDomains-1; ++step) {
00741 
00742     // initialisation des pointeurs
00743     // les sends
00744     atomSendBuf  = this->atomsBuffers[sendIdx];
00745     orderSendBuf = this->orderBuffers[sendIdx];
00746     ovrlpSendBuf = this->ovrlpBuffers[sendIdx];
00747     zoneSendBuf  = this->zoneBuffers [sendIdx];
00748     blockSendBuf = this->blockBuffers[sendIdx];
00749     densSendBuf  = this->densBuffers [sendIdx];
00750     
00751     
00752     // les recvs
00753     atomRecvBuf  = this->atomsBuffers[recvIdx];
00754     orderRecvBuf = this->orderBuffers[recvIdx];
00755     ovrlpRecvBuf = this->ovrlpBuffers[recvIdx];
00756     zoneRecvBuf  = this->zoneBuffers [recvIdx];
00757     blockRecvBuf = this->blockBuffers[recvIdx];
00758     densRecvBuf  = this->densBuffers [recvIdx];
00759 
00760 
00761     int fromId, toId;
00762 
00763     fromId = (myrank + nproc - 1)%nproc;
00764     toId   = (myrank + 1)%nproc;
00765     //
00766     //    A améliorer cela doit etre dans un seul message
00767     //
00768     MPI_Isend (atomSendBuf, maxNbAtom * QCAtoms::DIMENSION + ATOM_OFFSET,
00769                MPI_DOUBLE, toId, FCONTRIB_ATOM, QC_COMM_WORLD, &requests[0]);
00770     MPI_Isend (orderSendBuf, maxNbAtom * ORDR_CELL_SIZE,
00771                MPI_INT,   toId, FCONTRIB_ORDER, QC_COMM_WORLD, &requests[1]);
00772     MPI_Isend (ovrlpSendBuf, maxOvrlpSize,
00773                MPI_INT,   toId, FCONTRIB_OVRLP, QC_COMM_WORLD, &requests[2]);
00774     MPI_Isend (zoneSendBuf, maxNbAtom,
00775                MPI_INT,   toId, FCONTRIB_ZONE,  QC_COMM_WORLD, &requests[3]);
00776     MPI_Isend (blockSendBuf, maxNbAtom,
00777                MPI_INT,   toId, FCONTRIB_BLOCK, QC_COMM_WORLD, &requests[4]);
00778     MPI_Isend (densSendBuf, maxDensitySize,
00779                MPI_DOUBLE, toId, FCONTRIB_DENS, QC_COMM_WORLD, &requests[5]);
00781     MPI_Irecv (atomRecvBuf, maxNbAtom * QCAtoms::DIMENSION + ATOM_OFFSET,
00782                MPI_DOUBLE, fromId, FCONTRIB_ATOM, 
00783                QC_COMM_WORLD, &requests[F_NB_REQUESTS/2]);
00784     MPI_Irecv (orderRecvBuf, maxNbAtom * ORDR_CELL_SIZE,
00785                MPI_INT,   fromId, FCONTRIB_ORDER, 
00786                QC_COMM_WORLD, &requests[F_NB_REQUESTS/2 + 1]);
00787     MPI_Irecv (ovrlpRecvBuf, maxOvrlpSize,
00788                MPI_INT,   fromId, FCONTRIB_OVRLP, 
00789                QC_COMM_WORLD, &requests[F_NB_REQUESTS/2 + 2]);
00790     MPI_Irecv (zoneRecvBuf, maxNbAtom,
00791                MPI_INT,   fromId, FCONTRIB_ZONE,  
00792                QC_COMM_WORLD, &requests[F_NB_REQUESTS/2 + 3]);
00793     MPI_Irecv (blockRecvBuf, maxNbAtom,
00794                MPI_INT,   fromId, FCONTRIB_BLOCK, 
00795                QC_COMM_WORLD, &requests[F_NB_REQUESTS/2 + 4]);
00796     MPI_Irecv (densRecvBuf, maxDensitySize,
00797                MPI_DOUBLE, fromId, FCONTRIB_DENS, 
00798                QC_COMM_WORLD, &requests[F_NB_REQUESTS/2 + 5]);
00799 
00800     
00801     for (int i=1; i < QCMDSystem<TPSolver>::nbDomains; ++i) {
00802       int currIdx = (recvIdx+i)%(QCMDSystem<TPSolver>::nbDomains+1);
00803       pipeline[i].setAtomOverlpFrom(atomsBuffers[currIdx],
00804                                     orderBuffers[currIdx],
00805                                     ovrlpBuffers[currIdx],
00806                                     zoneBuffers [currIdx],
00807                                     blockBuffers[currIdx]);
00808       
00809 #ifdef QC_DEBUG_FOCK_CONTRIB_LEVEL2
00810       std::cout << "PROC " << myrank  << ": STEP " << step << " OF FOCK CONTRIB ( ): SD `"
00811            << QCMDSystem<TPSolver>::domains[i].getId()<< "' <= SD `" << pipeline[i].getId() 
00812            << "'" << std::endl;
00813 #endif
00814 
00815       
00816       QCMDSystem<TPSolver>::domains[i].getFockF().completeElems (// Les params
00817                                            params,
00818                                            // le domaine courant
00819                                            QCMDSystem<TPSolver>::domains[i],
00820                                            // le domaine contribuant
00821                                            pipeline[i],
00822                                            densBuffers[currIdx],
00823                                            // autres
00824                                            QCMDSystem<TPSolver>::domains[i].getRepInterElecIntegrals(),
00825                                            fockFAA, data, step, QCMDSystem<TPSolver>::totalNbDomains,
00826                                            data.getIntgAcquisitionMethod(), isFirstCall);
00827     }
00828     // On attend la reception des tableaux
00829     recvChrono.start();
00830     MPI_Waitall(F_NB_REQUESTS/2, requests + F_NB_REQUESTS/2, statuses + F_NB_REQUESTS/2);
00831     recvChrono.pause();
00832     //    
00833     int currIdx = (recvIdx)%(QCMDSystem<TPSolver>::nbDomains+1);
00834     pipeline[0].setAtomOverlpFrom(atomsBuffers[currIdx], orderBuffers[currIdx],
00835                                   ovrlpBuffers[currIdx], zoneBuffers [currIdx],
00836                                   blockBuffers[currIdx]);
00837     //
00838 #ifdef QC_DEBUG_FOCK_CONTRIB_LEVEL2
00839     std::cout << "PROC " << myrank << ": STEP " << step << " OF FOCK CONTRIB (*): SD `"
00840          << QCMDSystem<TPSolver>::domains[0].getId() << "' <= SD `" << pipeline[0].getId() 
00841          << "' WAIT=" << setprecision(6) << recvChrono.getvalsec() << " sec" << std::endl;
00842 #endif
00843     //
00844     QCMDSystem<TPSolver>::domains[0].getFockF().completeElems (// Les params
00845                                            params,
00846                                            // le domaine courant
00847                                            QCMDSystem<TPSolver>::domains[0],
00848                                            // le domaine contribuant
00849                                            pipeline[0], densBuffers[currIdx],
00850                                            // autres
00851                                            QCMDSystem<TPSolver>::domains[0].getRepInterElecIntegrals(),
00852                                            fockFAA, data, step, QCMDSystem<TPSolver>::totalNbDomains,  data.getIntgAcquisitionMethod(),
00853                                            isFirstCall);
00854 
00855 
00856     // On attend les envois
00857     MPI_Waitall(F_NB_REQUESTS/2, requests, statuses);
00858     //
00859     // mise a jour des indexes pour les buffers d'envoi et de reception
00860     //
00861     recvIdx = (recvIdx + QCMDSystem<TPSolver>::nbDomains) % (QCMDSystem<TPSolver>::nbDomains+1);
00862     sendIdx = (sendIdx + QCMDSystem<TPSolver>::nbDomains) % (QCMDSystem<TPSolver>::nbDomains+1);
00863   }
00864 #ifdef QC_VERBOSE
00865   std::cout << "PROC " << myrank  << ": * FOCK CONTRIBUTIONS (END RING)*" << std::endl;
00866 #endif
00867 #if defined(MPICL_TRACE) && defined(__QC_xlC__)
00868   /* finish task */
00869   work = iter;
00870   MPI_Pcontrol(TRACEEVENT, "exit", task, 1, &work);
00871 #endif
00872 #if defined (QC_TRACE_IN_FILE)  
00873   for (int i=0; i < QCMDSystem<TPSolver>::nbDomains; ++i) {
00874     qctrace[i].writeFockMatrix();
00875     }
00876 #endif
00877 
00878 #ifdef QC_OUTPUT_FOCK_CONTRIB
00879 
00880   for (int i=0; i < QCMDSystem<TPSolver>::nbDomains; ++i) {
00881     ostringstream osstr;
00882     osstr << "fock_dist_mine_" << QCMDSystem<TPSolver>::domains[i].getId();
00883     QCMDSystem<TPSolver>::domains[i].getFockF().getMatrix().printInFile(osstr.str().c_str());
00884   }
00885   
00886 #endif  
00887   //
00888   bool rewind = (data.getIntgAcquisitionMethod() == QC_INDIRECT_STORAGE &&  !isFirstCall);
00889   //
00890   // On prepare le calcul de l'energie
00891   //
00892   this->elecEnergy = QC_ZERO;
00893   //
00894   // Local contribution : 
00895   //
00896   for (int i=0; i < QCMDSystem<TPSolver>::nbDomains; ++i) {
00897     // Retour au debut des fichiers d'intégrales de repultion
00898     if (rewind) {
00899       QCMDSystem<TPSolver>::domains[i].getRepInterElecIntegrals().getIntgReader()->rewindFile();
00900     }
00901     // Calcul de l'energie electronique sur chaque domaine
00902     QCIterator iter (this, i);
00903     domainElecEnergy = model.computeElecEnergy(iter,
00904                                                *QCMDSystem<TPSolver>::domains[i].getWeightedDensityP());
00905 #ifdef QC_VERBOSE_ENERGY
00906     std::cout << setprecision(20) << "PROC " << myrank  << ": ELECTRONIC ENERGY (SD " << setw(3) << i << ") = "
00907          << domainElecEnergy << std::endl;
00908 #endif
00909 #if defined (QC_TRACE_IN_FILE)  
00910     qctrace[i].writeVal("* Local ELECTRONIC ENERGY of my domain  = " ,domainElecEnergy);
00911 #endif
00912     this->elecEnergy += domainElecEnergy;
00913   }
00914 #ifdef QC_VERBOSE_ENERGY
00915   std::cout << "PROC " << myrank  << ": ELECTRONIC ENERGY (ALL MY DOMAINS)  = " << this->elecEnergy << std::endl;
00916 #endif
00917   // 
00918 #if defined (QC_TRACE_IN_FILE)  
00919   for (int i=0; i < QCMDSystem<TPSolver>::nbDomains; ++i) {
00920     qctrace[i].writeVal("* Local ELECTRONIC ENERGY (ALL Domains in my DistDomain)  = " ,this->elecEnergy);
00921   }
00922 #endif
00923   //
00924   MPI_Allreduce(&(this->elecEnergy), &(this->totalElecEnergy), 1, MPI_DOUBLE, MPI_SUM, QC_COMM_WORLD);
00925 #ifdef QC_VERBOSE_ENERGY
00926   if (myrank == 0) {
00927     std::cout << std::endl  <<"* ELECTRONIC ENERGY (ALL PROCESSES)  = "  << setprecision(20)  << this->totalElecEnergy << " *" << std::endl << std::endl;
00928   }
00929 #endif  
00930 #if defined (QC_TRACE_IN_FILE)  
00931   for (int i=0; i < QCMDSystem<TPSolver>::nbDomains; ++i) {
00932     qctrace[i].writeVal("* Global ELECTRONIC ENERGY (all domains)  = " ,this->totalElecEnergy);
00933   }
00934 #endif
00935 
00936   QC_TRACE("END   QCDistMDSystem<TPSolver>::completeFockMatricesAndElecEnergy ");
00937   return this->totalElecEnergy;
00938 }
00939 
00940 /****************************************************************/
00941 /****************************************************************/
00945 template <class TPSolver>
00946 template <class TPManager>
00947 void 
00948 QCDistMDSystem<TPSolver>::completeDensityMatrices (TPManager& manager) {
00949   QC_TRACE_DENSITY("BEGIN QCDistMDSystem<TPSolver>::completeDensityMatrices");
00950   
00951   typedef typename TPManager::TModel         TModel;
00952   typedef typename TPManager::TModel::TParam TParam;
00957   TModel&         QCRestrict model  = manager.getModel();
00958   const TParam *  QCRestrict params = manager.getParameters();
00959   //
00960   // Les matrices densite
00961   QCDensityGtr<QCSymMatrix> * QCRestrict densityP;
00962   QCSymMatrix               * QCRestrict weightedDensity;
00963   QCSymMatrix                 remWeightedDensity;
00964   //
00965   // matrices de Hamilton locales a cette methode
00966   QCMatrix&    QCRestrict wDensityAB = model.getSpWorkingAB();
00967   QCMatrix&    QCRestrict interPAB   = model.getSpWorkingAB2();
00968   QCSymMatrix& QCRestrict wDensityAA = model.getSpWorkingAA();
00969   //  
00970   // le domaine apportant la contribution
00971   QCSubDomain  contributer;
00972   QCSymMatrix  remDensity;
00973   
00974 #ifdef QC_VERBOSE_DENSITY
00975   std::cout << "PROC " << myrank  << ": * DENSITY CONTRIBUTIONS *" << std::endl;
00976 #endif
00977   //
00978   // On copie la matrice densite ponderee dans la matrice densite
00979   //
00980   for (int i = 0; i < QCMDSystem<TPSolver>::nbDomains; ++i) {
00981     //
00982     densityP        = &QCMDSystem<TPSolver>::domains[i].getDensityP();
00983     weightedDensity = &QCMDSystem<TPSolver>::domains[i].getWeightedDensityP()->getMatrix();
00984     densityP->getMatrix().copy(*weightedDensity);
00985   }
00986   //
00987   MPI_Request sendReq [2][D_NB_REQUESTS];
00988   MPI_Request recvReq [2][D_NB_REQUESTS];
00989   MPI_Status  sendStat[2][D_NB_REQUESTS];
00990   MPI_Status  recvStat[2][D_NB_REQUESTS];
00991   //
00992   // INFO
00993   int  ierr, sendInfo[2][D_INFO_SIZE], recvInfo[2][D_INFO_SIZE];  
00994   int   currentStep, nbSteps, msgRead, msgSend, msgPost ;
00995   msgRead   = msgSend = currentStep = msgPost = 0 ; 
00996   nbSteps   = static_cast<int>(this->_sends.size()) ;
00997 
00998 #ifdef QC_DEBUG_DENSITY_SR
00999   std::cout << "                                RECVS                " << std::endl;
01000   printProcVect(_recvs, false);
01001   std::cout << "                                SEND                " << std::endl;
01002   printProcVect(_sends, true);
01003   std::cout << std::endl;
01004 #endif
01005   //
01006   //
01007   // On poste les premiers envoi/reception
01008   //    ++msgPost  ;
01009   int buffSend = -1 , buffRecv = -1 ;
01010   if ( _sends[currentStep].proc != -1 )  {
01011     buffSend                 = (buffSend+1)%2 ; // 0  first Buffer
01012     _sends[currentStep].step = buffSend ;
01013     postDensitySends(currentStep, sendInfo[buffSend], sendReq[buffSend]);   // 0  first Buffer
01014     //      std::cout << "             ISENDS done on buffer : " << buffSend << std::endl;
01015   }
01016   int recvIdx ;
01017   if ( _recvs[currentStep].proc != -1 ) {
01018     buffRecv                 = (buffRecv+1)%2 ; // 0  first Buffer
01019     _recvs[currentStep].step = buffRecv ;
01020     postDensityRecvs(currentStep, buffRecv, recvInfo[buffRecv], recvReq[buffRecv]); 
01021     //    std::cout << "             IRECVS done on buffer : " << buffRecv << std::endl;
01022 
01023   }
01024   //
01025   // On repercute les contributions locales
01026   // 
01027 #ifdef QC_DEBUG_DENSITY_CONTRIB_LEVEL2
01028   std::cout << "   BEGIN Local Contributers" << std::endl;
01029 #endif  
01030   for (int i = 0; i < _nbLocalContrib ; ++i) {
01031       int  SD1 =  _localContrib[2*i], SD2 =  _localContrib[2*i+1]  ;
01032 
01033 #ifdef QC_DEBUG_DENSITY_CONTRIB_LEVEL2
01034     std::cout << "         : (L) RECV" << ": CONTRIBUTER1 = `" << SD1    << ": CONTRIBUTER2 = `" << SD2 <<std::endl;
01035 #endif    
01036     this->mergeDensityAux(params,QCMDSystem<TPSolver>::domains[SD1],QCMDSystem<TPSolver>::domains[SD2],
01037                           QCMDSystem<TPSolver>::domains[SD2].getWeightedDensityP()->getMatrix(),
01038                           wDensityAB,interPAB, wDensityAA);
01039 // #ifdef QC_DEBUG_DENSITY_CONTRIB_LEVEL2
01040 //     std::cout << "         : (L) RECV" << ": CONTRIBUTER1 = `" << SD2    << ": CONTRIBUTER2 = `" << SD1 <<std::endl;
01041 // #endif    
01042 //     this->mergeDensityAux(params,QCMDSystem<TPSolver>::domains[SD2],QCMDSystem<TPSolver>::domains[SD1],
01043 //                        QCMDSystem<TPSolver>::domains[SD2].getWeightedDensityP()->getMatrix(),
01044 //                        wDensityAB,interPAB, wDensityAA);
01045   }
01046 #ifdef QC_DEBUG_DENSITY_CONTRIB_LEVEL2
01047   std::cout << "   END Local Contributers" << std::endl;
01048 #endif
01049   //
01050   // loop to read the remote density matrices
01051   //
01052 #ifdef QC_DEBUG_DENSITY_CONTRIB_LEVEL2
01053   std::cout << std::endl<< "               LOOPS ON REMOTE SD " << std::endl;
01054 #endif
01055   while ( currentStep <  nbSteps  ) { 
01056     //
01057 #ifdef QC_DEBUG_DENSITY_CONTRIB_LEVEL2
01058     std::cout << "         currentStep "<< currentStep << " buffSend : " << buffSend << "  buffRecv " << buffRecv << std::endl;
01059 #endif
01060     int sendIdx ;
01061     if( currentStep + 1 < nbSteps) {
01062      if (_sends[currentStep+1].proc != -1 )  {
01063         buffSend                   = (buffSend+1)%2 ; // 0  first Buffer
01064         _sends[currentStep+1].step = buffSend ;
01065         postDensitySends(currentStep+1, sendInfo[buffSend], sendReq[buffSend]);            //  Buffer
01066         //      std::cout << "             ISENDS done on buffer : " << buffSend << " for step : " << currentStep+1<< std::endl;
01067       }
01068       if ( _recvs[currentStep+1].proc != -1 )  {
01069         buffRecv                   = (buffRecv+1)%2 ; // 0  first Buffer
01070         _recvs[currentStep+1].step = buffRecv ;
01071         postDensityRecvs(currentStep+1, buffRecv, recvInfo[buffRecv], recvReq[buffRecv]);   //  Buffer
01072         //      std::cout << "             IRECVS done on buffer : " << buffRecv << " for step : " << currentStep+1<<std::endl;
01073       }
01074     }
01075     if (_sends[currentStep].proc != -1 ) { // on attend la fin de l'envoie du message sur le buffer 
01076       sendIdx = _sends[currentStep].step ;
01077 #ifdef QC_DEBUG_DENSITY_CONTRIB_LEVEL2
01078       std::cout << "      step : " << currentStep << " Wait end send to proc " 
01079                 << _sends[currentStep].proc << " on buffer : " <<sendIdx<<std::endl;
01080       std::cout << "                                SEND                " << std::endl;
01081       printProcVect(_sends, true);
01082       std::cout << std::endl;
01083 #endif      
01084 #ifdef  TRACE_MPI
01085       std::cout << "        MPI_Waitall (send) on index " << sendIdx <<std::endl ;
01086 #endif
01087       MPI_Waitall(D_NB_REQUESTS, sendReq[sendIdx], sendStat[sendIdx]);
01088       _sends[currentStep].step = currentStep ;
01089 #ifdef  TRACE_MPI
01090       std::cout << "                                   END MPI_Waitall SEND"  <<std::endl ;
01091  #endif
01092       ++msgSend ;
01093     }
01094     QCChrono recvChrono;
01095     if (_recvs[currentStep].proc != -1 ) {
01096       recvIdx = _recvs[currentStep].step;
01097 #ifdef QC_DEBUG_DENSITY_CONTRIB_LEVEL2
01098       std::cout << "   step : " << currentStep << "  Wait end receive to proc " << 
01099         _recvs[currentStep].proc << " on buffer : " <<recvIdx<<std::endl;
01100       printProcVect(_recvs, false);
01101 #endif      
01102       recvChrono.start();
01103 #ifdef  TRACE_MPI
01104       std::cout << "        MPI_Waitall receive on index  " << recvIdx <<std::endl ;
01105 #endif
01106       ierr = MPI_Waitall(D_NB_REQUESTS, recvReq[recvIdx], recvStat[recvIdx]);
01107       _recvs[currentStep].step = currentStep ;
01108 #ifdef  TRACE_MPI
01109       std::cout << "                                   END MPI_Waitall RCV" <<std::endl ;
01110 #endif
01111       recvChrono.pause();
01112       _recvs[currentStep].time = recvChrono.getvalsec();
01113       
01114 #ifdef QC_DEBUG_DENSITY_CONTRIB_LEVEL2
01115       std::cout << std::endl ;
01116       std::cout << "PROC " << myrank   << " step : " <<currentStep << ": (R) RECV "
01117                 << ": CONTRIBUTER=`(PROC:" << recvStat[recvIdx][0].MPI_SOURCE
01118                 << ", SD:" << recvInfo[recvIdx][D_INFO_SD]<< ")' EXPECTED=`(PROC:" << _recvs[currentStep].proc  
01119                 << ", SD:" << _recvs[currentStep].sbdId 
01120                 << ")' RECV IDX=" << recvIdx  << ", RECV Nb =" << currentStep << setprecision(6)
01121                 << ", WAIT=" << _recvs[currentStep].time << " s"   << std::endl;
01122       
01123       assert(recvInfo[recvIdx][D_INFO_SD]    == _recvs[currentStep].sbdId);
01124       assert(recvStat[recvIdx][0].MPI_SOURCE == _recvs[currentStep].proc);
01125 #endif
01126       //     
01127       // BEGIN contribution
01128       //
01129 #ifdef QC_DEBUG_DENSITY_CONTRIB_LEVEL2
01130       std::cout << "         BEGIN REMOTE CONTRIBUTION " << std::endl;
01131 #endif
01132       contributer.setAtomOverlpFrom(NULL,orderBuffers[recvIdx],ovrlpBuffers[recvIdx], NULL);      
01133       contributer.setId(_recvs[currentStep].sbdId);
01134       remWeightedDensity.setDimAndElems(recvInfo[recvIdx][D_INFO_DIM],wdensMatBuffs[recvIdx]);
01135       mergeDensityAux(params,contributer,remWeightedDensity, wDensityAB,interPAB,wDensityAA);
01136       contributer.unsetAtomOverlp();
01137       remWeightedDensity.unsetElems();
01138 #ifdef QC_DEBUG_DENSITY_CONTRIB_LEVEL2
01139       std::cout << "         END   REMOTE CONTRIBUTION " << std::endl;
01140 #endif
01141       //
01142       // END contribution
01143       //
01144     }
01145     ++currentStep; 
01146   }
01147 #ifdef OUTPUT_DENSITY_CONTRIB
01148   for (int i=0; i < QCMDSystem<TPSolver>::nbDomains; ++i) {
01149     ostringstream osstr;
01150     osstr << "merged_density_dist_mine_" << QCMDSystem<TPSolver>::domains[i].getId();
01151     QCMDSystem<TPSolver>::domains[i].getDensityP().getMatrix().printInFile(osstr.str().c_str());
01152   }
01153 #endif
01154   QC_TRACE_DENSITY("END   QCDistMDSystem<TPSolver>::completeDensityMatrices");
01155 }
01156 
01157 template <class TPSolver>
01158 void QCDistMDSystem<TPSolver>::postDensitySends (int currSend, int * sendInfo, MPI_Request *requests) {
01159   //
01160   QC_TRACE_DENSITY("BEGIN QCDistMDSystem<TPSolver>::postDensitySends NEW -- step : "<< currSend);
01161 
01162   // Pointeurs sur les buffers
01163   int       * orderSendBuf;
01164   int       * ovrlpSendBuf;
01165   QCFloat   * wdensSendMat;
01166 
01167 
01168   int sdtosend = _sends[currSend].sbdId;
01169   int toId     = _sends[currSend].proc;
01170   int step     = currSend ;
01171   int mpierr   = 0 ;
01172   int localNumberSD = sdtosend - QCMDSystem<TPSolver>::offset ; 
01173   //
01174   orderSendBuf = QCMDSystem<TPSolver>::domains[localNumberSD].getOrderingCells();
01175   ovrlpSendBuf = QCMDSystem<TPSolver>::domains[localNumberSD].getOverlapMap();
01176   wdensSendMat = QCMDSystem<TPSolver>::domains[localNumberSD].getWeightedDensityP()->getMatrix().getElems();
01177   //
01178   sendInfo[D_INFO_SD]  = sdtosend;
01179   sendInfo[D_INFO_DIM] = QCMDSystem<TPSolver>::domains[localNumberSD].getWeightedDensityP()->getMatrix().getDim();
01180   mpierr = MPI_Isend (sendInfo, D_INFO_SIZE, MPI_INT,    toId, (step * D_NB_REQUESTS) + DCONTRIB_INFO,  
01181                       QC_COMM_WORLD, &requests[0]);
01182   
01183   mpierr = MPI_Isend (orderSendBuf, maxNbAtom * ORDR_CELL_SIZE, MPI_INT, toId, (step * D_NB_REQUESTS) + DCONTRIB_ORDER, 
01184                       QC_COMM_WORLD, &requests[1]);
01185 
01186   mpierr = MPI_Isend (ovrlpSendBuf, maxOvrlpSize, MPI_INT, toId, (step * D_NB_REQUESTS) + DCONTRIB_OVRLP, 
01187                       QC_COMM_WORLD, &requests[2]);
01188 
01189   mpierr = MPI_Isend (wdensSendMat, (maxNbAO * (maxNbAO+1))/2, MPI_DOUBLE, toId, (step * D_NB_REQUESTS) + DCONTRIB_DENS,  
01190                       QC_COMM_WORLD, &requests[3]);
01191   
01192 #ifdef QC_DEBUG_DENSITY_CONTRIB_LEVEL2
01193   std::cout << "              PROC " << myrank     << ": DONE POST SEND `sd" << sdtosend
01194             << "' AT STEP `" << step << "' TO PROC `" << toId << "' on buffer `" << _sends[currSend].step << " with tags : " 
01195             << (step * D_NB_REQUESTS) + DCONTRIB_INFO          << " " 
01196             << (step * D_NB_REQUESTS) + DCONTRIB_ORDER         << " "  
01197             << (step * D_NB_REQUESTS) + DCONTRIB_OVRLP         << " " 
01198             << (step * D_NB_REQUESTS) + DCONTRIB_DENS          << " "  << std::endl
01199             << " Sizes of requests : " << D_INFO_SIZE << "  " << maxNbAtom * ORDR_CELL_SIZE 
01200             << "   " << maxOvrlpSize <<"  "<< (maxNbAO * (maxNbAO+1))/2 
01201             << std::endl;
01202 #endif
01203 
01204   QC_TRACE_DENSITY("END   QCDistMDSystem<TPSolver>::postDensitySends NEW");  
01205 }
01209 template <class TPSolver>
01210 void
01211 QCDistMDSystem<TPSolver>::postDensityRecvs (int currRecv, int recvIdx, int *recvInfo, MPI_Request *requests) {
01212   //
01213   QC_TRACE_DENSITY("BEGIN QCDistMDSystem<TPSolver>::postDensityRecvs NEW -- step : " << currRecv);
01214   // Pointeurs sur les buffers
01215   int       *orderRecvBuf;
01216   int       *ovrlpRecvBuf;
01217   QCFloat   *wdensRecvMat;
01218   int        mpierr = 0;
01219 
01220   int        fromId = _recvs[currRecv].proc;
01221   int        step   = currRecv ;
01222   //
01223   orderRecvBuf = orderBuffers [recvIdx];
01224   ovrlpRecvBuf = ovrlpBuffers [recvIdx];
01225   wdensRecvMat = wdensMatBuffs[recvIdx];
01226    
01227 //  #ifdef  TRACE_MPI
01228 //   std::cout << "            MPI_Irecv  from " << fromId << " with tags : " 
01229 //          << (step * D_NB_REQUESTS) + DCONTRIB_INFO          << " " 
01230 //          << (step * D_NB_REQUESTS) + DCONTRIB_ORDER         << " " 
01231 //          << (step * D_NB_REQUESTS) + DCONTRIB_OVRLP         << " " 
01232 //          << (step * D_NB_REQUESTS) + DCONTRIB_DENS          << " " 
01233 //          << std::endl ;
01234 //  #endif
01235   mpierr = MPI_Irecv (recvInfo, D_INFO_SIZE, MPI_INT, fromId, (step * D_NB_REQUESTS) + DCONTRIB_INFO,  
01236                       QC_COMM_WORLD, &requests[0]);
01237 //  #ifdef  TRACE_MPI
01238 //    std::cout << "                                 ierr : " << mpierr <<std::endl;
01239 //    std::cout << "  MPI_Irecv  from " << fromId << " request 1 " <<std::endl;
01240 //  #endif
01241   
01242   mpierr =  MPI_Irecv (orderRecvBuf, maxNbAtom * ORDR_CELL_SIZE,
01243              MPI_INT,    fromId, (step * D_NB_REQUESTS) + DCONTRIB_ORDER,  
01244              QC_COMM_WORLD, &requests[1]);
01245 //  #ifdef  TRACE_MPI
01246 //    std::cout << "                                 ierr : " << mpierr <<std::endl;
01247 //    std::cout << "  MPI_Irecv  from " << fromId << " request 2 " <<std::endl;
01248 //  #endif
01249   
01250   mpierr =  MPI_Irecv (ovrlpRecvBuf, maxOvrlpSize, MPI_INT, fromId, (step * D_NB_REQUESTS) + DCONTRIB_OVRLP,  
01251                        QC_COMM_WORLD, &requests[2]);
01252 //  #ifdef  TRACE_MPI
01253 //    std::cout << "                                 ierr : " << mpierr <<std::endl;
01254 //    std::cout << "  MPI_Irecv  from " << fromId << " request 3 " <<std::endl;
01255 //  #endif
01256   
01257   mpierr =  MPI_Irecv (wdensRecvMat, (maxNbAO * (maxNbAO+1))/2, MPI_DOUBLE, fromId, (step * D_NB_REQUESTS) + DCONTRIB_DENS,  
01258                        QC_COMM_WORLD, &requests[3]);
01259 
01260 #ifdef QC_DEBUG_DENSITY_CONTRIB_LEVEL2
01261   //  std::cout << "                                 ierr : " << mpierr <<std::endl;
01262   std::cout << "              PROC " << myrank << ": DONE POST RECV `sd" << _recvs[currRecv].sbdId
01263             << "' AT STEP `" << step << "' FROM PROC `" << fromId  << "' on buffer `" << _recvs[currRecv].step << "  with tags : " 
01264             << (step * D_NB_REQUESTS) + DCONTRIB_INFO          << " " 
01265             << (step * D_NB_REQUESTS) + DCONTRIB_ORDER         << " " 
01266             << (step * D_NB_REQUESTS) + DCONTRIB_OVRLP         << " " 
01267             << (step * D_NB_REQUESTS) + DCONTRIB_DENS          << " "  <<  std::endl
01268             << " Sizes of requests : " << D_INFO_SIZE << "  " << maxNbAtom * ORDR_CELL_SIZE 
01269             << "   " << maxOvrlpSize <<"  "<< (maxNbAO * (maxNbAO+1))/2 
01270             <<  std::endl;
01271 #endif
01272   //
01273   QC_TRACE_DENSITY("END   QCDistMDSystem<TPSolver>::postDensityRecvs NEW");
01274 }
01279 template <class TPSolver>
01280 template <class TPParam>
01281 QCFloat
01282 QCDistMDSystem<TPSolver>::applyOptimalDamping (const TPParam * QCRestrict params) {
01283 
01284   QC_TRACE_SCF("BEGIN QCDistMDSystem<TPSolver>::applyOptimalDamping");
01285 
01286   //Variables pour l'Optimal damping.
01287   QCFloat localTraces[4], globalTraces[4];
01288 
01289   QCFloat traceFP1_PP1OfSD, traceF_PP1OfSD;
01290   QCFloat traceFP1_POfSD, traceF_POfSD;
01291   
01292   QCFloat lambdaOptDamp;
01293 
01294 
01295   traceFP1_PP1OfSD = QC_ZERO;
01296   traceF_PP1OfSD   = QC_ZERO;
01297   traceFP1_POfSD   = QC_ZERO;
01298   traceF_POfSD     = QC_ZERO;
01299 
01300   for (int i=0; i < QCMDSystem<TPSolver>::nbDomains; ++i) {
01301     
01302     traceFP1_PP1OfSD +=  QCMDSystem<TPSolver>::domains[i].traceProductCoreShell1(params,
01303                                         QCMDSystem<TPSolver>::domains[i].getFockF().getMatrix(),
01304                                         QCMDSystem<TPSolver>::domains[i].getDensityP().getMatrix());
01305     
01306     traceF_PP1OfSD   +=  QCMDSystem<TPSolver>::domains[i].traceProductCoreShell1(params,
01307                                         QCMDSystem<TPSolver>::domains[i].getFockFTild()->getMatrix(),
01308                                         QCMDSystem<TPSolver>::domains[i].getDensityP().getMatrix());
01309     
01310     traceFP1_POfSD   +=  QCMDSystem<TPSolver>::domains[i].traceProductCoreShell1(params,
01311                                         QCMDSystem<TPSolver>::domains[i].getFockF().getMatrix(),
01312                                         QCMDSystem<TPSolver>::domains[i].getDensityPTild()->getMatrix());
01313     
01314     traceF_POfSD     +=  QCMDSystem<TPSolver>::domains[i].traceProductCoreShell1(params,
01315                                         QCMDSystem<TPSolver>::domains[i].getFockFTild()->getMatrix(),
01316                                         QCMDSystem<TPSolver>::domains[i].getDensityPTild()->getMatrix());
01317   }
01318   localTraces[0] = traceFP1_PP1OfSD;
01319   localTraces[1] = traceF_PP1OfSD;
01320   localTraces[2] = traceFP1_POfSD;
01321   localTraces[3] = traceF_POfSD;
01322   //
01323   //
01324   MPI_Allreduce(localTraces, globalTraces, 4, MPI_DOUBLE, MPI_SUM, QC_COMM_WORLD);
01325   
01326   lambdaOptDamp = QCMDSystem<TPSolver>::calculateLambdaOptDamp(globalTraces[0], globalTraces[1],
01327                                                                globalTraces[2], globalTraces[3]);
01328   
01329   QC_TRACE_SCF("            lambdaOptDamp : "<< lambdaOptDamp );
01330   for (int i=0; i < QCMDSystem<TPSolver>::nbDomains; ++i) {
01331     QCMDSystem<TPSolver>::domains[i].modifyMatricesAndEnergy(lambdaOptDamp); 
01332     
01333 #if defined (QC_TRACE_IN_FILE)  
01334     qctrace[i].writeMsg(" We apply Optimal Damping "); 
01335     qctrace[i].writeVal(" lambdaOptDamping : ",lambdaOptDamp); 
01336     qctrace[i].writeFockMatrix() ;
01337 #endif
01338   }
01339   QC_TRACE_SCF("END   QCDistMDSystem<TPSolver>::applyOptimalDamping");
01340   return lambdaOptDamp ;
01341 }
01342 
01346 template <class TPSolver>
01347 bool QCDistMDSystem<TPSolver>::testConvergence (QCFloat & error, QCFloat& threshold) {
01348   //  
01349   QC_TRACE_SCF("BEGIN QCDistMDSystem<TPSolver>::testConvergence");
01350   QCFloat localError; 
01351   bool convReached = QCMDSystem<TPSolver>::testConvergence(localError, threshold);
01352   //
01353   MPI_Allreduce(&localError, &error, 1, MPI_DOUBLE, MPI_MAX, QC_COMM_WORLD);
01354   //  
01355   convReached = convReached && (error <= threshold);
01356   QC_TRACE_SCF("END   QCDistMDSystem<TPSolver>::testConvergence");
01357 
01358   return convReached ;
01359 }
01360 
01364 template <class TPSolver>
01365 template <class TPParam>
01366 void
01367 QCDistMDSystem<TPSolver>::mergeDensityAux (const TPParam * QCRestrict params,
01368                                            QCSubDomain&    QCRestrict contributer,
01369                                            QCSymMatrix&    QCRestrict remWeightedDensity,
01370                                            QCMatrix&       QCRestrict wDensityAB,
01371                                            QCMatrix&       QCRestrict interPAB,
01372                                            QCSymMatrix&    QCRestrict wDensityAA) {
01373 
01374   QC_TRACE_DENSITY("BEGIN QCDistMDSystem<TPSolver>::mergeDensityAux ");
01375   QCDensityGtr<QCSymMatrix> * QCRestrict densityP;
01376 
01377   for (int i=0; i < QCMDSystem<TPSolver>::nbDomains; ++i) {
01378     
01379     densityP = &QCMDSystem<TPSolver>::domains[i].getDensityP();    
01380 #ifdef QC_DEBUG_DENSITY_CONTRIB_LEVEL2
01381     std::cout << "   BEGIN Contribution for domain "   << QCMDSystem<TPSolver>::domains[i].getId()   
01382               << " (local number " << i    << std::endl;
01383 #endif
01384     for (int j=0; j < QCMDSystem<TPSolver>::domains[i].getNbOverlpSD(); ++j) {
01385       
01386       int rId = QCMDSystem<TPSolver>::domains[i].getRemDomainId(j);
01387       if (rId == contributer.getId()) {
01388         
01389 #ifdef QC_DEBUG_DENSITY_CONTRIB_LEVEL2
01390         std::cout << "PROC " << myrank << ": DENSITY CONTRIB: SD `"  << QCMDSystem<TPSolver>::domains[i].getId()
01391                   << "' <= SD `" << rId  << "'" << std::endl;
01392 #endif
01393         
01394 #ifdef QC_DEBUG_DENSITY_CONTRIB_LEVEL2
01395         std::cout << "               Contribution from  sd " << QCMDSystem<TPSolver>::domains[i].getRemDomainId(j)   << std::endl;
01396 #endif
01397         densityP->mergeDensity(params, QCMDSystem<TPSolver>::domains[i],contributer, remWeightedDensity,
01398                                j, wDensityAB, interPAB, wDensityAA);
01399       }
01400     }
01401 #ifdef QC_DEBUG_DENSITY_CONTRIB_LEVEL2
01402     std::cout << "   END Contribution for domain "   << QCMDSystem<TPSolver>::domains[i].getId()   
01403               << " (local number " << i    << std::endl;
01404 #endif
01405   }
01406   QC_TRACE_DENSITY("END   QCDistMDSystem<TPSolver>::mergeDensityAux ");
01407 }
01408 
01412 template <class TPSolver>
01413 template <class TPParam>
01414 void
01415 QCDistMDSystem<TPSolver>::mergeDensityAux (const TPParam * QCRestrict params,
01416                                            QCSubDomain&    QCRestrict contributer1,
01417                                            QCSubDomain&    QCRestrict contributer2,
01418                                            QCSymMatrix&    QCRestrict remWeightedDensity,
01419                                            QCMatrix&       QCRestrict wDensityAB,
01420                                            QCMatrix&       QCRestrict interPAB,
01421                                            QCSymMatrix&    QCRestrict wDensityAA) {
01422 
01423   QC_TRACE_DENSITY("BEGIN QCDistMDSystem<TPSolver>::mergeDensityAux version 2");
01424   QCDensityGtr<QCSymMatrix> * QCRestrict densityP; 
01425   
01426   densityP = &contributer1.getDensityP();    
01427 #ifdef QC_DEBUG_DENSITY_CONTRIB_LEVEL2
01428   std::cout << "   BEGIN Contribution for domain "   << contributer1.getId()   <<std::endl;
01429 #endif
01430   for (int j=0; j <  contributer1.getNbOverlpSD(); ++j) {
01431     
01432     int rId =  contributer1.getRemDomainId(j);
01433     if (rId == contributer2.getId()) {
01434       
01435 #ifdef QC_DEBUG_DENSITY_CONTRIB_LEVEL2
01436       std::cout << "PROC " << myrank << ": DENSITY CONTRIB: SD `"  <<  contributer1.getId()
01437                 << "' <= SD `" << rId  << "'" << std::endl;
01438 #endif
01439       
01440 #ifdef QC_DEBUG_DENSITY_CONTRIB_LEVEL2
01441       std::cout << "               Contribution from  sd " <<  contributer1.getRemDomainId(j)   << std::endl;
01442 #endif
01443       densityP->mergeDensity(params,  contributer1,contributer2, remWeightedDensity,j, 
01444                              wDensityAB, interPAB, wDensityAA);
01445     }
01446   }
01447 #ifdef QC_DEBUG_DENSITY_CONTRIB_LEVEL2
01448   std::cout << "   END Contribution for domain "   << contributer1.getId()   << std::endl;
01449 #endif
01450   QC_TRACE_DENSITY("END   QCDistMDSystem<TPSolver>::mergeDensityAux version 2");
01451  }
01452 
01453 
01454 
01455 /************************************************
01456 *************************************************
01457 *               MPI Operation                   *
01458 *************************************************
01459 *************************************************/
01460 void compute_min_max (QCFloatFloat * in, QCFloatFloat * inout,
01461                  int          * len,  MPI_Datatype * /* dptr */) {
01462 
01463   int i;
01464   
01465   for (i = 0; i < *len; ++i, ++in, ++inout) {
01466     if (inout->frstval > in->frstval) {
01467       inout->frstval = in->frstval;
01468     }
01469     if (inout->scndval < in->scndval) {
01470       inout->scndval = in->scndval;
01471     }
01472   }
01473 }
01477 template <class TPSolver>
01478 void
01479 QCDistMDSystem<TPSolver>::getGlobalMinMax(QCFloat& min, QCFloat& max) {
01480 
01481   QCFloatFloat in, out;
01482 
01483   in.frstval = min;
01484   in.scndval = max;
01485   
01486   MPI_Allreduce(&in, &out, 1, QC_2DOUBLE, QC_MIN_MAX, QC_COMM_WORLD);
01487 
01488   min = out.frstval;
01489   max = out.scndval;
01490 }
01491 
01495 template <class TPSolver>
01496 void QCDistMDSystem<TPSolver>::getTotalSumNbElecs (QCFloat * nbelec, int size) {
01497   QCFloat * totalNbElec = new QCFloat [size];
01498   MPI_Allreduce(nbelec, totalNbElec, size, MPI_DOUBLE, MPI_SUM, QC_COMM_WORLD);
01499 
01500   for (int i = 0; i < size; ++i) {
01501     nbelec[i] = totalNbElec[i];
01502   }  
01503   delete [] totalNbElec;
01504 }
01508 template <class TPSolver>
01509 QCFloat
01510 QCDistMDSystem<TPSolver>::getTotalSumTrace (void) {
01511   QCFloat sumTrace = 0;
01512   QCFloat totalSumTrace;
01513 
01514   for (int i=0; i < QCMDSystem<TPSolver>::nbDomains; ++i) {
01515     sumTrace += QCMDSystem<TPSolver>::domains[i].getWeightedDensityP()->trace();
01516   }
01517 
01518   MPI_Reduce(&sumTrace, &totalSumTrace, 1,
01519              MPI_DOUBLE, MPI_SUM, 0, QC_COMM_WORLD);
01520   
01521   return totalSumTrace;
01522 }
01523 /********************************************************************/
01524 /********************************************************************/
01528 template <class TPSolver>
01529 void 
01530 QCDistMDSystem<TPSolver>::initFermiSeq (void) {
01531 
01532   int numSD;
01533   totalAO = 0;
01534 
01535   nbOM = new int [QCMDSystem<TPSolver>::nbDomains];
01536 
01537   for (numSD = 0; numSD < QCMDSystem<TPSolver>::nbDomains; ++numSD) {
01538     totalAO += QCMDSystem<TPSolver>::domains[numSD].getNbAtomicOrbitals();
01539   }
01540 
01541   sbuffer = new QCFloatFloat [totalAO];
01542   MPI_Type_contiguous(totalAO, QC_2DOUBLE, &stype);
01543   MPI_Type_commit(&stype);
01544   if (myrank == 0) {
01545     rcounts = new int [nproc];
01546     displs  = new int [nproc];
01547   }
01548 
01549   MPI_Gather(&totalAO, 1, MPI_INT, rcounts, 1, MPI_INT, 0, QC_COMM_WORLD);
01550   
01551   if (myrank == 0) {
01552     _totalAOAll = 0;
01553     for (int i = 0; i < nproc; ++i) {
01554       displs[i]   = _totalAOAll;  _totalAOAll += rcounts[i];
01555     }
01556     if(rbuffer != NULL)
01557       {delete [] rbuffer ;}
01558     rbuffer = new QCFloatFloat [_totalAOAll];
01559   }
01560 }
01564 template <class TPSolver>
01565 template <class TPMDSystem>
01566 void 
01567 QCDistMDSystem<TPSolver>::adjustFermiEnergySeq (TPMDSystem& /* mdsystem */,
01568                                                    QCMemory&   /* memory */) {
01569   //
01570   QC_TRACE_ENER("BEGIN QCDistMDSystem<TPSolver>::adjustFermiEnergySeq " );
01571   int numOM, numSD, gnumOM;
01572   QCFloat * eigenVal;
01573   QCFloat * bFactor;
01574   QCFloat * orbitalOccupN;
01575   QCFloat   sumNE, halfNBE = this->nbElectrons*QC_HALF;
01576   QCFloatFloat result;
01577 
01578 #ifdef QC_VERBOSE
01579   std::cout << "* ADJUST FERMI ENERGY (NEW SEQ VERSION) *" << std::endl;
01580 #endif
01581 
01582   // Copie dans le buffer d'envoie
01583   gnumOM = 0;
01584   for (numSD = 0; numSD < QCMDSystem<TPSolver>::nbDomains; ++numSD) {
01585     
01586     eigenVal      = QCMDSystem<TPSolver>::domains[numSD].getEigenVal();
01587     bFactor       = QCMDSystem<TPSolver>::domains[numSD].getBFactor();
01588     orbitalOccupN = QCMDSystem<TPSolver>::domains[numSD].getOrbitalOccupN();
01589     nbOM[numSD]   = QCMDSystem<TPSolver>::domains[numSD].getNbAtomicOrbitals();
01590 
01591     // On met le vecteur d'occupation a 0
01592     memset(orbitalOccupN, 0, nbOM[numSD] * sizeof(QCFloat));
01593 
01594     for (numOM = 0; numOM < nbOM[numSD]; ++numOM, ++gnumOM) {
01595       sbuffer[gnumOM].frstval = eigenVal[numOM];
01596       sbuffer[gnumOM].scndval = bFactor[numOM];
01597     }
01598   }
01599   MPI_Gatherv (sbuffer, 1, stype, rbuffer, rcounts, displs, QC_2DOUBLE, 0, QC_COMM_WORLD);
01600 
01601   if (myrank == 0) {
01602     stable_sort(rbuffer, rbuffer+_totalAOAll);
01603     sumNE = QC_ZERO;
01604     for (numOM = 0; numOM < _totalAOAll && sumNE < halfNBE; ++numOM) {
01605       sumNE += rbuffer[numOM].scndval; 
01606     }  
01607      
01608     result.frstval = rbuffer[numOM-1].frstval;
01609     result.scndval = sumNE * QC_TWO;
01610   } 
01611 
01612   MPI_Bcast(&result, 1, QC_2DOUBLE, 0, QC_COMM_WORLD); 
01613 
01614   this->fermiEnergy = result.frstval;
01615   sumNE             = result.scndval;
01616 
01617 #ifdef QC_VERBOSE_FERMI
01618   std::cout << "PROC " << myrank << ": Fermi Energy = " << this->fermiEnergy  << ", sumNE = " << sumNE << std::endl;
01619 #endif 
01620 
01621   // Set Orbital occupation
01622   for (numSD = 0; numSD < QCMDSystem<TPSolver>::nbDomains; ++numSD) {
01623     eigenVal      = QCMDSystem<TPSolver>::domains[numSD].getEigenVal();
01624     bFactor       = QCMDSystem<TPSolver>::domains[numSD].getBFactor();
01625     orbitalOccupN = QCMDSystem<TPSolver>::domains[numSD].getOrbitalOccupN();
01626     nbOM[numSD]   = QCMDSystem<TPSolver>::domains[numSD].getNbAtomicOrbitals();
01627     
01628     for (numOM = 0; numOM < nbOM[numSD] && eigenVal[numOM] < this->fermiEnergy; ++numOM) {
01629       orbitalOccupN[numOM] = QC_TWO;
01630     }
01631     QCMDSystem<TPSolver>::domains[numSD].setNbFilledOccupiedOM(numOM, numOM);
01632     
01633     if (numOM < nbOM[numSD] && eigenVal[numOM] == this->fermiEnergy) {
01634       orbitalOccupN[numOM] = QC_TWO + (this->nbElectrons - sumNE) / bFactor[numOM];
01635       QCMDSystem<TPSolver>::domains[numSD].setNbFilledOccupiedOM(numOM, numOM+1);
01636     }
01637   }
01638   QC_TRACE_ENER("END  QCDistMDSystem<TPSolver>::adjustFermiEnergySeq " );
01639 }
01640 
01641 
01642 
01643 /********************************************************************/
01644 /********************************************************************/
01645 
01646 
01647 
01648 
01652 template <class TPSolver>
01653 void
01654 QCDistMDSystem<TPSolver>::removeLevels (QCFloat& sumNbElecs, QCTopEnergyLevel * topEnergyLevels) {
01655 
01656   QC_TRACE_ENER("BEGIN QCDistMDSystem<TPSolver>::removeLevels -- Adjust Fermi Energy on proc "<<myrank);
01657   QC_TRACE_ENER("         nbDomains " << QCMDSystem<TPSolver>::nbDomains  );
01658 
01659   // On change toujours les elements de la derniere case
01660   // car tri par ordre croissant
01661   int& numSDToChange = topEnergyLevels[QCMDSystem<TPSolver>::nbDomains-1].numSD;
01662   int& numOMToChange = topEnergyLevels[QCMDSystem<TPSolver>::nbDomains-1].numOM;
01663   QCFloat& energyLevelToChange = topEnergyLevels[QCMDSystem<TPSolver>::nbDomains-1].energyLevel;
01664   int numSD, numOfHomoOfSD;
01665   
01666   // MPI vars
01667   QCFloatInt localpair, maxpair, resultpair;
01668 
01669   bool rightNumberOfElec = false;
01670   
01671     // On met dans topEnergyLevels les energies des homos.
01672   for (numSD = 0; numSD < QCMDSystem<TPSolver>::nbDomains; ++numSD) {
01673     numOfHomoOfSD = topEnergyLevels[numSD].numOM;
01674     if (numOfHomoOfSD >= 0) {
01675       topEnergyLevels[numSD].energyLevel = QCMDSystem<TPSolver>::domains[numSD].getEigenVal(numOfHomoOfSD);
01676     }
01677     else {
01678       // On ne peut plus enlever d'electrons sur ce sous-domaine.
01679       topEnergyLevels[numSD].energyLevel = -DBL_MAX;
01680     }
01681   }
01682   //
01683   int loop = 0, choice = -1;
01684   while (!rightNumberOfElec) {
01685   
01686     // On trie les niveaux par ordre croissant.
01687     stable_sort(topEnergyLevels, topEnergyLevels + QCMDSystem<TPSolver>::nbDomains);
01688 
01689     localpair.floatval = energyLevelToChange;
01690     localpair.intval   = myrank;
01691     
01692     MPI_Allreduce(&localpair, &maxpair, 1, MPI_DOUBLE_INT, MPI_MAXLOC, QC_COMM_WORLD);
01693 
01694     if (maxpair.intval == myrank) {
01695       
01696       // On enleve le nb de du niveau le + haut au niveau au total.
01697       if ( (sumNbElecs - QC_TWO * QCMDSystem<TPSolver>::domains[numSDToChange].getBFactor(numOMToChange)) >=
01698            this->nbElectrons) {
01699         choice = 0 ;
01700         sumNbElecs -= QC_TWO * QCMDSystem<TPSolver>::domains[numSDToChange].getBFactor(numOMToChange);
01701         QCMDSystem<TPSolver>::domains[numSDToChange].setOrbitalOccupN(numOMToChange, QC_ZERO);
01702         
01703         // On le remplace par le niveau le plus haut du meme sous domaine.
01704         --numOMToChange;
01705         if (numOMToChange >= 0) {
01706           choice = 1 ;
01707           energyLevelToChange = QCMDSystem<TPSolver>::domains[numSDToChange].getEigenVal(numOMToChange);
01708           
01709         } else {
01710           // On ne peut plus enlever d electrons sur ce sous-domaine.
01711           choice = 2 ;
01712           energyLevelToChange = -DBL_MAX;
01713         }
01714         
01715         
01716       } else {
01717         // On diminue le nombre d'occupation de l'homo.
01718         choice =3 ;
01719         QCMDSystem<TPSolver>::domains[numSDToChange].
01720           setOrbitalOccupN(numOMToChange, 
01721                            QC_TWO +
01722                            (this->nbElectrons - sumNbElecs) / QCMDSystem<TPSolver>::domains[numSDToChange].getBFactor(numOMToChange));
01723         rightNumberOfElec = true;
01724       }  // 
01725     }
01726   
01727     resultpair.floatval = sumNbElecs;
01728     resultpair.intval   = rightNumberOfElec;
01729     QC_TRACE_ENER("    root du bcast  : " << resultpair.intval <<"  sumNbElecs " << sumNbElecs 
01730                   <<"  maxpair.floatval " <<  maxpair.floatval  );
01731 
01732     MPI_Bcast(&resultpair, 1, MPI_DOUBLE_INT, maxpair.intval, QC_COMM_WORLD);
01733     
01734     if (maxpair.intval != myrank) {
01735       sumNbElecs = resultpair.floatval;
01736       rightNumberOfElec = static_cast<bool>(resultpair.intval);
01737     }
01738     QC_TRACE_ENER("     loop : " << loop <<"        choice " << choice  <<"  sumNbElecs " << sumNbElecs 
01739                   <<"  maxpair.floatval " <<  maxpair.floatval );
01740     
01741     ++loop;
01742   }
01743   //
01744   // On repositionne l'energie de Fermi. Ici l'energie de Fermi est l'energie 
01745   // de la Homo
01746   this->fermiEnergy = maxpair.floatval;
01747   //
01748   QC_TRACE_ENER("        (" << loop << " loops) fermi energy = " << this->fermiEnergy );
01749   //
01750   // On affecte nbFilledOM et nbOccupiedOM du domaine
01751   // Dans topEnergyLevels, on a les homos.
01752   for (numSD = 0; numSD < QCMDSystem<TPSolver>::nbDomains; ++numSD) {
01753     
01754     // Le nb d'orbitale occupe est egal au nombre d'orbitales remplies.
01755     
01756     // Le sous-domaine est vide.
01757     if (topEnergyLevels[numSD].energyLevel == -DBL_MAX) {
01758       QCMDSystem<TPSolver>::domains[topEnergyLevels[numSD].numSD].setNbFilledOccupiedOM(0, 0);
01759     }
01760     // L'homo du SD est remplie.
01761     else if (QCMDSystem<TPSolver>::domains[topEnergyLevels[numSD].numSD].getOrbitalOccupN(topEnergyLevels[numSD].numOM)
01762              == QC_TWO) {
01763       
01764       // La derniere orbitale est completement remplie.
01765       QCMDSystem<TPSolver>::domains[topEnergyLevels[numSD].numSD].setNbFilledOccupiedOM(topEnergyLevels[numSD].numOM + 1,
01766                                                                                         topEnergyLevels[numSD].numOM + 1);
01767     }
01768     
01769     // L'homo du SD est occupee.
01770     else {
01771 
01772       // La derniere orbitale n'est pas completement remplie.
01773       QCMDSystem<TPSolver>::domains[topEnergyLevels[numSD].numSD].
01774         setNbFilledOccupiedOM(topEnergyLevels[numSD].numOM,
01775                               topEnergyLevels[numSD].numOM + 1);
01776     }
01777   }
01778   QC_TRACE_ENER("END   QCDistMDSystem<TPSolver>::removeLevels -- Adjust Fermi Energy");
01779 }
01780 
01784 template <class TPSolver>
01785 void
01786 QCDistMDSystem<TPSolver>::addLevels (QCFloat& sumNbElecs, const int * nbOM,
01787                                      QCTopEnergyLevel * topEnergyLevels) {
01788 
01789   QC_TRACE_ENER("BEGIN QCDistMDSystem<TPSolver>::addLevels on proc " <<myrank);  
01790 
01791   // On change toujours les elements de la premiere case
01792   // car tri par ordre croissant
01793   int& numSDToChange = topEnergyLevels[0].numSD;
01794   int& numOMToChange = topEnergyLevels[0].numOM;
01795   QCFloat& energyLevelToChange = topEnergyLevels[0].energyLevel;
01796   
01797   // MPI vars
01798   QCFloatInt localpair, minpair, resultpair;
01799 
01800 
01801   int numSD;
01802   int numOfLumoOfSD;
01803 
01804   bool rightNumberOfElec = false;
01805 
01806   // On met dans topEnergyLevels les energies des lumos.
01807   for (numSD = 0; numSD < QCMDSystem<TPSolver>::nbDomains; ++numSD) {
01808     numOfLumoOfSD = topEnergyLevels[numSD].numOM + 1;
01809     if (numOfLumoOfSD < nbOM[numSD]) {
01810       topEnergyLevels[numSD].numOM = numOfLumoOfSD;
01811       topEnergyLevels[numSD].energyLevel = QCMDSystem<TPSolver>::domains[numSD].getEigenVal(numOfLumoOfSD);
01812       
01813     } else {
01814       // On ne peut plus enlever d'electrons sur ce sous-domaine.
01815       topEnergyLevels[numSD].energyLevel = DBL_MAX;
01816     }
01817   }
01818   
01819   int loop = 0;
01820   while (!rightNumberOfElec) {
01821     
01822     // On trie les niveaux par ordre croissant.
01823     stable_sort(topEnergyLevels, topEnergyLevels + QCMDSystem<TPSolver>::nbDomains);
01824     
01825     localpair.floatval = energyLevelToChange;
01826     localpair.intval   = myrank;
01827     
01828     MPI_Allreduce(&localpair, &minpair, 1, MPI_DOUBLE_INT, MPI_MINLOC, QC_COMM_WORLD);
01829 
01830     if (minpair.intval == myrank) {
01831 
01832       // On ajoute le nb d'e du niveau le + haut au niveau au total.
01833       if ( (sumNbElecs + QC_TWO * QCMDSystem<TPSolver>::domains[numSDToChange].getBFactor(numOMToChange)) <=
01834            this->nbElectrons) {
01835         
01836         sumNbElecs += QC_TWO * QCMDSystem<TPSolver>::domains[numSDToChange].getBFactor(numOMToChange);
01837         QCMDSystem<TPSolver>::domains[numSDToChange].setOrbitalOccupN(numOMToChange, QC_TWO);
01838         
01839         // On le remplace par le niveau le plus haut du meme sous domaine.
01840         numOMToChange++;
01841         if (numOMToChange < nbOM[numSDToChange]) {
01842           energyLevelToChange = QCMDSystem<TPSolver>::domains[numSDToChange].getEigenVal(numOMToChange);
01843           
01844         } else {
01845           // On ne peut plus ajouter d electrons sur ce sous-domaine.
01846           energyLevelToChange = DBL_MAX;
01847         }
01848         
01849 
01850       } else {
01851         // On augmente le nombre d'occupation de la lumo.
01852         QCMDSystem<TPSolver>::domains[numSDToChange].
01853           setOrbitalOccupN(numOMToChange, 
01854                            (this->nbElectrons - sumNbElecs) / QCMDSystem<TPSolver>::domains[numSDToChange].getBFactor(numOMToChange));
01855         
01856         // La lumo est la suivante
01857         numOMToChange++;
01858         rightNumberOfElec = true;
01859       }
01860 
01861     }
01862 
01863     resultpair.floatval = sumNbElecs;
01864     resultpair.intval   = rightNumberOfElec;
01865 
01866     MPI_Bcast(&resultpair, 1, MPI_DOUBLE_INT, minpair.intval, QC_COMM_WORLD);
01867     
01868     if (minpair.intval != myrank) {
01869       sumNbElecs = resultpair.floatval;
01870       rightNumberOfElec = static_cast<bool>(resultpair.intval);
01871     }
01872 
01873     ++loop;
01874   }
01875 
01876 
01877   // On repositionne l'energie de Fermi. Ici l'energie de Fermi est l'energie 
01878   // de la Lumo
01879   this->fermiEnergy = minpair.floatval;
01880   
01881   QC_TRACE_ENER("      proc "<<myrank<<"  (" << loop  << " loops) fermi energy = " << this->fermiEnergy );  
01882   // On affecte nbFilledOM et nbOccupiedOM du domaine
01883 
01884   // Dans topEnergyLevels, on a les lumos.
01885   for (numSD = 0; numSD < QCMDSystem<TPSolver>::nbDomains; ++numSD) {
01886 
01887     // Le nb d'orbitale occupe est egal au nombre d'orbitales remplies.
01888 
01889     // Le sous-domaine est rempli
01890     if (topEnergyLevels[numSD].energyLevel == DBL_MAX) {
01891       QCMDSystem<TPSolver>::domains[topEnergyLevels[numSD].numSD].
01892         setNbFilledOccupiedOM(nbOM[topEnergyLevels[numSD].numSD], 
01893                               nbOM[topEnergyLevels[numSD].numSD]);
01894     }
01895 
01896     // Le sous-domaine est vide
01897     else if (QCMDSystem<TPSolver>::domains[topEnergyLevels[numSD].numSD].getOrbitalOccupN(0)
01898              == QC_ZERO) {
01899       
01900       QCMDSystem<TPSolver>::domains[topEnergyLevels[numSD].numSD].
01901         setNbFilledOccupiedOM(0, 0);
01902     }
01903     
01904     // La lumo du SD est remplie
01905     else if (QCMDSystem<TPSolver>::domains[topEnergyLevels[numSD].numSD].getOrbitalOccupN(topEnergyLevels[numSD].numOM-1)
01906              == QC_TWO) {
01907 
01908       QCMDSystem<TPSolver>::domains[topEnergyLevels[numSD].numSD].
01909         setNbFilledOccupiedOM(topEnergyLevels[numSD].numOM,
01910                               topEnergyLevels[numSD].numOM);
01911     }
01912     
01913     // La lumo du SD est juste occupee
01914     else {
01915 
01916       // La derniere orbitale n'est pas completement remplie.
01917       QCMDSystem<TPSolver>::domains[topEnergyLevels[numSD].numSD].
01918         setNbFilledOccupiedOM(topEnergyLevels[numSD].numOM - 1,
01919                               topEnergyLevels[numSD].numOM);
01920     }
01921   }
01922   QC_TRACE_ENER("END   QCDistMDSystem<TPSolver>::addLevels on proc "<<myrank );    
01923 }
01924 
01928 template <class TPSolver>
01929 void
01930 QCDistMDSystem<TPSolver>::allocBuffers (bool density) {
01931   
01932   // Le pipeline
01933   pipeline  = new QCSubDomain [QCMDSystem<TPSolver>::nbDomains];
01934   
01935   atomsBuffers = new QCFloat * [QCMDSystem<TPSolver>::nbDomains + 1];
01936   orderBuffers = new int     * [QCMDSystem<TPSolver>::nbDomains + 1];
01937   ovrlpBuffers = new int     * [QCMDSystem<TPSolver>::nbDomains + 1];
01938   zoneBuffers  = new QCSubDomainZone     * [QCMDSystem<TPSolver>::nbDomains + 1];
01939   
01940   if (density) {
01941     blockBuffers = new int     * [QCMDSystem<TPSolver>::nbDomains + 1];
01942     densBuffers  = new QCFloat * [QCMDSystem<TPSolver>::nbDomains + 1];
01943   }
01944 
01945   for (int i=0; i <= QCMDSystem<TPSolver>::nbDomains; ++i) {
01946     atomsBuffers[i] = new QCFloat [maxNbAtom * QCAtoms::DIMENSION + ATOM_OFFSET];
01947     orderBuffers[i] = new int     [maxNbAtom * ORDR_CELL_SIZE];
01948     ovrlpBuffers[i] = new int     [maxOvrlpSize];
01949     zoneBuffers [i] = new QCSubDomainZone [maxNbAtom];
01950     
01951     if (density) {
01952       blockBuffers[i] = new int     [maxNbAtom];
01953       densBuffers [i] = new QCFloat [maxDensitySize];
01954       if (i < 2) {
01955         wdensMatBuffs[i] = new QCFloat [(maxNbAO * (maxNbAO + 1))/2 + 1];
01956       }
01957     }
01958   }
01959   
01960 }
01964 template <class TPSolver>
01965 void
01966 QCDistMDSystem<TPSolver>::freeBuffers (void) {
01967   QC_TRACE_END("BEGIN QCDistMDSystem<TPSolver>::freeBuffers");
01968 
01969   if (pipeline) {
01970     for (int i=0; i < QCMDSystem<TPSolver>::nbDomains; ++i) {
01971       pipeline[i].unsetAtomOverlp();
01972     }
01973     delete [] pipeline;
01974     pipeline = NULL;
01975   }
01976 
01977   for (int i=0; i <= QCMDSystem<TPSolver>::nbDomains; ++i) {
01978     if (atomsBuffers && atomsBuffers[i]) {
01979       delete [] atomsBuffers[i] ;
01980       atomsBuffers[i] = NULL;
01981     }
01982     if (orderBuffers && orderBuffers[i]) {
01983       delete [] orderBuffers[i] ;
01984       orderBuffers[i] = NULL;
01985     }
01986     if (ovrlpBuffers && ovrlpBuffers[i]) {
01987       delete [] ovrlpBuffers[i] ;
01988       ovrlpBuffers[i] = NULL;
01989     }
01990     if (zoneBuffers && zoneBuffers[i]) {
01991       delete [] zoneBuffers [i];
01992       zoneBuffers [i] = NULL;
01993     }
01994     if (blockBuffers && blockBuffers[i]) {
01995       delete [] blockBuffers[i];
01996       blockBuffers[i] = NULL;
01997     }
01998     if (densBuffers && densBuffers[i]) {
01999       delete [] densBuffers[i];
02000       densBuffers[i] = NULL;
02001     }
02002     if (i < 2 && wdensMatBuffs[i]) {
02003       delete [] wdensMatBuffs[i];
02004       wdensMatBuffs[i] = NULL;
02005     }
02006   }
02007 
02008   if (atomsBuffers) {
02009     delete [] atomsBuffers; atomsBuffers = NULL;
02010   }
02011   if (orderBuffers) {
02012     delete [] orderBuffers; orderBuffers = NULL;
02013   }
02014   if (ovrlpBuffers) {
02015     delete [] ovrlpBuffers; ovrlpBuffers = NULL;
02016   }
02017   if (zoneBuffers) {
02018     delete [] zoneBuffers; zoneBuffers = NULL;
02019   }
02020   if (blockBuffers) {
02021     delete [] blockBuffers; blockBuffers = NULL;
02022   }
02023   if (densBuffers) {
02024     delete [] densBuffers; densBuffers = NULL;
02025   }
02026 
02027   QC_TRACE_END("END   QCDistMDSystem<TPSolver>::freeBuffers");
02028 }
02032 template <class TPSolver>
02033 void  
02034 QCDistMDSystem<TPSolver>::fillCommunicationStruct(const int &maxOvrlpSDs) {
02035  
02036   // Allocation des buffers
02037   this->fillSDsCommArray();  // A fussionner avec fillSDsCommVect pour supprimer getProcRank(int&)
02038   // Echange globale des listes de reception
02039   //
02040   MPI_Alltoall(sdsFromProc, (maxOvrlpSDs + SDS_HEAD), MPI_INT,
02041                sdsIntoProc, (maxOvrlpSDs + SDS_HEAD), MPI_INT, QC_COMM_WORLD);
02042   //
02043   // Pour plus de facilite d'utilisation on construit
02044   // deux vecteurs correspondant aux echanges de messages 
02045   // a effectuer
02046   this->fillSDsCommVect();
02047 
02048 #ifdef QC_DEBUG_DENSITY_SR
02049   printProcArray(sdsFromProc, false);
02050   printProcArray(sdsIntoProc, true);
02051   std::cout << "                                RECVS                " << std::endl;
02052   printProcVect(_recvs, false);
02053   std::cout << "                                SEND                " << std::endl;
02054   printProcVect(_sends, true);
02055   std::cout << std::endl;
02056 #endif
02057 
02058 }
02062 template <class TPSolver>
02063 void  
02064 QCDistMDSystem<TPSolver>::fillSDsCommArray (void) {
02065 
02066   QC_TRACE_MPI("BEGIN QCDistMDSystem<TPSolver>::fillSDsCommArray ");
02067   sdsFromProc = new int [nproc * (maxOvrlpSDs + SDS_HEAD)];
02068   sdsIntoProc = new int [nproc * (maxOvrlpSDs + SDS_HEAD)];
02069   for (int i=0; i < nproc* (maxOvrlpSDs + SDS_HEAD); ++i) {
02070     sdsFromProc[i ] = -1;
02071     sdsIntoProc[i ] = -1;
02072   }
02073   for (int i=0; i < nproc; ++i) {
02074     sdsFromProc[i * (maxOvrlpSDs + SDS_HEAD) + NB_SDS] = 0;
02075     sdsIntoProc[i * (maxOvrlpSDs + SDS_HEAD) + NB_SDS] = 0;
02076   }
02077 
02078   
02079   for (int i=0; i < QCMDSystem<TPSolver>::nbDomains; ++i) {
02080 //     std::cout << "proc : "<< myrank <<"   SD  - local number " << i  
02081 //            << "  number of SD overlapping me : " <<  QCMDSystem<TPSolver>::domains[i].getNbOverlpSD() << std::endl;
02082    for (int j = 0 ; j < QCMDSystem<TPSolver>::domains[i].getNbOverlpSD(); ++j) {
02083       bool addfrom   = true;
02084       int  remSdId   = QCMDSystem<TPSolver>::domains[i].getRemDomainId(j);
02085       int  remProc   = this->getProcRank(remSdId);
02086       int  poffset   = remProc * (maxOvrlpSDs + SDS_HEAD);
02087       int& nbSDsFrom = sdsFromProc[poffset + NB_SDS];
02088 //       std::cout <<"                 communique avec " << remSdId  
02089 //              <<"  remProc " <<remProc<< " poffset "<< poffset
02090 //              <<"  nbSDsFrom "<<nbSDsFrom<<std::endl;
02091       for (int k=0; k < nbSDsFrom; ++k) {
02092         if (sdsFromProc[poffset + SDS_HEAD + k] == remSdId) {
02093           addfrom = false;
02094           break;
02095         }
02096       }
02097       if (addfrom) {
02098         sdsFromProc[poffset + SDS_HEAD + nbSDsFrom] = remSdId;
02099         nbSDsFrom++;
02100       }
02101     }
02102   }
02103   for (int i=0; i < nproc; ++i) {
02104     int  poffset   = i * (maxOvrlpSDs + SDS_HEAD);
02105     int  nbSDsFrom = sdsFromProc[poffset + NB_SDS];
02106     //int  nbSDsInto = sdsIntoProc[QCMDSystem<TPSolver>::poffset + NB_SDS];
02107     if (nbSDsFrom > 1) {
02108       stable_sort(sdsFromProc + poffset + SDS_HEAD, sdsFromProc + poffset + SDS_HEAD + nbSDsFrom);
02109     }
02110   }
02111   //
02112   //  Affichage
02113   //
02114 #ifdef TRACE_MPI
02115   std::cout << "#################################################################################" <<std::endl;
02116   std::cout << "############             sdsFromProc                   ##########################" <<std::endl;
02117   std::cout << "#################################################################################" <<std::endl;
02118   for (int i = 0 ; i <  nproc  ; ++i ){
02119     std::cout << "proc : " << i ;
02120     for (int j = 0 ; j<  (maxOvrlpSDs + SDS_HEAD) ; ++j ){
02121     
02122       std::cout << "  "<<sdsFromProc[i*(maxOvrlpSDs + SDS_HEAD)+j] ;
02123     }
02124     std::cout <<std::endl;
02125   }
02126 #endif
02127 
02128 
02129   QC_TRACE_MPI("END   QCDistMDSystem<TPSolver>::fillSDsCommArray ");
02130 
02131 }
02135 template <class TPSolver>
02136 void  
02137 QCDistMDSystem<TPSolver>::fillSDsCommVect (void) {
02138   //
02139   QC_TRACE_MPI("BEGIN QCDistMDSystem<TPSolver>::fillSDsCommVect ");
02140   //
02141   MPI_Datatype QCCOMM, type[2] ;
02142   type[0] = MPI_INT ; type[1] = MPI_DOUBLE ;
02143   int nbtype[2] ; nbtype[0]= 3 ;  nbtype[1]= 1 ; 
02144   MPI_Aint displacement[2] ; 
02145   displacement[0] = 0 ; displacement[1] = 3*sizeof(int) ; 
02146   MPI_Type_struct(2, nbtype, displacement,type, &QCCOMM);
02147   MPI_Type_commit(&QCCOMM);
02148   //
02149   int sizeMax = QCMDSystem<TPSolver>::totalNbDomains ;
02150   QCComm init(-1,-1,0) ;
02151   QCComm lineMat[QCMDSystem<TPSolver>::nbDomains][sizeMax] ;
02152   int sbdOnMyProc[sizeMax] ;  // correspondance entre le numero global et le numero local si le SD est sur le proc
02153   //
02154   for (int j = 0 ; j <  sizeMax  ; ++j ){
02155     sbdOnMyProc[j] = -1 ;
02156     for ( int i = 0; i < QCMDSystem<TPSolver>::nbDomains; ++i) {
02157       lineMat[i][j] = init ;
02158     }
02159   }
02160   //
02161   _nbLocalContrib = 0 ;
02162   for ( int i = 0; i < QCMDSystem<TPSolver>::nbDomains; ++i) {
02163     sbdOnMyProc[this->domains[i].getId()] = i ;
02164   }
02165   int globalIndexOfRemoteSD, myGlobalId ;
02166   for ( int i = 0; i < QCMDSystem<TPSolver>::nbDomains; ++i) {
02167     myGlobalId = this->domains[i].getId() ;
02168 //     std::cout << "proc : "<< myrank <<"   SD  - local number " << i << " global number " << myGlobalId 
02169 //            << "  number of SD overlapping me : " <<  QCMDSystem<TPSolver>::domains[i].getNbOverlpSD() << std::endl;
02170     for (int j = 0 ; j < QCMDSystem<TPSolver>::domains[i].getNbOverlpSD(); ++j) {
02171       globalIndexOfRemoteSD = QCMDSystem<TPSolver>::domains[i].getRemDomainId(j) ;
02172       //      std::cout <<"                 communique avec " << globalIndexOfRemoteSD  <<std::endl;
02173       lineMat[i][globalIndexOfRemoteSD].step = -1;
02174       if(sbdOnMyProc[globalIndexOfRemoteSD] != -1 ){    // SD not on my proc
02175         ++_nbLocalContrib ;
02176       }
02177     }
02178     lineMat[i][myGlobalId].proc   = myrank ;
02179     lineMat[i][myGlobalId].sbdId  = myGlobalId;
02180   }
02181   //  _nbLocalContrib = _nbLocalContrib  / 2
02182   _localContrib = new int[2*_nbLocalContrib] ;
02183 #ifdef QC_DEBUG_DENSITY_CONTRIB_LEVEL2
02184   print(sizeMax, sbdOnMyProc ,"   Correspondance globale locale ");
02185   print(QCMDSystem<TPSolver>::nbDomains,sizeMax,&(lineMat[0][0]),"  LINEMAT init ");
02186 #endif
02187  //
02188  //  suppression des communications locales 
02189  //
02190   int k = 0 ;
02191   for ( int i = 0; i < QCMDSystem<TPSolver>::nbDomains; ++i) {
02192     for (int j = 0 ; j < QCMDSystem<TPSolver>::domains[i].getNbOverlpSD(); ++j) {
02193       globalIndexOfRemoteSD = QCMDSystem<TPSolver>::domains[i].getRemDomainId(j) ;
02194       if(lineMat[i][globalIndexOfRemoteSD].step == -1 && sbdOnMyProc[globalIndexOfRemoteSD] !=-1 ) {     // not on my proc
02195         //      if( i < sbdOnMyProc[globalIndexOfRemoteSD] ) {
02196           _localContrib[2*k] = i ; _localContrib[2*k+1] = sbdOnMyProc[globalIndexOfRemoteSD] ; 
02197           ++k ;
02198           //    }
02199         lineMat[i][globalIndexOfRemoteSD].step = 0 ;
02200 #ifdef QC_DEBUG_DENSITY_MPI
02201         std::cout << " Pas de communication entre "<< i <<" (loc) "<<  this->domains[i].getId() 
02202                   <<" (num glob) et j " << sbdOnMyProc[globalIndexOfRemoteSD] <<" (loc) "
02203                   << j <<" ( num glob)   car sur le meme processeur" <<std::endl;
02204 #endif
02205        }
02206      }
02207   }
02208 #ifdef QC_DEBUG_DENSITY_MPI
02209   k = 2*_nbLocalContrib ;
02210   print(k, _localContrib,"  Contributions  locales ");
02211 #endif
02212   
02213   //
02214   //  suppression des communications en double lorsque l'on demande deux fois un sous domaine pour deux sousdomaines.
02215   //
02216   for (int j = 0 ; j < sizeMax; ++j) {
02217     bool firstUseSD = true ;
02218     for ( int i = 0; i < QCMDSystem<TPSolver>::nbDomains; ++i) {
02219       if( lineMat[i][j].step == -1){  // communication entre le SD[this->domains[i].getId()] et le SD[j]
02220        if( ! firstUseSD){ 
02221          lineMat[i][j].step = 0 ;
02222 #ifdef QC_DEBUG_DENSITY_MPI
02223          std::cout << " Suppression de la communication en double "<< i <<" (loc) "<<  this->domains[i].getId() 
02224                    <<" (num glob) et j " << j <<" ( num glob)   " <<std::endl;
02225 #endif
02226        }
02227        firstUseSD = false ;
02228       } 
02229     }
02230   }
02231   //  
02232 #ifdef QC_DEBUG_DENSITY_MPI
02233   print(QCMDSystem<TPSolver>::nbDomains,sizeMax,&(lineMat[0][0]),"  LINEMAT init ");
02234 #endif
02235   //
02236   //  MPI 
02237   int etape, msgSize = QCMDSystem<TPSolver>::nbDomains*sizeMax ; ;
02238   QCComm *matComm;
02239   int sdOnProc[sizeMax], sdOnline[sizeMax], lineOnSD[sizeMax] ;
02240   if (myrank == 0){
02241     int nbComm =0 ; 
02242     //
02243     matComm = new QCComm[sizeMax*sizeMax] ;
02244     MPI_Gather(lineMat,msgSize,QCCOMM,matComm,msgSize,QCCOMM,0,QC_COMM_WORLD);
02245     //
02246     // on repere l'affectation des sous domaines sur les procs
02247     //
02248     for (int i = 0 ; i <  sizeMax  ; ++i ){
02249       for (int j = 0 ; j <  sizeMax  ; ++j ){
02250         if( matComm[j+i*sizeMax].sbdId != -1){
02251           sdOnProc[ matComm[j+i*sizeMax].sbdId ] = matComm[j+i*sizeMax].proc;
02252           sdOnline[ matComm[j+i*sizeMax].sbdId ] = i ;
02253           lineOnSD[ i ]                          = matComm[j+i*sizeMax].sbdId  ;
02254           break ;
02255         }
02256       }
02257     }
02258     //
02259     for (int i = 0 ; i <  sizeMax*sizeMax  ; ++i ){
02260       nbComm +=  matComm[i].step ; 
02261     }   
02262     for (int i = 0 ; i <  sizeMax  ; ++i ){
02263       for (int j = 0 ; j<  sizeMax ; ++j ){
02264         matComm[j+i*sizeMax].sbdId = j ;
02265         matComm[j+i*sizeMax].proc  = sdOnProc[ lineOnSD[ j ] ] ;
02266       }   
02267     }
02268     nbComm *= -1;
02269 #ifdef QC_DEBUG_DENSITY_MPI
02270     print(sizeMax, sdOnProc," Affectation des SD sur les procs ");
02271     print(sizeMax, sdOnline," Affectation des SD sur les lignes de la matrice comm1 ");
02272     print(sizeMax, lineOnSD," Affectation des lignes de la matrice comm1 sur les SD");
02273 #endif
02274     //
02275     int tl[nproc], tc[nproc], sdPerProc[nproc] ;
02276     for (int i = 0 ; i <  nproc  ; ++i ){
02277       tc[i] = -1 ;    tl[i] = -1 ; sdPerProc[i] = 0 ;
02278     }
02279     for (int j = 0 ; j<  sizeMax ; ++j ){
02280       sdPerProc[sdOnProc[j]] += 1 ;
02281     }
02282 #ifdef QC_DEBUG_DENSITY_MPI
02283     print(nproc, sdPerProc," Number of sub domains per processor ");
02284 #endif
02285     //
02286     for (etape = 1 ; nbComm >0 ; ++etape){
02287       for ( int line = 0 ; line < sizeMax ; ++line){
02288         for ( int col = 0 ; col < sizeMax ; ++col){
02289           if( matComm[col+line*sizeMax].step < 0 && tc[sdOnProc[col]] < etape && tl[sdOnProc[line]] < etape ){
02290             matComm[col+line*sizeMax].step = etape ; 
02291             tc[sdOnProc[col]]              = etape ;
02292             tl[sdOnProc[line]]             = etape ; 
02293             --nbComm ; 
02294             break ;
02295           }
02296         }
02297       }
02298     }
02299     --etape;
02300  #ifdef QC_DEBUG_DENSITY_MPI
02301      std::cout << " Affichage de la matrice de communication " <<std::endl ;
02302      for (int i = 0 ; i <  sizeMax  ; ++i ){
02303        std::cout << "dom : " << i ;
02304        for (int j = 0 ; j<  sizeMax ; ++j ){
02305         std::cout << " ("<<matComm[j+i*sizeMax].proc<< ","<<matComm[j+i*sizeMax].sbdId 
02306                   << "," <<matComm[j+i*sizeMax].step <<") " ;
02307        }
02308        std::cout <<std::endl;
02309      }
02310 #endif
02311 
02312   }
02313   else {
02314     MPI_Gather(lineMat,msgSize,QCCOMM,matComm,msgSize,QCCOMM,0,QC_COMM_WORLD);
02315   }
02316   //  std::cout << " Diffusion et construction des etapes de communications " << std::endl;
02317   //
02318   //  Diffusion des etapes de communications
02319   //
02320   MPI_Bcast(&etape,1,MPI_INT,0,QC_COMM_WORLD);
02321   _recvs.resize(etape) ;   _sends.resize(etape) ; 
02322   //
02323   // La ligne = (proc + SD que l'on doir recevoir)
02324   MPI_Scatter(matComm, msgSize,QCCOMM,lineMat,msgSize,QCCOMM,0,QC_COMM_WORLD);
02325   QCComm empty(-1,-1,-1);
02326   //
02327 #ifdef QC_DEBUG_DENSITY_MPI
02328   std::cout <<"   Nombre max reel d'etapes :" << etape <<std::endl ;
02329   print(QCMDSystem<TPSolver>::nbDomains,sizeMax,&(lineMat[0][0]),"  LINEMAT after MPI_Scatter ");
02330 #endif
02331 
02332   for ( int i = 0; i < QCMDSystem<TPSolver>::nbDomains; ++i) {
02333     for (int j = 0 ; j <  sizeMax  ; ++j ){
02334       if(lineMat[i][j].step != 0) {
02335         _recvs[lineMat[i][j].step-1] =  lineMat[i][j];
02336         --_recvs[lineMat[i][j].step-1].step;
02337       }
02338     }
02339   }
02340   //  this->printProcVect(_recvs,false) ;
02341   //  std::cout << " ---------------  STEP 2  -------------------------------" <<std::endl ;
02342   if (myrank == 0){
02343     for (int i = 0 ; i <  sizeMax  ; ++i ){
02344       for (int j = 0 ; j<  sizeMax ; ++j ){
02345         matComm[j+i*sizeMax].sbdId = j ;
02346         matComm[j+i*sizeMax].proc  = sdOnProc[i ] ; 
02347      }   
02348     } 
02349     QCComm tmp;
02350     for (int i = 0 ; i <  sizeMax ; ++i ){
02351       for (int j = 0 ; j < i ; ++j ){
02352         tmp                   = matComm[j+i*sizeMax] ;
02353         matComm[j+i*sizeMax] = matComm[i+j*sizeMax];
02354         matComm[i+j*sizeMax] = tmp ;
02355       }
02356     }
02357   }
02358   //  std::cout << " Affichage de matcom1 STEP 2" <<std::endl ;
02359   MPI_Scatter(matComm, msgSize,QCCOMM,lineMat,msgSize,QCCOMM,0,QC_COMM_WORLD);
02360 //   printProcVect(lineMat1,true);
02361 //
02362   for ( int i = 0; i < QCMDSystem<TPSolver>::nbDomains; ++i) {
02363     for (int j = 0 ; j <  sizeMax  ; ++j ){
02364       if(lineMat[i][j].step != 0) {
02365         _sends[lineMat[i][j].step-1] =  lineMat[i][j];
02366         --_sends[lineMat[i][j].step-1].step;
02367       }
02368     }
02369   }
02370   MPI_Type_free(&QCCOMM); 
02371   //
02372   QC_TRACE_MPI("END   QCDistMDSystem<TPSolver>::fillSDsCommVect ");
02373 }
02374 
02378 template <class TPSolver>
02379 void  
02380 QCDistMDSystem<TPSolver>::printProcArray (int * array, bool issend) const {
02381 
02382   std::cout << "PROC " << myrank;
02383   
02384   if (issend) {
02385     std::cout << ": SENDS => {";
02386   } else {
02387     std::cout << ": RECVS <= {";
02388   }
02389   
02390   for (int i = 0; i < nproc; ++i) {
02391 
02392     int poffset = i * (maxOvrlpSDs + SDS_HEAD);
02393     
02394     if (array[poffset + NB_SDS] > 0) {
02395       std::cout << " p" << i << ":(";
02396     }
02397     
02398     for (int j=0; j < array[poffset + NB_SDS]; ++j) {
02399       std::cout << " " <<  array[poffset + SDS_HEAD + j];
02400     }
02401     
02402     if (array[poffset + NB_SDS] > 0) {
02403       std::cout << ")";
02404     }
02405   }
02406 
02407   std::cout << " }" << std::endl;
02408 
02409 }
02413 template <class TPSolver>
02414 void  
02415 QCDistMDSystem<TPSolver>::printProcVect (const vector<QCComm>& vect, 
02416                                          bool issend) const {
02417 
02418   std::cout << "PROC " << myrank;
02419   if (issend) {
02420     std::cout << ": SENDS => (PROC, SD, STEP) {";
02421   
02422   } else {
02423     std::cout << ": RECVS <= (PROC, SD, STEP) {";
02424   }
02425   
02426   for (int i=0; i < static_cast<int>(vect.size()); ++i) {
02427     std::cout << " (" << vect[i].proc << ", " << vect[i].sbdId
02428          << ", " << vect[i].step << ")";
02429   }
02430   std::cout << " }" << std::endl;
02431 }
02432 
02433 template <class TPSolver>
02434 QCFloat  
02435 QCDistMDSystem<TPSolver>::computeSumOfEigenValOccupied (){
02436   QCFloat t = 0.0, tloc=0.0;
02437   QCFloat * eigenVal = NULL;
02438   for (int i = 0 ; i < QCMDSystem<TPSolver>::nbDomains; ++i) {
02439         eigenVal      = QCMDSystem<TPSolver>::domains[i].getEigenVal();      
02440         for (int j = 0 ; j < QCMDSystem<TPSolver>::domains[i].getNbOccupiedOM(); ++j) {
02441           tloc += eigenVal[j] ;
02442         }
02443   }
02444   MPI_Allreduce(&tloc, &t, 1, MPI_DOUBLE, MPI_SUM, QC_COMM_WORLD);
02445 
02446   return t ;
02447 }
02448 
02449 
02453 template <class TPSolver>
02454 template <class TPManager>
02455 void 
02456 QCDistMDSystem<TPSolver>::writeDensityOnFileAscii(TPManager& manager,  const std::string& path){ 
02457  //
02458   QC_TRACE_OUT("BEGIN QCDistSystem<TPSolver>::writeDensityOnFileAscii   ");
02459   //
02460   std::string  fileNameOri ;
02461   fileNameOri = manager.getFiles().getResultFile() ;
02462   //  std::cout << " myrank : " << this->myrank <<std::endl;
02463   if(  this->myrank == 0){
02464     std::ofstream out;
02465     std::string fileName , 
02466     fileNameOri = manager.getFiles().getResultFile() ;
02467     fileName = path + "/" + manager.getFiles().getResultFile() +"-para"+ intToString( this->nproc) + "-density.ascii";
02468     //
02469     //    std::cout <<" FILENAME : " << fileName <<std::endl;
02470     out.open(fileName.c_str());
02471     out  << "BASIC-PARA     "<< this->nproc     << "  "<< manager.getPartitioner()->getNbPartitions() <<std::endl;
02472     for(int p = 0 ; p < this->nproc ; ++p ){
02473       out << "  "<< manager.files().resultFile() <<p  <<"-density.ascii"<<std::endl ;
02474     }
02475     //
02476     out.close();
02477   }
02478   manager.files().resultFile() +=  intToString(this->myrank) ;
02479   //
02480   //  std::cout <<"Fichier : " <<  manager.getFiles().resultFile()  <<"-density.ascii"<<std::endl;
02481   //
02482   // We store the local density matrices in the shared format
02483   //
02484   QCMDSystem<TPSolver>::writeDensityOnFileAscii(manager, path);
02485   //
02486   manager.files().resultFile() = fileNameOri ;
02487   //
02488   QC_TRACE_OUT("END   QCDistSystem<TPSolver>::writeDensityOnFileAscii   " ) ;
02489 }
02490 //
02491 //
02492 //
02493 template <class TPSolver>
02494 void 
02495 QCDistMDSystem<TPSolver>::readDensityFromFileAscii(const QCFiles & files,  const string& path){
02496   //
02497   QC_TRACE_INIT("BEGIN QCDistMDSystem<TPSolver>::readDensityOnFileAscii   "<<files.getDensityFile());
02498   //  
02499   std::ifstream data ;
02500   std::string fileName , type ;
02501   //  fileName = path + "/" + files.getDensityFile() ;
02502   fileName = files.getDensityFile() ;
02503   data.open(fileName.c_str());
02504   //
02505   if (!data){
02506     std::cerr << "Error to open density Matrix in File " << fileName <<std::endl; 
02507     exit(EXIT_FAILURE) ;
02508   }
02509   //
02510   int numberOfProc ;
02511   data >> type ; 
02512   if (type ==  "BASIC-PARA"){
02513     data >> numberOfProc ;
02514   }
02515   else {
02516     std::cerr << "Bad type ("<<type<<") of storage in density Matrix in File " 
02517               << fileName <<std::endl; 
02518     exit(EXIT_FAILURE) ;
02519   }
02520   std::cerr << " NOT YET IMPLEMENTED " <<std::endl;
02521   exit(EXIT_FAILURE);
02522   //
02523   //
02524   //
02525   QC_TRACE_INIT("END   QCDistMDSystem<TPSolver>::readDensityOnFileAscii   ");
02526 }
02527 //
02528 // Implicite 
02529 //
02530 template class QCDistMDSystem<QCDCAlgo>;
02531 //
02532 QCMANAGER_METH_EXPL_INST_DIST_MD_PARAM(void QCDistMDSystem<QCDCAlgo>::writeDensityOnFileAscii, 
02533                                        ONE_PARAM(const std::string&));
02534 
02535 QCMANAGER_METH_EXPL_INST_DIST_MD_PARAM(void QCDistMDSystem<QCDCAlgo>::init, ONE_PARAM(const std::string&));
02536 QCMANAGER_METH_EXPL_INST_DIST_MD_PARAM(QCFloat QCDistMDSystem<QCDCAlgo>::completeFockMatricesAndElecEnergy,
02537                                        TWO_PARAMS(bool, int));
02538 
02539 QCMANAGER_METH_EXPL_INST_DIST_MD(void QCDistMDSystem<QCDCAlgo>::completeHamiltonMatrices);
02540 QCMANAGER_METH_EXPL_INST_DIST_MD(void QCDistMDSystem<QCDCAlgo>::completeDensityMatrices);
02541 
02542 QCPARAMETER_METH_EXPL_INST(QCFloat QCDistMDSystem<QCDCAlgo>::applyOptimalDamping);
02543 
02544 
02545  template void
02546  QCDistMDSystem<QCDCAlgo>::adjustFermiEnergySeq(TQCDistMDSystem&, QCMemory&);
02547 
02548 
02549 
02550 
02551 
02552 

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