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()))
60 m_vec (gsl_vector_calloc(map.NumGlobalElements()))
80 m_vec (gsl_vector_calloc(map.NumGlobalElements()))
90 for (
unsigned int i = 0; i <
m_vec->size; ++i) {
91 double alpha = (double) i / ((
double)
m_vec->size - 1.);
92 (*this)[i] = (1.-alpha)*d1 + alpha*d2;
103 m_vec (gsl_vector_calloc(v.sizeLocal()))
113 for (
unsigned int i = 0; i <
m_vec->size; ++i) {
114 double alpha = (double) i / ((
double)
m_vec->size - 1.);
115 (*this)[i] = (1. - alpha) * start + alpha * end;
126 m_vec (gsl_vector_calloc(v.sizeLocal()))
169 iRC = gsl_vector_scale(
m_vec,a);
189 for (
unsigned int i = 0; i < size1; ++i) {
190 (*this)[i] *= rhs[i];
203 for (
unsigned int i = 0; i < size1; ++i) {
204 (*this)[i] /= rhs[i];
234 iRC = gsl_vector_memcpy(this->
m_vec, src.
m_vec);
283 return std::sqrt(this->
norm2Sq());
292 for (
unsigned int i = 0; i < size; ++i) {
293 result += fabs((*
this)[i]);
306 for (
unsigned int i = 0; i < size; ++i) {
307 aux = fabs((*
this)[i]);
308 if (aux > result) result = aux;
319 for (
unsigned int i = 0; i < size; ++i) {
320 result += (*this)[i];
330 for (
unsigned int i = 0; i < size; ++i) {
340 for (
unsigned int i = 0; i < this->
sizeLocal(); ++i) {
350 for (
unsigned int i = 0; i < this->
sizeLocal(); ++i) {
359 for (
unsigned int i = 0; i < this->
sizeLocal(); ++i) {
372 double tmpSample = 0.;
373 for (
unsigned int i = 0; i < this->
sizeLocal(); ++i) {
379 <<
", alpha[i] = " << alpha[i]
380 <<
", beta[i] = " << beta[i]
381 <<
", sample = " << tmpSample
384 if ((alpha[i] == 1. ) &&
386 if (tmpSample == 1.) {
391 <<
", alpha[i] = " << alpha[i]
392 <<
", beta[i] = " << beta[i]
393 <<
", sample = " << tmpSample
397 std::cerr <<
"Hitting 'sample = 1' in GslVector::cwSetBeta()"
400 <<
", alpha[i] = " << alpha[i]
401 <<
", beta[i] = " << beta[i]
402 <<
", sample = " << tmpSample
406 }
while (tmpSample == 1.);
407 std::cerr <<
"Code was able to get 'sample != 1' in GslVector::cwSetBeta()"
410 <<
", alpha[i] = " << alpha[i]
411 <<
", beta[i] = " << beta[i]
412 <<
", sample = " << tmpSample
417 (*this)[i] = tmpSample;
429 for (
unsigned int i = 0; i < this->
sizeLocal(); ++i) {
442 for (
unsigned int i = 0; i < this->
sizeLocal(); ++i) {
453 for (
unsigned int i = 0; i < v1.
sizeLocal(); ++i) {
457 for (
unsigned int i = 0; i < v2.
sizeLocal(); ++i) {
467 unsigned int cummulativeSize = 0;
468 for (
unsigned int i = 0; i < vecs.size(); ++i) {
470 for (
unsigned int j = 0; j < vecs[i]->sizeLocal(); ++j) {
471 (*this)[cummulativeSize+j] = tmpVec[j];
488 for (
unsigned int i = 0; i < vec.
sizeLocal(); ++i) {
489 (*this)[initialPos+i] = vec[i];
502 for (
unsigned int i = 0; i < vec.
sizeLocal(); ++i) {
503 vec[i] = (*this)[initialPos+i];
513 for (
unsigned int i = 0; i < size; ++i) {
514 (*this)[i] = 1./(*this)[i];
524 for (
unsigned int i = 0; i < size; ++i) {
525 (*this)[i] = sqrt((*
this)[i]);
533 unsigned int firstPositionToStoreDiff,
534 double valueForRemainderPosition,
543 for (
unsigned int i = 0; i < (size-1); ++i) {
544 outputVec[firstPositionToStoreDiff+i] = (*this)[i+1]-(*this)[i];
546 if (firstPositionToStoreDiff == 0) {
547 outputVec[size-1] = valueForRemainderPosition;
550 outputVec[0] = valueForRemainderPosition;
568 for (
unsigned int i = 1; i < x1Vec.
sizeLocal(); ++i) {
572 for (
unsigned int id2 = 0; id2 < x2Vec.
sizeLocal(); ++id2) {
573 double x2 = x2Vec[id2];
574 unsigned int id1 = 0;
576 for (id1 = 0; id1 < x1Vec.
sizeLocal(); ++id1) {
577 if (x2 <= x1Vec[id1]) {
582 bool makeLinearModel =
false;
587 if (x2 == x1Vec[id1]) {
588 (*this)[id2] = y1Vec[id1];
590 else if (x2 < x1Vec[0]) {
592 makeLinearModel =
true;
598 else if (found1 ==
true) {
600 makeLinearModel =
true;
608 makeLinearModel =
true;
615 if (makeLinearModel) {
616 double rate = (yb-ya)/(xb-xa);
617 (*this)[id2] = ya + (x2-xa)*rate;
629 gsl_sort_vector(
m_vec);
638 if (bcastComm.
MyPID() < 0)
return;
644 double localNumNodes = 1.;
645 double totalNumNodes = 0.;
647 "GslVector::mpiBcast()",
648 "failed MPI.Allreduce() for numNodes");
652 double localVectorSize = this->
sizeLocal();
653 double sumOfVectorSizes = 0.;
655 "GslVector::mpiBcast()",
656 "failed MPI.Allreduce() for vectorSize");
658 if ( ((
unsigned int) sumOfVectorSizes) != ((
unsigned int)(totalNumNodes*localVectorSize)) ) {
659 std::cerr <<
"rank " << bcastComm.
MyPID()
660 <<
": sumOfVectorSizes = " << sumOfVectorSizes
661 <<
", totalNumNodes = " << totalNumNodes
662 <<
", localVectorSize = " << localVectorSize
666 queso_require_equal_to_msg(((
unsigned int) sumOfVectorSizes), ((
unsigned int)(totalNumNodes*localVectorSize)),
"inconsistent vectorSize");
669 std::vector<double> dataBuffer((
unsigned int) localVectorSize, 0.);
670 if (bcastComm.
MyPID() == srcRank) {
671 for (
unsigned int i = 0; i < dataBuffer.size(); ++i) {
672 dataBuffer[i] = (*this)[i];
677 "GslVector::mpiBcast()",
678 "failed MPI.Bcast()");
680 if (bcastComm.
MyPID() != srcRank) {
681 for (
unsigned int i = 0; i < dataBuffer.size(); ++i) {
682 (*this)[i] = dataBuffer[i];
693 if (opComm.
MyPID() < 0)
return;
698 for (
unsigned int i = 0; i < size; ++i) {
699 double srcValue = (*this)[i];
700 double resultValue = 0.;
701 opComm.
Allreduce<
double>(&srcValue, &resultValue, (int) 1, mpiOperation,
702 "GslVector::mpiAllReduce()",
703 "failed MPI.Allreduce()");
704 resultVec[i] = resultValue;
714 if (opComm.
MyPID() < 0)
return;
721 for (
unsigned int i = 0; i < size; ++i) {
722 double auxDouble = (int) (*
this)[i];
723 std::vector<double> vecOfDoubles(opComm.
NumProc(),0.);
724 opComm.
Gather<
double>(&auxDouble, 1, &vecOfDoubles[0], (int) 1, 0,
725 "GslVector::mpiAllQuantile()",
726 "failed MPI.Gather()");
728 std::sort(vecOfDoubles.begin(), vecOfDoubles.end());
730 double result = vecOfDoubles[(
unsigned int)( probability*((
double)(vecOfDoubles.size()-1)) )];
733 "GslVector::mpiAllQuantile()",
734 "failed MPI.Bcast()");
736 resultVec[i] = result;
753 std::ostream::fmtflags curr_fmt = os.flags();
758 unsigned int savedPrecision = os.precision();
762 for (
unsigned int i = 0; i < size; ++i) {
763 os << std::scientific << (*this)[i]
768 for (
unsigned int i = 0; i < size; ++i) {
769 os << std::scientific << (*this)[i]
774 os.precision(savedPrecision);
780 for (
unsigned int i = 0; i < size; ++i) {
781 os << std::dec << (*this)[i]
786 for (
unsigned int i = 0; i < size; ++i) {
787 os << std::dec << (*this)[i]
804 const std::string& varNamePrefix,
805 const std::string& fileName,
806 const std::string& fileType,
807 const std::set<unsigned int>& allowedSubEnvIds)
const
829 *filePtrSet.
ofsVar << *
this;
834 *filePtrSet.
ofsVar <<
"];\n";
844 const std::string& fileName,
845 const std::string& fileType,
846 const std::set<unsigned int>& allowedSubEnvIds)
860 unsigned int idOfMyFirstLine = 1;
861 unsigned int idOfMyLastLine = this->
sizeLocal();
862 unsigned int numParams = 1;
866 std::string tmpString;
869 *filePtrSet.
ifsVar >> tmpString;
873 *filePtrSet.
ifsVar >> tmpString;
878 *filePtrSet.
ifsVar >> tmpString;
880 unsigned int posInTmpString = 6;
883 char nPositionsString[tmpString.size()-posInTmpString+1];
884 unsigned int posInPositionsString = 0;
887 nPositionsString[posInPositionsString++] = tmpString[posInTmpString++];
888 }
while (tmpString[posInTmpString] !=
',');
889 nPositionsString[posInPositionsString] =
'\0';
893 char nParamsString[tmpString.size()-posInTmpString+1];
894 unsigned int posInParamsString = 0;
897 nParamsString[posInParamsString++] = tmpString[posInTmpString++];
898 }
while (tmpString[posInTmpString] !=
')');
899 nParamsString[posInParamsString] =
'\0';
902 unsigned int sizeOfVecInFile = (
unsigned int) strtod(nPositionsString,NULL);
903 unsigned int numParamsInFile = (
unsigned int) strtod(nParamsString, NULL);
907 <<
", sizeOfVecInFile = " << sizeOfVecInFile
908 <<
", numParamsInFile = " << numParamsInFile
909 <<
", this->sizeLocal() = " << this->
sizeLocal()
917 queso_require_equal_to_msg(numParamsInFile, numParams,
"number of parameters of vec in file is different than number of parameters in this vec object");
920 unsigned int maxCharsPerLine = 64*numParams;
922 unsigned int lineId = 0;
923 while (lineId < idOfMyFirstLine) {
924 filePtrSet.
ifsVar->ignore(maxCharsPerLine,
'\n');
930 <<
": beginning to read input actual data"
937 *filePtrSet.
ifsVar >> tmpString;
941 *filePtrSet.
ifsVar >> tmpString;
946 std::streampos tmpPos = filePtrSet.
ifsVar->tellg();
947 filePtrSet.
ifsVar->seekg(tmpPos+(std::streampos)2);
951 <<
": beginning to read lines with numbers only"
952 <<
", lineId = " << lineId
953 <<
", idOfMyFirstLine = " << idOfMyFirstLine
954 <<
", idOfMyLastLine = " << idOfMyLastLine
958 while (lineId <= idOfMyLastLine) {
959 *filePtrSet.
ifsVar >> (*this)[lineId - idOfMyFirstLine];
983 while ((i < size) && (result ==
false)) {
984 result = ( (*this)[i] < rhs[i] );
999 while ((i < size) && (result ==
false)) {
1000 result = ( (*this)[i] > rhs[i] );
1012 bool result =
false;
1015 while ((i < size) && (result ==
false)) {
1016 result = ( (*this)[i] <= rhs[i] );
1028 bool result =
false;
1031 while ((i < size) && (result ==
false)) {
1032 result = ( (*this)[i] >= rhs[i] );
1042 return gsl_vector_max(
m_vec );
1048 return gsl_vector_min(
m_vec );
1054 return gsl_vector_max_index(
m_vec );
1060 return gsl_vector_min_index(
m_vec );
1086 unsigned int size = abs_of_this_vec.
sizeLocal();
1088 for(
unsigned int i = 0; i < size; ++i )
1090 abs_of_this_vec[i] = std::fabs( (*
this)[i] );
1093 return abs_of_this_vec;
1145 for (
unsigned int i = 0; i < size1; ++i) {
1146 result += x[i]*y[i];
1177 for (
unsigned int i = 0; i < size1; ++i) {
1178 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)