26 #include <queso/GslVector.h>
27 #include <queso/Defines.h>
28 #include <gsl/gsl_sort_vector.h>
36 m_vec (gsl_vector_calloc(map.NumGlobalElements()))
59 m_vec (gsl_vector_calloc(map.NumGlobalElements()))
79 m_vec (gsl_vector_calloc(map.NumGlobalElements()))
89 for (
unsigned int i = 0; i <
m_vec->size; ++i) {
90 double alpha = (double) i / ((
double)
m_vec->size - 1.);
91 (*this)[i] = (1.-alpha)*d1 + alpha*d2;
102 m_vec (gsl_vector_calloc(v.sizeLocal()))
112 for (
unsigned int i = 0; i <
m_vec->size; ++i) {
113 double alpha = (double) i / ((
double)
m_vec->size - 1.);
114 (*this)[i] = (1. - alpha) * start + alpha * end;
125 m_vec (gsl_vector_calloc(v.sizeLocal()))
168 iRC = gsl_vector_scale(
m_vec,a);
188 for (
unsigned int i = 0; i < size1; ++i) {
189 (*this)[i] *= rhs[i];
202 for (
unsigned int i = 0; i < size1; ++i) {
203 (*this)[i] /= rhs[i];
233 iRC = gsl_vector_memcpy(this->
m_vec, src.
m_vec);
282 return std::sqrt(this->
norm2Sq());
291 for (
unsigned int i = 0; i < size; ++i) {
292 result += fabs((*
this)[i]);
305 for (
unsigned int i = 0; i < size; ++i) {
306 aux = fabs((*
this)[i]);
307 if (aux > result) result = aux;
318 for (
unsigned int i = 0; i < size; ++i) {
319 result += (*this)[i];
329 for (
unsigned int i = 0; i < size; ++i) {
339 for (
unsigned int i = 0; i < this->
sizeLocal(); ++i) {
349 for (
unsigned int i = 0; i < this->
sizeLocal(); ++i) {
358 for (
unsigned int i = 0; i < this->
sizeLocal(); ++i) {
371 double tmpSample = 0.;
372 for (
unsigned int i = 0; i < this->
sizeLocal(); ++i) {
378 <<
", alpha[i] = " << alpha[i]
379 <<
", beta[i] = " << beta[i]
380 <<
", sample = " << tmpSample
383 if ((alpha[i] == 1. ) &&
385 if (tmpSample == 1.) {
390 <<
", alpha[i] = " << alpha[i]
391 <<
", beta[i] = " << beta[i]
392 <<
", sample = " << tmpSample
396 std::cerr <<
"Hitting 'sample = 1' in GslVector::cwSetBeta()"
399 <<
", alpha[i] = " << alpha[i]
400 <<
", beta[i] = " << beta[i]
401 <<
", sample = " << tmpSample
405 }
while (tmpSample == 1.);
406 std::cerr <<
"Code was able to get 'sample != 1' in GslVector::cwSetBeta()"
409 <<
", alpha[i] = " << alpha[i]
410 <<
", beta[i] = " << beta[i]
411 <<
", sample = " << tmpSample
416 (*this)[i] = tmpSample;
428 for (
unsigned int i = 0; i < this->
sizeLocal(); ++i) {
441 for (
unsigned int i = 0; i < this->
sizeLocal(); ++i) {
452 for (
unsigned int i = 0; i < v1.
sizeLocal(); ++i) {
456 for (
unsigned int i = 0; i < v2.
sizeLocal(); ++i) {
466 unsigned int cummulativeSize = 0;
467 for (
unsigned int i = 0; i < vecs.size(); ++i) {
469 for (
unsigned int j = 0; j < vecs[i]->sizeLocal(); ++j) {
470 (*this)[cummulativeSize+j] = tmpVec[j];
487 for (
unsigned int i = 0; i < vec.
sizeLocal(); ++i) {
488 (*this)[initialPos+i] = vec[i];
501 for (
unsigned int i = 0; i < vec.
sizeLocal(); ++i) {
502 vec[i] = (*this)[initialPos+i];
512 for (
unsigned int i = 0; i < size; ++i) {
513 (*this)[i] = 1./(*this)[i];
523 for (
unsigned int i = 0; i < size; ++i) {
524 (*this)[i] = sqrt((*
this)[i]);
532 unsigned int firstPositionToStoreDiff,
533 double valueForRemainderPosition,
542 for (
unsigned int i = 0; i < (size-1); ++i) {
543 outputVec[firstPositionToStoreDiff+i] = (*this)[i+1]-(*this)[i];
545 if (firstPositionToStoreDiff == 0) {
546 outputVec[size-1] = valueForRemainderPosition;
549 outputVec[0] = valueForRemainderPosition;
567 for (
unsigned int i = 1; i < x1Vec.
sizeLocal(); ++i) {
571 for (
unsigned int id2 = 0; id2 < x2Vec.
sizeLocal(); ++id2) {
572 double x2 = x2Vec[id2];
573 unsigned int id1 = 0;
575 for (id1 = 0; id1 < x1Vec.
sizeLocal(); ++id1) {
576 if (x2 <= x1Vec[id1]) {
581 bool makeLinearModel =
false;
586 if (x2 == x1Vec[id1]) {
587 (*this)[id2] = y1Vec[id1];
589 else if (x2 < x1Vec[0]) {
591 makeLinearModel =
true;
597 else if (found1 ==
true) {
599 makeLinearModel =
true;
607 makeLinearModel =
true;
614 if (makeLinearModel) {
615 double rate = (yb-ya)/(xb-xa);
616 (*this)[id2] = ya + (x2-xa)*rate;
628 gsl_sort_vector(
m_vec);
637 if (bcastComm.
MyPID() < 0)
return;
643 double localNumNodes = 1.;
644 double totalNumNodes = 0.;
646 "GslVector::mpiBcast()",
647 "failed MPI.Allreduce() for numNodes");
651 double localVectorSize = this->
sizeLocal();
652 double sumOfVectorSizes = 0.;
654 "GslVector::mpiBcast()",
655 "failed MPI.Allreduce() for vectorSize");
657 if ( ((
unsigned int) sumOfVectorSizes) != ((
unsigned int)(totalNumNodes*localVectorSize)) ) {
658 std::cerr <<
"rank " << bcastComm.
MyPID()
659 <<
": sumOfVectorSizes = " << sumOfVectorSizes
660 <<
", totalNumNodes = " << totalNumNodes
661 <<
", localVectorSize = " << localVectorSize
665 queso_require_equal_to_msg(((
unsigned int) sumOfVectorSizes), ((
unsigned int)(totalNumNodes*localVectorSize)),
"inconsistent vectorSize");
668 std::vector<double> dataBuffer((
unsigned int) localVectorSize, 0.);
669 if (bcastComm.
MyPID() == srcRank) {
670 for (
unsigned int i = 0; i < dataBuffer.size(); ++i) {
671 dataBuffer[i] = (*this)[i];
676 "GslVector::mpiBcast()",
677 "failed MPI.Bcast()");
679 if (bcastComm.
MyPID() != srcRank) {
680 for (
unsigned int i = 0; i < dataBuffer.size(); ++i) {
681 (*this)[i] = dataBuffer[i];
692 if (opComm.
MyPID() < 0)
return;
697 for (
unsigned int i = 0; i < size; ++i) {
698 double srcValue = (*this)[i];
699 double resultValue = 0.;
700 opComm.
Allreduce<
double>(&srcValue, &resultValue, (int) 1, mpiOperation,
701 "GslVector::mpiAllReduce()",
702 "failed MPI.Allreduce()");
703 resultVec[i] = resultValue;
713 if (opComm.
MyPID() < 0)
return;
720 for (
unsigned int i = 0; i < size; ++i) {
721 double auxDouble = (int) (*
this)[i];
722 std::vector<double> vecOfDoubles(opComm.
NumProc(),0.);
723 opComm.
Gather<
double>(&auxDouble, 1, &vecOfDoubles[0], (int) 1, 0,
724 "GslVector::mpiAllQuantile()",
725 "failed MPI.Gather()");
727 std::sort(vecOfDoubles.begin(), vecOfDoubles.end());
729 double result = vecOfDoubles[(
unsigned int)( probability*((
double)(vecOfDoubles.size()-1)) )];
732 "GslVector::mpiAllQuantile()",
733 "failed MPI.Bcast()");
735 resultVec[i] = result;
752 std::ostream::fmtflags curr_fmt = os.flags();
757 unsigned int savedPrecision = os.precision();
761 for (
unsigned int i = 0; i < size; ++i) {
762 os << std::scientific << (*this)[i]
767 for (
unsigned int i = 0; i < size; ++i) {
768 os << std::scientific << (*this)[i]
773 os.precision(savedPrecision);
779 for (
unsigned int i = 0; i < size; ++i) {
780 os << std::dec << (*this)[i]
785 for (
unsigned int i = 0; i < size; ++i) {
786 os << std::dec << (*this)[i]
803 const std::string& varNamePrefix,
804 const std::string& fileName,
805 const std::string& fileType,
806 const std::set<unsigned int>& allowedSubEnvIds)
const
828 *filePtrSet.
ofsVar << *
this;
833 *filePtrSet.
ofsVar <<
"];\n";
843 const std::string& fileName,
844 const std::string& fileType,
845 const std::set<unsigned int>& allowedSubEnvIds)
859 unsigned int idOfMyFirstLine = 1;
860 unsigned int idOfMyLastLine = this->
sizeLocal();
861 unsigned int numParams = 1;
865 std::string tmpString;
868 *filePtrSet.
ifsVar >> tmpString;
872 *filePtrSet.
ifsVar >> tmpString;
877 *filePtrSet.
ifsVar >> tmpString;
879 unsigned int posInTmpString = 6;
882 char nPositionsString[tmpString.size()-posInTmpString+1];
883 unsigned int posInPositionsString = 0;
886 nPositionsString[posInPositionsString++] = tmpString[posInTmpString++];
887 }
while (tmpString[posInTmpString] !=
',');
888 nPositionsString[posInPositionsString] =
'\0';
892 char nParamsString[tmpString.size()-posInTmpString+1];
893 unsigned int posInParamsString = 0;
896 nParamsString[posInParamsString++] = tmpString[posInTmpString++];
897 }
while (tmpString[posInTmpString] !=
')');
898 nParamsString[posInParamsString] =
'\0';
901 unsigned int sizeOfVecInFile = (
unsigned int) strtod(nPositionsString,NULL);
902 unsigned int numParamsInFile = (
unsigned int) strtod(nParamsString, NULL);
906 <<
", sizeOfVecInFile = " << sizeOfVecInFile
907 <<
", numParamsInFile = " << numParamsInFile
908 <<
", this->sizeLocal() = " << this->
sizeLocal()
916 queso_require_equal_to_msg(numParamsInFile, numParams,
"number of parameters of vec in file is different than number of parameters in this vec object");
919 unsigned int maxCharsPerLine = 64*numParams;
921 unsigned int lineId = 0;
922 while (lineId < idOfMyFirstLine) {
923 filePtrSet.
ifsVar->ignore(maxCharsPerLine,
'\n');
929 <<
": beginning to read input actual data"
936 *filePtrSet.
ifsVar >> tmpString;
940 *filePtrSet.
ifsVar >> tmpString;
945 std::streampos tmpPos = filePtrSet.
ifsVar->tellg();
946 filePtrSet.
ifsVar->seekg(tmpPos+(std::streampos)2);
950 <<
": beginning to read lines with numbers only"
951 <<
", lineId = " << lineId
952 <<
", idOfMyFirstLine = " << idOfMyFirstLine
953 <<
", idOfMyLastLine = " << idOfMyLastLine
957 while (lineId <= idOfMyLastLine) {
958 *filePtrSet.
ifsVar >> (*this)[lineId - idOfMyFirstLine];
982 while ((i < size) && (result ==
false)) {
983 result = ( (*this)[i] < rhs[i] );
998 while ((i < size) && (result ==
false)) {
999 result = ( (*this)[i] > rhs[i] );
1011 bool result =
false;
1014 while ((i < size) && (result ==
false)) {
1015 result = ( (*this)[i] <= rhs[i] );
1027 bool result =
false;
1030 while ((i < size) && (result ==
false)) {
1031 result = ( (*this)[i] >= rhs[i] );
1041 return gsl_vector_max(
m_vec );
1047 return gsl_vector_min(
m_vec );
1053 return gsl_vector_max_index(
m_vec );
1059 return gsl_vector_min_index(
m_vec );
1085 unsigned int size = abs_of_this_vec.
sizeLocal();
1087 for(
unsigned int i = 0; i < size; ++i )
1089 abs_of_this_vec[i] = std::fabs( (*
this)[i] );
1092 return abs_of_this_vec;
1143 for (
unsigned int i = 0; i < size1; ++i) {
1144 result += x[i]*y[i];
1175 for (
unsigned int i = 0; i < size1; ++i) {
1176 if (lhs[i] != rhs[i]) {
double norm2() const
Returns the 2-norm (Euclidean norm) of the vector.
const RngBase * rngObject() const
Access to the RNG object.
int NumProc() const
Returns total number of processes.
bool atLeastOneComponentSmallerOrEqualThan(const GslVector &rhs) const
This function returns true if at least one component of this is smaller than or equal to the respecti...
double norm2Sq() const
Returns the 2-norm squared of this vector.
GslMatrix operator+(const GslMatrix &m1, const GslMatrix &m2)
void cwExtract(unsigned int initialPos, GslVector &vec) const
This function sets the values of this starting at position initialPos ans saves them in vector vec...
int NumMyElements() const
Returns the number of elements owned by the calling processor.
const Map & m_map
Mapping variable.
GslMatrix operator-(const GslMatrix &m1, const GslMatrix &m2)
int getMaxValueIndex() const
This function returns the index of the maximum value in the vector this.
double getMaxValue() const
Returns the maximum value in the vector this.
void matlabDiff(unsigned int firstPositionToStoreDiff, double valueForRemainderPosition, GslVector &outputVec) const
void Barrier() const
Pause every process in *this communicator until all the processes reach this point.
void print(std::ostream &os) const
Print method. Defines the behavior of the std::ostream << operator inherited from the Object class...
unsigned int numOfProcsForStorage() const
void cwSetGamma(const GslVector &a, const GslVector &b)
This function returns a random variate from the gamma distribution with vector parameters a and b...
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
#define queso_require_greater_equal_msg(expr1, expr2, msg)
int getMinValueIndex() const
This function returns the index of the minimum value in the vector this.
int subRank() const
Access function for sub-rank.
std::ifstream * ifsVar
Provides a stream interface to read data from files.
void getMinValueAndIndex(double &value, int &index)
This function returns minimum value in the vector this and its the index.
This (virtual) class sets up the environment underlying the use of the QUESO library by an executable...
void copy(const GslVector &src)
This function copies the elements of the vector src into this.
int MyPID() const
Return my process ID.
void cwInvert()
This function inverts component-wise the element values of this.
void cwSetInverseGamma(const GslVector &alpha, const GslVector &beta)
This function returns a random variate from the inverse gamma distribution with vector parameters alp...
int NumGlobalElements() const
Returns the total number of elements across all processors.
void cwSetBeta(const GslVector &alpha, const GslVector &beta)
This function returns a random variate from the beta distribution, with vector parameters alpha and b...
void cwSetGaussian(double mean, double stdDev)
This function sets component-wise Gaussian random variates, with mean mean and standard deviation std...
GslMatrix operator*(double a, const GslMatrix &mat)
#define queso_require_equal_to_msg(expr1, expr2, msg)
double norm1() const
Returns the 1-norm of the vector.
bool openInputFile(const std::string &fileName, const std::string &fileType, const std::set< unsigned int > &allowedSubEnvIds, FilePtrSetStruct &filePtrSet) const
Opens an input file.
bool getPrintHorizontally() const
Checks if vector is printed horizontally.
#define queso_require_less_msg(expr1, expr2, msg)
GslVector & operator=(const GslVector &rhs)
Copies values from vector rhs to this.
Class for vector operations (virtual).
void subWriteContents(const std::string &varNamePrefix, const std::string &fileName, const std::string &fileType, const std::set< unsigned int > &allowedSubEnvIds) const
A class for partitioning vectors and matrices.
Struct for handling data input and output from files.
GslVector & operator*=(double a)
Stores in this the coordinate-wise multiplication of this and a.
Class for vector operations using GSL library.
GslVector & operator/=(double a)
Stores in this the coordinate-wise division of this by a.
std::ostream & operator<<(std::ostream &os, const BaseEnvironment &obj)
#define RawValue_MPI_DOUBLE
void cwSetUniform(const GslVector &aVec, const GslVector &bVec)
This function sets component-wise a number uniformly distributed in the range of elements of [aVec...
void cwSet(double value)
Component-wise sets all values to this with value.
GslVector & operator+=(const GslVector &rhs)
Stores in this the coordinate-wise addition of this and rhs.
virtual double uniformSample() const =0
Samples a value from a uniform distribution.
const BaseEnvironment & m_env
Environment variable.
const std::string & subIdString() const
Access to the attribute m_subIdString; which stores the string for the sub-environment, and it will be used, for instance, to create the output files for each sub-environment.
virtual double betaSample(double alpha, double beta) const =0
Samples a value from a Beta distribution.
#define queso_require_greater_msg(expr1, expr2, msg)
bool operator==(const GslVector &lhs, const GslVector &rhs)
void closeFile(FilePtrSetStruct &filePtrSet, const std::string &fileType) const
Closes the file.
int fullRank() const
Returns the process full rank.
void cwSetConcatenated(const GslVector &v1, const GslVector &v2)
This function concatenates GslVector v1 and GslVector v2 into this.
void cwSqrt()
This function returns component-wise the square-root of this.
void matlabLinearInterpExtrap(const GslVector &x1Vec, const GslVector &y1Vec, const GslVector &x2Vec)
std::ofstream * ofsVar
Provides a stream interface to write data to files.
bool m_printHorizontally
Flag for either or not print this matrix horizontally.
double scalarProduct(const GslVector &x, const GslVector &y)
bool openOutputFile(const std::string &fileName, const std::string &fileType, const std::set< unsigned int > &allowedSubEnvIds, bool writeOver, FilePtrSetStruct &filePtrSet) const
Opens an output file for each sub-environment that was chosen to send data to the file...
bool atLeastOneComponentSmallerThan(const GslVector &rhs) const
This function returns true if at least one component of this is smaller than the respective component...
#define queso_require_less_equal_msg(expr1, expr2, msg)
GslVector operator/(double a, const GslVector &x)
unsigned int displayVerbosity() const
virtual void base_copy(const Vector &src)
Copies base data from vector src to this vector.
void getMaxValueAndIndex(double &value, int &index)
This function returns maximum value in the vector this and its the index.
unsigned int sizeLocal() const
Returns the length of this vector.
virtual double gammaSample(double a, double b) const =0
Samples a value from a Gamma distribution.
#define queso_require_msg(asserted, msg)
The QUESO MPI Communicator Class.
gsl_vector * m_vec
GSL vector.
bool getPrintScientific() const
Checks if the vector should be printed in Scientific Notation.
void Gather(void *sendbuf, int sendcnt, RawType_MPI_Datatype sendtype, void *recvbuf, int recvcount, RawType_MPI_Datatype recvtype, int root, const char *whereMsg, const char *whatMsg) const
Gather values from each process to collect on all processes.
void sort()
This function sorts the elements of the vector this in ascending numerical order. ...
void setPrintHorizontally(bool value) const
Determines whether vector should be printed horizontally.
bool atLeastOneComponentBiggerOrEqualThan(const GslVector &rhs) const
This function returns true if at least one component of this is bigger than or equal to the respectiv...
double normInf() const
Returns the infinity-norm (maximum norm) of the vector.
void subReadContents(const std::string &fileName, const std::string &fileType, const std::set< unsigned int > &allowedSubEnvIds)
bool atLeastOneComponentBiggerThan(const GslVector &rhs) const
This function returns true if at least one component of this is bigger than the respective component ...
unsigned int sizeGlobal() const
Returns the global length of this vector.
virtual double gaussianSample(double stdDev) const =0
Samples a value from a Gaussian distribution with standard deviation given by stdDev.
double getMinValue() const
Returns minimum value in the vector this.
void Allreduce(void *sendbuf, void *recvbuf, int count, RawType_MPI_Datatype datatype, RawType_MPI_Op op, const char *whereMsg, const char *whatMsg) const
Combines values from all processes and distributes the result back to all processes.
void setPrintScientific(bool value) const
Determines whether vector should be printed in Scientific Notation.
void mpiAllReduce(RawType_MPI_Op mpiOperation, const MpiComm &opComm, GslVector &resultVec) const
void mpiAllQuantile(double probability, const MpiComm &opComm, GslVector &resultVec) const
bool m_printScientific
Flag for either or not print this matrix in scientific notation.
void mpiBcast(int srcRank, const MpiComm &bcastComm)
double sumOfComponents() const
Returns the sum of the components of the vector.
GslVector & operator-=(const GslVector &rhs)
Stores in this the coordinate-wise subtraction of this by rhs.
gsl_vector * data() const
GslVector()
Default Constructor.
GslVector abs() const
This function returns absolute value of elements in this.
void Bcast(void *buffer, int count, RawType_MPI_Datatype datatype, int root, const char *whereMsg, const char *whatMsg) const
Broadcast values from the root process to the slave processes.