Getting Started | Documentation | Glish | Learn More | Programming | Contact Us |
Version 1.9 Build 1367 |
|
August 5, 1992
This document is an introduction to the Rogue Wave Math Libraries Math.h++, Matrix.h++, and Linpack.h++ for C++. It contains information on how to compile a program that implements these libraries, as well as several examples of C++ programs. Examples have been taken from the three manuals Math.h++, Matrix.h++, and Linpack.h++. Specifically, the first eight are taken verbatim from section 14 of Math.h++. Some familarity with C++ programming is assumed; the beginning programmer may have difficulty with some of the terminology or code in this document.
Compilation:
To compile any C++ file that includes the RW math libraries, type:
The CC specifies the compiler, in this case, Sun C++; I briefly attempted to use the Gnu g++ compiler, but I could not make it the program compile correctly. The order of the ``-l'' libraries is very important. If some of these libraries are not used in the .cc file, the appropriate ``-l'' may be omitted. However, the ``-D ATT2 '' is absolutely mandatory; without it, the compiler will return strange error messages, and the file may not compile correctly. The -o <filename> option instructs the compiler to name the executable output file <filename>. Without this option, the compiler automatically creates an executable file named a.out.
I have used some of the examples in the manual to test the math libraries. Following is the source code of each and my impressions thereof.
Example 1: Vectors (Math.h++, §14.1, pg.81)
/* * Example program using class DoubleVec * to do simple statistics with double precision vectors. */ // Include the DoubleVec class header file. #include <rw/dvec.h> #include <rw/rstream.h> main() { // Declare the double precision vector V DoubleVec V; // Read the vector from standard input: cin >> V; /* * Output the length of the vector. * This illustrates the use of a DoubleVec * member function: */ cout << "The length of the vector is:\t" << V.length() << "\n"; /* * Calculate the mean of the vector. * This illustrates the use * of the global function mean(): */ double avg = mean(V); cout << "The mean of the vector is:\t" << avg << "\n"; /* * Find the index of the maximum element of * the vector and output the value: */ int maxV = maxIndex(V); cout << "The maximum value of the vector is:\t" << V(maxV) << "\n"; // Output the variance of the vector: cout << "The variance of the vector is:\t" << variance(V) << "\n"; /* * Here is how to demean and normalize a * vector and print out the result: */ V -= mean(V); // Remove the mean maxV = maxIndex(V); // Max element if ( V(maxV) != 0) { V /= fabs(V(maxV)); cout << "The normalized demeaned vector is:\n" << V; } return 0; }
Example 2: Matrices (Math.h++, §14.2, pg. 83)
/* * Example program showing the use of matrix * classes to multiply a double precision matrix by a vector. */ // Include the header files for double precision matrices // and vectors. #include <rw/dvec.h> #include <rw/dgenmat.h> #include <rw/rstream.h> main() { // Define the vector and the matrix: DoubleVec vector; DoubleGenMat matrix; // Read in the vector, then the matrix: cin >> vector >> matrix; // Matrix multiply (the inner product): DoubleVec answer = product(matrix, vector); // print out the results: cout << "answer = \n" << answer << "\n"; return 0; }
Example 3: FFT's (Math.h++, §14.3, pg. 84)
/* * Using the Complex FFT server class. */ // Include the header file for complex FFT server class: #include <rw/cfft.h> // Header file for complex vectors: #include <rw/cvec.h> #include <math.h> #include <rw/rstream.h> main() { // Set series length: unsigned npts = 12; // Allocate a server: DComplexFFTServer server; /* * Create two series, one a cosine series, the * other a sine series. This is done by first * using a constructor for a DoubleVec with * increasing elements, taking the cosine (or * sine) of this vector, then constructing a * complex vector from that. */ // one cycle: DComplexVec a = cos( DoubleVec(npts,0,2.0*M_PI/npts) ); // two cycles: DComplexVec a2 = sin( DoubleVec(npts,0,4.0*M_PI/npts) ); /* * Calculate the superposition of the two. Note that we * are adding two complex vectors here with one * simple statement: */ DComplexVec asum = a + a2; // Output the vectors: cout << "a:\n" << a << "\n"; cout << "a2:\n" << a2 << "\n"; cout << "asum:\n" << asum << "\n"; /* * Print out the transforms, normalized by the * number of points. The explicit DComplex constructors * are necessary for Zortech, which does not have * real to complex type conversion. */ cout << "\nTransform of a:\n"; cout << server.fourier(a)/DComplex(npts,0); cout << "\nTransform of a2:\n"; cout << server.fourier(a2)/DComplex(npts,0); DComplexVec asumFFT = server.fourier(asum); cout << "\nTransform of asum:\n"; cout << asumFFT/DComplex(npts,0); /* * Check Parseval's Theorem: First, calculate the variance of * the series asum, using the global function variance(): */ cout << "\nOriginal Variance: " << variance(asum) << "\n"; /* * Next calculate the spectral variance of the fourier * transformed series, using the global function * spectralVariance(): */ double var = spectralVariance(asumFFT/DComplex(npts,0)); cout << "Spectral Variance: " << var << "\n\n"; /* * Print out the normalized back transform of * asumFFT; This should equal the original asum: */ cout << "\nBack Transform of asum:\n" <<server.ifourier(asumFFT)/DComplex(npts,0); return 0; }
Example 4: Implicit Type Conversion (Math.h++, §14.4, pg. 86)
/* * Example program illustrating implicit type conversion. * C++ provides for implicit conversion between data types, * performed by the compiler without programmer intervention. * If the types of the arguments to a function do not match * the function prototype, a match is sought using * class-defined conversions. */ // Include the DoubleVec and IntVec class header files. #include <rw/dvec.h> #include <rw/ivec.h> #include <rw/rstream.h> // Initial data for the vectors double adata[] = {1., 3., -2., -6., 7.}; int idata[] = {2, 6, -4, 2, 1}; main() { // Construct the vectors V and iV from the // arrays defined above: DoubleVec V(adata, 5); IntVec iV(idata, 5); /* * Output the dot product of the two vectors. * Note that the function dot() has no prototype * dot(DoubleVec&, IntVec&). * However, the prototype * dot(DoubleVec&, DoubleVec&) * does exist. The conversion operator DoubleVec(), defined * for the class IntVec, provides for the implicit type * conversion: IntVec to DoubleVec. */ cout << dot (V, iV); return 0; }
Example 5: Persistence (saveOn / restoreFrom) (Math.h++, §14.5, pg. 87)
/* * Example program showing the use of the member functions * saveOn() and restoreFrom() for data storage and retrieval. */ #include <rw/ivec.h> #include <rw/pstream.h> #include <fstream.h> main() { /* * Construct an integer vector a, with 24 elements: * The first element has value 0; each succeeding element is * incremented by 1 */ IntVec a(24, 0, 1); // Store the vector to file "vec.dat" using class // RWpostream which saves in a portable ASCII format: { ofstream fstr("vec.dat", ios::out); RWpostream postr(fstr); a.saveOn(postr); } // A vector that has been saved using function saveOn() // may be restored using restoreFrom(): ifstream fstr("vec.dat", ios::in); // Construct a RWpistream from fstr RWpistream pistr(fstr); IntVec b; b.restoreFrom(pistr); // Restore from file "vec.dat" cout << a << NL; // Print the original 'a' cout << b << NL; // Print the restored 'b' return 0; }
Example 6: Class derivation (Math.h++, §14.6, pg. 88)
/* * This example shows how a CosVec class may be derived from * the class DoubleVec. Class CosVec objects are double * precision vectors that have been initialized with cosines * with a specified number of integral cycles. */ // Include the DoubleVec class header file. #include <rw/dvec.h> #include <rw/rstream.h> /* * Start derived class declarations: * DoubleVec is declared as a public base class for CosVec: */ class CosVec : public DoubleVec { public: // unadorned constructor: CosVec(); /* This is the useful constructor. We first use the * DoubleVec constructor * DoubleVec(n, val, by); * to create a vector of phases. Taking the cosine of this * vector will yield the desired initial values for the base * class of CosVec: */ CosVec (unsigned N, int cycles = 1) : DoubleVec( cos(DoubleVec(N,0,2*M_PI*cycles/N)) ) { } // Copy constructor: CosVec (const CosVec& v) : DoubleVec(v) { } }; main() { /* * Initialize a CosVec object with 12 elements and one * cycle: */ CosVec c(12,1); /* * Output the vector: note that class CosVec inherits * operator<< from the base class DoubleVec: */ cout << c; return 0; }
Example 7: Linear Algebra (Math.h++, §15.7, pg. 90)
/* * This example shows how to use the class DoubleGenFact to do * linear algebra with double precision matrices. */ // Include the DoubleGenFact class header file. #include <rw/dgenfct.h> // Include the double precision matrix header file: #include <rw/dgenmat.h> #include <rw/rstream.h> // Initial data for the vectors const double adata[] = {-3.0, 2.0, 1.0, 8.0, -7.0, 9.0, 5.0, 4.0, -6.0}; const double arhs[] = {6.0, 9.0, 1.0}; main() { // Construct a test matrix and print it out: DoubleGenMat testmat(adata,3,3); cout << "test matrix:\n" << testmat << "\n"; /* * Calculate and print the inverse. * Note that a type conversion occurs: * testmat is converted to type DoubleGenFact * before the inverse is computed: */ cout << "inverse:\n" << inverse(testmat); // Now construct a DoubleGenFact from the matrix: DoubleGenFact LUtest(testmat); // Once constructed, LUtest may be reused as required: // Find the determinant: cout << "\ndeterminant\n" << determinant(LUtest) << NL; // Solve the linear system testmat*x = rhs: DoubleVec rhs(arhs, 3); DoubleVec x = solve(LUtest, rhs); cout << "solution:\n" << x << "\n"; return 0; }
Example 8: Multiple inheritance (Math.h++, § 14.8, pg. 91)
/* * Example to illustrate the use of multiple inheritance to do * run time binding. Although the Rogue Wave Math.h++ Class * Library does not use any virtual functions, it is easy * enough to implement them in a derived class by using * Multiple Inheritance. */ #include <rw/cvec.h> #include <rw/dvec.h> #include <rw/rstream.h> /* * This is the abstract base class where the virtual functions * that we plan to use are declared (but not defined). */ class Vector { public: virtual void print() = 0; }; /* * Here are two classes that inherit not only the abstract * base class above, but also the properties of a Math.h++ * class. The first class will inherit DComplexVec, the * second DoubleVec. */ class CVector : public DComplexVec, public Vector { public: CVector(unsigned n, DComplex start, DComplex inc) : DComplexVec(n, start, inc) {} virtual void print(); }; class DVector : public DoubleVec, public Vector { public: DVector(unsigned n, double start, double inc) : DoubleVec(n, start, inc) {} virtual void print(); }; /* * A little test program to illustrate. */ main() { int type; cout << "Demonstration of runtime binding.\n"; cout << "Type 1 for a vector of complexes;\n"; cout << "Type 2 for a vector of doubles: "; cin >> type; Vector* avec; switch (type) { case 1 : cout << "Constructing a CVector\n"; avec = new CVector(10, DComplex(0,0), DComplex(1,0) ); break; default : cout << "Constructing a DVector\n"; avec = new DVector(10, 0, 1); } /* * Note that "avec" points to the *base* class "Vector". * The actual type of what it points to is unknown at * compile time --- it depends on what the user typed above. * Runtime binding happens for the virtual function print(): */ avec->print(); return 0; } void CVector::print() { cout <<"Vector of complex:\n"; // Inherit the operator<<() function of DComplexVec: cout << *this << "\n"; } void DVector::print() { cout <<"Vector of doubles:\n"; // Inherit the operator<<() function of DoubleVec: cout << *this << "\n"; }
LU-Decomposition: (Linpack.h++, §5.2, pg.LU-Decomposition: (Linpack.h++, pg. 50)
(Linpack.h++, §4.2, pg. 19)
(Linpack.h++, §5.3, pg. 25)
\* LU Decomposition: lsinc.cc *\ \* Used to test LU Decomposition classes *\ #include <rw/rstream.h> #include <rw/dlsinc.h> main(){ DoubleGenMat A; DoubleVec b; cin >> A >> b; DoubleLeastSqInc decomp(A.cols()); for (int i =0; i <A.rows(); i++) decomp.addEquation(A.row(i), b[i]); if (decomp.fail()) cout << "Can't find a solution" << NL; else { DoubleVec x = decomp.solve(); cout << "the solution is: " << NL << x << NL; } return 0; } \* LU-Decomposition - inc.cc *\ \* Used to test the LU incremental decomposition classes. *\ #include <rw/rstream.h> #include <rw/dlsinc.h> main(){ unsigned numUnknowns; cin >> numUnknowns; DoubleLeastSqInc S(numUnknowns); DoubleVec a; double b; while (cin >> a >> b){ S.addEquation(a,b); if (S.good()){ cout << "Solution:" << NL; cout << S.solve() << NL; cout << "Error: " << S.residualNorm() << NL; } } } \* LU Decomposition: several.cc *\ \* Used to test Lu incremental decomposion classes *\ #include <rw/rstream.h> #include <rw/dgenfct.h> main(){ DoubleGenMat A; cin >> A; DoubleGenFact fact(A); if (fact.good()){ DoubleVec b; while (cin >> b){ DoubleVec x = solve(fact,b); cout << x << NL; } } else { cout << "Not 'good'\n"; } }
Least Square Fit: (Linpack.h++, §5.1, pg. 24)
\* Least Square Fit: lsch.cc *\ \* Used to test linear least square fit class *\ #include <rw/rstream.h> #include <rw/flsch.h> main(){ FloatGenMat A; FloatVec b; cin >> A >> b; FloatLeastSqCh ls(A); FloatVec x = solve (ls,b); cout << x << NL; }
QR decomposition: (Linpack.h++, pg. 59)
\* QR decomposition: *\ \* Used to test QR decomposition class *\ #include <rw/rstream.h> #include <rw/dqr.h> #include <rw/dutrimat.h> main(){ DoubleGenMat A; cin >> A; DoubleQRDecomp QR; QR.doPivoting(FALSE); QR.factor(A); DoubleGenMat Rsq = QR; if (QR.rows() > QR. cols()) Rsq.resize(QR.cols(), QR.cols()); else Rsq.resize(QR.rows(), QR.rows()); DoubleUpperTriMat R = toUpperTriMat(Rsq); DoubleGenMat Q(QR.rows(),QR.rows()); DoubleVec e(QR.rows(),0); for (int i = 0; i < QR.rows();i++){ e[i]=1; Q.col(i) = QR.computeQy(e); e[i] =0; } cout << "Q is " << NL << Q << NL; cout << "R is " << NL << R << NL; }
Cholesky decomposition: (There is no code in the manuals for this example.)
/* Cholesky Decompositions: chdec.cc */ /* This file tests the RW Cholesky decomposition class. */ #include <rw/dch.h> #include <rw/dgenmat.h> #include <rw/rstream.h> #include <rw/dutrimat.h> main(){ DoubleGenMat A; cin >> A; cout << "Original:\n" << A << NL; DoubleChDecomp d; d.doPivoting(FALSE); d.factor(A); DoubleGenMat D(d); for (int i = -1; i > -(D.rows()); i--) D.diagonal(i) =0; if (d.isPD()==TRUE) cout << ".isPD sucessful\n"; if (d.fail()) cout << "Decomposition returned \"fail\"!\n"; DoubleGenMat Dt = transpose(D); cout << "Result:\n" << product(Dt,D); }
Singular Value Decomposition: (Linpack.h++, pg. 64)
/* Singular Value Decomposition: svd.cc */ /* Used to test singular value decomposion class */ #include <rw/rstream.h> #include <rw/dsv.h> main(){ DoubleGenMat A; cin >> A; DoubleSVDecomp SVD; SVD.computeLeftVecs(FALSE); SVD.computeRightVecs(FALSE); SVD.factor(A); if (SVD.fail()) cout << "Couldn't calculate singular values"; else { cout << "Singular values are: " <<NL; cout << SVD.singularValues() << NL; } return 0; }
Matrix Arithmetic: (Matrix.h++, §4.2, pg. 16)
\* Matrix Arithmatic: matarith.cc *\ \* used to test simple matrix arithmatic *\ #include <rw/ftrdgmat.h> #include <rw/rstream.h> main(){ FloatTriDiagMat T(6,6); T.diagonal(-1)=1; T.diagonal(0)=2; T.diagonal(1)=3; FloatVec x(6,1,1); FloatVec Tx = product(T,x); cout << "T = " << T << NL; cout << "x = " << x << NL; cout << "Tx = " << Tx << NL; }
FloatTriDiagMat T(6,6)and not
FloatTriDiagMat T(6,6,0)as is in the manual.
Non-Linear Least Squares Fit: (There is no code in the manual for this example.)
/* RWFIT.CC -- performs leas - squares - fit on data from file or */ /* entered by keyboard. Also does some error analysis on the */ /* calculated values. This version of the program uses the Rogue */ /* Wave matrix libraries. */ #include <fstream.h> #include <rw/dgenfct.h> #include <rw/dgenmat.h> #include <rw/rstream.h> #include <rw/pstream.h> double power(double base, double exp){ if (exp == 0) return 1; else if (base == 0) return 0; else return pow(base,exp); } class LeastSquareFit{ // class to do a polynomial least-squares-fit on public: // supplied data matrices LeastSquareFit(); void inputxy(); // for keyboard input of data void calcFit(); // calculates actual LU fit from existing arrays void errorAnal(); // performs error analysis from calculated values void fileLoad(); // loads data into matrices from given file private: int n; // number of data points DoubleGenMat *y; // pointer to array of y data DoubleGenMat *a; // pointer to x matrix, size depends on order of fit DoubleGenMat *calcValues; // pointer to calculated values int order; // order of polynomial to fit DoubleGenMat *X; // matrix of polynomial coefficients double errorSqd(); // finds sum of error squared double rms(); // finds rms of errors double std(); // finds standard deviation of errors double corrCoeff(); // finds correlation coeff. ("linearness") of data }; LeastSquareFit::LeastSquareFit(){ cout << "What order should the fit be?"; cin >> order; X = new DoubleGenMat(order+1, 1); char answer; cout << "Load from file? (y/n)"; cin >> answer; if ((answer == 'Y') || (answer == 'y')){ fileLoad(); } else inputxy(); } void LeastSquareFit::inputxy(){ cout << "How many sets of data are there? "; cin >> n; a = new DoubleGenMat(n, order+1); y = new DoubleGenMat(n,1); double xtemp, ytemp; for (int i = 0;i < n; i++){ cout << "#" << i+1; cout << "\tX:"; cin >> xtemp; cout << "\tY:"; cin >> ytemp; (*y)(i,0) = ytemp; for (int j = 0; j <= order; j++) (*a)(i,j) = power(xtemp,j); } } void LeastSquareFit::calcFit(){ DoubleGenMat At = transpose(*a); DoubleGenMat G(At.rows(), a->cols()); G = product(At,*a); if(DoubleGenFact(G).fail()){ cout << "Matrix constructed was singular,\n"; cout << "perhaps you need more data points.\n"; exit (1); } DoubleGenMat ginv = inverse(G); *X = product(ginv, product(At,*y)); printf ("\n\nThe coeff's. are, in increasing order:\n"); for (int i = 0; i <= order; i++) printf("x^%i = %.2f\n", i, (*X)(i,0)); printf("\n\n"); } void LeastSquareFit::errorAnal(){ DoubleGenMat yCalc(n,1); for (int i = 0; i < n; i++){ double temp = 0; for (int j = order; j >= 0; j--) temp += (*X)(j,0) * power ((*a)(i,1),j); yCalc(i,0) = temp; } calcValues = &yCalc; printf("Sum of Errors Squared = %.3f\n", errorSqd()); printf("Rms of Error = %.3f\n", rms()); printf("Standard Deviation = %.3f\n", std()); printf("Correlation Coeff. = %f\n", corrCoeff()); } double LeastSquareFit::errorSqd(){ double sum = 0; for (int i = 0; i < n; i++) sum += power ((*calcValues)(i,0) - (*y)(i,0), 2); return sum; } double LeastSquareFit::rms(){ return sqrt( errorSqd() / n ); } double LeastSquareFit::std(){ // find standard dev. of errors; double std = 0; // find mean of errors; double mean = 0; for (int i = 0;i < n; i++) mean += (*calcValues)(i,0) - (*y)(i,0); mean /= n; // now calculate variance; for (i = 0; i < n; i++) std += power ( ((*calcValues)(i,0) - (*y)(i,0)) - mean, 2); return sqrt ( std / n); } double LeastSquareFit::corrCoeff(){ double ymean = 0; double xmean = 0; double calcMean = 0; for (int i = 0; i < n; i++){ ymean += (*y)(i,0); xmean += (*a)(i,1); calcMean += (*calcValues)(i,0); } xmean /= n; ymean /=n; calcMean /=n; double xy = 0, xx = 0, yy = 0; for (i = 0; i < n; i++){ xy += ((*a)(i,1) - xmean) * ((*y)(i,0)-ymean); xx += power ((*a)(i,1) - xmean,2); yy += power ((*y)(i,0) - ymean,2); } return (xy/ sqrt (xx*yy)); } void LeastSquareFit::fileLoad(){ char *filename; filename = new char[100]; cout << "What is the filename? "; cin >> filename; FILE *input; input = fopen(filename, "r"); if (!input) while (!input){ cout << "File \"" << filename << "\" does not exist\n"; cout << "Load from what file? "; cin >> filename; input = fopen(filename, "r"); } delete filename; float tempx, tempy; n = 0; while (1){ (void) fscanf(input, "%f%f", &tempx, &tempy); if feof(input) break; n++; } double xtemp, ytemp; a = new DoubleGenMat (n, order+1); y = new DoubleGenMat (n, 1); rewind (input); for (int i = 0; i < n; i++){ (void)fscanf(input, "%lf%lf", &xtemp, &ytemp); (*y)(i,0) = ytemp; for (int j = 0; j <= order; j++) (*a)(i,j) = power(xtemp,j); } fclose (input); } int main(){ LeastSquareFit x; x.calcFit(); x.errorAnal(); }
Summary: