_WHY C++ WILL REPLACE FORTRANB_ by Thomas Keffer [LISTING ONE] // This class holds a reference count and a pointer to the vector data class DataBlock { private: unsigned short refs; // Number of references void* array; // Pointer to raw, untyped data public: DataBlock(unsigned nelem, size_t elemsize) {array = new char[nelem*elemsize]; refs = 1; } ~DataBlock() { delete array; } void addReference() {refs++;} unsigned short removeReference() {return --refs;} void* data() {return array;} }; [LISTING TWO] class DoubleVec { DataBlock* block; // Pointer to the DataBlock double* begin; // Start of data unsigned npts; // Length of the vector int step; // Stride length DoubleVec(const DoubleVec&, int, unsigned, int); // For slices public: DoubleVec(unsigned n, double val); // n elements, initialized to val DoubleVec(const DoubleVec& a); // Copy constructor ~DoubleVec(); DoubleVec slice(int start, int stride, unsigned lgt) const; unsigned length() const {return npts;} int stride() const {return step;} double& operator()(int i) const {return begin[i*step];} DoubleVec& operator=(const DoubleVec& v); }; // Copy constructor: DoubleVec::DoubleVec(const DoubleVec& a) { block = a.block; block->addReference(); npts = a.npts; begin = a.begin; step = a.step; } // Destructor: DoubleVec::~DoubleVec() { if (block->removeReference()==0) delete block; } [LISTING THREE] // Construct a vector that is a slice of an existing vector: DoubleVec::DoubleVec(const DoubleVec& v, int start, unsigned n, int str) { block = v.block; block->addReference(); npts = n; begin = v.begin + start*v.stride; // Add start points stride = str*v.stride; // Accumulate strides } [LISTING FOUR] class DoubleMatrix : private DoubleVec { unsigned ncols; // Number of columns unsigned nrows; // Number of rows public: DoubleMatrix(); DoubleMatrix(unsigned rows, unsigned cols, double initval); DoubleMatrix(const DoubleMatrix& m); // Reference to m will be made unsigned cols() const {return ncols;} unsigned rows() const {return nrows;} DoubleVec col(unsigned j) const // Return col as slice { return DoubleVec::slice(j*nrows, 1, nrows); } DoubleVec row(unsigned i) const; // Return row as slice { return DoubleVec::slice(i, nrows, ncols); } DoubleVec diagonal() const; // Return diagonal as slice double& operator()(int i, int j) const; // Subscripting ... // Other operators }; [LISTING FIVE] class LUDecomp : private DoubleMatrix { int* permute; // Row permutations public: LUDecomp(const DoubleMatrix&); ~LUDecomp() { delete permute; } int isSingular() const; friend double determinant(const LUDecomp&); friend DoubleMatrix inverse(const LUDecomp&); friend DoubleVec solve(const LUDecomp&, const DoubleVec&); }; [LISTING SIX] class FFTServer { unsigned npts; ComplexVec theRoots; void calculateRoots(); public: FFTServer(unsigned n = 0) {setOrder(n);} void setOrder(unsigned n) { npts = n; calculateRoots(); } ComplexVec fourier(const ComplexVec& v){ if(v.length() != npts) setOrder(v.length()); ComplexVec tran; // ... (calculate transform, put it in tran) return tran; } }; Example 1 ComplexVec b; // Declare a complex vector cin >> b; // Read it in FFTServer s; // Allocate an FFT server // Calculate the transform of b: ComplexVec theTransform = s.fourier(b); cout << theTransform; // Print it out Example 2 (a) DoubleVec a; // Null vector - no elements at all, but can be resized DoubleVec b(8); // 8 elements long, uninitialized DoubleVec c(8,1); // 8 elements, initialized to 1.0 DoubleVec d(8,1,2); // 8 elements, initialized to 1, 3, 5, ... (b) b = 2.0; // Set all elements of b to 2 b = c + d; // Set b to the sum of c and d b *= 2; // Multiply each element in b by 2. (c) b[2] = 4.0; // Set the 3'rd element of b to 4.0. c[1] = b[3]; // Set the 2'nd element of c to the 4'th element of b (d) DoubleVec a(10, 0); // 10 element vector for(int i = 0; i<10; i+=2) a[i] = 1; Example 3 (a) a.slice(0, 2) = 1; // Starting with element 0; set every other element to 1 (b) class DoubleVec { ... public: ... operator DoubleSlice(); // For type conversion DoubleSlice slice(int start, int step, unsigned N) const; }; // The "helper class": class DoubleSlice { DoubleVec* theVector; int startElement; unsigned sliceLength; int step; public: ... friend DoubleVec operator+(const DoubleSlice&, const DoubleSlice&); ... }; (c) DoubleVec b(10, 0), c(10, 1); DoubleVec d = b.slice(0, 2) + c.slice(1, 2); (d) DoubleVec g = b + c; // DoubleVec to DoubleSlice type conversion occurs Example 4 (a) DoubleVec real(const ComplexVec&); (b) ComplexVec a(10, 0); // (0,0), (0,0), (0,0), ... real(a) = 1.0; // (1,0), (1,0), (1,0), ... Example 5 (a) DoubleMatrix a(10, 10, 0); // 10 by 10 matrix, initialized to zero a.row(3) = 1; // Set row 3 to 1 a.col(2) = a.col(4); // Copy column 4 to column 2 (b) DoubleMatrix I(10, 10, 0); // 10x10 initialized to 0 I.diagonal() = 1; // Create an identity matrix Example 6 Example 6 (a) DoubleMatrix a(10, 10); // ... (initialize a somehow) // Construct the LU decomposition of a: LUDecomp aLU(a); // Now use it: double det = determinant(aLU); DoubleMatrix aInverse = inverse(aLU); (b) // 5 different sets of linear equations to be solved: DoubleVec b[5], x[5]; // ... (set up the 5 vectors b and the 5 vectors x, each // with 10 elements as per the matrix a above) for(int i = 0; i < 5; i++) x[i] = solve(aLU, b[i]); (c) [SET THIS IN HELV, NOT PRESTIGE] a xi = bi; i = 0, ..., 4. (d) DoubleMatrix a(10, 10); // ... (initialize a) // Calculate the inverse directly from a. // A DoubleMatrix to LUDecomp type conversion takes place automatically: DoubleMatrix aInverse = inverse(a); (e) DoubleMatrix aInverse = inverse(LUDecomp(a)); Example 7 ComplexVec timeVector(30); FFTServer aServer; // Allocate a server // Will automagically reconfigure for a vector of length 30: ComplexVec spectrum = aServer.fourier(timeVector); Example 8 class Force; // 1 class String { // 2 public: String(double length, double tension, double density); // 3 void setPoints(int nx); // 4 void timeStep(double dt, const Force& force); // 5 DoubleVec displacement() const; // 6 private: DoubleVec u; // 7 double cSquared; double length; }; Example 9 (a) class Force { public: virtual DoubleVec value() = 0; // 1 }; (b) class WindForceString : public Force { public: WindForceString( String& string, double dragCoeff ); void setVelocity(double windspeed); virtual DoubleVec value(); // 1 private: String& myString; // The string we are tracking double wind; // Present wind speed. double drag; }; (c) DoubleVec WindForceString::value() { // Get present string displacement: DoubleVec d = string.displacement(); return -drag*wind*wind*d; // Some (bogus) calculation } Figure 1: C++ Fortran #include parameter (M=8000) #include double precision a(M), b(M) #include double precision c(M) main(int argc, char* argv[]) { read(5,*) N, Niter unsigned N = atoi(argv[1]); write(6,1000)Niter, N unsigned Niter = atoi(argv[2]); 1000 format(i8,' iterations of ', * i8, ' points each.') cout << Niter << "iterations of " << N << " points each.\n"; do 50 j=1,N a(j) = 3 DoubleVec a(N, 3), b(N, 5); 50 b(j) = 5 while(Niter--) do 100 i=1,Niter DoubleVec c = a * b; do 100 j=1,N } 100 c(j) = a(j) * b(j) stop end Figure 2 (a) #include #include class DTestMatrix : public DoubleGenMat { public: DTestMatrix(unsigned order); }; class DTestRHS : public DoubleVec { public: DTestRHS(const DTestMatrix&); }; double second(); double epslon(double); const unsigned N = 90; const unsigned long ops = 2.0*N*N*N/3.0 + 2.0*N*N; void main(){ DTestMatrix a(N); double norma = maxVal(abs(a)); DTestRHS b(a); double t1 = second(); // Construct the LU Factorization: DoubleGenFact fact(a); t1 = second() - t1; double t2 = second(); DoubleVec x = solve(fact, b); t2 = second() - t2; double total = t1+t2; double mflops = ops / (1.0e6*total); DoubleVec tol = a.product(x) - b; double resid = maxVal(abs(tol)); double normx = maxVal(abs(x)); double eps = epslon(1.0); double residn= resid / (N*norma*normx*eps); cout << "Normalized residual = " << residn << NL; cout << "Residual = " << resid << NL; cout << "Machine precision = " << eps << NL; cout << "Factorization time = " << t1 << NL; cout << "Solution time = " << t2 << NL; cout << "Total time = " << total << NL; cout << "MFLOPS = " << mflops << NL; } DTestMatrix::DTestMatrix(unsigned n) : DoubleGenMat(n,n) { long init = 1325; for(int j=0; j