/************************************************************************** * * Subroutine P11 performs the matrix multiplication C = A x B * * Parameters: * * Provided by calling routine: * A = Matrix of multiple precision numbers * B = Matrix of multiple precision numbers * MODULUS = Large integer modulus, expressed as an array of integers * Q = The size of the grid of PEs, Q x Q * M = The size of the sub-matrices on each PE * NUMSIZE = The number of words used to store an MP integer * * Returned by this routine: * C = Matrix of multiple precision numbers * *************************************************************************** * * This version for non-shared-memory systems uses Fox`s algorithm for * matrix multiplication. * **************************************************************************/ #include "bench11.h" /* CVS info */ /* \$Date: 2005/01/10 20:57:47 \$ */ /* \$Revision: 1.2 \$ */ /* \$RCSfile: p11.c,v \$ */ /* \$Name: rel_5 \$ */ static char cvs_info[] = "BMARKGRP \$Date: 2005/01/10 20:57:47 \$ \$Revision: 1.2 \$ \$RCSfile: p11.c,v \$ \$Name: rel_5 \$"; #define TIMEIT 0 /* use macros to address 3-dim arrays with variable dimension */ #define c(I,J,K) c[ (K) + numsize * ( (J) + n * (I) ) ] #define a(I,J,K) a[ (K) + numsize * ( (J) + n * (I) ) ] #define b(I,J,K) b[ (K) + numsize * ( (J) + n * (I) ) ] #define c(I,J,K) c[ (K) + numsize * ( (J) + n * (I) ) ] void p11 ( mplong a[], mplong b[], mplong c[], mplong modulus[], int n, int numsize ) { /* Multiply into tempmp, accumulate sums in tempc */ mplong tempmp [ 2*maxnumsize ]; mplong tempc [ 2*maxnumsize ]; int i, j, k, abcsize, tempsize; double wall(), T0, multtime = 0; /* used if TIMEIT is 1 */ void mult(); abcsize = numsize - MPHL; tempsize = 2*maxnumsize - MPHL; /* initialize c to hold product a*b */ for ( i = 0; i < n; i++ ) { for ( j = 0; j < n; j++ ) { c(i,j,MPALLOC) = abcsize; #if (GMP32 == 0) c(i,j,MPSIZE) = 0; #endif #ifdef GMP c(i,j,MPPTR) = (mplong) & c(i,j,MPDATA); #endif } } /* initialize temporaries */ tempmp[MPALLOC] = tempsize; #if (GMP32 == 0) tempmp[MPSIZE] = 0; #endif tempc[MPALLOC] = tempsize; #ifdef GMP tempmp[MPPTR] = (mplong) & tempmp[MPDATA]; tempc[MPPTR] = (mplong) & tempc[MPDATA]; #endif /* do the matrix multiplication */ /* extern int VERBOSE; */ for ( i = 0; i < n; i++ ) for ( j = 0; j < n; j++ ) { /* zero tempc - accumulate sum during loop */ #if (GMP32 == 0) tempc[MPSIZE] = 0; #else /* erase upper 32 bits of dual-use word 0 */ tempc[MPALLOC] = tempc[MPALLOC] & BitMask32; #endif /* VERBOSE = 2; */ for ( k = 0; k < n; k++ ) { mpmult ( tempmp, &a(i,k,0), &b(k,j,0) ); mpadd ( tempc, tempc, tempmp ); } mpmod ( tempc, tempc, modulus ); #ifdef GMP mpset ( &c(i,j,0), tempc ); #else c(i,j,MPSIZE) = tempc[MPSIZE]; for ( k = MPDATA; k < MPDATA + tempc[MPSIZE]; k++ ) c(i,j,k) = tempc[k]; #endif } }