/*
Copyright (C) 2000-2004
Code contributed by Greg Collecutt, Joseph Hope and Paul Cochrane
This file is part of xmds.
This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
*/
/*
$Id: xmdssimulation.cc,v 1.52 2005/10/27 00:58:39 joehope Exp $
*/
/*! @file xmdssimulation.cc
@brief Simulation element parsing classes and methods
More detailed explanation...
*/
#include<xmlbasics.h>
#include<dom3.h>
#include<xmdsutils.h>
#include<xmdsclasses.h>
#include<config.h>
// ******************************************************************************
// ******************************************************************************
// xmdsSimulation public
// ******************************************************************************
// ******************************************************************************
extern bool debugFlag;
long nxmdsSimulations=0; //!< The number of xmds simulation objects
// ******************************************************************************
xmdsSimulation::xmdsSimulation(
const char *const yourRawFileName,
const bool& yourVerboseMode,
const bool& mpiAvailable) :
xmdsElement(this,yourVerboseMode){
if(debugFlag) {
nxmdsSimulations++;
printf("xmdsSimulation::xmdsSimulation\n");
printf("nxmdsSimulations=%li\n",nxmdsSimulations);
}
myParameters.rawFileName = yourRawFileName;
myParameters.stochastic = 0;
myParameters.nThreads = 1;
myParameters.nPaths = 1;
myParameters.mpiAvailable = mpiAvailable;
myParameters.usempi = 0;
myParameters.bing = 0;
myParameters.seed[0] = 1;
myParameters.seed[1] = 2;
myParameters.nNoises = 0;
myParameters.noiseKind = "gaussian";
myParameters.noiseMean = 5.0;
myParameters.errorCheck = 1;
myParameters.useWisdom = 0;
myParameters.binaryOutput = 0;
myParameters.useDouble = 1;
myParameters.benchmark = 0;
myParameters.version = VERSION;
myParameters.release = RELEASE;
myCurrentSegmentNumber = 0;
};
// ******************************************************************************
xmdsSimulation::~xmdsSimulation() {
if(debugFlag) {
nxmdsSimulations--;
printf("xmdsSimulation::~xmdsSimulation\n");
printf("nxmdsSimulations=%li\n",nxmdsSimulations);
}
};
// ******************************************************************************
void xmdsSimulation::processElement(
const Element *const yourElement) {
if(debugFlag) {
printf("xmdsSimulation::processElement\n");
}
// ************************************
// process parameters
// ************************************
list<XMLString> myXMLStringList;
list<bool> myBoolList;
list<unsigned long> myULongList;
// ************************************
// find 'name'
getAssignmentStrings(yourElement,"name",0,1,myXMLStringList);
if(myXMLStringList.size()==1) {
myParameters.simulationName=*myXMLStringList.begin();
if(verbose()) {
printf("simulation name = '%s'\n",myParameters.simulationName.c_str());
}
}
else {
// Create default simulation file name
unsigned long lastdot=0;
const XMLString tempName = myParameters.rawFileName;
unsigned long i=0;
while(i<tempName.length()) {
if(tempName.data(i)=='.') {
lastdot=i;
}
i++;
}
if(lastdot>0) {
tempName.subString(myParameters.simulationName,0,lastdot);
}
else {
myParameters.simulationName=tempName;
}
if(verbose()) {
printf("simulation name defaulting to '%s'\n",myParameters.simulationName.c_str());
}
}
// ************************************
// find 'prop_dim'
getAssignmentStrings(yourElement,"prop_dim",1,1,myXMLStringList);
myParameters.propDimName=*myXMLStringList.begin();
if(verbose()) {
printf("simulation prop_dim = '%s'\n",myParameters.propDimName.c_str());
}
// ************************************
// find 'author'
getAssignmentStrings(yourElement,"author",NOT_REQD,0,myXMLStringList);
if (myXMLStringList.size() > 0) {
myParameters.authorName = *myXMLStringList.begin();
myXMLStringList.pop_front();
for(list<XMLString>::const_iterator pXMLString = myXMLStringList.begin(); pXMLString != myXMLStringList.end(); pXMLString++) {
myParameters.authorName += " ";
myParameters.authorName += *pXMLString;
}
if(verbose()) {
printf("simulation author = '%s'\n",myParameters.authorName.c_str());
}
}
else {
// this warning may need to be taken out somehow, but I sort of want people
// to be nice little coders and document their code nicely, and this is
// one way to do it... (PTC)
printf("No <author> tag found. It's not required, but it's a Good Idea.\n");
}
// ************************************
// find 'description'
getAssignmentStrings(yourElement,"description",NOT_REQD,0,myXMLStringList);
// Storing the description can cause overflow errors if it is too long, so we'll
// comment out the actual loading of the description, and put into the description
// variable a note that it actually exists.
if (myXMLStringList.size() > 0) {
myParameters.description += "Description found. See xmds file for the rest of it.";
/*
myParameters.description = *myXMLStringList.begin();
myXMLStringList.pop_front();
for(list<XMLString>::const_iterator pXMLString = myXMLStringList.begin(); pXMLString != myXMLStringList.end(); pXMLString++) {
myParameters.description += " ";
myParameters.description += *pXMLString;
}
*/
if(verbose()) {
printf("simulation description = '%s'\n",myParameters.description.c_str());
}
}
else {
// this warning may need to be taken out somehow, but I sort of want people
// to be nice little coders and document their code nicely, and this is
// one way to do it... (PTC)
printf("No <description> tag found. It's not required, but it's a Good Idea.\n");
}
// ************************************
// find 'error_check'
getAssignmentBools(yourElement,"error_check",0,1,myBoolList);
if(myBoolList.size()==1) {
myParameters.errorCheck=*myBoolList.begin();
if(verbose()) {
if(myParameters.errorCheck) {
printf("simulation error_check is 'yes'\n");
}
else {
printf("simulation error_check is 'no'\n");
}
}
}
else {
if(verbose()) {
printf("simulation error_check defaulting to 'yes'\n");
}
myParameters.errorCheck=1;
}
// ************************************
// find 'use_wisdom'
getAssignmentBools(yourElement,"use_wisdom",0,1,myBoolList);
if(myBoolList.size()==1) {
myParameters.useWisdom=*myBoolList.begin();
if(verbose()) {
if(myParameters.useWisdom) {
printf("simulation use_wisdom is 'yes'\n");
}
else {
printf("simulation use_wisdom is 'no'\n");
}
}
}
else {
if(verbose()) {
printf("simulation use_wisdom defaulting to 'no'\n");
}
myParameters.useWisdom=0;
}
// ************************************
// find 'use_prefs'
getAssignmentBools(yourElement,"use_prefs",0,1,myBoolList);
if(myBoolList.size()==1) {
myParameters.usePrefs=*myBoolList.begin();
if(verbose()) {
if(myParameters.usePrefs) {
printf("simulation use_prefs is 'yes'\n");
}
else {
printf("simulation use_prefs is 'no'\n");
}
}
}
else {
if(verbose()) {
printf("simulation use_prefs defaulting to 'yes'\n");
}
myParameters.usePrefs=1;
}
// ***********************************
// find binary output flag
getAssignmentBools(yourElement,"binary_output",0,1,myBoolList);
if(myBoolList.size()==1) {
// now we need to print a warning and tell the user what they can do to fix the problem
sprintf(errorMessage(),"Using the <binary_output> tag is obsolete.\n"
" Please replace this tag and use the format=\"ascii|binary\"\n"
" attribute of the <output> tag instead.");
throw xmdsException(errorMessage());
// this is the old code
/*
myParameters.binaryOutput=*myBoolList.begin();
if(verbose()) {
if(myParameters.binaryOutput) {
printf("simulation binary_output is 'yes'\n");
}
else
printf("simulation binary_output is 'no'\n");
}
*/
}
else {
if(verbose()) {
printf("simulation binary_output defaulting to 'no'\n");
}
myParameters.binaryOutput=0;
}
// ***********************************
// find use_double element
getAssignmentBools(yourElement,"use_double",0,1,myBoolList);
if(myBoolList.size()==1) {
// now we need to print a warning and tell the user what they can do to fix the problem
sprintf(errorMessage(),"Using the <use_double> tag is obsolete.\n"
" Please replace this tag and use the precision=\"double|single\"\n"
" attribute of the <output> tag instead.");
throw xmdsException(errorMessage());
// this is the old code
/*
myParameters.useDouble=*myBoolList.begin();
if(verbose()) {
if(myParameters.useDouble) {
printf("simulation use_double is 'yes'\n");
}
else
printf("simulation use_double is 'no'\n");
}
*/
}
else {
if(verbose()) {
printf("simulation use_double defaulting to 'yes'\n");
}
myParameters.useDouble=1;
}
// ***********************************
// find benchmark element
getAssignmentBools(yourElement,"benchmark",0,1,myBoolList);
if(myBoolList.size()==1) {
myParameters.benchmark=*myBoolList.begin();
if(verbose()) {
if(myParameters.benchmark) {
printf("simulation benchmark is 'yes'\n");
}
else {
printf("simulation benchmark is 'no'\n");
}
}
}
else {
if(verbose()) {
printf("simulation benchmark defaulting to 'no'\n");
}
myParameters.benchmark=0;
}
// ************************************
// find 'threads'
getAssignmentULongs(yourElement,"threads",0,1,myULongList);
if(myULongList.size()==1) {
myParameters.nThreads = *myULongList.begin();
if(myParameters.nThreads < 1) {
throw xmdsException(yourElement,"number of threads must be > 0");
}
if(verbose()) {
printf("number of threads = %li\n",myParameters.nThreads);
}
}
else {
if(verbose()) {
printf("number of threads defaulting to 1\n");
}
myParameters.nThreads=1;
}
// ************************************
// find 'stochastic'
getAssignmentBools(yourElement,"stochastic",0,1,myBoolList);
if(myBoolList.size()==1) {
myParameters.stochastic=*myBoolList.begin();
if(verbose()) {
if(myParameters.stochastic) {
printf("simulation stochastic is 'yes'\n");
}
else {
printf("simulation stochastic is 'no'\n");
}
}
}
else {
if(verbose()) {
printf("simulation stochastic defaulting to 'no'\n");
}
myParameters.stochastic=0;
}
// ************************************
// find 'use_mpi' assignment
getAssignmentBools(yourElement,"use_mpi",0,1,myBoolList);
if(myBoolList.size()==1) {
myParameters.usempi=*myBoolList.begin();
if(myParameters.usempi&!myParameters.mpiAvailable) {
throw xmdsException(yourElement,"MPI was not available when xmds was installed");
}
if(verbose()) {
if(myParameters.usempi) {
printf("use_mpi is 'yes'\n");
}
else {
printf("use_mpi is 'no'\n");
}
}
}
else {
if(verbose()) {
printf("use_mpi defaulting to 'no'\n");
}
myParameters.usempi=0;
}
if((myParameters.usempi)&(myParameters.nThreads>1)&(!myParameters.stochastic)) {
throw xmdsException(yourElement,"Non-determinstic MPI and threaded FFTs are not compatible");
}
// ************************************
// find 'bing' assignment
getAssignmentBools(yourElement,"bing",0,1,myBoolList);
if(myBoolList.size()==1) {
myParameters.bing=*myBoolList.begin();
if(verbose()) {
if(myParameters.bing) {
printf("bing is 'yes'\n");
}
else {
printf("bing is 'no'\n");
}
}
}
else {
if(verbose()) {
printf("bing defaulting to 'no'\n");
}
myParameters.bing=0;
}
// ************************************
// if stochastic get paths, and seed assignments
if(myParameters.stochastic) {
// ************************************
// get paths
getAssignmentULongs(yourElement,"paths",1,1,myULongList);
myParameters.nPaths = *myULongList.begin();
if(myParameters.nPaths < 1) {
throw xmdsException(yourElement,"number of paths must be > 0");
}
if(verbose()) {
printf("number of paths = %li\n",myParameters.nPaths);
}
// ************************************
// get seeds
getAssignmentULongs(yourElement,"seed",1,2,myULongList);
list<unsigned long>::const_iterator pULong = myULongList.begin();
myParameters.seed[0] = *pULong;
pULong++;
myParameters.seed[1] = *pULong;
if((myParameters.seed[0] < 0)|(myParameters.seed[1] < 0)) {
throw xmdsException(yourElement,"seeds must be positive");
}
if(myParameters.seed[0]==myParameters.seed[1]) {
throw xmdsException(yourElement,"seeds are indentical");
}
if(verbose()) {
printf("seeds are %li and %li\n",myParameters.seed[0],myParameters.seed[1]);
}
// ************************************
// get noises
getAssignmentULongs(yourElement,"noises",1,1,myULongList);
myParameters.nNoises = 2*((*myULongList.begin()+1)/2);
if(myParameters.nNoises < 2) {
throw xmdsException(yourElement,"number of noises must be > 0");
}
if(verbose()) {
printf("number of noises = %li\n",myParameters.nNoises);
}
// find out what kind of noise is specified, if any
getAttributeStrings(yourElement, "noises", "kind", NOT_REQD, myParameters.noiseKind);
// the kind has several values, current possible values are:
// (this list may be augmented in time)
// "gaussian"
// "gaussFast"
// "poissonian"
// "uniform"
if (debugFlag) {
printf("After getAttributeStrings(), noiseKind = %s\n",myParameters.noiseKind.c_str());
}
if (myParameters.noiseKind == "") {
if (verbose()) {
printf("Setting the default noise kind (i.e. \"gaussian\")\n");
}
myParameters.noiseKind = "gaussian"; // this is the default option
}
else if (myParameters.noiseKind == "gaussian") {
if (verbose()) {
printf("Noise kind set to \"gaussian\"\n");
}
}
else if (myParameters.noiseKind == "gaussFast") {
if (verbose()) {
printf("Noise kind set to \"gaussFast\"\n");
}
}
else if (myParameters.noiseKind == "poissonian") {
if (verbose()) {
printf("Noise kind set to \"poissonian\"\n");
}
// if we have poissonian specified, we must have a (nonzero) mean for the
// distribution specified as well
XMLString noiseMeanString;
getAttributeStrings(yourElement, "noises", "mean", REQD, noiseMeanString);
if (debugFlag) {
printf("After getAttributeStrings(), noiseMean = %s\n",noiseMeanString.c_str());
}
// now convert the string to its relevant floating point value
myParameters.noiseMean = atof(noiseMeanString.c_str());
if (verbose()) {
printf("Mean for Poissonian noise set to %g\n",myParameters.noiseMean);
}
}
else if (myParameters.noiseKind == "uniform") {
if (verbose()) {
printf("Noise kind set to \"uniform\"\n");
}
}
else {
throw xmdsException(yourElement,
"You must specify a correct noise kind. Current choices are:\n gaussian\n gaussFast\n poissonian\n uniform");
}
// find out which scheduling method to use
if(myParameters.usempi) {
getAssignmentStrings(yourElement,"MPI_Method",0,1,myXMLStringList);
if(myXMLStringList.size()==0) {
// Default mpi method is "scheduling"
myXMLStringList.push_back("Scheduling");
myParameters.mpiMethod=*myXMLStringList.begin();
if(verbose()) {
printf("MPI method defaulting to '%s'\n",myParameters.mpiMethod.c_str());
}
}
else {
if(!((*myXMLStringList.begin()=="Scheduling")||(*myXMLStringList.begin()=="Uniform"))) {
sprintf(errorMessage(),"MPI method '%s' unknown.\nCurrent methods are 'Scheduling' and 'Uniform'.",myXMLStringList.begin()->c_str());
throw xmdsException(yourElement,errorMessage());
}
myParameters.mpiMethod=*myXMLStringList.begin();
if(verbose()) {
printf("MPI method = '%s'\n",myParameters.mpiMethod.c_str());
}
}
}
}
// ************************************
// find and process children
// ************************************
const NodeList* candidateElements;
// ************************************
// find argv element and process
candidateElements = yourElement->getElementsByTagName("argv",0);
if(candidateElements->length()>1) {
throw xmdsException(yourElement,"Multiple <argv> elements!");
}
else if(candidateElements->length()==1) {
myArgv = createxmdsArgv();
const Element* yourElement = dynamic_cast<const Element*> (candidateElements->item(0));
myArgv->processElement(yourElement);
}
// ************************************
// find globals element and process if present
candidateElements = yourElement->getElementsByTagName("globals",0);
if(candidateElements->length()>1) {
throw xmdsException(yourElement,"Multiple <globals> elements!");
}
else if(candidateElements->length()==1) {
xmdsElement* newGlobals = createxmdsGlobals();
const Element* yourElement = dynamic_cast<const Element*> (candidateElements->item(0));
newGlobals->processElement(yourElement);
}
// ************************************
// find field element and process
candidateElements = yourElement->getElementsByTagName("field",0);
if(candidateElements->length()>1) {
throw xmdsException(yourElement,"Multiple <field> elements!");
}
else if(candidateElements->length()==1) {
myField = createxmdsField();
const Element* yourElement = dynamic_cast<const Element*> (candidateElements->item(0));
myField->processElement(yourElement);
}
else
throw xmdsException(yourElement,
"Cannot find <field> element");
// ************************************
// find output element and process
candidateElements = yourElement->getElementsByTagName("output",0);
if(candidateElements->length()>1) {
throw xmdsException(yourElement,"Multiple <output> elements!");
}
else if(candidateElements->length()==1) {
// now grab the format attribute if available
XMLString outputFormat = "ascii"; // set it to the default while declaring it
const Node *thisElement = candidateElements->item(0);
getAttributeStrings(yourElement, "output", "format", NOT_REQD, outputFormat);
if (outputFormat == "") {
outputFormat = "ascii"; // the default option
}
else {
if (outputFormat == "ascii") {
myParameters.binaryOutput = false;
if (verbose()) {
printf("Output format set to 'ascii' in output element\n");
}
}
else if (outputFormat == "binary") {
myParameters.binaryOutput = true;
if (verbose()) {
printf("Output format set to 'binary' in output element\n");
}
}
else {
throw xmdsException(thisElement, "Output format should be either 'ascii' or 'binary'");
}
}
// now grab the precision attribute if available
XMLString outputPrecision = "double"; // set it to the default while declaring it
getAttributeStrings(yourElement, "output", "precision", NOT_REQD, outputPrecision);
if (outputPrecision == "") {
outputPrecision = "double"; // the default option
}
else {
if (outputPrecision == "double") {
myParameters.useDouble = true;
if (verbose()) {
printf("Output precision set to 'double' in output element\n");
}
}
else if (outputPrecision == "single") {
myParameters.useDouble = false;
if (verbose()) {
printf("Output precision set to 'single' in output element\n");
}
}
else {
throw xmdsException(thisElement, "Output precision should be either 'double' or 'single'");
}
}
// now process the output element
myOutput = createxmdsOutput();
const Element* yourElement = dynamic_cast<const Element*> (candidateElements->item(0));
myOutput->processElement(yourElement);
}
else
throw xmdsException(yourElement,"Cannot find <output> element");
// ************************************
// find sequence element and process
candidateElements = yourElement->getElementsByTagName("sequence",0);
if(candidateElements->length()>1) {
throw xmdsException(yourElement,"Multiple <sequence> elements at top level");
}
else if(candidateElements->length()==1) {
mySequence = createxmdsSequence();
const Element* yourElement = dynamic_cast<const Element*> (candidateElements->item(0));
mySequence->processElement(yourElement);
}
else
throw xmdsException(yourElement,"Cannot find <sequence> element");
myField->outputSampleCount();
mySequence->outputSampleCount();
myOutput->finaliseGeometry();
// ************************************
// Allow non-stochastic mpi?
if(myParameters.usempi) {
if((myField->geometry()->nDims()<2)&(!myParameters.stochastic)) {
myParameters.usempi=0;
printf("MPI disabled. MPI cannot be used for deterministic problems with transverse dimensions < 2\n");
}
if((myField->geometry()->nDims()>1)&(!myParameters.stochastic)){
printf("MPI used for a deterministic problem. Ain't life great? \n");
}
}
};
// ******************************************************************************
const xmdsSimulation::simulationParametersStruct*
xmdsSimulation::parameters() const {
return &myParameters;
}
// ******************************************************************************
xmdsField* xmdsSimulation::field() const {
return myField;
};
// ******************************************************************************
xmdsSimulation::argvStruct*
xmdsSimulation::argStruct() const {
return &myArgStruct;
};
// ******************************************************************************
const xmdsOutput* xmdsSimulation::output() const {
return myOutput;
};
// ******************************************************************************
const xmdsSequence* xmdsSimulation::sequence() const {
return mySequence;
};
// ******************************************************************************
unsigned long xmdsSimulation::nextSegmentNumber() const {
if(debugFlag) {
printf("xmdsSimulation::nextSegmentNumber\n");
}
myCurrentSegmentNumber++;
return myCurrentSegmentNumber-1;
};
// ******************************************************************************
void xmdsSimulation::makeCode(
const unsigned long& inFileSplitPoint) const {
if(debugFlag) {
printf("xmdsSimulation::makeCode\n");
}
if(myOutput != 0) {
myOutput->setInFileSplitPoint(inFileSplitPoint);
}
XMLString fileName = myParameters.simulationName;
fileName += ".cc";
FILE* outfile = fopen(fileName.c_str(),"w");
if(outfile==0) {
sprintf(errorMessage(),"Unable to open file '%s' for write access",fileName.c_str());
throw xmdsException(errorMessage());
}
if(verbose()) {
printf("Beginning to write code ...\n");
}
try {
writeIncludes(outfile);
writeDefines(outfile);
writeGlobals(outfile);
writePrototypes(outfile);
writeRoutines(outfile);
}
catch (xmdsException xmdsExceptionErr) {
printf("Error: simulation failed to write output code\n");
printf("due to the following xmdsException:\n");
printf("%s",xmdsExceptionErr.getError());
fclose(outfile);
throw xmdsException("Internal Error");
}
fclose(outfile);
};
// ******************************************************************************
// ******************************************************************************
// xmdsSimulation private
// ******************************************************************************
// ******************************************************************************
// ******************************************************************************
void xmdsSimulation::writeIncludes(
FILE *const outfile) const {
if(debugFlag) {
printf("xmdsSimulation::writeIncludes\n");
}
if(verbose()) {
printf("Writing simulation includes ...\n");
}
fprintf(outfile,
"// ********************************************************\n"
"// simulation includes\n"
"\n"
"#include <stdio.h>\n"
"#include <math.h>\n"
"#include <string>\n"
"#include <fstream>\n"
"#include <iostream>\n"
"#include <sstream>\n"
"#include <cstdlib>\n"
"#include <unistd.h>\n"
"#include <sys/time.h>\n"
"#include <time.h>\n");
if(simulation()->argStruct()->nameList.size() != 0) {
fprintf(outfile,"#include <getopt_xmds.h>\n");
}
fprintf(outfile,"#include <vector>\n");
if(myParameters.usempi) {
fprintf(outfile,"#include <mpi.h>\n");
}
if(myParameters.nThreads>1) {
fprintf(outfile,"#include <fftw_threads.h>\n");
}
else if((!myParameters.stochastic)&(myParameters.usempi)) {
fprintf(outfile,
"#include <fftw_mpi.h>\n"
"#include <fftw.h>\n");
}
else
fprintf(outfile,"#include <fftw.h>\n");
fprintf(outfile,"#include <xmdscomplex.h>\n");
if (myParameters.binaryOutput) {
fprintf(outfile,"#include <xmdsconfig.h>\n");
}
fprintf(outfile,
"using namespace std;\n"
"\n");
};
// ******************************************************************************
void xmdsSimulation::writeDefines(
FILE *const outfile) const {
if(debugFlag) {
printf("xmdsSimulation::writeDefines\n");
}
if(verbose()) {
printf("Writing simulation defines ...\n");
}
if(myParameters.mpiMethod=="Scheduling") {
fprintf(outfile,
"// ********************************************************\n"
"//Adaptive Scheduler defines\n"
"\n"
"#define MIN(a,b) (((a)<(b))?(a):(b))\n"
"#define MAX(a,b) (((a)<(b))?(b):(a))\n"
"#define DEBUG_MPI_Scheduling 1\n\n");
}
fprintf(outfile,
"// ********************************************************\n"
"// simulation defines\n"
"\n");
if(myParameters.nThreads>1) {
fprintf(outfile,"#define _num_threads %li\n\n",myParameters.nThreads);
}
if(myParameters.stochastic) {
fprintf(outfile,"#define _n_noises %li\n",myParameters.nNoises);
for(long i=0;i<myParameters.nNoises;i++) {
fprintf(outfile,"#define n_%li _noises[%li]\n",i+1,i);
}
if(myParameters.nPaths>1) {
fprintf(outfile,"#define _n_paths %li\n",myParameters.nPaths);
}
}
fprintf(outfile,"#define d%s _step\n",myParameters.propDimName.c_str());
fprintf(outfile,"\n");
xmdsElement::writeDefines(outfile);
fprintf(outfile,"\n");
};
// ******************************************************************************
void xmdsSimulation::writeGlobals(
FILE *const outfile) const {
if(debugFlag) {
printf("xmdsSimulation::writeGlobals\n");
}
if(verbose()) {
printf("Writing simulation globals ...\n");
}
if(myParameters.mpiMethod=="Scheduling") {
fprintf(outfile,
"// ********************************************************\n"
"//Adaptive Scheduler globals\n"
"\n"
"pthread_mutex_t tasklock; /*Ensures mutual exclusion when assigning tasks*/\n"
"pthread_mutex_t finlock; /*Lock to synchronize completion of thread and master*/\n"
"\n"
"int total_paths; /*total number of paths of each type to be completed*/\n"
"int *slave_stat; /*Array of slaves indicating their present status IDLE(0) or BUSY(1)*/\n"
"\n"
"int *fschedcount; /*Array of num full paths scheduled for each slave*/\n"
"double *fcommave; /*Array of full path communication latency for each slave*/\n"
"double *fthroughput; /*Array of full path throughput (task/sec) for each slave*/\n"
"int fullpaths_assigned=0; /*number of full paths assigned or completed*/\n"
"static int fpathnum = 0; /*for debugging only*/\n"
"\n");
if(myParameters.errorCheck) {
fprintf(outfile,
"int *hschedcount; /*Array of num half paths scheduled for each slave*/\n"
"double *hcommave; /*Array of half path communication latency for each slave*/\n"
"double *hthroughput; /*Array of half path throughput (task/sec) for each slave*/\n"
"int halfpaths_assigned=0; /*number of half paths assigned or completed*/\n"
"static int hpathnum = 0; /*for debugging only*/\n\n");
}
}
fprintf(outfile,
"// ********************************************************\n"
"// simulation globals\n"
"\n");
fprintf(outfile,"double %s;\n",myParameters.propDimName.c_str());
fprintf(outfile,"\n");
if(myParameters.errorCheck) {
fprintf(outfile,"bool _half_step;\n\n");
}
if(myParameters.usempi){
fprintf(outfile,"int rank, size;\n\n");
}
if(myParameters.stochastic) {
if (myParameters.noiseKind == "poissonian") {
fprintf(outfile,
"double _noiseMean = %g;\n",myParameters.noiseMean);
}
if(myParameters.errorCheck) {
fprintf(outfile,
"unsigned short _gen1[3];\n"
"unsigned short _gen2[3];\n");
}
else
fprintf(outfile,"unsigned short _gen[3];\n");
fprintf(outfile,"\n");
}
else if(myParameters.usempi){
fprintf(outfile,
"int local_nx, local_x_start, local_ny_after_transpose, local_y_start_after_transpose, total_local_size;\n\n");
}
xmdsElement::writeGlobals(outfile);
fprintf(outfile,"\n");
};
// ******************************************************************************
void xmdsSimulation::writePrototypes(
FILE *const outfile) const {
if(debugFlag) {
printf("xmdsSimulation::writePrototypes\n");
}
if(verbose()) {
printf("Writing simulation prototypes ...\n");
}
if(myParameters.stochastic) {
if(myParameters.mpiMethod=="Scheduling") {
fprintf(outfile,
"// ********************************************************\n"
"//Adaptive Scheduler prototypes\n"
"\n"
"int print_array2(int *array, int size); /*prints int arrays*/\n"
"int print_array(double *array, int size); /*prints double arrays*/\n"
"void *slave(void *ptr); /*slave routine*/\n"
"int master(); /*master routine*/\n"
"void *mslave(void *ptr); /*thread slave*/\n\n"
);
}
if (debugFlag) {
printf("xmdsSimulation::writePrototypes noiseKind = %s\n",myParameters.noiseKind.c_str());
}
if (myParameters.noiseKind == "gaussian") {
fprintf(outfile,
"// ********************************************************\n"
"// simulation prototypes\n"
"\n"
"void _make_noises(\n"
" unsigned short *const _generator,\n"
" const double& _var,\n"
" double *const _noise_vector,\n"
" const unsigned long& _n);\n"
"\n");
}
else if (myParameters.noiseKind == "gaussFast") {
fprintf(outfile,
"// ********************************************************\n"
"// simulation prototypes\n"
"\n"
"void _make_noises(\n"
" unsigned short *const _generator,\n"
" const double& _var,\n"
" double *const _noise_vector,\n"
" const unsigned long& _n);\n"
"\n");
}
else if (myParameters.noiseKind == "uniform") {
fprintf(outfile,
"// ********************************************************\n"
"// simulation prototypes\n"
"\n"
"void _make_noises(\n"
" unsigned short *const _generator,\n"
" const double& _var,\n"
" double *const _noise_vector,\n"
" const unsigned long& _n);\n"
"\n");
}
else if (myParameters.noiseKind == "poissonian") {
fprintf(outfile,
"// ********************************************************\n"
"// simulation prototypes\n"
"\n"
"void _make_noises(\n"
" unsigned short *const _generator,\n"
" const double& _var,\n"
" double *const _noise_vector,\n"
" const unsigned long& _n);\n"
"\n"
"double gammln(double _xx);\n"
"\n"
"double poidev(double xm, unsigned short *const _generator);\n"
"\n");
}
else {
sprintf(errorMessage(),"Incorrectly specified noise kind (%s) in xmdsSimulation::writePrototypes",myParameters.noiseKind.c_str());
throw xmdsException(errorMessage());
}
}
xmdsElement::writePrototypes(outfile);
};
// ******************************************************************************
void xmdsSimulation::writeRoutines(
FILE *const outfile) const {
if(debugFlag) {
printf("xmdsSimulation::writeRoutines\n");
}
if(verbose()) {
printf("Writing simulation routines ...\n");
}
// NOTE: it turns out the erand48() is more likely to be thread safe, and generate
// decent random numbers. rand() isn't necessarily thread safe, and the thread safe
// version (rand_r()) isn't as good a random number generator.
if(myParameters.stochastic) {
if (debugFlag) {
printf("xmdsSimulation::writeRoutines noiseKind = %s\n",myParameters.noiseKind.c_str());
}
if (myParameters.noiseKind == "gaussian") {
fprintf(outfile,
"// ********************************************************\n"
"// simulation routines\n"
"// ********************************************************\n"
"\n"
"// *************************\n"
"void _make_noises(\n"
" unsigned short *const _generator,\n"
" const double& _var,\n"
" double *const _noise_vector,\n"
" const unsigned long& _n) {\n"
"\n"
"for(unsigned long _i0=0; _i0<_n; _i0+=2) {\n"
" const double _mag = sqrt(-2*log(erand48(_generator))*_var);\n"
" const double _arg = 2*M_PI*erand48(_generator);\n"
" _noise_vector[_i0 + 0] = _mag*cos(_arg);\n"
" _noise_vector[_i0 + 1] = _mag*sin(_arg);\n"
" }\n"
"}\n"
"\n");
}
else if (myParameters.noiseKind == "gaussFast") {
fprintf(outfile,
"// ********************************************************\n"
"// simulation routines\n"
"// ********************************************************\n"
"\n"
"// *************************\n"
"void _make_noises(\n"
" unsigned short *const _generator,\n"
" const double& _var,\n"
" double *const _noise_vector,\n"
" const unsigned long& _n) {\n"
"\n"
"for(unsigned long _i0=0; _i0<_n; _i0+=2) {\n"
" double _v1, _v2, _rsq;\n"
" do {\n"
" _v1 = 2.0*erand48(_generator)-1.0;\n"
" _v2 = 2.0*erand48(_generator)-1.0;\n"
" _rsq = _v1*_v1 + _v2*_v2;\n"
" } while(_rsq >= 1.0 || _rsq == 0.0);\n"
" const double _fac = sqrt(-2.0*_var*log(_rsq)/_rsq);\n"
" _noise_vector[_i0 + 0] = _v1*_fac;\n"
" _noise_vector[_i0 + 1] = _v2*_fac;\n"
" }\n"
"}\n"
"\n");
}
else if (myParameters.noiseKind == "uniform") {
fprintf(outfile,
"// ********************************************************\n"
"// simulation routines\n"
"// ********************************************************\n"
"\n"
"// *************************\n"
"void _make_noises(\n"
" unsigned short *const _generator,\n"
" const double& _var,\n"
" double *const _noise_vector,\n"
" const unsigned long& _n) {\n"
"\n"
"for(unsigned long _i0=0; _i0<_n; _i0+=2) {\n"
" const double _fac = sqrt(_var);\n"
" _noise_vector[_i0 + 0] = _fac*erand48(_generator);\n"
" _noise_vector[_i0 + 1] = _fac*erand48(_generator);\n"
" }\n"
"}\n"
"\n");
}
else if (myParameters.noiseKind == "poissonian") {
fprintf(outfile,
"// ********************************************************\n"
"// simulation routines\n"
"// ********************************************************\n"
"\n"
"// *************************\n"
"double gammln(double _xx) {\n"
" double _x, _y, _tmp, _ser;\n"
" const double _cof[6] = {76.18009172947146,-86.50532032941677,\n"
" 24.01409824083091,-1.231739572450155,\n"
" 0.1208650973866179e-2,-0.5395239384953e-5};\n"
" int _j;\n"
" _y = _x = _xx;\n"
" _tmp = _x + 5.5;\n"
" _tmp -= (_x+0.5)*log(_tmp);\n"
" _ser = 1.000000000190015;\n"
" for (_j=0; _j<=5; _j++) _ser += _cof[_j]/++_y;\n"
" return -_tmp+log(2.5066282746310005*_ser/_x);\n"
"}\n"
"\n"
"double poidev(double xm, unsigned short *const _generator) {\n"
" static double sq, alxm, g;\n"
" double em,t,y;\n"
" if (xm < 12.0) { // Use direct method \n"
" g=exp(-xm);\n"
" em = -1;\n"
" t=1.0;\n"
" // Instead of adding exponential deviates it is equivalent\n"
" // to multiply uniform deviates. We never actually have to\n"
" // take the log, merely compare to the pre-computed exponential\n"
" do {\n"
" ++em;\n"
" t *= erand48(_generator);\n"
" } while (t > g);\n"
" } else { // Use rejection method\n"
" sq=sqrt(2.0*xm);\n"
" alxm=log(xm);\n"
" g=xm*alxm-gammln(xm+1.0);\n"
" do {\n"
" do { // y is a deviate from a Lorenzian comparison function\n"
" y=tan(M_PI*erand48(_generator));\n"
" em=sq*y+xm; // em is y, shifted and scaled\n"
" } while (em < 0.0); // Reject if in regime of zero probability\n"
" em=floor(em); // The trick for integer-valued distributions\n"
" t=0.9*(1.0+y*y)*exp(em*alxm-gammln(em+1.0)-g);\n"
" // The ratio of the desired distribution to the comparison\n"
" // function; we reject by comparing it to another uniform\n"
" // deviate. The factor 0.9 so that t never exceeds 1.\n"
" } while (erand48(_generator) > t);\n"
" }\n"
" return em;\n"
" }\n"
"\n"
"void _make_noises(\n"
" unsigned short *const _generator,\n"
" const double& _var,\n"
" double *const _noise_vector,\n"
" const unsigned long& _n) {\n"
"\n"
"for(unsigned long _i0=0; _i0<_n; _i0+=2) {\n"
" const double _fac = sqrt(_var);\n"
" _noise_vector[_i0 + 0] = _fac*poidev(_noiseMean, _generator);\n"
" _noise_vector[_i0 + 1] = _fac*poidev(_noiseMean, _generator);\n"
" }\n"
"}\n"
"\n");
}
else {
sprintf(errorMessage(),"Incorrectly specified noise kind (%s) in xmdsSimulation::writeRoutines",myParameters.noiseKind.c_str());
throw xmdsException(errorMessage());
}
}
xmdsElement::writeRoutines(outfile);
};
// ******************************************************************************
xmdsGlobals* xmdsSimulation::createxmdsGlobals() {
if(debugFlag) {
printf("xmdsSimulation::createxmdsGlobals\n");
}
xmdsGlobals* newxmdsGlobals =
new xmdsGlobals(this,verbose());
addChild((xmdsElement*) newxmdsGlobals);
return newxmdsGlobals;
};
// ******************************************************************************
xmdsField* xmdsSimulation::createxmdsField() {
if(debugFlag) {
printf("xmdsSimulation::createxmdsField\n");
}
xmdsField* newxmdsField =
new xmdsField(this,verbose());
addChild((xmdsElement*) newxmdsField);
return newxmdsField;
};
// ******************************************************************************
xmdsArgv* xmdsSimulation::createxmdsArgv() {
if(debugFlag) {
printf("xmdsSimulation::createxmdsArgv\n");
}
xmdsArgv* newxmdsArgv = new xmdsArgv(this,verbose());
addChild((xmdsElement*) newxmdsArgv);
return newxmdsArgv;
};
// ******************************************************************************
xmdsOutput* xmdsSimulation::createxmdsOutput() {
if(debugFlag) {
printf("xmdsSimulation::createxmdsOutput\n");
}
xmdsOutput* newxmdsOutput =
new xmdsOutput(this,verbose());
addChild((xmdsElement*) newxmdsOutput);
return newxmdsOutput;
};
// ******************************************************************************
xmdsSequence* xmdsSimulation::createxmdsSequence() {
if(debugFlag) {
printf("xmdsSimulation::createxmdsSequence\n");
}
xmdsSequence* newxmdsSequence =
new xmdsSequence(this,verbose());
addChild((xmdsElement*) newxmdsSequence);
return newxmdsSequence;
};
syntax highlighted by Code2HTML, v. 0.9.1