QCLinearPartitioner.cpp

Go to the documentation of this file.
00001 //*****************************************************************************//
00002 //                                                                             //
00003 //   Copyright (c) 2001                                                        //
00004 //      INRIA                                                                  //
00005 //      54600 VILLERS LES NANCY                                                //
00006 //      France                                                                 //
00007 //                                                                             //
00008 //*****************************************************************************//
00009 //                                                                             //
00010 //               *** NOTICE OF PROPRIETARY INFORMATION ***                     //
00011 //                                                                             //
00012 // The information contained in this file is considered proprietary and the    //
00013 // exclusive property of  INRIA. This information may not be disclosed,        //
00014 // duplicated or used, in whole or in part, for  any purpose  whatsoever       //
00015 // without express written authorization from INRIA                            //
00016 //                                                                             //
00017 //*****************************************************************************//
00018 
00019 
00020 
00021 #include "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 *          /*load*/,
00080                                         const int            /*sizeLoad*/) {
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   }//for i
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   // * Boucle pour l'equilibrage du nombre d'OA par SD *
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     //  Core
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     //  Overlap 
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     //  Cut
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   // Init domain number 
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   //  std::cout <<"   System original " << std::endl << system <<std::endl;
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   //  std::cout <<"   System Projete " << std::endl << system <<std::endl;
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   // Set Global Number
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   // Les voisins
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   //  Allocation
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     //    std::cout << "   Subdomain  " << nsd <<   " Nombre de voisins : "<<domains[nsd].getNbNeighbors()<<std::endl;
00792     //    domains[nsd].write();
00793     for (k = 0; k < domains[nsd].getNbNeighbors(); ++k) {
00794       neighbor = domains[nsd].getNeighborIdx(k);
00795       if (neighbor < nsd) {continue; }
00796       //      std::cout << "              neighbor : " << neighbor << std::endl;
00797       domains[nsd].addAtomsInBuffers(system.getAtoms(), domains[neighbor], k,
00798                                      this->radius1, this->radius2, this->useFrag);
00799     }
00800     //    domains[nsd].write();
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   // gamma = facteur de réduction :
00871   gamma       = (totalWideSD - removeLength ) / totalTmpWideSD;
00872   totalWideSD = 0.0 ; 
00873   diff = 0.0; 
00874   errorLoad = 0 ;
00875 
00876   // Les nouveaux wideSD[i] :
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   // Free memory    
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]) { // on permute
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   //  std::cout <<"   System Renumerote " << std::endl << system <<std::endl;
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

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