Index: configure.in =================================================================== RCS file: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/configure.in,v retrieving revision 1.121.2.7 retrieving revision 1.132 diff -u -r1.121.2.7 -r1.132 --- configure.in 16 Aug 2004 10:54:32 -0000 1.121.2.7 +++ configure.in 19 Aug 2004 12:17:40 -0000 1.132 @@ -164,8 +164,8 @@ OL_ARG_WITH(metis,[ --with-metis with Metis model partitioning support],auto,[auto yes no]) dnl temporarily disabled until we figure out how to make chaco work... -dnl OL_ARG_WITH(chaco,[ --with-chaco with Chaco model partitioning support],auto,[auto yes no]) -ol_with_chaco=no +OL_ARG_WITH(chaco,[ --with-chaco with Chaco model partitioning support],auto,[auto yes no]) +dnl ol_with_chaco=no OL_ARG_WITH(threads,[ --with-threads with threads],auto,[auto yes no]) OL_ARG_WITH(rtai,[ --with-rtai with RTAI support],no,[auto yes no]) @@ -360,7 +361,7 @@ ;; esac if test "$ol_with_metis" = no -a "$ol_with_chaco" = no ; then - AC_MSG_ERROR([Schur parallel solver needs either Metis or Chaco]) + AC_MSG_WARN([No partitioning library available; Schur parallel solver will allow manual partitioning only...]) fi fi Index: mbdyn/base/chaco_interface.cc =================================================================== RCS file: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/base/chaco_interface.cc,v retrieving revision 1.2 retrieving revision 1.3 diff -u -r1.2 -r1.3 --- mbdyn/base/chaco_interface.cc 12 Oct 2003 21:14:23 -0000 1.2 +++ mbdyn/base/chaco_interface.cc 19 Aug 2004 12:17:54 -0000 1.3 @@ -20,68 +20,147 @@ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ +#ifdef HAVE_CONFIG_H +#ifdef HAVE_CONFIG_H +#include /* This goes first in every *.c,*.cc file */ +#endif /* HAVE_CONFIG_H */ + +#include /* This goes first in every *.c,*.cc file */ +#endif /* HAVE_CONFIG_H */ + #include #include -using namespace std; - extern "C" { int interface( -int nvtxs, /* number of vertices in full graph */ -int *start, /* start of edge list for each vertex */ -int *adjacency, /* edge list data */ -int *vwgts, /* weights for all vertices */ -float *ewgts, /* weights for all edges */ -float *x, float *y, float *z,/* coordinates for inertial method */ -char *outassignname, /* name of assignment output file */ -char *outfilename, /* output file name */ -short *assignment, /* set number of each vtx (length n) */ -int architecture, /* 0 => hypercube, d => d-dimensional mesh */ -int ndims_tot, /* total number of cube dimensions to divide */ -int mesh_dims[3], /* dimensions of mesh of processors */ -double *goal, /* desired set sizes for each set */ -int global_method, /* global partitioning algorithm */ -int local_method, /* local partitioning algorithm */ -int rqi_flag, /* should I use RQI/Symmlq eigensolver? */ -int vmax, /* how many vertices to coarsen down to? */ -int ndims, /* number of eigenvectors (2^d sets) */ -double eigtol, /* tolerance on eigenvectors */ -long seed /* for random graph mutations */ -); + int nvtxs, /* number of vertices in full graph */ + int *start, /* start of edge list for each vertex */ + int *adjacency, /* edge list data */ + int *vwgts, /* weights for all vertices */ + float *ewgts, /* weights for all edges */ + float *x, float *y, float *z,/* coordinates for inertial method */ + char *outassignname, /* name of assignment output file */ + char *outfilename, /* output file name */ + short *assignment, /* set number of each vtx (length n) */ + int architecture, /* 0 => hypercube, d => d-dimensional mesh */ + int ndims_tot, /* total number of cube dimensions to divide */ + int mesh_dims[3], /* dimensions of mesh of processors */ + double *goal, /* desired set sizes for each set */ + int global_method, /* global partitioning algorithm */ + int local_method, /* local partitioning algorithm */ + int rqi_flag, /* should I use RQI/Symmlq eigensolver? */ + int vmax, /* how many vertices to coarsen down to? */ + int ndims, /* number of eigenvectors (2^d sets) */ + double eigtol, /* tolerance on eigenvectors */ + long seed /* for random graph mutations */ + ); } -void chaco_interface(const int iTotVertices, int *start, int *adjacency, - int *vertex_weights, float *edge_weights, - const int num_processors, int *pParAmgProcs) +extern "C" int FREE_GRAPH; + +void +chaco_interface( + const int iTotVertices, + int *start, + int *adjacency, + int *vertex_weights, + int *comm_weights, /* unsupported */ + int *edge_weights, + const int num_processors, + int *pParAmgProcs + ) { #ifdef USE_CHACO - vector set_assignment(iTotVertices); + std::vector set_assignment(iTotVertices); - const int architecture(1), /* The computers are connected as a - simple one-dimensional mesh */ - global_method(2), /* Use Spectral decomposition */ - local_method(1), /* With KL local refinement */ - rqi_flag(0), /* Use Lanczos for solving the - eigenvalue problem for spectral - decomposition */ - vmax(100), /* When to stop coarsening. Not used */ - seed(7654321); /* A random number seed */ - int mesh_dims[]={1,1,1}; - - double eigtol(1e-3); - - - if(interface(iTotVertices,start,adjacency,vertex_weights,edge_weights, - NULL,NULL,NULL,NULL,NULL,&set_assignment[0],architecture, - num_processors,mesh_dims,NULL,global_method,local_method,rqi_flag, - vmax,1,eigtol,seed)!=0) - { - cerr << "Partition failed" << endl; - } - - for(int i=0;i hypercube, d => d-dimensional mesh */ + global_method(2), /* Use Spectral decomposition */ + local_method(1), /* With KL local refinement */ + rqi_flag(0), /* Use Lanczos for solving the + * eigenvalue problem for spectral + * decomposition */ + vmax(100), /* When to stop coarsening. Not used */ + seed(7654321); /* A random number seed */ + int mesh_dims[] = { /* if architecture > 0, mesh size goes here */ + num_processors, 1, 1 + }; + double eigtol(1.e-3); /* */ + + /* Chaco vertices are base 1 */ + for (int i = 0; i < start[iTotVertices]; i++) { + adjacency[i]++; + } + + /* weights = 0 are not allowed */ + std::vector ivertex_weights(iTotVertices); + if (vertex_weights) { + for (int i = 0; i < iTotVertices; i++) { + /* trying to emulate communication weights? */ + if (comm_weights) { + vertex_weights[i] += comm_weights[i]; + } + + if (vertex_weights[i] == 0) { + ivertex_weights[i] = 1; + } else { + ivertex_weights[i] = vertex_weights[i]; + } + } + } + + /* Chaco uses floats for the communication weights */ + std::vector fedge_weights(2*start[iTotVertices]); + if (edge_weights) { + for (int i = 0; i < iTotVertices; i++) { + fedge_weights[i] = edge_weights[i]; + fedge_weights[start[iTotVertices] + i] = edge_weights[i]; + } + } + + /* NOTE: don't free() adjacency !!! (default in Chaco 2.2) */ + if (FREE_GRAPH) { + FREE_GRAPH = 0; + } + + int rc = interface( + iTotVertices, /* number of vertices in full graph */ + start, /* start of edge list for each vertex */ + adjacency, /* edge list data */ + vertex_weights ? &ivertex_weights[0] : 0, /* weights for all vertices */ + edge_weights ? &fedge_weights[0] : 0, /* weights for all edges */ + NULL, /* coordinates for inertial method */ + NULL, /* ... */ + NULL, /* ... */ + NULL, /* name of assignment output file */ + NULL, /* output file name */ + &set_assignment[0], /* set number of each vtx (length n) */ + architecture, /* 0 => hypercube, d => d-dimensional mesh */ + num_processors, /* total number of cube dimensions to divide */ + mesh_dims, /* dimensions of mesh of processors */ + NULL, /* desired set sizes for each set */ + global_method, /* global partitioning algorithm */ + local_method, /* local partitioning algorithm */ + rqi_flag, /* should I use RQI/Symmlq eigensolver? */ + vmax, /* how many vertices to coarsen down to? */ + 1, /* number of eigenvectors (2^d sets) */ + eigtol, /* tolerance on eigenvectors */ + seed /* for random graph mutations */ + ); + + /* Chaco vertices are base 1 */ + for (int i = 0; i < start[iTotVertices]; i++) { + adjacency[i]--; + } + + if (rc != 0) { + silent_cerr("chaco_interface(): partition failed" << std::endl); + } + + for (int i = 0; i < iTotVertices; i++) { + pParAmgProcs[i] = set_assignment[i]; + } #endif /* USE_CHACO */ } + Index: mbdyn/base/chaco_interface.h =================================================================== RCS file: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/base/chaco_interface.h,v retrieving revision 1.2 retrieving revision 1.3 diff -u -r1.2 -r1.3 --- mbdyn/base/chaco_interface.h 12 Oct 2003 21:14:23 -0000 1.2 +++ mbdyn/base/chaco_interface.h 19 Aug 2004 12:17:54 -0000 1.3 @@ -23,9 +23,17 @@ #ifndef CHACO_INTERFACE #define CHACO_INTERFACE -void chaco_interface(const int iTotVertices, int *start, int *adjacency, - int *vertex_weights, float *edge_weights, - const int num_processors, int *pParAmgProcs); +extern void +chaco_interface( + const int iTotVertices, + int *start, + int *adjacency, + int *vertex_weights, + int *comm_weights, + int *edge_weights, + const int num_processors, + int *pParAmgProcs + ); #endif /* CHACO_INTERFACE */ Index: mbdyn/base/dataman.cc =================================================================== RCS file: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/base/dataman.cc,v retrieving revision 1.42.2.5 retrieving revision 1.55 diff -u -r1.42.2.5 -r1.55 --- mbdyn/base/dataman.cc 12 Aug 2004 13:17:45 -0000 1.42.2.5 +++ mbdyn/base/dataman.cc 19 Aug 2004 12:17:54 -0000 1.55 @@ -370,9 +370,13 @@ } ReadElems(HP); +#if 0 + /* FIXME: we don't check extra statements after "end: elements;" + * because there might be more... */ try { CurrDesc = KeyWords(HP.GetDescription()); } catch (EndOfFile) {} +#endif } else { DEBUGCERR(""); silent_cerr("warning, no elements are defined" << std::endl); Index: mbdyn/base/metiswrap.cc =================================================================== RCS file: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/base/metiswrap.cc,v retrieving revision 1.8 retrieving revision 1.9 diff -u -r1.8 -r1.9 --- mbdyn/base/metiswrap.cc 3 Aug 2004 18:12:28 -0000 1.8 +++ mbdyn/base/metiswrap.cc 19 Aug 2004 12:17:54 -0000 1.9 @@ -39,6 +39,8 @@ } #undef ASSERT /* kill annoying redefiniton message */ +#include + #endif /* !USE_METIS */ int @@ -47,7 +49,7 @@ int *pAdjncy, int *pVertexWgts, int *pCommWgts, - int wgtflag, + int *pEdgeWgts, int DataCommSize, int *pParAmgProcs) { @@ -57,25 +59,104 @@ int numflag = 0; /* if options[0] == 0, the rest is ignored */ int options[5] = { 0 }; - /* number of edges cut */ - int edgecut = 0; + /* total communication volume */ + int volume = 0; + /* weiht flags */ + int wgtflag; + + idxtype *xadj = 0, + *adjncy = 0, + *vwgt = 0, + *adjwgt = 0, + *part = 0; + + if (sizeof(idxtype) == sizeof(int)) { + xadj = pXadj; + adjncy = pAdjncy; + vwgt = pVertexWgts; + adjwgt = pCommWgts; + part = pParAmgProcs; + + } else { + SAFENEWARR(xadj, idxtype, iTotVertices); + SAFENEWARR(adjncy, idxtype, pXadj[iTotVertices]); + if (pVertexWgts) { + SAFENEWARR(vwgt, idxtype, iTotVertices); + } + if (pCommWgts) { + SAFENEWARR(adjwgt, idxtype, pXadj[iTotVertices]); + } + SAFENEWARR(part, idxtype, iTotVertices); + + for (int i = 0; i < iTotVertices; i++) { + xadj[i] = pXadj[i]; + part[i] = pParAmgProcs[i]; + } + + if (pVertexWgts) { + for (int i = 0; i < iTotVertices; i++) { + vwgt[i] = pVertexWgts[i]; + } + } + + for (int i = 0; i < pXadj[iTotVertices]; i++) { + adjncy[i] = pAdjncy[i]; + } + + if (pCommWgts) { + for (int i = 0; i < pXadj[iTotVertices]; i++) { + adjwgt[i] = pCommWgts[i]; + } + } + } + + if (pVertexWgts && pCommWgts) { + wgtflag = 3; + } else if (pVertexWgts) { + wgtflag = 2; + } else if (pCommWgts) { + wgtflag = 1; + } else { + wgtflag = 0; + } METIS_PartGraphVKway(&iTotVertices, - pXadj, - pAdjncy, - pVertexWgts, - pCommWgts, + xadj, + adjncy, + vwgt, + adjwgt, &wgtflag, &numflag, &DataCommSize, options, - &edgecut, - pParAmgProcs); + &volume, + part); + + + if (sizeof(idxtype) != sizeof(int)) { + for (int i = 0; i < iTotVertices; i++) { + pParAmgProcs[i] = part[i]; + } + + SAFEDELETEARR(xadj); + SAFEDELETEARR(adjncy); + + if (pVertexWgts) { + SAFEDELETEARR(vwgt); + } + + if (pCommWgts) { + SAFEDELETEARR(adjwgt); + } + + SAFEDELETEARR(part); + } /* NOTE: the manual suggests to use * METIS_PartGraphRecursive if DataCommSize < 8 */ -#endif /* !USE_METIS */ +#endif /* USE_METIS */ + return 0; } Index: mbdyn/base/metiswrap.h =================================================================== RCS file: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/base/metiswrap.h,v retrieving revision 1.4 retrieving revision 1.5 diff -u -r1.4 -r1.5 --- mbdyn/base/metiswrap.h 3 Jan 2004 14:01:12 -0000 1.4 +++ mbdyn/base/metiswrap.h 19 Aug 2004 12:17:54 -0000 1.5 @@ -37,7 +37,7 @@ int *pAdjncy, int *pVertexWgts, int *pCommWgts, - int wgtflag, + int *pEdgeWgts, int DataCommSize, int *pParAmgProcs); Index: mbdyn/base/schurdataman.cc =================================================================== RCS file: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/base/schurdataman.cc,v retrieving revision 1.40.2.1 retrieving revision 1.42 diff -u -r1.40.2.1 -r1.42 --- mbdyn/base/schurdataman.cc 16 Aug 2004 10:54:33 -0000 1.40.2.1 +++ mbdyn/base/schurdataman.cc 19 Aug 2004 12:17:54 -0000 1.42 @@ -51,6 +51,7 @@ #undef min #undef max +#include #include #include @@ -188,7 +189,7 @@ iNumMyInt(0), pMyIntDofs(NULL), pLabelsList(NULL), -wgtflag(1), +wgtflag(WEIGHT_VERTICES), pParAmgProcs(NULL), Partitioner(PARTITIONER_DEFAULT), pRotorComm(NULL), @@ -206,7 +207,7 @@ DEBUGCOUT("Communicator Size: " << DataCommSize << std::endl); iTotVertices = iTotNodes + iTotElem; - DEBUGCOUT("iTotVertexes: " << iTotVertices << std::endl); + DEBUGCOUT("iTotVertices: " << iTotVertices << std::endl); /* parole chiave del blocco parallelizzazione */ const char* sKeyWords[] = { @@ -282,13 +283,13 @@ try { if (KeyWords(HP.GetDescription()) != BEGIN) { - pedantic_cerr("no explicit connection declared " + pedantic_cerr("no explicit connections declared " "for this input file" << std::endl); return; } } catch (EndOfFile) { - pedantic_cerr("no explicit connection declared " + pedantic_cerr("no explicit connections declared " "for this input file" << std::endl); return; } @@ -296,8 +297,8 @@ int iNumElems = 0; int iNumNodes = 0; if (KeyWords(HP.GetWord()) != PARALLEL) { - silent_cerr("Error: expected at line " - << HP.GetLineData() << "; aborting..." << std::endl); + silent_cerr("Error: \"begin: parallel;\" expected at line " + << HP.GetLineData() << std::endl); throw ErrGeneric(); } @@ -307,23 +308,28 @@ if (!HP.IsArg()) { silent_cerr("Error: Weight flag expected " "at line " << HP.GetLineData() - << "; aborting ..." << std::endl); + << std::endl); throw ErrGeneric(); } - if (HP.IsKeyWord("no" "weights")) { - wgtflag = 0; - } else if (HP.IsKeyWord("edges" "only")) { - wgtflag = 1; - } else if (HP.IsKeyWord("vertices" "only")) { - wgtflag = 2; - } else if (HP.IsKeyWord("vertices" "and" "edges")) { - wgtflag = 3; - } else { - wgtflag = HP.GetInt(); - if (wgtflag < 0 || wgtflag > 3) { - silent_cerr("invalid weights " << wgtflag - << " at line " << HP.GetLineData() + while (HP.IsArg()) { + if (HP.IsKeyWord("none")) { + wgtflag = WEIGHT_NONE; + + } else if (HP.IsKeyWord("vertices") + || HP.IsKeyWord("computation")) + { + wgtflag |= WEIGHT_VERTICES; + + } else if (HP.IsKeyWord("communication")) { + wgtflag |= WEIGHT_COMM; + + } else if (HP.IsKeyWord("edges")) { + wgtflag |= WEIGHT_EDGES; + + } else { + silent_cerr("invalid weight " + " at line " << HP.GetLineData() << std::endl); throw ErrGeneric(); } @@ -334,16 +340,27 @@ SAFENEWARR(pParAmgProcs, int, iTotVertices); for (int i = 0; i < iTotVertices; i++) { if (!HP.IsArg()) { - silent_cerr("Error: the partition " - "assignment is not complete, " + silent_cerr("the partition " + "assignment is not complete; " "only " << i << " vertices " "input so far " "at line " << HP.GetLineData() - << "; aborting ..." << std::endl); throw ErrGeneric(); } pParAmgProcs[i] = HP.GetInt(); + if (pParAmgProcs[i] < 0 || + pParAmgProcs[i] >= DataCommSize) + { + silent_cerr("illegal value " + << pParAmgProcs[i] + << " for partition " + "assignment[" << i << "] " + "at line " << HP.GetLineData() + << "; must be between 0 and " + << DataCommSize - 1 << std::endl); + throw ErrGeneric(); + } } break; @@ -363,10 +380,13 @@ silent_cerr("CHACO partitioner not available; aborting..." << std::endl); throw ErrGeneric(); #endif /* ! USE_CHACO */ + } else if (HP.IsKeyWord("manual")) { + Partitioner = PARTITIONER_MANUAL; } else { - silent_cerr("unknown partitioner at line " << HP.GetLineData() - << "; aborting..." << std::endl); + silent_cerr("unknown partitioner " + "at line " << HP.GetLineData() + << std::endl); throw ErrGeneric(); } break; @@ -672,12 +692,13 @@ { DEBUGCOUT("Entering SchurDataManager::CreatePartition()" << std::endl); + std::vector > adj(iTotVertices); /* adjacency */ + std::vector CommWgts(iTotVertices); /* communication weights */ + Adjacency Vertices; /* Struttura contenente le connessioni fra i vertici */ - int* pVertexWgts = 0; /* Pesi dei vertici = dofs x ogni v. utile per METIS */ - int* pCommWgts = 0; + std::vector VertexWgts(iTotVertices); /* Pesi dei vertici = dofs x ogni v. utile per METIS */ Vertices.pXadj = 0; Vertices.pAdjncy = 0; - pVertexWgts = 0; integer iMax = 0; integer iRMax = 0; int iCount = 0; @@ -688,11 +709,7 @@ /* Costruisco e inizializzo le due strutture */ SAFENEWARR(Vertices.pXadj, int, iTotVertices + 1); - SAFENEWARR(pVertexWgts, int, iTotVertices*2); - SAFENEWARR(pCommWgts, int, iTotVertices); memset(Vertices.pXadj, 0, (iTotVertices + 1)*sizeof(int)); - memset(pVertexWgts, 0, iTotVertices*2*sizeof(int)); - memset(pCommWgts, 0, iTotVertices*sizeof(int)); /* Ciclo per la scrittura degli array delle connessioni. * Il ciclo viene ripetuto se per un vertice si ha un numero @@ -721,28 +738,26 @@ memset(pLabelsList, 0, iTotNodes*sizeof(unsigned int)); while (true) { - InitList(Vertices.pXadj, iTotVertices+1, 0); - InitList(pVertexWgts, iTotVertices*2, 0); - InitList(pCommWgts, iTotVertices, 0); + InitList(Vertices.pXadj, iTotVertices + 1, 0); SAFENEWARR(Vertices.pAdjncy, int, iTotVertices*iMaxConnectionsPerVertex); InitList(Vertices.pAdjncy, iTotVertices*iMaxConnectionsPerVertex, ADJ_UNDEFINED); ASSERT(ppElems != NULL); + for (unsigned int i = 0; i < iTotNodes; i++) { + pLabelsList[i] = ppNodes[i]->GetLabel(); + } + /* ciclo sugli elementi per assemblare la struttura delle connessioni */ Node** ppCurrNode = NULL; /* per numerare i nodi prendo la posizione del puntatore * al nodo nell'array ppNodes */ int position; iCount = iTotNodes; - for (unsigned int i = 0; i < iTotNodes; i++) { - pLabelsList[i] = ppNodes[i]->GetLabel(); - } - iNumRt = 0; for (Elem** ppTmpEl = ppElems; - ppTmpEl < ppElems+iTotElem; + ppTmpEl < ppElems + iTotElem; ppTmpEl++, iCount++) { if ((*ppTmpEl)->GetElemType() == Elem::GRAVITY) { @@ -764,13 +779,12 @@ /* peso dell'elemento */ integer dimA, dimB; (*ppTmpEl)->WorkSpaceDim(&dimA, &dimB); - pVertexWgts[iCount] = dimA * dimB; + VertexWgts[iCount] = dimA * dimB; - pCommWgts[iCount] = - (*ppTmpEl)->iGetNumDof()*(*ppTmpEl)->iGetNumDof(); + CommWgts[iCount] = (*ppTmpEl)->iGetNumDof()*(*ppTmpEl)->iGetNumDof(); - for (int i = 0; i <= iNumberOfNodes-1; i++) { - Vertices.pXadj[iCount+1] += 1; + for (int i = 0; i < iNumberOfNodes; i++) { + Vertices.pXadj[iCount + 1]++; /* trovo la pos. del nodo nella lista dei puntatori ai nodi */ ppCurrNode = SearchNode(NodeData[pMyTypes[i]].ppFirstNode, @@ -782,24 +796,24 @@ * di ciascun nodo connesso */ #if 0 /* FIXME: i nodi hanno peso comp. nullo */ - pVertexWgts[iCount] += (*ppCurrNode)->iGetNumDof(); + VertexWgts[iCount] += (*ppCurrNode)->iGetNumDof(); #endif /* aggiungo fra le connessioni dell'elemento il nodo attuale */ - if ((iCount*iMaxConnectionsPerVertex) + Vertices.pXadj[iCount+1] - 1 < iTotVertices*iMaxConnectionsPerVertex) { + if ((iCount*iMaxConnectionsPerVertex) + Vertices.pXadj[iCount + 1] - 1 < iTotVertices*iMaxConnectionsPerVertex) { Vertices.pAdjncy[(iCount*iMaxConnectionsPerVertex) - + Vertices.pXadj[iCount+1] - 1] = position; + + Vertices.pXadj[iCount + 1] - 1] = position; } /* aggiungo alle connessioni del nodo l'elemento attuale */ - Vertices.pXadj[position+1] += 1; - if ((position*iMaxConnectionsPerVertex) + Vertices.pXadj[position+1] - 1 < iTotVertices*iMaxConnectionsPerVertex) { + Vertices.pXadj[position + 1]++; + if ((position*iMaxConnectionsPerVertex) + Vertices.pXadj[position + 1] - 1 < iTotVertices*iMaxConnectionsPerVertex) { Vertices.pAdjncy[(position*iMaxConnectionsPerVertex) - + Vertices.pXadj[position+1] - 1] = iCount; + + Vertices.pXadj[position + 1] - 1] = iCount; } /* peso (di comunicazione) del nodo */ - pCommWgts[position] = (*ppCurrNode)->iGetNumDof(); + CommWgts[position] = (*ppCurrNode)->iGetNumDof(); } } @@ -818,12 +832,12 @@ } iElPos = j; - Vertices.pXadj[iNdPos+1] += 1; - Vertices.pXadj[iTotNodes + iElPos+1] += 1; + Vertices.pXadj[iNdPos + 1]++; + Vertices.pXadj[iTotNodes + iElPos + 1]++; Vertices.pAdjncy[(iNdPos*iMaxConnectionsPerVertex) - + Vertices.pXadj[iNdPos+1] - 1] = iElPos+iTotNodes; - Vertices.pAdjncy[((iTotNodes+iElPos)*iMaxConnectionsPerVertex) - + Vertices.pXadj[iTotNodes+iElPos+1] - 1] = iNdPos; + + Vertices.pXadj[iNdPos + 1] - 1] = iElPos + iTotNodes; + Vertices.pAdjncy[((iTotNodes + iElPos)*iMaxConnectionsPerVertex) + + Vertices.pXadj[iTotNodes + iElPos + 1] - 1] = iNdPos; } } @@ -844,7 +858,7 @@ } for (int i = 1; i <= iTotVertices; i++) { - Vertices.pXadj[i] += Vertices.pXadj[i-1]; + Vertices.pXadj[i] += Vertices.pXadj[i - 1]; } /* Compatta il vettore delle adiacenze */ @@ -856,85 +870,111 @@ SAFENEWARR(pParAmgProcs, int, iTotVertices); #ifdef DEBUG - ofstream ofMetis; if (MyRank == 0) { - ofMetis.open("metis_conn.debug"); - ofMetis << "# METIS Input File" << std::endl - << "Column 1 is for Computational weights " << std::endl - << "Column 2 is for Comunicational weight " << std::endl - << "Total Vertexes: " << iTotVertices << std::endl + std::ofstream ofPartition; + + ofPartition.open("partition.debug"); + ofPartition << "# METIS-like Input File" << std::endl + << "# Column 1: computational weight" << std::endl + << "# Column 2: communication weight" << std::endl + << "# Total Vertices: " << iTotVertices << std::endl << "# Nodes" << std::endl; - for (int i = 0; i < iTotNodes; i++) { - ofMetis << "# " << i << " Node Type: " + for (unsigned int i = 0; i < iTotNodes; i++) { + ofPartition << "# " << i << " Node Type: " << "(" << psNodeNames[ppNodes[i]->GetNodeType()] << ")" << " Label: " << ppNodes[i]->GetLabel() << std::endl - << pVertexWgts[i] << " " << pCommWgts[i]; - for (int j = Vertices.pXadj[i]; j < Vertices.pXadj[i+1]; j++) { - ofMetis << " " << Vertices.pAdjncy[j]; - } - ofMetis << std::endl; - } - ofMetis << "# Elements" << std::endl; - for (int i = 0; i < iTotElem; i++) { - ofMetis << "# " << i+iTotNodes << " Element Type: " + << VertexWgts[i] << " " << CommWgts[i]; + for (int j = Vertices.pXadj[i]; j < Vertices.pXadj[i + 1]; j++) { + ofPartition << " " << Vertices.pAdjncy[j]; + } + ofPartition << std::endl; + } + ofPartition << "# Elements" << std::endl; + for (unsigned int i = 0; i < iTotElem; i++) { + ofPartition << "# " << i + iTotNodes << " Element Type: " << "(" << psElemNames[ppElems[i]->GetElemType()] << ")" << " Label: " << ppElems[i]->GetLabel() << std::endl - << pVertexWgts[i+iTotNodes] - << " " << pCommWgts[i+iTotNodes]; - for (int j = Vertices.pXadj[i+iTotNodes]; - j < Vertices.pXadj[i+iTotNodes+1]; + << VertexWgts[i + iTotNodes] + << " " << CommWgts[i + iTotNodes]; + for (int j = Vertices.pXadj[i + iTotNodes]; + j < Vertices.pXadj[i + iTotNodes + 1]; j++) { - ofMetis << " " << Vertices.pAdjncy[j]; + ofPartition << " " << Vertices.pAdjncy[j]; } - ofMetis << std::endl; + ofPartition << std::endl; } - ofMetis.close(); + ofPartition.close(); } #endif /* DEBUG */ switch (Partitioner) { + case PARTITIONER_MANUAL: + ASSERT(0); + break; + case PARTITIONER_CHACO: { - silent_cerr("CHACO partitioning algorithm not supported yet" << std::endl); - throw ErrGeneric(); + int *vwgt = 0; + int *cwgt = 0; + int *ewgt = 0; - /* Chaco uses floats for the communication weights */ - float* pfCommWgts = NULL; + if (wgtflag & WEIGHT_VERTICES) { + vwgt = &VertexWgts[0]; + } - SAFENEWARR(pfCommWgts, float, iTotVertices); + if (wgtflag & WEIGHT_COMM) { + /* unsupported */ + silent_cout("communication weights currently unsupported by CHACO; trying to emulate..." << std::endl); + cwgt = &CommWgts[0]; + } - for (int i = 0; i < iTotVertices; i++) { - pfCommWgts[i] = pCommWgts[i]; + if (wgtflag & WEIGHT_EDGES) { + /* unsupported ... */ + silent_cout("edges weights currently unsupported by MBDyn" << std::endl); } chaco_interface(iTotVertices, Vertices.pXadj, Vertices.pAdjncy, - pVertexWgts, - pfCommWgts, - /* wgtflag */ + vwgt, + cwgt, + ewgt, DataCommSize, pParAmgProcs); - - SAFEDELETEARR(pfCommWgts); break; } - case PARTITIONER_METIS: + case PARTITIONER_METIS: { + int *vwgt = 0; + int *cwgt = 0; + + if (wgtflag & WEIGHT_VERTICES) { + vwgt = &VertexWgts[0]; + } + + if (wgtflag & WEIGHT_EDGES) { + /* unsupported ... */ + silent_cout("edges weights currently unsupported by MBDyn" << std::endl); + } + + if (wgtflag & WEIGHT_COMM) { + cwgt = &CommWgts[0]; + } + mbdyn_METIS_PartGraph(iTotVertices, Vertices.pXadj, Vertices.pAdjncy, - pVertexWgts, - pCommWgts, - wgtflag, + vwgt, + cwgt, + NULL, /* ewgt */ DataCommSize, pParAmgProcs); break; + } default: - silent_cerr("Sorry. You need to compile with -DUSE_METIS " // "or -DUSE_CHACO." - << std::endl - << "No other partition library is implemented yet." - " Aborting ..." << std::endl); + silent_cerr("no partition library is available; " + "partition assignments must be provided " + "manually" << std::endl); throw ErrGeneric(); } } @@ -974,7 +1014,7 @@ if (pParAmgProcs[i] == MyRank) { /* se uno dei nodi e' connesso ad un elemento non appartenente * a questo processo e' un nodo di interfaccia */ - for (int j = Vertices.pXadj[i]; j < Vertices.pXadj[i+1]; j++) { + for (int j = Vertices.pXadj[i]; j < Vertices.pXadj[i + 1]; j++) { TmpPrc = pParAmgProcs[Vertices.pAdjncy[j]]; if (TmpPrc != MyRank) { InterfNodes.pXadj[TmpPrc] += 1; @@ -1050,7 +1090,7 @@ /* FIXME: there's a better way to find GravityPos and so on... */ ppMyElems[iNumLocElems] = ppElems[GravityPos]; iNumLocElems++; - pParAmgProcs[GravityPos+iTotNodes] = -1; + pParAmgProcs[GravityPos + iTotNodes] = -1; move++; } @@ -1058,7 +1098,7 @@ if (ElemData[Elem::AIRPROPERTIES].iNum != 0) { ppMyElems[iNumLocElems] = ppElems[AirPropPos]; iNumLocElems++; - pParAmgProcs[AirPropPos+iTotNodes] = -1; + pParAmgProcs[AirPropPos + iTotNodes] = -1; move++; } @@ -1108,14 +1148,14 @@ key = MyRank; } else { color = 1; - if (pParAmgProcs[pMyRot[i]+iTotNodes] == MyRank) { + if (pParAmgProcs[pMyRot[i] + iTotNodes] == MyRank) { silent_cout("Rotor " << ppElems[pRotPos[i]]->GetLabel() << " assigned to process " << MyRank << std::endl); iRotorIsMine = 1; key = 0; } else { - key = MyRank+1; + key = MyRank + 1; } } pRotorComm[i] = MBDynComm.Split(color, key); @@ -1128,13 +1168,13 @@ for (int i = 0; i < iMyTotRot; i++) { ppMyElems[iNumLocElems] = ppElems[pMyRot[i]]; - if (pParAmgProcs[pMyRot[i]+iTotNodes] == MyRank) { + if (pParAmgProcs[pMyRot[i] + iTotNodes] == MyRank) { iNumLocDofs += ppMyElems[iNumLocElems]->iGetNumDof(); } else { move++; } iNumLocElems++; - pParAmgProcs[pMyRot[i]+iTotNodes] = -1; + pParAmgProcs[pMyRot[i] + iTotNodes] = -1; } } @@ -1220,8 +1260,8 @@ ppMyElems[i]->GetConnectedNodes(iNumberOfNodes, pMyTypes, pMyLabels); for (int j = 0; j < iNumberOfNodes; j++) { unsigned int* p = - std::find(llabels, llabels+iNumIntNodes, pMyLabels[j]); - if (p != llabels+iNumIntNodes) { + std::find(llabels, llabels + iNumIntNodes, pMyLabels[j]); + if (p != llabels + iNumIntNodes) { ppMyIntElems[iNumIntElems] = ppMyElems[i]; iNumIntDofs += ppMyElems[i]->iGetNumDof(); iNumLocDofs -= ppMyElems[i]->iGetNumDof(); @@ -1243,8 +1283,8 @@ pMyTypes, pMyLabels); for (int j = 0; j < iNumberOfNodes; j++) { unsigned int* p = - std::find(llabels, llabels+iNumIntNodes, pMyLabels[j]); - if (p != llabels+iNumIntNodes) { + std::find(llabels, llabels + iNumIntNodes, pMyLabels[j]); + if (p != llabels + iNumIntNodes) { ppMyIntElems[iNumIntElems] = ppMyElems[i]; iNumIntDofs += ppMyElems[i]->iGetNumDof(); iNumLocDofs -= ppMyElems[i]->iGetNumDof(); @@ -1267,7 +1307,7 @@ if (ppMyNodes[i]->iGetNumDof() != 0) { integer First = (ppMyNodes[i])->iGetFirstIndex(); - pLocalDofs[iCount] = First+1; + pLocalDofs[iCount] = First + 1; iCount++; for (unsigned int j = 1; j < ppMyNodes[i]->iGetNumDof(); j++) { pLocalDofs[iCount] = First + j + 1; @@ -1287,7 +1327,7 @@ ElemWithDofs* pWithDofs = ppMyElems[i]->pGetElemWithDofs(); integer First = (pWithDofs)->iGetFirstIndex(); - pLocalDofs[iCount] = First+1; + pLocalDofs[iCount] = First + 1; iCount++; for (int j = 1; j < TmpDofNum; j++) { pLocalDofs[iCount] = First + 1 + j; @@ -1383,14 +1423,6 @@ SAFEDELETEARR(pLabelsList); } - if ( pVertexWgts != NULL) { - SAFEDELETEARR(pVertexWgts); - } - - if ( pCommWgts != NULL) { - SAFEDELETEARR(pCommWgts); - } - if ( pPosIntElems != NULL) { SAFEDELETEARR(pPosIntElems); } @@ -1417,22 +1449,47 @@ time_t tCurrTime(time(NULL)); OutHdl.Partition() << "# Partition file for MBDyn. Time: " << ctime(&tCurrTime) << std::endl - << "# Partition produced with METIS" << std::endl << std::endl << std::endl; + << "# Partition produced "; + switch (Partitioner) { + case PARTITIONER_MANUAL: + OutHdl.Partition() + << "manually"; + break; + + case PARTITIONER_METIS: + OutHdl.Partition() + << "with METIS"; + break; + + case PARTITIONER_CHACO: + OutHdl.Partition() + << "with CHACO"; + break; + + default: + ASSERT(0); + } + OutHdl.Partition() + << std::endl + << std::endl + << std::endl; /* Dati */ OutHdl.Partition() << "# Control data useful to verify the partition " - << std::endl << std::endl + << std::endl + << std::endl << "Total number of processes: " << DataCommSize << ";" << std::endl - << " Process #: " << MyRank << ";" << std::endl + << "Process #: " << MyRank << ";" << std::endl << "Total number of Nodes: " << iTotNodes << ";" << std::endl << "Total number of Elements: " << iTotElem << ";" << std::endl << "Total number of Dofs: " << iTotDofs << ";" << std::endl - << std::endl << std::endl; + << std::endl + << std::endl; - OutHdl.Partition() << " Local Dofs number: " << iNumLocDofs << std::endl; + OutHdl.Partition() << "Local Dofs number: " << iNumLocDofs << std::endl; - OutHdl.Partition() << " Local Nodes: " << iNumLocNodes << std::endl; + OutHdl.Partition() << "Local Nodes: " << iNumLocNodes << std::endl; for (int i = 0; i < iNumLocNodes; i++) { ASSERT(ppMyNodes[i] != NULL); OutHdl.Partition() @@ -1444,7 +1501,10 @@ } OutHdl.Partition() - << std::endl << std::endl << "Local Elements: "<< iNumLocElems << std::endl << std::endl; + << std::endl + << std::endl + << "Local Elements: "<< iNumLocElems << std::endl + << std::endl; for (int i = 0; i < iNumLocElems; i++) { ASSERT(ppMyElems[i] != NULL); OutHdl.Partition() @@ -1456,9 +1516,12 @@ } OutHdl.Partition() - << std::endl << std::endl << " Interface Dofs number: " << iNumIntDofs << std::endl; + << std::endl + << std::endl + << "Interface Dofs number: " << iNumIntDofs << std::endl; OutHdl.Partition() - << std::endl << "Interface Nodes: " << iNumIntNodes << std::endl; + << std::endl + << "Interface Nodes: " << iNumIntNodes << std::endl; for (int i = 0; i < iNumIntNodes; i++) { ASSERT(ppIntNodes[i] != NULL); OutHdl.Partition() @@ -1470,8 +1533,11 @@ } OutHdl.Partition() - << std::endl << std::endl << "Elements whose internal dofs are interface dofs: " - << iNumIntElems << std::endl << std::endl; + << std::endl + << std::endl + << "Elements whose internal dofs are interface dofs: " + << iNumIntElems << std::endl + << std::endl; for (int i = 0; i < iNumIntElems; i++) { OutHdl.Partition() << "Element Type: " @@ -1482,7 +1548,9 @@ } #ifdef DEBUG - OutHdl.Partition() << std::endl << "Local Dofs List:" << std::endl; + OutHdl.Partition() + << std::endl + << "Local Dofs List:" << std::endl; int j = 0; for (int i = 0; i < iNumLocDofs; i++) { @@ -1493,9 +1561,12 @@ j = 0; } } - OutHdl.Partition() << std::endl; + OutHdl.Partition() + << std::endl; - OutHdl.Partition() << std::endl << "Local Interface Dofs List:" <