[gnu.g++] C++ vrs Fortran for linear algebra.

grunwald@flute.cs.uiuc.edu (Dirk Grunwald) (07/26/89)

I modified your example slightly and tried it on a Sun-3/60 w/68881
running SunOS 4.0.3, using g++ 1.35.1+

I got the following results:

Floating Point		Compiled by:				Time
--------------------------------------------------------------------------
real			f77 -O					84.7
real			f77 -O4					85.4
double precision	f77 -O4					88.1

float			g++ -O -fstrength-reduce		85.5
double			g++ -O -fstrength-reduce		95.9
double, global vars	g++ -O -fstrength-reduce		95.3
double, no dope vector	g++ -O -fstrength-reduce		97.3

Conclusions:
	+ I can't wait to get a DEC-3100

	+ G++ isn't so bad for numeric codes

	+ G++/Gcc needs better optimization for `doubles'

	+ G++ needs to improve optimization for induction variables
	  that are class member functions [posted to gnu.g++.bug]

Changes to the code & why:

	+ g++ doesn't grok inlines if they're in the class defn, so I
	  pulled it out

	+ G++ has dynamic local vectors, so you don't need the mallocs,
	  but you do need a macro [curse the language].

	  I think this saves MUCHO time. You could pull a similar trick
	  using ``alloca'' on a vanilla-flavoured C++ backend.
	  You could also make the variables be global, since FORTRAN
	  does that (no saving in of itself)

	+ used inline math library (<math.h> includes <math-68881.h>)

//--------------- C++ header to handle arrays like fortran.
// by dave nichols, hack by d. grunwald
class Matrix {
private:
    int *index;
    double *data;
public:
    Matrix( int n1, int n2, int *indx, double *dat);
    double& operator()( int i, int j );
};

static inline double&
Matrix::operator()( int i, int j ) {
    return ( data[i + index[j] ] );
}

static inline
Matrix::Matrix( int n1, int n2, int *indx, double *dat) {
    index = indx;
    data = dat;
    for( int i =0; i<n2; i++ ) {
	index[i] = i*n1;
    }
}

#include <stdio.h>
#include <math.h>
#include <generic.h>

#define MATRIX(name,n1,n2)\
    int name2(name,index)[n1];\
	double name2(name,data)[n1*n2]; \
	    Matrix name(n1, n2, name2(name,index), name2(name,data));

main () {
    int	i, j, k;
    double	total=0.;

    const int n=200;
    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",
	   (double) total );
}
--
Dirk Grunwald -- Univ. of Illinois 		  (grunwald@flute.cs.uiuc.edu)