00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include "config.h"
00022 #include "QCLinearPartitioner.hpp"
00023 #include "QCGeneralData.hpp"
00024 #include "QCSystem.hpp"
00025 #include "QCLapack.hpp"
00026 #include "QCTable.hpp"
00027 #include "QCMDSystem.hpp"
00028 #include "QCDistMDSystem.hpp"
00029 #include "QCGlobalSystem.hpp"
00030 #include "QCDCAlgo.hpp"
00031 #include "QCSCFAlgo.hpp"
00032
00033
00034
00035 template <class TPSystem>
00036 const int QCLinearPartitioner<TPSystem>::DIM;
00040 template <class TPSystem>
00041 QCLinearPartitioner<TPSystem>::QCLinearPartitioner (void) :
00042 QCPartitioner<TPSystem>(),
00043 gravityCenter(),
00044 xAxis(),
00045 yAxis(),
00046 zAxis(),
00047 domains(NULL),
00048 atomsInZone(NULL)
00049 {
00050 for(int i= 0 ; i < DIM ; ++i){
00051 partitionGrid[i] = 1 ;
00052 for(int j= 0 ; j < DIM ; ++j){
00053 inertiaMatrix[i][j] = QC_ZERO;
00054 }
00055 }
00056 }
00060 template <class TPSystem>
00061 QCLinearPartitioner<TPSystem>::~QCLinearPartitioner (void) {
00062
00063 if (domains) {
00064 delete [] domains;
00065 domains = NULL;
00066 }
00067 if (atomsInZone) {
00068 delete [] atomsInZone;
00069 atomsInZone = NULL;
00070 }
00071 }
00072
00073
00074
00075 template <class TPSystem>
00076 void
00077 QCLinearPartitioner<TPSystem>::initData(const QCGeneralData& data,
00078 const QCFiles& files,
00079 const int * ,
00080 const int ) {
00081 QC_TRACE_INIT("BEGIN QCLinearPartitioner<TPSystem>::init");
00082
00083 QCPartitioner<TPSystem>::initData(data, files);
00084
00085 if (!this->readPartFile){
00086 this->domains = new QCPartSubDomain [this->nbPartitions];
00087 partitionGrid[COORDX] = this->nbPartitions;
00088 partitionGrid[COORDY] = partitionGrid[COORDZ] = 1;
00089 }
00090 QC_TRACE_INIT("END QCLinearPartitioner<TPSystem>::init");
00091 }
00095 template <class TPSystem>
00096 void
00097 QCLinearPartitioner<TPSystem>::partitioning (TPSystem& system) {
00098
00099 QC_TRACE_PART("BEGIN QCLinearPartitioner<TPSystem>::partitioning");
00100 QCSystem& root = system.getRootSystem();
00101
00102
00104 this->_ordering = new int [root.getNbAtoms()] ;
00105 for(int i = 0 ; i <root.getNbAtoms() ; ++i) {
00106 this->_ordering[i] = i;
00107 }
00111 this->computeInnerReference(root);
00112 this->projectOnInnerReference(root);
00116 this->sortFirstAxis(root);
00117
00118 if(this->domains == NULL){
00119 this->domains = new QCPartSubDomain [this->nbPartitions];
00120 partitionGrid[COORDX] = this->nbPartitions;
00121 partitionGrid[COORDY] = partitionGrid[COORDZ] = 1;
00122 }
00126 this->setNeighbors();
00130 this->setXMinXMax(root);
00131
00132
00133
00134
00135 const int MAX_LOOPS = 1;
00136 for (int cpt = 0; cpt < MAX_LOOPS; ++cpt) {
00137
00138 #ifndef QC_NO_DEBUG
00139 std::cout << std::endl << "** Equilibrate: loop " << cpt+1 << "/" << MAX_LOOPS << std::endl;
00140 #endif
00141
00144 this->setAtomsInKernel(root);
00148 this->setAtomsInBuffers(root);
00149
00150 for (int nsd = 0 ; nsd < this->nbPartitions ; ++nsd){
00151 domains[nsd].computeNumbers();
00152 }
00156 this->setLoad(root);
00160 bool cvReached = equilibrateAOs(this->radius1, this->radius2);
00161 if(cvReached || cpt + 1 == MAX_LOOPS) {
00162 break;
00163 }
00167 this->resetData();
00168 }
00169 this->setAtomZones(root.getNbAtoms());
00171 this->extractDomains(system);
00172 this->updateInternalStructure(system);
00173
00174
00175 QC_TRACE_PART("END QCLinearPartitioner<TPSystem>::partitioning");
00176
00177 }
00178
00182 template <class TPSystem>
00183 void QCLinearPartitioner<TPSystem>::setAtomZones (int nbAtoms) {
00184
00185 QC_TRACE_PART("BEGIN QCLinearPartitioner<TPSystem>::setAtomZones");
00186
00187 QCAtomIn atomNSD ;
00188 int num ;
00189
00190 if (!atomsInZone) {
00191 atomsInZone = new vector<QCAtomIn > [nbAtoms];
00192 }
00193 else {
00194 for (int i = 0 ; i < nbAtoms ; ++i){
00195 atomsInZone[i].clear() ;
00196 }
00197 }
00198 for(int i = 0; i < this->nbPartitions; ++i){
00199 num = 0 ;
00200 atomNSD.numDomain = i ;
00201
00202
00203
00204 vector<int> & indexCORE = domains[i].getAtomsInKernel() ;
00205 atomNSD.typeZone = QC_CORE;
00206
00207 for(unsigned int j = 0; j < indexCORE.size() ; ++j){
00208 atomNSD.localNum = j ;
00209 atomsInZone[indexCORE[j]].push_back(atomNSD);
00210 }
00211
00212
00213
00214 num = indexCORE.size() ;
00215 atomNSD.typeZone = QC_SHELL1;
00216 for ( int n = 0 ; n < domains[i].getNbNeighbors() ; ++n){
00217 vector<int> & indexOVERLAP = domains[i].getAtomsInShell1(n) ;
00218 for(unsigned int j = 0; j < indexOVERLAP.size() ; ++j){
00219 atomNSD.localNum = num + j ;
00220 atomsInZone[indexOVERLAP[j]].push_back(atomNSD);
00221 }
00222 num += indexOVERLAP.size() ;
00223 }
00224
00225
00226
00227 atomNSD.typeZone = QC_SHELL2 ;
00228 for ( int n = 0 ; n < domains[i].getNbNeighbors() ; ++n){
00229 vector<int> & indexCUT = domains[i].getAtomsInShell2(n) ;
00230 for(unsigned int j = 0; j < indexCUT.size() ; ++j){
00231 atomNSD.localNum = num + j ;
00232 atomsInZone[indexCUT[j]].push_back(atomNSD);
00233 }
00234 num += indexCUT.size() ;
00235 }
00236 }
00237 #ifndef QC_NO_DEBUG
00238 std::cout <<"#######################################################################################################"<<std::endl;
00239 std::cout <<"#######################################################################################################"<<std::endl;
00240 for (int i = 0; i < nbAtoms; ++i) {
00241 vector<QCAtomIn>& v = atomsInZone[i];
00242 std::cout << "Atom : "<< i << " " ;
00243 for(unsigned int j = 0; j < v.size(); ++j) {
00244 std::cout << "( "<< v[j].numDomain << ", " << v[j].localNum << ", "<< v[j].typeZone << ") ";
00245 }
00246 std::cout <<std::endl ;
00247 }
00248 std::cout <<"#######################################################################################################"<<std::endl;
00249 std::cout <<"#######################################################################################################"<<std::endl;
00250 #endif
00251
00252 for (int i = 0; i < this->nbPartitions; ++i) {
00253
00254 domains[i].getSharedAtomsVect().resize(this->nbPartitions);
00255 for (int j = 0; j < this->nbPartitions; ++j) {
00256 domains[i].getSharedAtomsVect()[j] = 0;
00257 }
00258 }
00259
00260 int sd1, sd2;
00261 for (int i = 0; i < nbAtoms; ++i) {
00262 vector<QCAtomIn>& v = atomsInZone[i];
00263
00264 for (unsigned int j = 0; j < v.size()-1; ++j) {
00265 sd1 = v[j].numDomain;
00266
00267 for (unsigned int k = j+1; k < v.size(); ++k) {
00268 sd2 = v[k].numDomain;
00269
00270 if (v[j].typeZone != QC_CORE && v[k].typeZone != QC_SHELL2) {
00271 domains[sd1].getSharedAtomsVect()[sd2]++;
00272 }
00273 if (v[k].typeZone != QC_CORE && v[j].typeZone != QC_SHELL2) {
00274 domains[sd2].getSharedAtomsVect()[sd1]++;
00275 }
00276 }
00277 }
00278 }
00279
00280 for (int i = 0; i < this->nbPartitions; ++i) {
00281 for (int j = 0; j < this->nbPartitions; ++j) {
00282 if (domains[i].getSharedAtomsVect()[j] != 0) {
00283 domains[i].getNbOverlp()++;
00284 }
00285 }
00286 }
00287 QC_TRACE_PART("END QCLinearPartitioner<TPSystem>::setAtomZones");
00288 }
00289
00293 template <class TPSystem>
00294 void QCLinearPartitioner<TPSystem>::setLoad (const QCSystem& system) {
00295
00296 QC_TRACE_PART("BEGIN QCLinearPartitioner<TPSystem>::setLoad");
00297 int nsd, k, loadOverlap, loadCut,loadKernel;
00298 unsigned int at;
00299
00300
00301 #ifndef QC_NO_DEBUG
00302 for (int i = 0; i < system.getNbAtoms(); ++i) {
00303 std::cout << "(" << table[system.getType(i)] << "," << this->charge[i] << ") ";
00304 }
00305 std::cout << std::endl;
00306 #endif
00307
00308 for(nsd = 0 ; nsd < this->nbPartitions ; ++nsd) {
00309
00310 #ifndef QC_NO_DEBUG
00311 std::cout << "-> Domain " << nsd << std::endl;
00312 #endif
00313 loadOverlap = loadCut = loadKernel = 0 ;
00314 vector<int> & indexK = domains[nsd].getAtomsInKernel() ;
00315
00316 for (k = 0 ; k < domains[nsd].getNbNeighbors() ; ++k){
00317 vector<int> & indexO = domains[nsd].getAtomsInShell1(k) ;
00318 vector<int> & indexC = domains[nsd].getAtomsInShell2(k) ;
00319 #ifndef QC_NO_DEBUG
00320 std::cout << "overlap(" << k << "): ";
00321 #endif
00322 for(at = 0 ; at < indexO.size() ; ++at){
00323 loadOverlap += this->charge[indexO[at]] ;
00324 #ifndef QC_NO_DEBUG
00325 std::cout << " (" << indexO[at] << "," << this->charge[indexO[at]] << ")";
00326 #endif
00327 }
00328 #ifndef QC_NO_DEBUG
00329 std::cout << std::endl << "cut(" << k << "): ";
00330 #endif
00331 for(at = 0 ; at < indexC.size() ; ++at){
00332 loadCut += this->charge[indexC[at]] ;
00333 #ifndef QC_NO_DEBUG
00334 std::cout << " (" << indexC[at] << "," << this->charge[indexC[at]] << ")";
00335 #endif
00336 }
00337 std::cout << std::endl;
00338 }
00339
00340
00341 #ifndef QC_NO_DEBUG
00342 std::cout << "kernel: ";
00343 #endif
00344
00345
00346 for(at = 0 ; at < indexK.size() ; ++at){
00347 loadKernel += this->charge[indexK[at]] ;
00348
00349 #ifndef QC_NO_DEBUG
00350 std::cout << " (" << indexK[at] << "," << this->charge[indexK[at]] << ")";
00351 #endif
00352 }
00353
00354
00355 domains[nsd].setLoadInKernel(loadKernel);
00356 domains[nsd].setLoadInShell1(loadOverlap);
00357 domains[nsd].setLoadInShell2(loadCut);
00358
00359
00360 #ifndef QC_NO_DEBUG
00361 std::cout << std::endl << std::endl << "Domain " << nsd
00362 << " (AOs): core = " << loadKernel
00363 << " overlap = " << loadOverlap
00364 << " cut = " << loadCut
00365 << " total = " << loadKernel + loadOverlap + loadCut
00366 << std::endl;
00367 #endif
00368 }
00369 QC_TRACE_PART("END QCLinearPartitioner<TPSystem>::setLoad");
00370 }
00371
00375 template <class TPSystem>
00376 void QCLinearPartitioner<TPSystem>::extractDomains(TPSystem& mdSystem) {
00377 QC_TRACE_PART("BEGIN QCLinearPartitioner<TPSystem>::extractDomains");
00378
00379
00380
00381 mdSystem.allocDomains(this->nbPartitions);
00382 for (int i=0; i < this->nbPartitions; ++i) {
00383 QCPartSubDomain& domain = this->getDomain(i);
00384 mdSystem.getSubDomain(i).allocStructures(domain);
00385 std::cout <<"****"<<std::endl << " SD : "<< i <<std::endl;
00386 mdSystem.getSubDomain(i).writeOverlapStructure();
00387 }
00388
00389 for (int i=0; i < mdSystem.getNbAtoms(); ++i) {
00390
00391 QCPoint3D current;
00392 mdSystem.getRootSystem().getPointAt(i, current);
00393
00394 int atomIdx, nbsbd;
00395 const vector<QCAtomIn>& atomInZone = getAtomsInZone(i);
00396
00397 atomIdx = i;
00398 nbsbd = atomInZone.size();
00399
00400 int * domainIds = new int [nbsbd];
00401 QCSubDomainZone * zones = new QCSubDomainZone [nbsbd];
00402 int * atomIdxs = new int [nbsbd];
00403
00404 for (int j=0; j < nbsbd; ++j) {
00405
00406 domainIds[j] = atomInZone[j].numDomain;
00407 atomIdxs [j] = atomInZone[j].localNum;
00408 zones [j] = atomInZone[j].typeZone;
00409
00410 mdSystem.addAtom(current, mdSystem.getRootSystem().getType(i),
00411 i, domainIds[j], zones[j], atomIdxs[j]);
00412 }
00413
00414 mdSystem.fillMaps(nbsbd, domainIds, zones, atomIdxs);
00415
00416 delete [] domainIds;
00417 delete [] zones;
00418 delete [] atomIdxs;
00419 }
00420 for (int i=0; i < this->nbPartitions; ++i) {
00421 std::cout <<"****"<<std::endl << " SD : "<< i <<std::endl;
00422 mdSystem.getSubDomain(i).writeOverlapStructure();
00423 }
00424 QC_TRACE_PART("END QCLinearPartitioner<TPSystem>::extractDomains");
00425 }
00426
00430 template <class TPSystem>
00431 void
00432 QCLinearPartitioner<TPSystem>::setToZero (void) {
00433 for (int i=0; i < DIM; ++i) {
00434 gravityCenter [i] = QC_ZERO;
00435 xAxis [i] = QC_ZERO;
00436 yAxis [i] = QC_ZERO;
00437 zAxis [i] = QC_ZERO;
00438 inertiaMatrix [COORDX][i] = QC_ZERO;
00439 inertiaMatrix [COORDY][i] = QC_ZERO;
00440 inertiaMatrix [COORDZ][i] = QC_ZERO;
00441 }
00442 }
00443
00444
00445
00449 template <class TPSystem>
00450 void
00451 QCLinearPartitioner<TPSystem>::resetData (void) {
00452
00453 for (int i = 0; i < this->nbPartitions; ++i) {
00454 domains[i].resetData();
00455 }
00456 }
00457
00458
00459
00460
00461
00462
00466 template <class TPSystem>
00467 void
00468 QCLinearPartitioner<TPSystem>::computeInnerReference (const QCSystem& system) {
00469
00470 QC_TRACE_PART("BEGIN QCLinearPartitioner<TPSystem>::computeInnerReference");
00471 system.buildGravityCenter(gravityCenter);
00472
00473 QCFloat Ixx, Iyy, Izz, Ixy, Ixz, Iyz ;
00474 QCFloat * xyz;
00475
00476 Ixx = Iyy = Izz = Ixy = Ixz = Iyz = 0.0 ;
00477
00478 for(int i = 0 ; i < system.getNbAtoms() ; ++i) {
00479
00480 xyz = system.getCoordsAt(i);
00481
00482
00483 Ixx += (QCPow<2>(xyz[COORDY] - gravityCenter[COORDY]) +
00484 QCPow<2>(xyz[COORDZ] - gravityCenter[COORDZ]));
00485 Iyy += (QCPow<2>(xyz[COORDX] - gravityCenter[COORDX]) +
00486 QCPow<2>(xyz[COORDZ] - gravityCenter[COORDZ]));
00487 Izz += (QCPow<2>(xyz[COORDX] - gravityCenter[COORDX]) +
00488 QCPow<2>(xyz[COORDY] - gravityCenter[COORDY]));
00489
00490 Ixy -= ((xyz[COORDX] - gravityCenter[COORDX]) *
00491 (xyz[COORDY] - gravityCenter[COORDY]));
00492 Ixz -= ((xyz[COORDX] - gravityCenter[COORDX]) *
00493 (xyz[COORDZ] - gravityCenter[COORDZ]));
00494 Iyz -= ((xyz[COORDY] - gravityCenter[COORDY]) *
00495 (xyz[COORDZ] - gravityCenter[COORDZ]));
00496
00497 }
00498
00499
00500
00501 inertiaMatrix[COORDX][COORDX] = Ixx;
00502 inertiaMatrix[COORDY][COORDY] = Iyy;
00503 inertiaMatrix[COORDZ][COORDZ] = Izz;
00504 inertiaMatrix[COORDX][COORDY] = inertiaMatrix[COORDY][COORDX] = Ixy;
00505 inertiaMatrix[COORDX][COORDZ] = inertiaMatrix[COORDZ][COORDX] = Ixz;
00506 inertiaMatrix[COORDY][COORDZ] = inertiaMatrix[COORDZ][COORDY] = Iyz;
00507
00508
00509
00510
00511
00516 const int dim = DIM;
00517 const int lda = DIM;
00518 const int lwork = 3*DIM-1;
00519
00520 int info;
00521
00522 QCFloat a [DIM*DIM];
00523 QCFloat work [lwork];
00524
00525 QCFloat eigenVals [DIM];
00526 QCFloat * eigenVects[DIM];
00527
00531 for (int i = 0; i < DIM; ++i) {
00532 for (int j = 0; j < DIM; ++j) {
00533 a[i+j*DIM] = inertiaMatrix[i][j];
00534 }
00535 }
00536 QC_TRACE_PART("QCLapack(dsyev)");
00537 QCLapack(dsyev)("V", "U", &dim, a, &lda, eigenVals, work, &lwork, &info);
00538
00539
00543 for (int i = 0; i < DIM; ++i) {
00544 eigenVects[i] = new QCFloat [DIM];
00545
00546 for (int j = 0; j < DIM; ++j) {
00547 eigenVects[i][j] = a[i+j*DIM];
00548 }
00549 }
00550
00551
00555 sortEigenVects(eigenVals, eigenVects, DIM);
00559 for (int j = 0; j < DIM; ++j) {
00560
00561 xAxis[j] = eigenVects[j][COORDX];
00562 yAxis[j] = eigenVects[j][COORDY];
00563 zAxis[j] = eigenVects[j][COORDZ];
00564 }
00568 for (int i = 0; i < DIM; ++i) {
00569 delete [] eigenVects[i];
00570 }
00574 #ifndef QC_NO_DEBUG
00575 std::cout << "*** System Inner Reference: " << std::endl
00576 << " \tGravity center: " << gravityCenter << std::endl
00577 << " \tX axis : " << xAxis << std::endl
00578 << " \tY axis : " << yAxis << std::endl
00579 << " \tZ axis : " << zAxis << std::endl;
00580 std::cout << " System Orientation "<< std::endl
00581 << " \te1^e2 : " << xAxis.vect(yAxis)<< std::endl
00582 << " \te2^e3 : " << yAxis.vect(zAxis)<< std::endl
00583 << " \te3^e1 : " << zAxis.vect(xAxis)<< std::endl;
00584 #endif
00585 QC_TRACE_PART("END QCLinearPartitioner<TPSystem>::computeInnerReference");
00586 }
00587
00588
00589
00590
00591
00595 template <class TPSystem>
00596 void
00597 QCLinearPartitioner<TPSystem>::projectOnInnerReference (QCSystem& system) {
00598
00599 QC_TRACE_PART("BEGIN QCLinearPartitioner<TPSystem>::projectOnInnerReference");
00600
00601 QCFloat s[DIM];
00602 QCPoint3D xyz ;
00603
00604 for(int i = 0 ; i < system.getNbAtoms() ; ++i) {
00605
00606 system.getPointAt(i,xyz);
00607
00608 xyz[COORDX] -= gravityCenter[COORDX];
00609 xyz[COORDY] -= gravityCenter[COORDY];
00610 xyz[COORDZ] -= gravityCenter[COORDZ];
00611
00612 s[COORDX] = xyz.dot(xAxis) ;
00613 s[COORDY] = xyz.dot(yAxis) ;
00614 s[COORDZ] = xyz.dot(zAxis) ;
00615
00616 system.setPointAt(i, s);
00617 }
00618
00619 QC_TRACE_PART("END QCLinearPartitioner<TPSystem>::projectOnInnerReference");
00620 }
00621
00622
00623
00624
00628 template <class TPSystem>
00629 void QCLinearPartitioner<TPSystem>::setNeighbors (void) {
00630
00631 QC_TRACE_PART("BEGIN QCLinearPartitioner<TPSystem>::setNeighbors "<<domains<< " "<<this->nbPartitions );
00632
00633 int nsd ;
00634
00635
00636
00637 for (nsd = 0 ; nsd < this->nbPartitions ; ++nsd) {
00638 this->domains[nsd].setIds(nsd, nsd) ;
00639 }
00640 QC_TRACE_PART(" 1 QCLinearPartitioner<TPSystem>::setNeighbors");
00641
00642
00643
00644 domains[0].addNeighbor(1) ;
00645 for (nsd = 1 ; nsd < this->nbPartitions-1; ++nsd){
00646 domains[nsd].addNeighbor(nsd-1) ;
00647 domains[nsd].addNeighbor(nsd+1) ;
00648 }
00649 domains[this->nbPartitions-1].addNeighbor(this->nbPartitions - 2) ;
00650
00651
00652
00653
00654 QC_TRACE_PART(" 2 QCLinearPartitioner<TPSystem>::setNeighbors");
00655 for (nsd = 0 ; nsd < this->nbPartitions ; ++nsd){
00656 domains[nsd].initNbOverlapDomain() ;
00657 }
00658 QC_TRACE_PART("END QCLinearPartitioner<TPSystem>::setNeighbors");
00659 }
00660
00661
00662
00663
00667 template <class TPSystem>
00668 void
00669 QCLinearPartitioner<TPSystem>::setXMinXMax (const QCSystem& root) {
00670
00671 QC_TRACE_PART("BEGIN QCLinearPartitioner<TPSystem>::setXMinXMax");
00672 QCFloat moleculSize[DIM], kernelSize[DIM];
00673 QCPoint3D coinMin, coinMax, coin1, coin2, s1, s2;
00674
00675 root.findMinMax(coinMin, coinMax);
00676
00677 for ( int j = 0 ; j < DIM ; ++j){
00678 coinMax[j] += 0.1 ;
00679 coinMin[j] -= 0.1 ;
00680 }
00681
00682 moleculSize[COORDX] = coinMax[COORDX] - coinMin[COORDX];
00683 moleculSize[COORDY] = coinMax[COORDY] - coinMin[COORDY];
00684 moleculSize[COORDZ] = coinMax[COORDZ] - coinMin[COORDZ];
00685
00686
00687 kernelSize[COORDX] = moleculSize[COORDX]
00688 / static_cast<QCFloat>(partitionGrid[COORDX]);
00689
00690 kernelSize[COORDY] = moleculSize[COORDY]
00691 / static_cast<QCFloat>(partitionGrid[COORDY]);
00692
00693 kernelSize[COORDZ] = moleculSize[COORDZ]
00694 / static_cast<QCFloat>(partitionGrid[COORDZ]);
00695
00696 for(int i = 0 ; i < DIM ; ++i){
00697 coin1[i] = coinMin[i];
00698 coin2[i] = coinMax[i];
00699 s1[i] = coinMin[i];
00700 s2[i] = coinMax[i];
00701 }
00702
00703 int nyz, nsd;
00704
00705
00706 for (int nsdz = 0; nsdz < partitionGrid[COORDZ]; ++nsdz){
00707 coin1[COORDZ] =
00708 coinMin[COORDZ] + static_cast<QCFloat>(nsdz) * kernelSize[COORDZ];
00709 coin2[COORDZ] = coin1[COORDZ] + kernelSize[COORDZ];
00710 s1[COORDZ] = coin1[COORDZ] ;
00711 s2[COORDZ] = coin2[COORDZ] ;
00712
00713 for (int nsdy = 0; nsdy < partitionGrid[COORDY]; ++nsdy){
00714 nyz = partitionGrid[COORDX]*(nsdy + partitionGrid[COORDY] * nsdz);
00715 s1[COORDY] = coin1[COORDY] + static_cast<double>(nsdy) * kernelSize[COORDY];
00716 s2[COORDY] = s1[COORDY] + kernelSize[COORDY];
00717
00718 s1[COORDX] = coinMin[COORDX] ;
00719
00720 s2[COORDX] = coinMin[COORDX]
00721 + static_cast<double>(partitionGrid[COORDX]) * kernelSize[COORDX];
00722
00723 domains[ nyz].setXMin(s1) ;
00724 domains[partitionGrid[COORDX]-1 + nyz].setXMax(s2) ;
00725
00726 for (int nsdx = 1; nsdx < partitionGrid[COORDX]; ++nsdx){
00727 nsd = nsdx + nyz ;
00728 s1[COORDX] = coin1[COORDX]
00729 + static_cast<double>( nsdx) * kernelSize[COORDX];
00730
00731 s2[COORDX] = s1[COORDX] ;
00732
00733 domains[nsd-1].setXMax( s2) ;
00734 domains[nsd ].setXMin( s1);
00735 }
00736 }
00737 }
00738 QC_TRACE_PART("END QCLinearPartitioner<TPSystem>::setXMinXMax");
00739
00740 }
00741
00742
00743
00744
00745
00749 template <class TPSystem>
00750 void
00751 QCLinearPartitioner<TPSystem>::setAtomsInKernel (const QCSystem& system) {
00752
00753 QC_TRACE_PART("BEGIN QCLinearPartitioner<TPSystem>::setAtomsInKernel");
00754
00755 int nsd;
00756
00757 for (nsd = 0; nsd < this->nbPartitions; ++nsd) {
00758 domains[nsd].addAtomsInKernel(system.getAtoms(), this->useFrag);
00759 }
00760
00761 #ifndef QC_NO_DEBUG
00762 int nb, nbg = 0;
00763 for (nsd = 0; nsd < this->nbPartitions; ++nsd) {
00764 nb = domains[nsd].getNbAtomsKernel();
00765 nbg += nb;
00766
00767 if (nb == 0) {
00768 cerr << "*** WARNING: no atoms in kernel of domain " << nsd << std::endl;
00769 }
00770 }
00771 int diff = system.getNbAtoms() - nbg;
00772 if (diff) {
00773 cerr << "*** WARNING: "
00774 << diff << "/" << system.getNbAtoms() << " atoms missing" << std::endl;
00775 }
00776 #endif
00777 QC_TRACE_PART("END QCLinearPartitioner<TPSystem>::setAtomsInKernel");
00778
00779 }
00783 template <class TPSystem>
00784 void
00785 QCLinearPartitioner<TPSystem>::setAtomsInBuffers (const QCSystem& system) {
00786
00787 QC_TRACE_PART("BEGIN QCLinearPartitioner<TPSystem>::setAtomsInBuffers");
00788 int nsd, k, neighbor;
00789
00790 for (nsd = 0; nsd < this->nbPartitions; ++nsd) {
00791
00792
00793 for (k = 0; k < domains[nsd].getNbNeighbors(); ++k) {
00794 neighbor = domains[nsd].getNeighborIdx(k);
00795 if (neighbor < nsd) {continue; }
00796
00797 domains[nsd].addAtomsInBuffers(system.getAtoms(), domains[neighbor], k,
00798 this->radius1, this->radius2, this->useFrag);
00799 }
00800
00801 }
00802 QC_TRACE_PART("END QCLinearPartitioner<TPSystem>::setAtomsInBuffers");
00803 }
00807 template <class TPSystem>
00808 bool
00809 QCLinearPartitioner<TPSystem>::equilibrateAOs (QCFloat buffer1, QCFloat buffer2) {
00810
00811 QC_TRACE_PART("BEGIN QCLinearPartitioner<TPSystem>::equilibrateAOs");
00812 int nsd, meanCharge, nb, errorLoad;
00813 QCFloat totalWideSD, totalTmpWideSD, removeLength, gamma, diff;
00814 QCFloat RhalfOverlapAndCut = 0.5*buffer1 + buffer2 ;
00815 QCFloat RhalfOverlap = 0.5*buffer1 , tmp , error ;
00816
00817 QCFloat * tmpWideSD = new double [this->nbPartitions];
00818 QCFloat * wideSD = new double [this->nbPartitions];
00819 QCFloat * c = new double [this->nbPartitions];
00820
00821
00822 bool cvReached;
00823
00824 if( !(partitionGrid[COORDY] == 1 && partitionGrid[COORDZ] == 1) ) {
00825 std::cerr << "** equilibrateAOs: not implemented for 3D partitionning"
00826 << std::endl;
00827 return true;
00828 }
00829
00830 cvReached = false;
00831
00832
00833 meanCharge = 0;
00834 totalWideSD = 0.0;
00835
00836 for(nsd = 0 ; nsd < this->nbPartitions; ++nsd) {
00837 wideSD[nsd] = domains[nsd].getXLength() ;
00838 meanCharge += domains[nsd].getLoad() ;
00839 totalWideSD += wideSD[nsd] ;
00840 c[nsd] = static_cast<double>(domains[nsd].getNbNeighbors()) ;
00841 }
00842 meanCharge /= this->nbPartitions ;
00843
00844
00850 totalTmpWideSD = 0;
00851 nb = 0 ;
00852 removeLength = 0.0 ;
00853
00854 for(nsd = 0 ; nsd < this->nbPartitions ; ++nsd) {
00855 tmp = c[nsd]* RhalfOverlapAndCut ;
00856 tmpWideSD[nsd] =
00857 ( (meanCharge /
00858 static_cast<QCFloat>(domains[nsd].getLoad()) )
00859 *(wideSD[nsd] + tmp) - tmp );
00860
00861 if(tmpWideSD[nsd] < c[nsd] * RhalfOverlap) {
00862 removeLength += c[nsd] * RhalfOverlap ;
00863 } else {
00864 ++nb ;
00865 totalTmpWideSD += tmpWideSD[nsd] ;
00866 }
00867 }
00868
00869
00870
00871 gamma = (totalWideSD - removeLength ) / totalTmpWideSD;
00872 totalWideSD = 0.0 ;
00873 diff = 0.0;
00874 errorLoad = 0 ;
00875
00876
00877 for(nsd = 0 ; nsd < this->nbPartitions ; ++nsd) {
00878 if(tmpWideSD[nsd] < c[nsd] * RhalfOverlap) {
00879 tmpWideSD[nsd] = c[nsd] * RhalfOverlap ;
00880 } else {
00881 tmpWideSD[nsd] = gamma * tmpWideSD[nsd];
00882 }
00883 errorLoad = FQCMax(errorLoad , meanCharge - domains[nsd].getLoad()) ;
00884 diff = FQCMax( diff , FQCAbs(wideSD[nsd] - tmpWideSD[nsd]) );
00885 wideSD[nsd] = tmpWideSD[nsd];
00886 totalWideSD += wideSD[nsd] ;
00887 }
00888
00889 error = static_cast<QCFloat>(errorLoad) / static_cast<QCFloat>(meanCharge) ;
00890
00891 #ifndef QC_NO_DEBUG
00892 std::cout << std::endl << std::endl
00893 << "=========================================================" << std::endl
00894 << "Equilibrate : " << meanCharge << std::endl
00895 << " absolute error on Load : "
00896 << errorLoad << " error in % " << error*100 << std::endl
00897 << " absolute error on distance : "
00898 << diff << std::endl << std::endl;
00899 #endif
00900
00904 if( error < 1.0e-2 ) {
00905 cvReached = true ;
00906 #ifndef QC_NO_DEBUG
00907 std::cout << " La methode a convergee." << std::endl;
00908 #endif
00909 }
00910
00911
00915 QCPoint3D s1, s2;
00916 domains[0].getXMin(s1);
00917 domains[0].getXMax(s2);
00918
00919 for (int nsd = 0 ; nsd < this->nbPartitions-1; ++nsd) {
00920 s1[COORDX] += wideSD[nsd] ;
00921 s2[COORDX] = s1[COORDX] ;
00922 domains[nsd].setXMax( s2);
00923 domains[nsd+1].setXMin( s1);
00924 }
00925
00926
00927 delete[] wideSD ; wideSD = NULL;
00928 delete[] tmpWideSD; tmpWideSD = NULL;
00929
00930 QC_TRACE_PART("END QCLinearPartitioner<TPSystem>::equilibrateAOs");
00931 return cvReached;}
00932
00936 template <class TPSystem>
00937 void QCLinearPartitioner<TPSystem>::sortFirstAxis (QCSystem& system) {
00938 QC_TRACE_PART("BEGIN QCLinearPartitioner<TPSystem>::sortFirstAxis");
00939
00940 QCFloat * xyz1, * xyz2;
00941 int i, j, k;
00942
00943 for (j = 0; j < system.getNbAtoms(); ++j) {
00944 for (i = 0; i < system.getNbAtoms()-1; ++i) {
00945
00946 xyz1 = system.getCoordsAt(i);
00947 xyz2 = system.getCoordsAt(i+1);
00948
00949 if (xyz1[COORDX] > xyz2[COORDX]) {
00950 system.swap(i,i+1) ;
00951 k = this->_ordering[i] ;
00952 this->_ordering[i] = this->_ordering[i+1] ;
00953 this->_ordering[i+1] = k ;
00954
00955 k = this->charge[i ];
00956 this->charge[i ] = this->charge[i+1];
00957 this->charge[i+1] = k;
00958 }
00959 }
00960 }
00961
00962
00963 QC_TRACE_PART("END QCLinearPartitioner<TPSystem>::sortFirstAxis");
00964 }
00965
00969 template <class TPSystem>
00970 void
00971 QCLinearPartitioner<TPSystem>::sortEigenVects (QCFloat *eigenVals,
00972 QCFloat ** eigenVects,
00973 int dim) {
00974
00975 QC_TRACE_PART("BEGIN QCLinearPartitioner<TPSystem>::sortEigenVects");
00976 QCFloat temp;
00977 int i, j, k;
00978
00979 for (j = 0; j < dim; ++j) {
00980 std::cout << "Val "<<j<<" : "<< eigenVals[j]<< " ";
00981 for (i = 0; i < dim; ++i) {
00982 std::cout << eigenVects [j][i] <<" ";
00983 }
00984 std::cout <<std::endl;
00985 }
00986
00987 for (j = 0; j < dim; ++j) {
00988 for (i = 0; i < dim-1; ++i) {
00989
00990 if (FQCAbs(eigenVals[i]) > FQCAbs(eigenVals[i+1])) {
00991
00992 temp = eigenVals [i ];
00993 eigenVals [i ] = eigenVals [i+1];
00994 eigenVals [i+1] = temp;
00995
00996 for (k = 0; k < dim; ++k) {
00997
00998 temp = eigenVects [k][i ];
00999 eigenVects [k][i ] = eigenVects [k][i+1];
01000 eigenVects [k][i+1] = temp;
01001 }
01002 }
01003 }
01004 }
01005 for (j = 0; j < dim; ++j) {
01006 std::cout << "Val "<<j<<" : "<< eigenVals[j]<< " ";
01007 for (i = 0; i < dim; ++i) {
01008 std::cout << eigenVects [j][i] <<" ";
01009 }
01010 std::cout <<std::endl;
01011 }
01012 QC_TRACE("END QCLinearPartitioner<TPSystem>::sortEigenVects");
01013 }
01014
01015
01016
01017
01018
01022 template class QCLinearPartitioner<TQCMDSystem>;
01023
01024
01025
01026 #if defined(HAVE_MPI) && defined(WITH_MPI_SUPPORT)
01027 template class QCLinearPartitioner<TQCDistMDSystem>;
01028 #endif