shootout

["E. Robert Tisdale" <edwin at netwood dot net>]

Can numerical application programmers use
the GNU Scientific Library (GSL) to write clear, concise code
that other numerical programmers can read, understand and maintain?

The best way to test this is to compare implementations using the GSL
with implementations using other numerical libraries such as
the Vector, Signal and Image Processing Library (VSIPL) for example.
Please find attached, a copy of the shootout program
which was designed to test the efficiency of implementing applications
using various numerical class libraries.
The problem is well documented in

	Comparison of APIs on VSIP ADBF Example

which can be obtained from the VSIPL website

	http://www.vsipl.org/

Just click on the EXAMPLES radio button on the left hand side of the page
then click on Comparison of APIs on VSIPL ADBF Example (August 1997).

Both Microsoft Word 6.0 and Adobe portable document files are available.

Numerical application programs implemented in ANSI C
using numerical class libraries are about three times larger
than the same numerical application program implemented in C++
using numerical class libraries.
An implementation of the shootout program implemented in C++ using
The C++ Scalar, Vector, Matrix and Tensor class library

	http://www.netwood.net/~edwin/svmt/

is attached below as well.

---------------------------------------------------------------------------

/* shootout.c Version 0.000 February 20, 1997
//
// This file contains VSIP Memory Hierarchy Test Application #1.
//
// Copyright February 1997, Hughes Research Laboratories, Inc.
// This is an unpublished work. All rights reserved.
//
// Authors: David Schwartz and E. Robert Tisdale
//
// Code is an adaptive beamformer,
// based on the MountainTop code by Yaron Seliktar, GA Tech.
// Comments may contain the equivalent MATLAB code.
//
// Revision History:
// Feb. 28, 1997 Added copyright notice.
// Mar. 20, 1997 Modified VSIP class and function names.
// July 10, 1997 RJudd modified to reflect current function names.
// July 11, 1997 MARichards modified to create complex weight vectors
*/

#include <vsip_cmview_d.h>

void adaptive(
  const vsip_cmview_d *X,	/* x[MN][L]		*/
  const double d_n,		/* Dn			*/
  const double f_n,		/* Fn			*/
  const vsip_vview_d *w_s,	/* spatial_window[N]	*/
  const vsip_vview_d *w_t,	/* temporal_window[M]	*/
  const vsip_cvview_d *a) {	/* weight[MN]		*/

  vsip_length N = vsip_vgetlength_d(w_s);
  vsip_length M = vsip_vgetlength_d(w_t);
  vsip_length MN = vsip_cmgetlength2_d(X);

  int singular;			/* for testing singularity of LU decomp	*/

  if (MN == M*N) {			/* conformant inputs		*/
    if (MN == vsip_cvgetlength_d(a)) {	/* conformant output		*/
      vsip_hint		 h = 0;		/* no optimization hints	*/
      vsip_cvview_d	*v = vsip_cvcreate_d(MN, h);
					/* steer[MN]			*/
      vsip_cmview_d	*R = vsip_cmcreate_d(MN, MN, h);
					/* R[MN][MN]			*/
      vsip_length	 L = vsip_cmlength1_d(X);
      vsip_couter_d(X, X, R);		/* R = x*x'			*/
      vsip_rcsmmul_d(1.0/L,R, R);	/* R = R/L			*/
      {					/* block statement		*/
	vsip_cvview_d	*d_R = vsip_cmdiagview_d(R);
					/* d_R = diag(R)		*/
	/* The diagonal elements of complex matrix R are real!		*/
	vsip_vview_d	*rd_R = vsip_cvrealview_d(d_R);
	vsip_vview_d	*mag_rd_R = vsip_vcreate_d(MN, h);
	vsip_vmag_d(rd_R, mag_rd_R);
	{				/* block statement		*/
	  double	 sigma = vsip_vmax_d(mag_rd_R)*1.0e-12;
	  /* noisef = 1e-12*max(mag(diag(R)))				*/
	  vsip_blockdestroy_d(vsip_vdestroy_d(mag_rd_R));
	  vsip_svadd_d(sigma, rd_R, rd_R);
	  /* R = R + noisef*eye(size(R))				*/
	  }
	vsip_vdestroy_d(rd_R);
	vsip_cvdestroy_d(d_R);
	}
      {					/* block statement		*/
	vsip_cvview_d	*v_t = vsip_cvcreate_d(M, h);
					/* temporal[M]			*/
	vsip_cvview_d	*v_s = vsip_cvcreate_d(N, h);
					/* spatial[N]			*/

	vsip_vview_d	*ramp = vsip_vcreate_d(N,h);
	vsip_vview_d	*routput = vsip_vcreate_d(N,h);
	vsip_vview_d	*ioutput = vsip_v vcreate_d(N,h);
	vsip_vramp_d(0.0, (2.0 * PI * d_n) , ramp);
	vsip_vcos_d(ramp, routput);
	vsip_vsin_d(ramp, ioutput);
	vsip_vcmplx_d( routput, ioutput, v_s );
	vsip_rcvmul_d( w_s, v_s, v_s );
	/* spatial = spatial_window .* exp(j*2*pi*Dn*[0:N-1]')		*/

	vsip_vramp_d(0.0, (2.0 * PI * f_n) , ramp);
	vsip_vcos_d(ramp, routput);
	vsip_vsin_d(ramp, ioutput);
	vsip_vcmplx_d( routput, ioutput, v_s );
	vsip_rcvmul_d( w_t, v_t, v_t );
	/* temporal = temporal_window .* exp(j*2*pi*Fn*[0:M-1]')	*/

	{				/* block statement		*/
	  vsip_vview_d	*t = vsip_vcreate_d(MN, h);
	  vsip_ckronprod_d(v_t, v_s, v);
	  /* steer = kron(temporal, spatial)				*/
	  vsip_blockdestroy_d(vsip_vdestroy_d(routput));
	  vsip_blockdestroy_d(vsip_vdestroy_d(ioutput));
	  vsip_blockdestroy_d(vsip_vdestroy_d(ramp));
	  vsip_blockdestroy_d(vsip_cvdestroy_d(v_s));
	  vsip_blockdestroy_d(vsip_cvdestroy_d(v_t));
	  {				/* block statement		*/
	    vsip_cvmag_d(v, t);
	    vsip_scalar_d	mean = vsip_vmeanval_d(t);
					/* avg = mean(mag(steer))	*/
	    vsip_blockdestroy_d(vsip_vdestroy_d(t));
	    vsip_rcsvmul_d(1.0/mean, v, v);
					/* steer = steer/avg		*/
	    }
	  }
	}
      {					/* block statement		*/
	vsip_vview_i	*pivot = vsip_vcreate_i(MN,h);
	int error = vsip_clud_d( R, singular, pivot);
	vsip_cvsolvelu( R, pivot, v, a);
					/* weight = R\steer		*/
	vsip_blockdestroy_d(vsip_cmdestroy_d(R));
	vsip_blockdestroy_d(vsip_cvdestroy_d(v));
	vsip_blockdestroy_i(vsip_vdestory_i(pivot));
	if(error) { /* Handle the error. */
	  }
	}
      }
    else { /* Handle the error. */
      }
    }
  return;
  }

---------------------------------------------------------------------------

// The C++ Scalar, Vector, Matrix and Tensor classes.
// Copyright (C) 1998 E. Robert Tisdale
// 
// This file is part of The C++ Scalar, Vector, Matrix and Tensor classes.
// This library is free software which you can redistribute and/or modify
// under the terms of the GNU Library General Public License
// as published by the Free Software Foundation;
// either version 2, or (at your option) any later version.
// 
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty
// of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU General Public License for more details.
// 
// You should have received a copy of the GNU Library General Public License
// along with this library.  If not, write to the Free Software Foundation,
// Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
// 
// Written by E. Robert Tisdale

#include <svmt.h>

doubleComplexVector
        XCS(const doubleSubVector& v, const double s, const double t = 0) {
  doubleComplexVector   result(v.length());
  result.imag().ramp() *= s;            // s*[0: v.length() - 1]
  result.real()  = cos(result.imag());
  result.imag()  = sin(result.imag());
  return v*result; }
inline
double  mean(const doubleVector& v) { return v.sum()/v.length(); }

inline
doubleComplexVector                     // solve b = xA^T
        solve(const doubleComplexVector& b, doubleComplexMatrix& A) {
  return b.pl(A.lud(), A).du(A); }

// The adaptive beamformer is part of the MountainTop code
// written in MATLAB by Yaron Seliktar at Georgia Tech. 
// The original MATLAB code appears in C++ comments
// immediately after the equivalent SVMT code.

doubleComplexVector adaptive(
  const doubleComplexSubMatrix& X,      // x[MN][L]
  const double& d,                      // Dn
  const double& f,                      // Fn
  const doubleSubVector& sw,            //  spatial_window[N]
  const doubleSubVector& tw) {          // temporal_window[M]
  if (X.extent2() != tw.length()*sw.length()) {
                                        // Handle the error.
    }
  doubleComplexVector     t = XCS(tw, 2.0*M_PI*f);
                // temporal = temporal_window .* exp(j*2*pi*Fn*[0: M-1]')
  doubleComplexVector     s = XCS(sw, 2.0*M_PI*d);
                //  spatial =  spatial_window .* exp(j*2*pi*Dn*[0: N-1]')
  doubleComplexVector     v = t.kron(s);// steer = kron(temporal, spatial)
  doubleComplexMatrix     R = X.dot(conj(X))/X.extent1();
                                        // R = x*x'/L
  R.diag() += 1.0e-12*max(abs(R.diag().real()));
                                        // R = R + noisef*eye(size(R))
  return solve(v/mean(abs(v)), R);      // weight = R\(steer/mean(abs(steer)))
  }