/*
|
|	example.c
|	22 Feb 94
|	
|	Example of how to use the DSP matrix multiplication functions for the Atari Falcon030.
|	Below the product of two matrices AB is computed three ways: using the matmlflt()
|	function, using the matmlfxp() function, and using normal floating point calculations.
|	Version for Turbo C.
|
*/
 
#include <stdio.h>
#include <math.h>
#include <tos.h>
#include <errno.h>
#include <time.h>
#include "dspmath.h"

float   			A[1600],B[1600],C[1600];   
unsigned long    	Ad[1600],Bd[1600],Cd[1600];

main() {
int				 retval;
unsigned long    m=3,n=3,r=2;	/* Matrix A is of size m x n.
								   Matrix B is of size n x r.
								   Result AB is thus of size m x r. */ 				
unsigned long    i,ii,j,k,maxexp,starttm,endtm,*lptr;
float   		 maxaval,maxbval;

/* 
 * A is a 3 x 3 matrix
 */
A[0] = 0.10; A[1]=  4.24;  A[2]=  0.56; 
A[3]= -1.75; A[4]=  0.001; A[5]=  0.2421; 
A[6]=  0.25; A[7]=  2.5;   A[8]=  0.95;

/*
 * B is a 3 x 2 matrix
 */
B[0] = 0.5;  B[1]= 0.10; 
B[2] = 0.94; B[3]= 10.1; 
B[4] = -6.43; B[5]= 0.0;

retval= Grab_DSP();   /* Initialize the DSP */

if(retval == UH_OH) { 
	printf("Uh oh. DSP initialization failed.\r\n");
	printf("Was DSPMATH.LOD in the current directory?\r\n");
	printf("Hit enter to continue...\r\n");
	Bconin(2);
	return(-1);
} else if(retval == DSP_BUSY) {
	printf("Some other application already using the DSP...\r\n");
	printf("Hit enter to continue...\r\n");
	Bconin(2);
	return(-1);
}



/*
+-------------------------------------------------------------------------------+
| Example of Matrix multiplication using matmlflt(), which takes two            |
| floating point matrices, converts them to fixed point, does the               |
| matrix multiplication, and converts it back.                                  |
+-------------------------------------------------------------------------------+
*/

/*
 * See the header file DSPMATH.H for details.
 */
starttm = clock() * 5;
for(ii=0; ii<100; ii++){  			/* Do 100 AB matrix multiplies */
	matmlflt(A,B,C,m,n,r,Ad,Bd,Cd);
}
endtm = clock() * 5;


printf("\r\nHere's the result AB computed using matmlflt():\r\n");
k=0;

for(i=0; i< m; i++) { 
    for(j=0;j<r;j++,k++) printf(" %9.2f ",C[k]);
    printf("\r\n");
}

printf("Execution time == %ld ms\r\n", (long) endtm - starttm);




/*
+-------------------------------------------------------------------------------+
| Example of Matrix multiplication using plain old floating point arithmetic.   |
|                                                                               |
|                                                                               |
+-------------------------------------------------------------------------------+
*/

printf("starting floating point matrix multiply...\r\n");
starttm = clock() * 5;
for(ii=0; ii<100; ii++){
	for(i=0;i<m;i++) {
   		for(j=0;j<r;j++) {
        	C[i*r+j] = 0;
        	for(k=0;k<n;k++) C[i*r+j] += A[i*n+k]*B[k*r+j];
    	}
	}
}
endtm = clock() * 5;

printf("Here's the result AB computed using\r\nthe floating point calculations:\r\n");
k=0;

for(i=0; i< m; i++) { 
    for(j=0;j<r;j++,k++) printf(" %9.2f ",C[k]);
    printf("\r\n");
}

printf("Execution time == %ld ms\r\n", (long) endtm - starttm);



/*
+-------------------------------------------------------------------------------+
| Example of Matrix multiplication using matmlfxp(), which takes two            |
| fixed point matrices. Because it does not have the conversion overhead of     |
| of matmlflt(), it is much faster.                                             |
+-------------------------------------------------------------------------------+
*/

/*
 * If the numbers are stored as fixed point numbers (scaled to  >= -1.0 and < 1.0), 
 * then we can use the fixed point multiplication which doesn't have the float-to-fixed
 * conversion overhead, and is therefore much faster. Here's the example:
 */

/* Get the largest value in the A matrix, then scale each value to be between 1.0 and -1.0 
 * By the way, this isn't the way matmlflt() does the conversion; it doesn't use floating
 * point operations, rather it shifts the mantissa and keeps track of the exponent and
 * signs, but I wanted to do it this way here so it was clearer what is going on.
 */
for(i=0,maxaval= -HUGE_VAL;i<m*n;i++) if(fabs(A[i]) > maxaval) maxaval = 1.1*A[i];
for(i=0;i<m*n;i++) Ad[i] = f_to_DSP((float) A[i]/maxaval ); 

/* Get the largest value in the B matrix, then scale each value to be between 1.0 and -1.0 */
for(i=0,maxbval= -HUGE_VAL;i<n*r;i++) if(fabs(B[i]) >maxbval) maxbval = 1.1*B[i];
for(i=0;i<n*r;i++) Bd[i] = f_to_DSP((float) B[i]/maxbval ); 



/* OK, pretend like we did the above conversion along the way as user data
 * was entered, and our main storage is to keep the numbers stored in the fixed
 * point arrays Ad and Bd. Then when it comes time to do the arithmetic, we do this: 
 */  
starttm = clock() * 5;
for(ii=0; ii<100; ii++){  			/* Do 100 AB matrix multiplies */
	matmlfxp(m,n,r,Ad,Bd,Cd);
}
endtm = clock() * 5;


/* Before the data is used in some way, we may need to convert it back to 
 * floating point numbers, like this:
 */
for(i=0;i<m*r;i++) C[i] = maxaval*maxbval*DSP_to_f( (unsigned long) Cd[i]);


printf("\r\nHere's the result AB computed using matmlfxp():\r\n");
k=0;

for(i=0; i< m; i++) { 
    for(j=0;j<r;j++,k++) printf(" %9.2f ",C[k]);
    printf("\r\n");
}

printf("Execution time == %ld ms\r\n", (long) endtm - starttm);

   
printf("Hit enter to continue...\r\n");
Release_DSP(); /* Release the DSP prior to exiting so other applications can use it. */
Bconin(2);
return(0);
}