Scientific Computing: C++ Versus Fortran
by Todd Veldhuizen 


Listing One
Matrix<float> A(3,3);
Vector<float> b(3), c(3);
// ...
c = product(A,b);

Listing Two
TinyMatrix<float,3,3> A;
TinyVector<float,3> b, c;
c = product(A,b);

Listing Three
template<class T, int N>
class TinyVector {
private:
    float data[N];
};

Listing Four
c[0] = A[0][0]*b[0] + A[0][1]*b[1] + A[0][2]*b[2];
c[1] = A[1][0]*b[0] + A[1][1]*b[1] + A[1][2]*b[2];
c[2] = A[2][0]*b[0] + A[2][1]*b[1] + A[2][2]*b[2];


Listing Five
typedef TinyVector<double,3> ray;
inline void reflect(ray& y, const ray& x,
  const ray& n)
{
  y = x - 2 * dot(x,n) * n;
}

Listing Six
double _t1 = x[0]*n[0] + x[1]*n[1]
    + x[2]*n[2];
double _t2 = _t1 + _t1;
y[0] = x[0] - _t2 * n[0];
y[1] = x[1] - _t2 * n[1];
y[2] = x[2] - _t2 * n[2];

Listing Seven
void test(ray& y)
{
  ray x, n;
  x = 1.00,  0.40,  -1.00;
  n = 0.31,  0.20,   0.93;
  reflect(y, x, n);
}

Listing Eight
void test(ray& y)
{
  y[0] = 1.3348;
  y[1] = 0.6160;
  y[2] = 0.0044;
}


Listing Nine
Vector<double> x(N), y(N);
double a;
   ... 
y += a * x;

Listing Ten
DO k=2, N-1
 DO j=2, N-1
  DO i=2, N-1
   A(i,j,k) = (B(i,j,k)   + B(i+1,j,k) + B(i-1,j,k)
   .         + B(i,j+1,k) + B(i,j-1,k) + B(i,j,k+1)
   .         + B(i,j,k-1)) / 7.0;
  END DO
 END DO
END DO

Listing Eleven
const int N = 64;
Array<double,3> A(N,N,N), B(N,N,N);
 ...
Range I(1,N-2), J(1,N-2), K(1,N-2);
A(I,J,K) = (B(I,J,K) + B(I+1,J,K) + B(I-1,J,K)
  + B(I,J-1,K) + B(I,J+1,K) + B(I,J,K-1)
  + B(I,J,K+1)) / 7.0;

Listing Twelve
typedef TinyMatrix<complex<double>,3,3> gaugeElement;
typedef TinyMatrix<complex<double>,3,2> spinor;
void qcdCalc(Vector<spinor>& res, Vector<gaugeElement>& M,
  Vector<spinor>& src)
{
  for (int i=0; i < res.length(); ++i)
    res[i] = product(M[i], src[i]);
}

Listing Thirteen
subroutine qcdf(M, res, src, V)
       integer V, iters
       double complex M(V,3,3), res(V,3,2), src(V,3,2)

       DO site=1,V
         DO spin=1,2
           DO col=1,3
              res(site,col,spin) = M(site,col,1) * src(site,1,spin)
     .          + M(site,col,2) * src(site,2,spin)
     .          + M(site,col,3) * src(site,3,spin)
           ENDDO
         ENDDO
       ENDDO

       return
       end

Listing Fourteen
class latticeUnit {
 ...
protected:
  spinor one;
  gaugeElement ge;
  spinor two;
};

3