librobotcontrol-sys 0.4.0

/**
 * @file rc_benchmark_algebra.c
 * @example    rc_benchmark_algebra
 *
 * @brief      benchmarks the linear algebra functions to test floating point
 *             speed
 *
 *             This example prints the time to execute all functions and reports
 *             the speed of basic matrix multiplication in MFLOPS. It can be
 *             used as a test to see if the compiled math library is using the
 *             CPU hardware vectorized floating point units.
 *
 *
 * @author     James Strawson
 * @date       1/29/2018
 */

#include <stdio.h>
#include <stdint.h>
#include <getopt.h>
#include <stdlib.h> // for atoi
#include <rc/time.h>
#include <rc/math.h>

#define DEFAULT_DIM	140
#define MIN_DIM		1
#define MAX_DIM		500

#define TIMER rc_nanos_thread_time()

// ns consumed just by reading the thread time
#define TIMER_DELAY 2100


static void __print_usage(void)
{
	printf("\n");
	printf("-d         use default matrix size (%dx%d)\n",DEFAULT_DIM,DEFAULT_DIM);
	printf("-s {size}  use custom matrix size\n");
	printf("-h         print this help message\n");
	printf("\n");
}


int main(int argc, char *argv[])
{
	int dim = 0;
	int c,diff,mflops;
	uint64_t t1, t2, flops;
	rc_vector_t b = RC_VECTOR_INITIALIZER;
	rc_vector_t x = RC_VECTOR_INITIALIZER;
	rc_matrix_t A =  RC_MATRIX_INITIALIZER;
	rc_matrix_t AA =  RC_MATRIX_INITIALIZER;
	rc_matrix_t B =  RC_MATRIX_INITIALIZER;
	rc_matrix_t L =  RC_MATRIX_INITIALIZER;
	rc_matrix_t U =  RC_MATRIX_INITIALIZER;
	rc_matrix_t P =  RC_MATRIX_INITIALIZER;
	rc_matrix_t Q =  RC_MATRIX_INITIALIZER;
	rc_matrix_t R =  RC_MATRIX_INITIALIZER;

	// make sure user gave an argument
	if(argc>3){
		printf("Too many arguments given.\n");
		__print_usage();
		return -1;
	}
	if(argc<2){
		printf("Not enough arguments given.\n");
		__print_usage();
		return -1;
	}
	// parse arguments
	opterr = 0;
	while ((c = getopt(argc, argv, "ds:h")) != -1){
		switch (c){
		case 'd': // default size option
			if(dim!=0){
				printf("invalid combination of arguments\n");
				__print_usage();
				return -1;
			}
			dim = DEFAULT_DIM;
			break;
		case 's': // custom size option
			if(dim!=0){
				printf("invalid combination of arguments\n");
				__print_usage();
				return -1;
			}
			dim = atoi(optarg);
			if(dim>MAX_DIM || dim<MIN_DIM){
				printf("requested size out of bounds\n");
				__print_usage();
				return -1;
			}
			break;
		case 'h':
			__print_usage();
			return 0;
		default:
			printf("inavlid argument\n");
			__print_usage();
			return -1;
		}
	}

	// set clock speed to 1000mhz to make sure scaling doesn't effect results
	//rc_set_cpu_freq(FREQ_1000MHZ);
	printf("Starting\n");

	// create a random nxn matrix for later use
	t1 = TIMER;
	rc_matrix_random(&A,dim,dim);
	rc_vector_zeros(&b,dim);
	t2 = TIMER;
	diff = (int)((t2-t1-TIMER_DELAY)/(uint64_t)1000);
	printf("%10dus Time to make random matrix & vector\n", diff);

	// duplicate matrix
	t1 = TIMER;
	rc_matrix_duplicate(A,&AA);
	t2 = TIMER;
	diff = (int)((t2-t1-TIMER_DELAY)/(uint64_t)1000);
	printf("%10dus Time to duplicate matrix\n", diff);

	// Multiply matrices
	rc_matrix_alloc(&B,dim,dim);
	t1 = TIMER;
	rc_matrix_multiply(A, AA, &B);
	t2 = TIMER;
	diff = (int)((t2-t1-TIMER_DELAY)/(uint64_t)1000);
	printf("%10dus Time to multiply matrices\n", diff);

	// calculate floating pointer operations per second, both multiplication
	// and addition count as operations, hence multiply by 2
	flops = ((uint64_t)2*dim*dim*dim*1000000000)/(diff);
	mflops = flops/(uint64_t)1000000;
	printf("%10d MFLOPS multiplying matrices\n", mflops);

	// find determinant
	t1 = TIMER;
	rc_matrix_determinant(A);
	t2 = TIMER;
	diff = (int)((t2-t1-TIMER_DELAY)/(uint64_t)1000);
	printf("%10dus Time to find matrix determinant\n", diff);

	// find inverse
	t1 = TIMER;
	rc_algebra_invert_matrix(A, &AA);
	t2 = TIMER;
	diff = (int)((t2-t1-TIMER_DELAY)/(uint64_t)1000);
	printf("%10dus Time to invert matrix\n", diff);

	// LUP
	rc_matrix_alloc(&L,dim,dim);
	rc_matrix_alloc(&U,dim,dim);
	rc_matrix_alloc(&P,dim,dim);
	t1 = TIMER;
	rc_algebra_lup_decomp(A,&L,&U,&P);
	t2 = TIMER;
	diff = (int)((t2-t1-TIMER_DELAY)/(uint64_t)1000);
	printf("%10dus Time to do LUP decomposition\n", diff);

	// do a QR decomposition on A
	rc_matrix_alloc(&Q,dim,dim);
	rc_matrix_alloc(&R,dim,dim);
	t1 = TIMER;
	rc_algebra_qr_decomp(A,&Q,&R);
	t2 = TIMER;
	diff = (int)((t2-t1-TIMER_DELAY)/(uint64_t)1000);
	printf("%10dus Time to do QR decomposition\n", diff);

	// do a QR decomposition on A
	rc_vector_alloc(&x,dim);
	t1 = TIMER;
	rc_algebra_lin_system_solve(A,b,&x);
	t2 = TIMER;
	diff = (int)((t2-t1-TIMER_DELAY)/(uint64_t)1000);
	printf("%10dus Time to solve linear system\n", diff);

	printf("DONE\n");
	//rc_set_cpu_freq(FREQ_ONDEMAND);
	return 0;
}