Xref: utzoo comp.lang.c++:4133 gnu.g++:303
Path: utzoo!attcan!utgpu!jarvis.csri.toronto.edu!mailrus!csd4.milw.wisc.edu!bionet!apple!brutus.cs.uiuc.edu!flute.cs.uiuc.edu!grunwald
From: grunwald@flute.cs.uiuc.edu (Dirk Grunwald)
Newsgroups: comp.lang.c++,gnu.g++
Subject: Re: C++ vrs Fortran for linear algebra.
Message-ID: <GRUNWALD.89Jul26104549@flute.cs.uiuc.edu>
Date: 26 Jul 89 15:47:57 GMT
References: <3967@portia.Stanford.EDU>
Sender: news@brutus.cs.uiuc.edu
Reply-To: grunwald@flute.cs.uiuc.edu
Organization: University of Illinois, Urbana-Champaign
Lines: 96
In-reply-to: jon@hanauma's message of 26 Jul 89 05:05:42 GMT


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)