ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/devel_omp/src/parallel/ForceMatrixDecomposition.cpp
(Generate patch)

Comparing branches/devel_omp/src/parallel/ForceMatrixDecomposition.cpp (file contents):
Revision 1598 by mciznick, Wed Jul 27 14:26:53 2011 UTC vs.
Revision 1614 by mciznick, Tue Aug 23 20:55:51 2011 UTC

# Line 50 | Line 50 | ForceMatrixDecomposition::ForceMatrixDecomposition(Sim
50  
51   ForceMatrixDecomposition::ForceMatrixDecomposition(SimInfo* info, InteractionManager* iMan) :
52          ForceDecomposition(info, iMan) {
53 –
53          // In a parallel computation, row and colum scans must visit all
54          // surrounding cells (not just the 14 upper triangular blocks that
55          // are used when the processor can see all pairs)
# Line 325 | Line 324 | void ForceMatrixDecomposition::distributeInitialData()
324                  }
325          }
326  
327 +        Globals* simParams_ = info_->getSimParams();
328 +        if (simParams_->haveNeighborListReorderFreq())
329 +        {
330 +                neighborListReorderFreq = simParams_->getNeighborListReorderFreq();
331 +        } else
332 +        {
333 +                neighborListReorderFreq = 0;
334 +        }
335 +        reorderFreqCounter = 1;
336 +
337          createGtypeCutoffMap();
338  
339   }
# Line 796 | Line 805 | void ForceMatrixDecomposition::collectData() {
805          pairwisePot += pot_temp[ii];
806   #endif
807  
808 <        cerr << "pairwisePot = " << pairwisePot << "\n";
808 >        //      cerr << "pairwisePot = " << pairwisePot << "\n";
809   }
810  
811   int ForceMatrixDecomposition::getNAtomsInRow() {
# Line 847 | Line 856 | Vector3d ForceMatrixDecomposition::getIntergroupVector
856          Vector3d d;
857  
858          d = snap_->cgData.position[cg2->getLocalIndex()] - snap_->cgData.position[cg1->getLocalIndex()];
859 < /*      cerr << "cg1_gid = " << cg1->getGlobalIndex() << "\tcg1_lid = " << cg1->getLocalIndex() << "\tcg1p = "
860 <                        << snap_->cgData.position[cg1->getLocalIndex()] << "\n";
861 <        cerr << "cg2_gid = " << cg2->getGlobalIndex() << "\tcg2_lid = " << cg2->getLocalIndex() << "\tcg2p = "
862 <                        << snap_->cgData.position[cg2->getLocalIndex()] << "\n";*/
859 >        /*      cerr << "cg1_gid = " << cg1->getGlobalIndex() << "\tcg1_lid = " << cg1->getLocalIndex() << "\tcg1p = "
860 >         << snap_->cgData.position[cg1->getLocalIndex()] << "\n";
861 >         cerr << "cg2_gid = " << cg2->getGlobalIndex() << "\tcg2_lid = " << cg2->getLocalIndex() << "\tcg2p = "
862 >         << snap_->cgData.position[cg2->getLocalIndex()] << "\n";*/
863  
864          snap_->wrapVector(d);
865          return d;
# Line 924 | Line 933 | bool ForceMatrixDecomposition::skipAtomPair(int atom1,
933   bool ForceMatrixDecomposition::skipAtomPair(int atom1, int atom2) {
934          int unique_id_1, unique_id_2;
935  
936 < //      cerr << "sap with atom1, atom2 =\t" << atom1 << "\t" << atom2 << "\n";
936 >        //      cerr << "sap with atom1, atom2 =\t" << atom1 << "\t" << atom2 << "\n";
937   #ifdef IS_MPI
938          // in MPI, we have to look up the unique IDs for each atom
939          unique_id_1 = AtomRowToGlobal[atom1];
# Line 992 | Line 1001 | void ForceMatrixDecomposition::fillInteractionData(Int
1001  
1002   // filling interaction blocks with pointers
1003   void ForceMatrixDecomposition::fillInteractionData(InteractionData &idat, int atom1, int atom2) {
1004 +
1005 +        idat.excluded = excludeAtomPair(atom1, atom2);
1006 +
1007 + #ifdef IS_MPI
1008 +        idat.atypes = make_pair( atypesRow[atom1], atypesCol[atom2]);
1009 +        //idat.atypes = make_pair( ff_->getAtomType(identsRow[atom1]),
1010 +        //                         ff_->getAtomType(identsCol[atom2]) );
1011 +
1012 +        if (storageLayout_ & DataStorage::dslAmat)
1013 +        {
1014 +                idat.A1 = &(atomRowData.aMat[atom1]);
1015 +                idat.A2 = &(atomColData.aMat[atom2]);
1016 +        }
1017 +
1018 +        if (storageLayout_ & DataStorage::dslElectroFrame)
1019 +        {
1020 +                idat.eFrame1 = &(atomRowData.electroFrame[atom1]);
1021 +                idat.eFrame2 = &(atomColData.electroFrame[atom2]);
1022 +        }
1023 +
1024 +        if (storageLayout_ & DataStorage::dslTorque)
1025 +        {
1026 +                idat.t1 = &(atomRowData.torque[atom1]);
1027 +                idat.t2 = &(atomColData.torque[atom2]);
1028 +        }
1029 +
1030 +        if (storageLayout_ & DataStorage::dslDensity)
1031 +        {
1032 +                idat.rho1 = &(atomRowData.density[atom1]);
1033 +                idat.rho2 = &(atomColData.density[atom2]);
1034 +        }
1035 +
1036 +        if (storageLayout_ & DataStorage::dslFunctional)
1037 +        {
1038 +                idat.frho1 = &(atomRowData.functional[atom1]);
1039 +                idat.frho2 = &(atomColData.functional[atom2]);
1040 +        }
1041 +
1042 +        if (storageLayout_ & DataStorage::dslFunctionalDerivative)
1043 +        {
1044 +                idat.dfrho1 = &(atomRowData.functionalDerivative[atom1]);
1045 +                idat.dfrho2 = &(atomColData.functionalDerivative[atom2]);
1046 +        }
1047 +
1048 +        if (storageLayout_ & DataStorage::dslParticlePot)
1049 +        {
1050 +                idat.particlePot1 = &(atomRowData.particlePot[atom1]);
1051 +                idat.particlePot2 = &(atomColData.particlePot[atom2]);
1052 +        }
1053 +
1054 +        if (storageLayout_ & DataStorage::dslSkippedCharge)
1055 +        {
1056 +                idat.skippedCharge1 = &(atomRowData.skippedCharge[atom1]);
1057 +                idat.skippedCharge2 = &(atomColData.skippedCharge[atom2]);
1058 +        }
1059 +
1060 + #else
1061 +
1062 +        idat.atypes = make_pair(atypesLocal[atom1], atypesLocal[atom2]);
1063 +        //idat.atypes = make_pair( ff_->getAtomType(idents[atom1]),
1064 +        //                         ff_->getAtomType(idents[atom2]) );
1065 +
1066 +        if (storageLayout_ & DataStorage::dslAmat)
1067 +        {
1068 +                idat.A1 = &(snap_->atomData.aMat[atom1]);
1069 +                idat.A2 = &(snap_->atomData.aMat[atom2]);
1070 +        }
1071 +
1072 +        if (storageLayout_ & DataStorage::dslElectroFrame)
1073 +        {
1074 +                idat.eFrame1 = &(snap_->atomData.electroFrame[atom1]);
1075 +                idat.eFrame2 = &(snap_->atomData.electroFrame[atom2]);
1076 +        }
1077 +
1078 +        if (storageLayout_ & DataStorage::dslTorque)
1079 +        {
1080 +                idat.t1 = &(snap_->atomData.torque[atom1]);
1081 +                idat.t2 = &(snap_->atomData.torque[atom2]);
1082 +        }
1083 +
1084 +        if (storageLayout_ & DataStorage::dslDensity)
1085 +        {
1086 +                idat.rho1 = &(snap_->atomData.density[atom1]);
1087 +                idat.rho2 = &(snap_->atomData.density[atom2]);
1088 +        }
1089  
1090 +        if (storageLayout_ & DataStorage::dslFunctional)
1091 +        {
1092 +                idat.frho1 = &(snap_->atomData.functional[atom1]);
1093 +                idat.frho2 = &(snap_->atomData.functional[atom2]);
1094 +        }
1095 +
1096 +        if (storageLayout_ & DataStorage::dslFunctionalDerivative)
1097 +        {
1098 +                idat.dfrho1 = &(snap_->atomData.functionalDerivative[atom1]);
1099 +                idat.dfrho2 = &(snap_->atomData.functionalDerivative[atom2]);
1100 +        }
1101 +
1102 +        if (storageLayout_ & DataStorage::dslParticlePot)
1103 +        {
1104 +                idat.particlePot1 = &(snap_->atomData.particlePot[atom1]);
1105 +                idat.particlePot2 = &(snap_->atomData.particlePot[atom2]);
1106 +        }
1107 +
1108 +        if (storageLayout_ & DataStorage::dslSkippedCharge)
1109 +        {
1110 +                idat.skippedCharge1 = &(snap_->atomData.skippedCharge[atom1]);
1111 +                idat.skippedCharge2 = &(snap_->atomData.skippedCharge[atom2]);
1112 +        }
1113 + #endif
1114 + }
1115 +
1116 + // filling interaction blocks with pointers
1117 + void ForceMatrixDecomposition::fillInteractionDataOMP(InteractionDataPrv &idat, int atom1, int atom2) {
1118 +
1119          idat.excluded = excludeAtomPair(atom1, atom2);
1120  
1121   #ifdef IS_MPI
# Line 1120 | Line 1243 | void ForceMatrixDecomposition::unpackInteractionData(I
1243  
1244   }
1245  
1246 + void ForceMatrixDecomposition::unpackInteractionDataOMP(InteractionDataPrv &idat, int atom1, int atom2) {
1247 +                pairwisePot += idat.pot;
1248 +
1249 +                snap_->atomData.force[atom1] += idat.f1;
1250 +                snap_->atomData.force[atom2] -= idat.f1;
1251 + }
1252 +
1253   void ForceMatrixDecomposition::reorderGroupCutoffs(vector<int> &order) {
1254          vector<int> tmp = vector<int> (groupToGtype.size());
1255  
# Line 1144 | Line 1274 | void ForceMatrixDecomposition::reorderPosition(vector<
1274                  tmp[i] = snap_->cgData.position[i];
1275          }
1276  
1277 <        vector<int> mapPos = vector<int>(nGroups_);
1277 >        vector<int> mapPos = vector<int> (nGroups_);
1278          for (int i = 0; i < nGroups_; ++i)
1279          {
1280                  snap_->cgData.position[i] = tmp[order[i]];
# Line 1163 | Line 1293 | void ForceMatrixDecomposition::reorderPosition(vector<
1293                          cg->setLocalIndex(mapPos[cg->getLocalIndex()]);
1294                  }
1295          }
1166 –
1167 – /*      if (info_->getNCutoffGroups() > 0)
1168 –        {
1169 –                for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi))
1170 –                {
1171 –                        for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci))
1172 –                        {
1173 –                                printf("gbI:%d locI:%d x:%f y:%f z:%f\n", cg->getGlobalIndex(), cg->getLocalIndex(),
1174 –                                                cgConfig->position[cg->getLocalIndex()].x(), cgConfig->position[cg->getLocalIndex()].y(),
1175 –                                                cgConfig->position[cg->getLocalIndex()].z());
1176 –                        }
1177 –                }
1178 –        } else
1179 –        {
1180 –                // center of mass of the group is the same as position of the atom
1181 –                // if cutoff group does not exist
1182 –                printf("ERROR!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
1183 –                //                      cgConfig->position = config->position;
1184 –        }*/
1296   }
1297  
1298   void ForceMatrixDecomposition::reorderGroupList(vector<int> &order) {
# Line 1200 | Line 1311 | void ForceMatrixDecomposition::reorderMemory(vector<ve
1311  
1312   void ForceMatrixDecomposition::reorderMemory(vector<vector<CutoffGroup *> > &H_c_l) {
1313          int n = 0;
1203 – //      printf("Reorder memory time:%f!!!!!!!!!!!!!!!!!!!!!!!!!!!\n",
1204 – //                      info_->getSnapshotManager()->getCurrentSnapshot()->getTime());
1314  
1315          /* record the reordered atom indices */
1316          vector<int> k = vector<int> (nGroups_);
# Line 1222 | Line 1331 | vector<vector<CutoffGroup *> > ForceMatrixDecompositio
1331   }
1332  
1333   vector<vector<CutoffGroup *> > ForceMatrixDecomposition::buildLayerBasedNeighborList() {
1225 – //      printf("buildLayerBasedNeighborList; nGroups:%d\n", nGroups_);
1334          // Na = nGroups_
1335          /* cell occupancy counter */
1336 < //      vector<int> k_c;
1336 >        //      vector<int> k_c;
1337          /* c_i - has cell containing atom i (size Na) */
1338          vector<int> c = vector<int> (nGroups_);
1339          /* l_i - layer containing atom i (size Na) */
1340 < //      vector<int> l;
1340 >        //      vector<int> l;
1341  
1342          RealType rList_ = (largestRcut_ + skinThickness_);
1343          Snapshot* snap_ = sman_->getCurrentSnapshot();
# Line 1248 | Line 1356 | vector<vector<CutoffGroup *> > ForceMatrixDecompositio
1356          int cellIndex;
1357          int nCtot = nCells_.x() * nCells_.y() * nCells_.z();
1358  
1359 < //      k_c = vector<int> (nCtot, 0);
1359 >        //      k_c = vector<int> (nCtot, 0);
1360  
1361          SimInfo::MoleculeIterator mi;
1362          Molecule* mol;
# Line 1277 | Line 1385 | vector<vector<CutoffGroup *> > ForceMatrixDecompositio
1385                          whichCell.y() = nCells_.y() * scaled.y();
1386                          whichCell.z() = nCells_.z() * scaled.z();
1387  
1388 < //                      printf("pos x:%f y:%f z:%f cell x:%d y:%d z:%d\n", rs.x(), rs.y(), rs.z(), whichCell.x(), whichCell.y(),
1389 < //                                      whichCell.z());
1388 >                        //                      printf("pos x:%f y:%f z:%f cell x:%d y:%d z:%d\n", rs.x(), rs.y(), rs.z(), whichCell.x(), whichCell.y(),
1389 >                        //                                      whichCell.z());
1390  
1391                          // find single index of this cell:
1392                          cellIndex = Vlinear(whichCell, nCells_);
# Line 1287 | Line 1395 | vector<vector<CutoffGroup *> > ForceMatrixDecompositio
1395                  }
1396          }
1397  
1398 < //      int k_c_curr;
1399 < //      int k_c_max = 0;
1398 >        //      int k_c_curr;
1399 >        //      int k_c_max = 0;
1400          /* the cell-layer occupancy matrix */
1401          vector<vector<CutoffGroup *> > H_c_l = vector<vector<CutoffGroup *> > (nCtot);
1402  
# Line 1305 | Line 1413 | vector<vector<CutoffGroup *> > ForceMatrixDecompositio
1413                          //                      {
1414                          //                              k_c_max = k_c_curr;
1415                          //                      }
1308 –
1416                          H_c_l[c[cg->getGlobalIndex()]].push_back(/*l[*/cg/*]*/);
1417                  }
1418          }
1419  
1420 <        reorderMemory(H_c_l);
1420 >        /* Frequency of reordering the memory */
1421 >        if (neighborListReorderFreq != 0)
1422 >        {
1423 >                if (reorderFreqCounter == neighborListReorderFreq)
1424 >                {
1425 >                        reorderMemory(H_c_l);
1426 >                        reorderFreqCounter = 1;
1427 >                } else
1428 >                {
1429 >                        reorderFreqCounter++;
1430 >                }
1431 >        }
1432  
1433          int m;
1434 <        /* the neighbor matrix */
1434 >        /* The neighbor matrix */
1435          vector<vector<CutoffGroup *> > neighborMatW = vector<vector<CutoffGroup *> > (nGroups_);
1436  
1437          groupCutoffs cuts;
1438          CutoffGroup *cg1;
1439  
1440 <        /* loops over objects(atoms, rigidBodies, cutoffGroups, etc.) */
1440 >        /* Loops over objects(atoms, rigidBodies, cutoffGroups, etc.) */
1441          for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi))
1442          {
1443                  for (cg1 = mol->beginCutoffGroup(ci); cg1 != NULL; cg1 = mol->nextCutoffGroup(ci))
# Line 1358 | Line 1476 | vector<vector<CutoffGroup *> > ForceMatrixDecompositio
1476                                  }
1477  
1478                                  int c2 = Vlinear(c2v, nCells_);
1479 <                                /* loops over layers l to access the neighbor atoms */
1479 >                                /* Loops over layers l to access the neighbor atoms */
1480                                  for (vector<CutoffGroup *>::iterator cg2 = H_c_l[c2].begin(); cg2 != H_c_l[c2].end(); ++cg2)
1481                                  {
1482 <                                        //                              if i'' = 0 then break // doesn't apply to vector implementation of matrix
1482 >                                        //                              if i'' = 0 then break // doesn't apply to vector implementation of the matrix
1483                                          //                              if(i != *j)
1484                                          if (c2 != c1 || (*cg2)->getGlobalIndex() < cg1->getGlobalIndex())
1485                                          {
# Line 1370 | Line 1488 | vector<vector<CutoffGroup *> > ForceMatrixDecompositio
1488                                                  cuts = getGroupCutoffs(cg1->getGlobalIndex(), (*cg2)->getGlobalIndex());
1489                                                  if (dr.lengthSquare() < cuts.third)
1490                                                  {
1491 <                                                        /* transposed version of Rapaport W mat, to occupy successive memory locations on CPU */
1491 >                                                        /* Transposed version of Rapaport W mat, to occupy successive memory locations on CPU */
1492                                                          neighborMatW[cg1->getGlobalIndex()].push_back((*cg2));
1493                                                  }
1494                                          }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines