Path: utzoo!utgpu!jarvis.csri.toronto.edu!rutgers!apple!sun-barr!decwrl!shelby!portia!hanauma!jon
From: jon@hanauma (Jon Claerbout)
Newsgroups: comp.lang.fortran
Subject: C++ vrs Fortran for linear algebra.
Message-ID: <3969@portia.Stanford.EDU>
Date: 26 Jul 89 05:08:43 GMT
Sender: USENET News System <news@Portia.stanford.edu>
Reply-To: jon@hanauma (Jon Claerbout)
Organization: Stanford University, Dept. of Geophysics
Lines: 70

ABSTRACT:	C++ can be customized to look like Fortran in its
		array subscription, but it has an advantage
		over Fortran that variable dimensioned arrays
		can be dynamically allocated.
		A benchmark shows C++ about half as fast on a DEC 3100.

ANALYSIS:	fortran			C
		real x(n,m)		x = alloc( n*m* sizeof(float))
		x(i,j) = etc		x[ i + index[j] ] = etc

COMMENT:	This is a great example for C++ textbook writers.

COMPILER:	AT&T C++ version 1.2.1

DEMONSTRATION BELOW -------------------------------------------------
/*	Benchmark Fortran equivalent (actually ratfor)
	real total, x(200,200)
	integer i,j,n
	n = 200
	total = 0.
	do i=1,n
		do j=1,n
			x(i,j) = sqrt(1.+i+j)
	do i=1,n
		do j=1,n
			do k=1,n
				total = total + x(i,k) * x(k,j)
	write(6,10) total
10      format('    1576677724. if double precision '/f15.0,' here in single')
	call exit(0);	end
*/
//------------------------------------ C++ user program equivalent to Fortran
#include <CC/stdio.h>
main () {
	double sqrt(float);
	int	i, j, k, n=200;
	float	total=0.;
	Matrix	x(n,n);
	for( i=0; i<n; i++ )
		for( j=0; j<n; j++ )
			x(i,j) = sqrt(1.+(i+1)+(j+1));
	for( i=0; i<n; i++ )
		for( j=0; j<n; j++ )
			for( k=0; k<n; k++ )
				total = total + x(i,k) * x(k,j);
	printf("total=%f\n    1576677724. if double precision\n",total );
}
//   7.0sec	1575556480. 	DEC 3100  (f77 -O3)
//  13.1sec	1575559552	DEC 3100  (C++ -O2)
//
//--------------- C++ header to handle arrays like fortran.
class Matrix {				// by dave nichols
	private:
		int *index;
		float *data;
	public:
		Matrix( int, int);
		inline float& operator()( int i, int j ) {
			return ( data[i + index[j] ] );
		}
		~Matrix();
	};
Matrix::Matrix( int n1, int n2) {
	index = new int[n2];
	for( int i =0; i<n2; i++ ) {
		index[i] = i*n1;
	}
	data = new float[n1*n2];
}
Matrix::~Matrix() { delete index; delete data; }