Basic operators and adjoints

Spray and sum : scatter and gather

Perhaps the most common operation is the summing of many values to get one value. Its adjoint operation takes a single input value and throws it out to a space of many values. The summation operator is a row vector of ones. Its adjoint is a column vector of ones. In one dimension, this operator is almost too easy for us to bother showing a routine. But it is more interesting in three dimensions, in which we could be summing or spraying on any of three subscripts, or even summing on some and spraying on others. In module spraysum, both input and output are taken to be three-dimensional arrays. Externally, however, either could be a scalar, vector, plane, or cube. For example, the internal array xx(n1,1,n3) could be externally the matrix map(n1,n3). When module spraysum is given the input dimensions and output dimensions stated in the following, the operations stated alongside are implied.

	(n1,n2,n3)	(1,1,1)	Sum a cube into a value.
	(1,1,1)	(n1,n2,n3)	Spray a value into a cube.
	(n1,1,1)	(n1,n2,1)	Spray a column into a matrix.
	(1,n2,1)	(n1,n2,1)	Spray a row into a matrix.
	(n1,n2,1)	(n1,n2,n3)	Spray a plane into a cube.
	(n1,n2,1)	(n1,1,1)	Sum rows of a matrix into a column.
	(n1,n2,1)	(1,n2,1)	Sum columns of a matrix into a row.
	(n1,n2,n3)	(n1,n2,n3)	Copy and add the whole cube.

If an axis is not of unit length on either input or output, then both lengths must be the same; otherwise, there is an error. Normally, after (possibly) erasing the output, we simply loop over all points on each axis, adding the input to the output. Either a copy or an add is done, depending on the add parameter. It is either a spray, a sum, or a copy, according to the specified axis lengths.

user/gee/spraysum.c

    for (i3=0; i3 < SF_MAX(n3,m3); i3++) {  
	x = SF_MIN(i3,n3-1); 
	y = SF_MIN(i3,m3-1);
	for (i2=0; i2 < SF_MAX(n2,m2); i2++) {   
	    x = x*n2 + SF_MIN(i2,n2-1);   
	    y = y*m2 + SF_MIN(i2,m2-1); 
	    for (i1=0; i1 < SF_MAX(n1,m1); i1++) {   
		x = x*n1 + SF_MIN(i1,n1-1);   
		y = y*m1 + SF_MIN(i1,m1-1); 

		if( adj)  xx[x]  +=  yy[y];
		else      yy[y]  +=  xx[x];
	    }
	}
    }

Basic operators and adjoints