00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
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
00057
00058 #if defined(MPICL_TRACE) && defined(__QC_xlC__)
00059
00060 # include "pcontrol.h"
00061
00062
00063
00064
00065
00066
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
00166
00167
00168
00169
00170
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
00240 if (isMine(domainIdx)) {
00241 int localIdx = _globalToLocalNum[domainIdx] ;
00242 QCMDSystem<TPSolver>::addAtom(coords,type, globalIdx,localIdx, zone, atomIdx);
00243 }
00244
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
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
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
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
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
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
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445
00446
00447
00448
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
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
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
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
00528 for (int step=0; step < QCMDSystem<TPSolver>::totalNbDomains-1; ++step) {
00529
00530
00531 atomSendBuf = atomsBuffers[sendIdx];
00532 orderSendBuf = orderBuffers[sendIdx];
00533 ovrlpSendBuf = ovrlpBuffers[sendIdx];
00534 zoneSendBuf = zoneBuffers [sendIdx];
00535
00536
00537
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 (
00589 params,
00590
00591 QCMDSystem<TPSolver>::domains[i],
00592
00593 pipeline[i],
00594
00595 QCMDSystem<TPSolver>::domains[i].getRepInterElecIntegrals(),
00596 hamiltonHAA,
00597 manager.getGeneralData(),
00598 step,
00599 QCMDSystem<TPSolver>::totalNbDomains);
00600
00601 }
00602
00603
00604
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 (
00624 params,
00625
00626 QCMDSystem<TPSolver>::domains[0],
00627
00628 pipeline[0],
00629
00630 QCMDSystem<TPSolver>::domains[0].getRepInterElecIntegrals(),
00631 hamiltonHAA,
00632 manager.getGeneralData(),
00633 step,
00634 QCMDSystem<TPSolver>::totalNbDomains);
00635
00636
00637
00638 MPI_Waitall(H_NB_REQUESTS/2, requests, statuses);
00639
00640
00641
00642
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 ) {
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
00690 const QCGeneralData& data = manager.getGeneralData();
00691
00692
00693 QCSymMatrix& QCRestrict fockFAA = model.getSpWorkingAA();
00694
00695
00696 QCFloat domainElecEnergy;
00697
00698
00699 #ifdef QC_VERBOSE
00700 std::cout << "PROC " << myrank << ": * FOCK CONTRIBUTIONS *" << std::endl;
00701 #endif
00702
00703
00704
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
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
00736
00737
00738
00739
00740 for (int step=0; step < QCMDSystem<TPSolver>::totalNbDomains-1; ++step) {
00741
00742
00743
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
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
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 (
00817 params,
00818
00819 QCMDSystem<TPSolver>::domains[i],
00820
00821 pipeline[i],
00822 densBuffers[currIdx],
00823
00824 QCMDSystem<TPSolver>::domains[i].getRepInterElecIntegrals(),
00825 fockFAA, data, step, QCMDSystem<TPSolver>::totalNbDomains,
00826 data.getIntgAcquisitionMethod(), isFirstCall);
00827 }
00828
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 (
00845 params,
00846
00847 QCMDSystem<TPSolver>::domains[0],
00848
00849 pipeline[0], densBuffers[currIdx],
00850
00851 QCMDSystem<TPSolver>::domains[0].getRepInterElecIntegrals(),
00852 fockFAA, data, step, QCMDSystem<TPSolver>::totalNbDomains, data.getIntgAcquisitionMethod(),
00853 isFirstCall);
00854
00855
00856
00857 MPI_Waitall(F_NB_REQUESTS/2, requests, statuses);
00858
00859
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
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
00891
00892 this->elecEnergy = QC_ZERO;
00893
00894
00895
00896 for (int i=0; i < QCMDSystem<TPSolver>::nbDomains; ++i) {
00897
00898 if (rewind) {
00899 QCMDSystem<TPSolver>::domains[i].getRepInterElecIntegrals().getIntgReader()->rewindFile();
00900 }
00901
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
00961 QCDensityGtr<QCSymMatrix> * QCRestrict densityP;
00962 QCSymMatrix * QCRestrict weightedDensity;
00963 QCSymMatrix remWeightedDensity;
00964
00965
00966 QCMatrix& QCRestrict wDensityAB = model.getSpWorkingAB();
00967 QCMatrix& QCRestrict interPAB = model.getSpWorkingAB2();
00968 QCSymMatrix& QCRestrict wDensityAA = model.getSpWorkingAA();
00969
00970
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
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
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
01008
01009 int buffSend = -1 , buffRecv = -1 ;
01010 if ( _sends[currentStep].proc != -1 ) {
01011 buffSend = (buffSend+1)%2 ;
01012 _sends[currentStep].step = buffSend ;
01013 postDensitySends(currentStep, sendInfo[buffSend], sendReq[buffSend]);
01014
01015 }
01016 int recvIdx ;
01017 if ( _recvs[currentStep].proc != -1 ) {
01018 buffRecv = (buffRecv+1)%2 ;
01019 _recvs[currentStep].step = buffRecv ;
01020 postDensityRecvs(currentStep, buffRecv, recvInfo[buffRecv], recvReq[buffRecv]);
01021
01022
01023 }
01024
01025
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
01040
01041
01042
01043
01044
01045 }
01046 #ifdef QC_DEBUG_DENSITY_CONTRIB_LEVEL2
01047 std::cout << " END Local Contributers" << std::endl;
01048 #endif
01049
01050
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 ;
01064 _sends[currentStep+1].step = buffSend ;
01065 postDensitySends(currentStep+1, sendInfo[buffSend], sendReq[buffSend]);
01066
01067 }
01068 if ( _recvs[currentStep+1].proc != -1 ) {
01069 buffRecv = (buffRecv+1)%2 ;
01070 _recvs[currentStep+1].step = buffRecv ;
01071 postDensityRecvs(currentStep+1, buffRecv, recvInfo[buffRecv], recvReq[buffRecv]);
01072
01073 }
01074 }
01075 if (_sends[currentStep].proc != -1 ) {
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
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
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
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
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
01228
01229
01230
01231
01232
01233
01234
01235 mpierr = MPI_Irecv (recvInfo, D_INFO_SIZE, MPI_INT, fromId, (step * D_NB_REQUESTS) + DCONTRIB_INFO,
01236 QC_COMM_WORLD, &requests[0]);
01237
01238
01239
01240
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
01246
01247
01248
01249
01250 mpierr = MPI_Irecv (ovrlpRecvBuf, maxOvrlpSize, MPI_INT, fromId, (step * D_NB_REQUESTS) + DCONTRIB_OVRLP,
01251 QC_COMM_WORLD, &requests[2]);
01252
01253
01254
01255
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
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
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
01458
01459
01460 void compute_min_max (QCFloatFloat * in, QCFloatFloat * inout,
01461 int * len, MPI_Datatype * ) {
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& ,
01568 QCMemory& ) {
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
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
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
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
01660
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
01667 QCFloatInt localpair, maxpair, resultpair;
01668
01669 bool rightNumberOfElec = false;
01670
01671
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
01679 topEnergyLevels[numSD].energyLevel = -DBL_MAX;
01680 }
01681 }
01682
01683 int loop = 0, choice = -1;
01684 while (!rightNumberOfElec) {
01685
01686
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
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
01704 --numOMToChange;
01705 if (numOMToChange >= 0) {
01706 choice = 1 ;
01707 energyLevelToChange = QCMDSystem<TPSolver>::domains[numSDToChange].getEigenVal(numOMToChange);
01708
01709 } else {
01710
01711 choice = 2 ;
01712 energyLevelToChange = -DBL_MAX;
01713 }
01714
01715
01716 } else {
01717
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
01745
01746 this->fermiEnergy = maxpair.floatval;
01747
01748 QC_TRACE_ENER(" (" << loop << " loops) fermi energy = " << this->fermiEnergy );
01749
01750
01751
01752 for (numSD = 0; numSD < QCMDSystem<TPSolver>::nbDomains; ++numSD) {
01753
01754
01755
01756
01757 if (topEnergyLevels[numSD].energyLevel == -DBL_MAX) {
01758 QCMDSystem<TPSolver>::domains[topEnergyLevels[numSD].numSD].setNbFilledOccupiedOM(0, 0);
01759 }
01760
01761 else if (QCMDSystem<TPSolver>::domains[topEnergyLevels[numSD].numSD].getOrbitalOccupN(topEnergyLevels[numSD].numOM)
01762 == QC_TWO) {
01763
01764
01765 QCMDSystem<TPSolver>::domains[topEnergyLevels[numSD].numSD].setNbFilledOccupiedOM(topEnergyLevels[numSD].numOM + 1,
01766 topEnergyLevels[numSD].numOM + 1);
01767 }
01768
01769
01770 else {
01771
01772
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
01792
01793 int& numSDToChange = topEnergyLevels[0].numSD;
01794 int& numOMToChange = topEnergyLevels[0].numOM;
01795 QCFloat& energyLevelToChange = topEnergyLevels[0].energyLevel;
01796
01797
01798 QCFloatInt localpair, minpair, resultpair;
01799
01800
01801 int numSD;
01802 int numOfLumoOfSD;
01803
01804 bool rightNumberOfElec = false;
01805
01806
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
01815 topEnergyLevels[numSD].energyLevel = DBL_MAX;
01816 }
01817 }
01818
01819 int loop = 0;
01820 while (!rightNumberOfElec) {
01821
01822
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
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
01840 numOMToChange++;
01841 if (numOMToChange < nbOM[numSDToChange]) {
01842 energyLevelToChange = QCMDSystem<TPSolver>::domains[numSDToChange].getEigenVal(numOMToChange);
01843
01844 } else {
01845
01846 energyLevelToChange = DBL_MAX;
01847 }
01848
01849
01850 } else {
01851
01852 QCMDSystem<TPSolver>::domains[numSDToChange].
01853 setOrbitalOccupN(numOMToChange,
01854 (this->nbElectrons - sumNbElecs) / QCMDSystem<TPSolver>::domains[numSDToChange].getBFactor(numOMToChange));
01855
01856
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
01878
01879 this->fermiEnergy = minpair.floatval;
01880
01881 QC_TRACE_ENER(" proc "<<myrank<<" (" << loop << " loops) fermi energy = " << this->fermiEnergy );
01882
01883
01884
01885 for (numSD = 0; numSD < QCMDSystem<TPSolver>::nbDomains; ++numSD) {
01886
01887
01888
01889
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
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
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
01914 else {
01915
01916
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
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
02037 this->fillSDsCommArray();
02038
02039
02040 MPI_Alltoall(sdsFromProc, (maxOvrlpSDs + SDS_HEAD), MPI_INT,
02041 sdsIntoProc, (maxOvrlpSDs + SDS_HEAD), MPI_INT, QC_COMM_WORLD);
02042
02043
02044
02045
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
02081
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
02089
02090
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
02107 if (nbSDsFrom > 1) {
02108 stable_sort(sdsFromProc + poffset + SDS_HEAD, sdsFromProc + poffset + SDS_HEAD + nbSDsFrom);
02109 }
02110 }
02111
02112
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] ;
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
02169
02170 for (int j = 0 ; j < QCMDSystem<TPSolver>::domains[i].getNbOverlpSD(); ++j) {
02171 globalIndexOfRemoteSD = QCMDSystem<TPSolver>::domains[i].getRemDomainId(j) ;
02172
02173 lineMat[i][globalIndexOfRemoteSD].step = -1;
02174 if(sbdOnMyProc[globalIndexOfRemoteSD] != -1 ){
02175 ++_nbLocalContrib ;
02176 }
02177 }
02178 lineMat[i][myGlobalId].proc = myrank ;
02179 lineMat[i][myGlobalId].sbdId = myGlobalId;
02180 }
02181
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
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 ) {
02195
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
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){
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
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
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
02317
02318
02319
02320 MPI_Bcast(&etape,1,MPI_INT,0,QC_COMM_WORLD);
02321 _recvs.resize(etape) ; _sends.resize(etape) ;
02322
02323
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
02341
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
02359 MPI_Scatter(matComm, msgSize,QCCOMM,lineMat,msgSize,QCCOMM,0,QC_COMM_WORLD);
02360
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
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
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
02481
02482
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
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
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