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]) {
unsigned int displayVerbosity() const
bool getPrintScientific() const
Checks if the vector should be printed in Scientific Notation.
void subReadContents(const std::string &fileName, const std::string &fileType, const std::set< unsigned int > &allowedSubEnvIds)
bool getPrintHorizontally() const
Checks if vector is printed horizontally.
std::ofstream * ofsVar
Provides a stream interface to write data to files.
int NumGlobalElements() const
Returns the total number of elements across all processors.
void cwSetConcatenated(const GslVector &v1, const GslVector &v2)
This function concatenates GslVector v1 and GslVector v2 into this.
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.
Class for vector operations (virtual).
int subRank() const
Access function for sub-rank.
void cwSqrt()
This function returns component-wise the square-root of this.
void sort()
This function sorts the elements of the vector this in ascending numerical order. ...
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.
void cwSetGamma(const GslVector &a, const GslVector &b)
This function returns a random variate from the gamma distribution with vector parameters a and b...
bool atLeastOneComponentSmallerThan(const GslVector &rhs) const
This function returns true if at least one component of this is smaller than the respective component...
GslVector & operator+=(const GslVector &rhs)
Stores in this the coordinate-wise addition of this and rhs.
virtual double betaSample(double alpha, double beta) const =0
Samples a value from a Beta distribution.
double getMinValue() const
Returns minimum value in the vector this.
void cwSetInverseGamma(const GslVector &alpha, const GslVector &beta)
This function returns a random variate from the inverse gamma distribution with vector parameters alp...
virtual double gaussianSample(double stdDev) const =0
Samples a value from a Gaussian distribution with standard deviation given by stdDev.
void getMaxValueAndIndex(double &value, int &index)
This function returns maximum value in the vector this and its the index.
void cwSetBeta(const GslVector &alpha, const GslVector &beta)
This function returns a random variate from the beta distribution, with vector parameters alpha and b...
A class for partitioning vectors and matrices.
#define queso_require_less_msg(expr1, expr2, msg)
const BaseEnvironment & m_env
Environment variable.
int fullRank() const
Returns the process full rank.
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...
GslVector & operator-=(const GslVector &rhs)
Stores in this the coordinate-wise subtraction of this by rhs.
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...
void matlabLinearInterpExtrap(const GslVector &x1Vec, const GslVector &y1Vec, const GslVector &x2Vec)
GslVector()
Default Constructor.
GslMatrix operator*(double a, const GslMatrix &mat)
unsigned int sizeGlobal() const
Returns the global length of this vector.
#define queso_require_less_equal_msg(expr1, expr2, msg)
void setPrintScientific(bool value) const
Determines whether vector should be printed in Scientific Notation.
virtual double gammaSample(double a, double b) const =0
Samples a value from a Gamma distribution.
#define queso_require_msg(asserted, msg)
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.
double getMaxValue() const
Returns the maximum value in the vector this.
#define queso_require_equal_to_msg(expr1, expr2, msg)
void Barrier() const
Pause every process in *this communicator until all the processes reach this point.
void mpiAllQuantile(double probability, const MpiComm &opComm, GslVector &resultVec) const
#define queso_require_greater_equal_msg(expr1, expr2, msg)
void cwSet(double value)
Component-wise sets all values to this with value.
virtual void base_copy(const Vector &src)
Copies base data from vector src to this vector.
GslVector operator/(double a, const GslVector &x)
void closeFile(FilePtrSetStruct &filePtrSet, const std::string &fileType) const
Closes the file.
std::ostream & operator<<(std::ostream &os, const BaseEnvironment &obj)
GslVector abs() const
This function returns absolute value of elements in this.
double norm2Sq() const
Returns the 2-norm squared of this vector.
const Map & m_map
Mapping variable.
void setPrintHorizontally(bool value) const
Determines whether vector should be printed horizontally.
Struct for handling data input and output from files.
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
GslVector & operator=(const GslVector &rhs)
Copies values from vector rhs to this.
int getMaxValueIndex() const
This function returns the index of the maximum value in the vector this.
void mpiAllReduce(RawType_MPI_Op mpiOperation, const MpiComm &opComm, GslVector &resultVec) const
This (virtual) class sets up the environment underlying the use of the QUESO library by an executable...
void print(std::ostream &os) const
Print method. Defines the behavior of the std::ostream << operator inherited from the Object class...
int MyPID() const
Return my process ID.
double sumOfComponents() const
Returns the sum of the components of the vector.
int getMinValueIndex() const
This function returns the index of the minimum value in the vector this.
bool m_printScientific
Flag for either or not print this matrix in scientific notation.
gsl_vector * data() const
GslVector & operator/=(double a)
Stores in this the coordinate-wise division of this by a.
unsigned int sizeLocal() const
Returns the length of this vector.
int NumProc() const
Returns total number of processes.
std::ifstream * ifsVar
Provides a stream interface to read data from files.
GslVector & operator*=(double a)
Stores in this the coordinate-wise multiplication of this and a.
GslMatrix operator-(const GslMatrix &m1, const GslMatrix &m2)
void copy(const GslVector &src)
This function copies the elements of the vector src into this.
void cwInvert()
This function inverts component-wise the element values of this.
double normInf() const
Returns the infinity-norm (maximum norm) of the vector.
GslMatrix operator+(const GslMatrix &m1, const GslMatrix &m2)
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.
virtual double uniformSample() const =0
Samples a value from a uniform distribution.
bool atLeastOneComponentBiggerThan(const GslVector &rhs) const
This function returns true if at least one component of this is bigger than the respective component ...
const RngBase * rngObject() const
Access to the RNG object.
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...
gsl_vector * m_vec
GSL vector.
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...
Class for vector operations using GSL library.
double norm1() const
Returns the 1-norm of the vector.
bool operator==(const GslVector &lhs, const GslVector &rhs)
void subWriteContents(const std::string &varNamePrefix, const std::string &fileName, const std::string &fileType, const std::set< unsigned int > &allowedSubEnvIds) const
The QUESO MPI Communicator Class.
int NumMyElements() const
Returns the number of elements owned by the calling processor.
unsigned int numOfProcsForStorage() const
void matlabDiff(unsigned int firstPositionToStoreDiff, double valueForRemainderPosition, GslVector &outputVec) const
bool m_printHorizontally
Flag for either or not print this matrix horizontally.
void mpiBcast(int srcRank, const MpiComm &bcastComm)
#define queso_require_greater_msg(expr1, expr2, msg)
void getMinValueAndIndex(double &value, int &index)
This function returns minimum value in the vector this and its the index.
void cwSetUniform(const GslVector &aVec, const GslVector &bVec)
This function sets component-wise a number uniformly distributed in the range of elements of [aVec...
double norm2() const
Returns the 2-norm (Euclidean norm) of the vector.
#define RawValue_MPI_DOUBLE
void cwSetGaussian(double mean, double stdDev)
This function sets component-wise Gaussian random variates, with mean mean and standard deviation std...
bool openInputFile(const std::string &fileName, const std::string &fileType, const std::set< unsigned int > &allowedSubEnvIds, FilePtrSetStruct &filePtrSet) const
Opens an input file.
double scalarProduct(const GslVector &x, const GslVector &y)