Chris@29: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
Chris@29: 
Chris@29: #include "bqvec/VectorOpsComplex.h"
Chris@29: 
Chris@29: #include <iostream>
Chris@29: #include <cstdlib>
Chris@29: 
Chris@29: #include <time.h>
Chris@29: 
Chris@29: using namespace std;
Chris@29: 
Chris@29: namespace breakfastquay {
Chris@29: 
Chris@29: namespace Test {
Chris@29: 
Chris@29: #ifdef _WIN32
Chris@29: #define drand48() (-1+2*((float)rand())/RAND_MAX)
Chris@29: #endif
Chris@29: 
Chris@29: bool
Chris@29: testMultiply()
Chris@29: {
Chris@29:     cerr << "testVectorOps: testing v_multiply complex" << endl;
Chris@29: 
Chris@29:     const int N = 1024;
Chris@29:     bq_complex_t target[N];
Chris@29:     bq_complex_t src1[N];
Chris@29:     bq_complex_t src2[N];
Chris@29: 
Chris@29:     for (int i = 0; i < N; ++i) {
Chris@29: 	src1[i].re = drand48();
Chris@29: 	src1[i].im = drand48();
Chris@29: 	src2[i].re = drand48();
Chris@29: 	src2[i].im = drand48();
Chris@29:     }
Chris@29: 
Chris@29:     double mean, first, last, total = 0;
Chris@29:     for (int i = 0; i < N; ++i) {
Chris@29:         bq_complex_t result;
Chris@29:         c_multiply(result, src1[i], src2[i]);
Chris@29: 	if (i == 0) first = result.re;
Chris@29: 	if (i == N-1) last = result.im;
Chris@29: 	total += result.re;
Chris@29: 	total += result.im;
Chris@29:     }
Chris@29:     mean = total / (N*2);
Chris@29:     cerr << "Naive method: mean = " << mean << ", first = " << first
Chris@29: 	 << ", last = " << last << endl;
Chris@29: 
Chris@29:     v_multiply(target, src1, src2, N);
Chris@29:     total = 0;
Chris@29: 
Chris@29:     for (int i = 0; i < N; ++i) {
Chris@29: 	if (i == 0) first = target[i].re;
Chris@29: 	if (i == N-1) last = target[i].im;
Chris@29: 	total += target[i].re;
Chris@29: 	total += target[i].im;
Chris@29:     }
Chris@29:     mean = total / (N*2);
Chris@29:     cerr << "v_multiply: mean = " << mean << ", first = " << first
Chris@29: 	 << ", last = " << last << endl;
Chris@29: 
Chris@29:     int iterations = 50000;
Chris@29:     cerr << "Iterations: " << iterations << endl;
Chris@29: 	
Chris@29:     cerr << "CLOCKS_PER_SEC = " << CLOCKS_PER_SEC << endl;
Chris@29:     float divisor = float(CLOCKS_PER_SEC) / 1000.f;
Chris@29: 
Chris@29:     clock_t start = clock();
Chris@29: 
Chris@29:     for (int j = 0; j < iterations; ++j) {
Chris@29: 	for (int i = 0; i < N; ++i) {
Chris@29:             c_multiply(target[i], src1[i], src2[i]);
Chris@29: 	}
Chris@29:     }
Chris@29:     
Chris@29:     clock_t end = clock();
Chris@29: 
Chris@29:     cerr << "Time for naive method: " << float(end - start)/divisor << endl;
Chris@29: 
Chris@29:     start = clock();
Chris@29: 
Chris@29:     for (int j = 0; j < iterations; ++j) {
Chris@29:         v_multiply(target, src1, src2, N);
Chris@29:     }
Chris@29:     
Chris@29:     end = clock();
Chris@29: 
Chris@29:     cerr << "Time for v_multiply: " << float(end - start)/divisor << endl;
Chris@29: 
Chris@29:     return true;
Chris@29: }
Chris@29: 
Chris@29: bool
Chris@29: testPolarToCart()
Chris@29: {
Chris@29:     cerr << "testVectorOps: testing v_polar_to_cartesian" << endl;
Chris@29: 
Chris@29:     const int N = 1024;
Chris@29:     bq_complex_t target[N];
Chris@29:     double mag[N];
Chris@29:     double phase[N];
Chris@29: 
Chris@29:     for (int i = 0; i < N; ++i) {
Chris@29: 	mag[i] = drand48();
Chris@29: 	phase[i] = (drand48() * M_PI * 2) - M_PI;
Chris@29:     }
Chris@29: 
Chris@29:     double mean, first, last, total = 0;
Chris@29:     for (int i = 0; i < N; ++i) {
Chris@29: 	double real = mag[i] * cos(phase[i]);
Chris@29: 	double imag = mag[i] * sin(phase[i]);
Chris@29: 	if (i == 0) first = real;
Chris@29: 	if (i == N-1) last = imag;
Chris@29: 	total += real;
Chris@29: 	total += imag;
Chris@29:     }
Chris@29:     mean = total / (N*2);
Chris@29:     cerr << "Naive method: mean = " << mean << ", first = " << first
Chris@29: 	 << ", last = " << last << endl;
Chris@29: 
Chris@29:     v_polar_to_cartesian(target, mag, phase, N);
Chris@29: 
Chris@29:     total = 0;
Chris@29: 
Chris@29:     for (int i = 0; i < N; ++i) {
Chris@29: 	if (i == 0) first = target[i].re;
Chris@29: 	if (i == N-1) last = target[i].im;
Chris@29: 	total += target[i].re;
Chris@29: 	total += target[i].im;
Chris@29:     }
Chris@29:     mean = total / (N*2);
Chris@29:     cerr << "v_polar_to_cartesian: mean = " << mean << ", first = " << first
Chris@29: 	 << ", last = " << last << endl;
Chris@29: 
Chris@29:     int iterations = 10000;
Chris@29:     cerr << "Iterations: " << iterations << endl;
Chris@29: 	
Chris@29:     cerr << "CLOCKS_PER_SEC = " << CLOCKS_PER_SEC << endl;
Chris@29:     float divisor = float(CLOCKS_PER_SEC) / 1000.f;
Chris@29: 
Chris@29:     clock_t start = clock();
Chris@29: 
Chris@29:     for (int j = 0; j < iterations; ++j) {
Chris@29: 	for (int i = 0; i < N; ++i) {
Chris@29: 	    target[i].re = mag[i] * cos(phase[i]);
Chris@29: 	    target[i].im = mag[i] * sin(phase[i]);
Chris@29: 	}
Chris@29:     }
Chris@29:     
Chris@29:     clock_t end = clock();
Chris@29: 
Chris@29:     cerr << "Time for naive method: " << float(end - start)/divisor << endl;
Chris@29: 
Chris@29:     start = clock();
Chris@29: 
Chris@29:     for (int j = 0; j < iterations; ++j) {
Chris@29: 	v_polar_to_cartesian(target, mag, phase, N);
Chris@29:     }
Chris@29:     
Chris@29:     end = clock();
Chris@29: 
Chris@29:     cerr << "Time for v_polar_to_cartesian: " << float(end - start)/divisor << endl;
Chris@29: 
Chris@29:     return true;
Chris@29: }
Chris@29: 
Chris@29: bool
Chris@29: testPolarToCartInterleaved()
Chris@29: {
Chris@29:     cerr << "testVectorOps: testing v_polar_interleaved_to_cartesian" << endl;
Chris@29: 
Chris@29:     const int N = 1024;
Chris@29:     bq_complex_t target[N];
Chris@29:     double source[N*2];
Chris@29: 
Chris@29:     for (int i = 0; i < N; ++i) {
Chris@29: 	source[i*2] = drand48();
Chris@29: 	source[i*2+1] = (drand48() * M_PI * 2) - M_PI;
Chris@29:     }
Chris@29: 
Chris@29:     double mean, first, last, total = 0;
Chris@29:     for (int i = 0; i < N; ++i) {
Chris@29: 	double real = source[i*2] * cos(source[i*2+1]);
Chris@29: 	double imag = source[i*2] * sin(source[i*2+1]);
Chris@29: 	if (i == 0) first = real;
Chris@29: 	if (i == N-1) last = imag;
Chris@29: 	total += real;
Chris@29: 	total += imag;
Chris@29:     }
Chris@29:     mean = total / (N*2);
Chris@29:     cerr << "Naive method: mean = " << mean << ", first = " << first
Chris@29: 	 << ", last = " << last << endl;
Chris@29: 
Chris@29:     v_polar_interleaved_to_cartesian(target, source, N);
Chris@29: 
Chris@29:     total = 0;
Chris@29: 
Chris@29:     for (int i = 0; i < N; ++i) {
Chris@29: 	if (i == 0) first = target[i].re;
Chris@29: 	if (i == N-1) last = target[i].im;
Chris@29: 	total += target[i].re;
Chris@29: 	total += target[i].im;
Chris@29:     }
Chris@29:     mean = total / (N*2);
Chris@29:     cerr << "v_polar_interleaved_to_cartesian: mean = " << mean << ", first = " << first
Chris@29: 	 << ", last = " << last << endl;
Chris@29: 
Chris@29:     int iterations = 10000;
Chris@29:     cerr << "Iterations: " << iterations << endl;
Chris@29: 	
Chris@29:     cerr << "CLOCKS_PER_SEC = " << CLOCKS_PER_SEC << endl;
Chris@29:     float divisor = float(CLOCKS_PER_SEC) / 1000.f;
Chris@29: 
Chris@29:     clock_t start = clock();
Chris@29: 
Chris@29:     for (int j = 0; j < iterations; ++j) {
Chris@29: 	for (int i = 0; i < N; ++i) {
Chris@29: 	    target[i].re = source[i*2] * cos(source[i*2+1]);
Chris@29: 	    target[i].im = source[i*2] * sin(source[i*2+1]);
Chris@29: 	}
Chris@29:     }
Chris@29:     
Chris@29:     clock_t end = clock();
Chris@29: 
Chris@29:     cerr << "Time for naive method: " << float(end - start)/divisor << endl;
Chris@29: 
Chris@29:     start = clock();
Chris@29: 
Chris@29:     for (int j = 0; j < iterations; ++j) {
Chris@29: 	v_polar_interleaved_to_cartesian(target, source, N);
Chris@29:     }
Chris@29:     
Chris@29:     end = clock();
Chris@29: 
Chris@29:     cerr << "Time for v_polar_interleaved_to_cartesian: " << float(end - start)/divisor << endl;
Chris@29: 
Chris@29:     return true;
Chris@29: }
Chris@29: 
Chris@29: bool
Chris@29: testVectorOps()
Chris@29: {
Chris@29:     if (!testMultiply()) return false;
Chris@29:     if (!testPolarToCart()) return false;
Chris@29:     if (!testPolarToCartInterleaved()) return false;
Chris@29:     
Chris@29:     return true;
Chris@29: }
Chris@29: 
Chris@29: }
Chris@29: 
Chris@29: }
Chris@29: