Vector operations

[Chris Torek]

[long quote]
In article <270 at sdrc.UUCP> scjones at sdrc.UUCP (Larry Jones) writes:
>most users of vector processors are surprised that the following function:
>
>	void addvec(double a[], double b[], double c[], int n) {
>	   int i;
>
>	   for (i = 0; i < n; i++)
>	      a[i] = b[i] + c[i];
>	   }
>
>does not result in a single vector addition (especially since writing the
>"same" function in FORTRAN does).  Actually, "angry" is probably a better
>term that "surprised" in this case.  The reason for this is, of course, that
>in C there is nothing wrong with calling this function like:
>
>	addvec(a+1, a, a+1, sizeof a / sizeof a[0] - 1);
>
>which sums the vector into its last element.  This is invalid in FORTRAN
>which is why the FORTRAN version can vectorize but the C version can't.

I would like to mention that while something like

	CALL ADDVEC(A(2), A(1), A(2), N - 1)

is technically illegal given a SUBROUTINE ADDVEC as defined above,
transliterated to F77, most compilers would never notice this.  The
call would behave on scalar machines and do something odd on vector
machines.  This sort of thing does happen, and I doubt not that (say)
Convex (who have vector-reducing F77 and C compilers) get bug reports
from irate users demanding that the compiler be fixed, when in fact it
is the users' programs that are broken.  One might think `I wrote a
loop, not a vector instruction'---and in a language like C, this claim
has far more weight than in a higher level language such as FP or
Prolog.

Indeed, given _builtin functions, one might note that a reasonably
portable vector application could be written without `noalias'.
Consider the following:

	#include "vector.h"

		...
		vec_addD3(a, b, c, VEC_N_ELEM(c));
		vec_mulD2(c, b, VEC_N_ELEM(c));
		...
	$ cc main.c vecsupport.c

Now look at `vector.h':

	/* get the configuration information */
	#include "vec_config.h"

	#define	VEC_N_ELEM(x) (sizeof(x) / sizeof(*(x)))

	#ifdef VEC_NONE
	/* use simulation library */
	#define VEC_3ADDRESS
	#endif

	#if defined(VEC_3ADDRESS)
	/* all built-in vectors are three-address-code versions */
	#define vec_addF2(a, b, n) vec_addF3(a, b, b, n)
	#define	vec_addD2(a, b, n) vec_addD3(a, b, b, n)
	#define vec_mulF2(a, b, n) vec_mulF3(a, b, b, n)
	#define vec_mulD2(a, b, n) vec_mulD3(a, b, b, n)
	...
	#else
	void vec_addF2(float a[], float b[], size_t n);
	void vec_addD2(double a[], double b[], size_t n);
	...
	#endif
	...
	#ifdef VEC_NONE
	void vec_addF3(float a[], float b[], float c[], size_t n)
	void vec_addD3(double a[], double b[], double c[], size_t n)
	void vec_mulF3(float a[], float b[], float c[], size_t n)
	...
	#else
	/* builtin versions figure type for themselves */
	#define vec_addF3(a, b, c, n) _builtin_vec_add(a, b, c, n)
	#define vec_addD3(a, b, c, n) _builtin_vec_add(a, b, c, n)
	#define vec_mulF3(a, b, c, n) _builtin_vec_mul(a, b, c, n)
	...

The greatest problem with this system is the necessity for applying the
proper type to the vector sum.  The C++ language can do this
automatically; indeed, with the proper class definitions, and a bit of
operator overloading, one can simply write

		c = a + b;		// vector sum

Perhaps, now that everyone is jumping on the C bandwagon, we should aim
the number crunching folk at C++ instead.  Tell them `It is just like
C, except that you can write matrix operations directly.'  Perhaps that
will win them over. :-)
-- 
In-Real-Life: Chris Torek, Univ of MD Comp Sci Dept (+1 301 454 7163)
Domain:	chris at mimsy.umd.edu	Path:	uunet!mimsy!chris