/*************************************************************************** * * Benchmark #11 -- Multiple Precision Arithmetic * * Given two NxN matrices with entries in the ring of integers modulo * some number MODULUS, compute the product of the matrices. * * Default Parameters: N = 225 MODULUS = (7**183) - 6 * **************************************************************************** * * Main Program for Benchmark #11 * * Call time parameters: * * -v = Verbose mode, invokes printing of intermediate information. * -V = Very verbose mode. * * N = The size of the matrices * Min allowed = 4 * Max allowed = 512 * Default = 225 * * MODULUS = The modulus * Max allowed = 1024 bits * Default = 7^183-6 (514 bits) * A "-" indicates reading the modulus from standard-in. * * Usage: * bench11 32 123456789 (32x32 matrices modulo 123456789) * bench11 N MODULUS (NxN matrices modulo MODULUS) * bench11 N (use the default value for MODULUS) * bench11 N 0 (use the default value for MODULUS) * bench11 0 0 (use the default values for N, MODULUS) * bench11 (use the default values for N, MODULUS) * * bench11 N i (i=1,2,...,10) uses the modulus Mi from the list below * * echo "2^512-1" | bc | bench11 100 - (100x100 matrices mod 2^512-1) * **************************************************************************** * * The check routine is designed to verify the results for any * combination of N and a modulus from the following list: * N = 5, 25, 100, 225, 256, 512 * MODULUS: * M1 = 7^23-6 (65 bits) * M2 = 2^89-1 (89 bits) * M3 = 2^127-1 (127 bits) * M4 = 3^121-2 (192 bits) * M5 = 9^81-8 (257 bits) * M6 = 5^166-4 (386 bits) * M7 = 2^512-1 (512 bits) * M8 = 7^183-6 (514 bits) (the default) * M9 = 3^477-2 (757 bits) * M10 = 5^430-4 (999 bits) * * Other combinations cannot be verified. * **************************************************************************** */ /* global variables - set here, extern in other modules */ int verbose = 0, debug = 0; /* This includes an architecture-specific .h file */ #include "bench11.h" /* on Cray, need to convert fcd returned for modulus into C string */ #ifdef _MPUTIL # include # ifdef _CRAY # include # endif #endif /* CVS info */ /* \$Date: 2005/08/12 19:46:38 \$ */ /* \$Revision: 1.3 \$ */ /* \$RCSfile: m11.c,v \$ */ /* \$Name: rel_5 \$ */ static char cvs_info[] = "BMARKGRP \$Date: 2005/08/12 19:46:38 \$ \$Revision: 1.3 \$ \$RCSfile: m11.c,v \$ \$Name: rel_5 \$"; /*----------------------------------------------------------------------------*/ void main ( int argc, char ** argv ) { /* Arrays to hold the matrices */ /* We will malloc them to be the size needed by each PE - m x m x numsize */ mplong * a, * b, * c; /* An array to hold the modulus */ mplong modulus[maxmodsize]; /* A 2-dim array with the 10 checkable moduli, defined and initialized according to number of bits per MP word */ # if BITSPERWORD == 32 /* Data for a machine with 32 bits per multiple precision word */ static mplong defmods[10][33] = { { 3, 0X30F0A771, 0X7BD152B3, 1 }, { 3, 0XFFFFFFFF, 0XFFFFFFFF, 0X01FFFFFF }, { 4, 0XFFFFFFFF, 0XFFFFFFFF, 0XFFFFFFFF, 0X7FFFFFFF }, { 6, 0X12CACD21, 0X22D2F6C3, 0X92AB1311, 0X9D135AC8, 0XCD3D4156, 0XDBDCFD23 }, { 9, 0XF8DB7C81, 0X786C0065, 0X95D05DC0, 0X5CC1C941, 0XA7987CBA, 0XD6FD8182, 0X278B4D09, 0XB2B6F77A, 1 }, { 13, 0X5BB45685, 0X2F03B89A, 0X0C238E65, 0XD0979D30, 0X7D1B506D, 0X334B46E8, 0XF63B8E74, 0XDC82E031, 0X6842B9A8, 0XBD586A9E, 0X80152BA6, 0XB69CB330, 2 }, { 16, 0XFFFFFFFF, 0XFFFFFFFF, 0XFFFFFFFF, 0XFFFFFFFF, 0XFFFFFFFF, 0XFFFFFFFF, 0XFFFFFFFF, 0XFFFFFFFF, 0XFFFFFFFF, 0XFFFFFFFF, 0XFFFFFFFF, 0XFFFFFFFF, 0XFFFFFFFF, 0XFFFFFFFF, 0XFFFFFFFF, 0XFFFFFFFF }, { 17, 0XF123C871, 0X2B7BAFD2, 0X51D886BB, 0X62E76C16, 0XAD76E211, 0XF2B2DAC1, 0X0EE8263D, 0X97DD013D, 0X49FE574B, 0X590FEBF5, 0XF4018B98, 0X91B9142B, 0XA0C310FE, 0X33304958, 0XF94260BA, 0X5AAA1CA6, 3 }, { 24, 0XBFB72E91, 0X3D890F36, 0XA692A7BB, 0XFCA7523A, 0XBD1C5E4D, 0X88CEDE98, 0X19C47C10, 0X302879A8, 0X98D73E3D, 0X1C54BA18, 0X7EFE1607, 0X5D42F516, 0X5290AE47, 0X4D35BA75, 0XC2A13902, 0XF567F6AB, 0X5C5D0361, 0X66153C40, 0XF9DC5FB4, 0X88E1811F, 0X7923A048, 0XDDCAF0C6, 0X6CCF51C0, 0X104DB4 }, { 32, 0X22E30765, 0XC81109B4, 0XC07755A8, 0XBAE3EC77, 0XAB13676F, 0X4E263E53, 0X3D1E4D7B, 0XD2E5D070, 0X5E140CBB, 0X79C1E1FB, 0X361F35F8, 0X826C9A87, 0X53A3D469, 0XDA1DFDE4, 0XFFCDA2F1, 0X9BC1690A, 0XD6468B5D, 0X856A8777, 0X79E33E50, 0X83412B3C, 0X3391E087, 0XD6C198C0, 0X04D38D4C, 0X34AF758C, 0XB09D292C, 0XE3E3BECA, 0X132CF817, 0X8BF9D3BB, 0X837B535A, 0X1C6AFDFA, 0X2B040960, 0X56 } }; #endif # if BITSPERWORD == 64 /* Data for a machine with 64 bits per multiple precision word */ static mplong defmods[10][17] = { { 2, 0X7BD152B330F0A771, 1 }, { 2, 0XFFFFFFFFFFFFFFFF, 0X01FFFFFF }, { 2, 0XFFFFFFFFFFFFFFFF, 0X7FFFFFFFFFFFFFFF }, { 3, 0X22D2F6C312CACD21, 0X9D135AC892AB1311, 0XDBDCFD23CD3D4156 }, { 5, 0X786C0065F8DB7C81, 0X5CC1C94195D05DC0, 0XD6FD8182A7987CBA, 0XB2B6F77A278B4D09, 1 }, { 7, 0X2F03B89A5BB45685, 0XD0979D300C238E65, 0X334B46E87D1B506D, 0XDC82E031F63B8E74, 0XBD586A9E6842B9A8, 0XB69CB33080152BA6, 2 }, { 8, 0XFFFFFFFFFFFFFFFF, 0XFFFFFFFFFFFFFFFF, 0XFFFFFFFFFFFFFFFF, 0XFFFFFFFFFFFFFFFF, 0XFFFFFFFFFFFFFFFF, 0XFFFFFFFFFFFFFFFF, 0XFFFFFFFFFFFFFFFF, 0XFFFFFFFFFFFFFFFF }, { 9, 0X2B7BAFD2F123C871, 0X62E76C1651D886BB, 0XF2B2DAC1AD76E211, 0X97DD013D0EE8263D, 0X590FEBF549FE574B, 0X91B9142BF4018B98, 0X33304958A0C310FE, 0X5AAA1CA6F94260BA, 3 }, { 12, 0X3D890F36BFB72E91, 0XFCA7523AA692A7BB, 0X88CEDE98BD1C5E4D, 0X302879A819C47C10, 0X1C54BA1898D73E3D, 0X5D42F5167EFE1607, 0X4D35BA755290AE47, 0XF567F6ABC2A13902, 0X66153C405C5D0361, 0X88E1811FF9DC5FB4, 0XDDCAF0C67923A048, 0X00104DB46CCF51C0 }, { 16, 0XC81109B422E30765, 0XBAE3EC77C07755A8, 0X4E263E53AB13676F, 0XD2E5D0703D1E4D7B, 0X79C1E1FB5E140CBB, 0X826C9A87361F35F8, 0XDA1DFDE453A3D469, 0X9BC1690AFFCDA2F1, 0X856A8777D6468B5D, 0X83412B3C79E33E50, 0XD6C198C03391E087, 0X34AF758C04D38D4C, 0XE3E3BECAB09D292C, 0X8BF9D3BB132CF817, 0X1C6AFDFA837B535A, 0X000000562B040960 } }; #endif /* user selects one of above moduli by using 1,2...,10 as 2nd arg. */ int defmodnum; /* Variables for timing purposes */ double cputime(), wall(), x0, y0; double cset, wset, crun, wrun, ccheck, wcheck; /* Other variables */ int n, numsize, modbits, modsize; int o, argind = 1, i, err = 0; double sqrt (); #ifdef _MPUTIL char tmpmod[310]; int tmpmodlen; /* # ifdef _CRAY */ #if defined _CRAY1||_CRAYMPP _fcd FCD; # endif #endif /* * If bench11 is run on a machine without a leading-zero intrinsic, * must use the software version of leadz(n), found in util.c. */ #ifdef SW_LEADZ initlz(); #endif /* Read the command line arguments */ if (argc > 1) { if (argv[1][0] == '-') { if (argv[1][1] == 'v') verbose = 1; if (argv[1][1] == 'V') verbose = 2; argind = 2; } } #ifdef _MPUTIL /* pass flags to Fortran MPUTIL */ mputil(&verbose, &debug); #endif /* Read the dimension of the matrices, check in range */ if (argc > argind) n = atoi(argv[argind]); else n = 0; if(n <= 0) n = 225; /* default */ if (n < 4) { printf("Minimum allowed value for N is 4, N = %d \n", n); exit(0); } if (n > maxmatsize) { printf("Maximum allowed value for N is %d, N = %d \n", maxmatsize, n); exit(0); } /* Prepare the structure for the modulus */ modulus[MPALLOC] = maxmodsize - MPHL; #if (GMP32 == 0) modulus[MPSIZE] = 0; #endif #ifdef GMP modulus[MPPTR] = (mplong) (& modulus[MPDATA]); #endif /* Read the modulus */ if (argc > argind+1) { if (argv[argind+1][0] == '-') { if (verbose > 0) printf( "Reading MODULUS from stdin: \n" ); mpread(modulus); } else { #ifdef _MPUTIL /* must pass Fortran MPFROMSTRING a string ending with a blank */ memset(tmpmod, 0, 310); strcpy(tmpmod, argv[argind+1]); tmpmodlen = (int)strlen(tmpmod); tmpmod[tmpmodlen] = ' '; tmpmod[tmpmodlen+1] = 0; /* # ifdef _CRAY */ #if defined _CRAY1||_CRAYMPP /* Cray needs Fortran Char Descriptor (FCD) */ FCD = _cptofcd( tmpmod, tmpmodlen+1 ); mpfromstring(modulus, FCD); # else /* DEC can just get C string */ mpfromstring(modulus, tmpmod); # endif #else /* call GMP's mpz_set_str with usual C string */ mpfromstring(modulus, argv[argind+1]); #endif } } #if (GMP32 == 0) modsize = modulus[MPSIZE]; #else modsize = modulus[MPALLOC] >> 32; #endif /* Check for default modulus or one of the 10 checkable moduli */ defmodnum = 0; if (modsize == 0) { defmodnum = 8; } else if ( (modsize == 1) && (modulus[MPDATA] <= 10) ) { defmodnum = modulus[MPDATA]; if (defmodnum == 0) defmodnum = 8; } if (defmodnum != 0) { modsize = defmods[defmodnum-1][0]; for (i = 0; i < modsize; i++ ) modulus[MPDATA+i] = defmods[defmodnum-1][i+1]; } numsize = modsize + MPHL; #if (GMP32 == 0) modulus[MPALLOC] = modsize; modulus[MPSIZE] = modsize; #else modulus[MPALLOC] = modsize + ( (mplong)modsize << 32 ); #endif /* Calculate number of bits in the modulus; start with highest limb */ modbits = WORDLEN - leadz( modulus[numsize-1] ); /* Now add total bits in lower-order limbs */ modbits = modbits + BITSPERWORD * (modsize-1); if (verbose) { printf( "Using a modulus of: \n" ); mpwrite( modulus ); printf( " of %d bits \n", modbits ); printf( "We are multiplying two %d x %d matrices \n", n, n ); } cputime(); wall(); /* now we're ready to call s11, p11, c11, r11 to do the work */ /* s11 builds the "random" matrices A and B */ x0 = cputime(); y0 = wall(); /* malloc the matrices A and B */ a = (mplong *) malloc ( 8 * n * n * numsize ); b = (mplong *) malloc ( 8 * n * n * numsize ); if (verbose) printf( "Calling subroutine S11 \n" ); s11 ( a, b, n, numsize, modbits ); if (verbose) printf( "Returned from subroutine S11 \n" ); /* malloc the matrix C used from now on */ c = (mplong *) malloc ( 8 * n * n * numsize ); cset = cputime() - x0; wset = wall() - y0; /* p11 does all the real work */ x0 = cputime(); y0 = wall(); if (verbose) printf( "Calling subroutine P11 \n" ); p11 ( a, b, c, modulus, n, numsize ); if (verbose) printf( "Returned from subroutine P11 \n" ); crun = cputime() - x0; wrun = wall() - y0; /* c11 checks the results */ x0 = cputime(); y0 = wall(); if (verbose) printf( "Calling subroutine C11 \n" ); err = c11 ( c, modulus, n, numsize, modbits ); if (verbose) printf( "Returned from subroutine C11 %d \n", err ); ccheck = cputime() - x0; wcheck = wall() - y0; /* r11 outputs the results */ if (verbose) printf( "Calling subroutine R11 \n" ); r11 ( c, modulus, n, numsize, err, cset, wset, crun, wrun, ccheck, wcheck ); if (verbose) printf( "Returned from subroutine R11 \n" ); }