clinpack.c

来自「开放源码的编译器open watcom 1.6.0版的源代码」· C语言 代码 · 共 1,405 行 · 第 1/3 页

/*
Translated to C by Bonnie Toy 5/88

You MUST specify one of -DSP   or -DDP     to compile correctly.
You MUST specify one of -DROLL or -DUNROLL to compile correctly.
You MUST specify a timer option(see below) to compile correctly.

To compile double precision version for Sun-4:
   cc -DUNIX -DDP -DROLL -O4 clinpack.c

To compile single precision version for Sun-4:
   cc -DUNIX -DSP -DROLL -O4 -fsingle -fsingle2 clinpack.c

To obtain   rolled source BLAS, add -DROLL   to the command lines.
To obtain unrolled source BLAS, add -DUNROLL to the command lines.

PLEASE NOTE: You can also just 'uncomment' one of the options below.
*/

/* #define SP     */
/* #define DP     */
/* #define ROLL   */
/* #define UNROLL */

/***************************************************************/
/* Timer options. You MUST uncomment one of the options below  */
/* or compile, for example, with the '-DUNIX' option.          */
/***************************************************************/
/* #define Amiga       */
/* #define UNIX        */
/* #define UNIX_Old    */
/* #define VMS         */
/* #define BORLAND_C   */
/* #define MSC         */
/* #define MAC         */
/* #define IPSC        */
/* #define FORTRAN_SEC */
/* #define GTODay      */
/* #define CTimer      */
/* #define UXPM        */
/* #define MAC_TMgr    */
/* #define PARIX       */
/* #define POSIX       */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "timer.h"
#include "report.h"

#ifdef SP
#define REAL float
#define ZERO 0.0
#define ONE  1.0
#define PREC "Single "
#endif

#ifdef DP
#define REAL double
#define ZERO 0.0e0
#define ONE  1.0e0
#define PREC "Double "
#endif

#define NTIMES 500

#ifdef ROLL
#define ROLLING "Rolled "
#endif

#ifdef UNROLL
#define ROLLING "Unrolled "
#endif

static double st[8][6];

extern void print_time(), matgen(), dgefa(), dgesl(), daxpy(), dscal(), dmxpy();

void main ()
{
   static REAL aa[200][200],a[200][201],b[200],x[200];
   REAL cray,ops,total,norma,normx;
   REAL resid,residn,eps;
   REAL epslon(),kf;
#if 0
   double t1;
   double tm;
#endif
   double tm2;
   double dtime();
   static int ipvt[200],n,i,ntimes,info,lda,ldaa,kflops;
   static user_timer second_timer;

   lda = 201;
   ldaa = 200;
   cray = .056; 
   n = 100;

   printf(ROLLING); printf(PREC);
   printf("Precision Linpack\n\n");

	ops = (2.0e0*(n*n*n))/3.0 + 2.0*(n*n);

	matgen(a,lda,n,b,&norma);
	TimerOn();
	dgefa(a,lda,n,ipvt,&info);
	TimerOff();
	st[0][0] = TimerElapsed();
	Report( "clinpack(dgefa#1)", st[0][0] );
	
	TimerOn();
	dgesl(a,lda,n,ipvt,b,0);
	TimerOff();
	st[1][0] = TimerElapsed();
	Report( "clinpack(dgesl#1)", st[1][0] );
	total = st[0][0] + st[1][0];

/*     compute a residual to verify results.  */ 

	for (i = 0; i < n; i++)
	   {
		 x[i] = b[i];
	   }
	matgen(a,lda,n,b,&norma);
	for (i = 0; i < n; i++) 
	   {
		 b[i] = -b[i];
	   }
	dmxpy(n,b,n,lda,x,a);
	resid = 0.0;
	normx = 0.0;
	for (i = 0; i < n; i++)
	 {
		 resid = (resid > fabs((double)b[i])) 
	 ? resid : fabs((double)b[i]);
		 normx = (normx > fabs((double)x[i])) 
	 ? normx : fabs((double)x[i]);
	 }
	eps = epslon((REAL)ONE);
	residn = resid/( n*norma*normx*eps );
   
   printf("   norm. resid      resid           machep");
   printf("         x[0]-1        x[n-1]-1\n");
   printf("%8.1f      %16.8e%16.8e%16.8e%16.8e\n",
	  (double)residn, (double)resid, (double)eps, 
		 (double)x[0]-1, (double)x[n-1]-1);

printf(" times are reported for matrices of order %5d\n",n);
printf("      dgefa      dgesl      total       kflops     unit");
printf("      ratio\n");

	st[2][0] = total;
	st[3][0] = ops/(1.0e3*total);
	st[4][0] = 2.0e3/st[3][0];
	st[5][0] = total/cray;

   printf(" times for array with leading dimension of%5d\n",lda);
   print_time(0);

	matgen(a,lda,n,b,&norma);
	TimerOn();
	dgefa(a,lda,n,ipvt,&info);
	TimerOff();
	st[0][1] = TimerElapsed();
	Report( "clinpack(dgefa#2)", st[0][1] );
	
	TimerOn();
	dgesl(a,lda,n,ipvt,b,0);
	TimerOff();
	st[1][1] = TimerElapsed();
	Report( "clinpack(dgesl#2)", st[1][1] );
	total = st[0][1] + st[1][1];
	
	st[2][1] = total;
	st[3][1] = ops/(1.0e3*total);
	st[4][1] = 2.0e3/st[3][1];
	st[5][1] = total/cray;

	matgen(a,lda,n,b,&norma);
	
	TimerOn();
	dgefa(a,lda,n,ipvt,&info);
	TimerOff();
	st[0][2] = TimerElapsed();
	Report( "clinpack(dgefa#3)", st[0][2] );
	
	TimerOn();
	dgesl(a,lda,n,ipvt,b,0);
	TimerOff();
	st[1][2] = TimerElapsed();
	Report( "clinpack(dgesl#3)", st[1][2] );
	
	total = st[0][2] + st[1][2];
	st[2][2] = total;
	st[3][2] = ops/(1.0e3*total);
	st[4][2] = 2.0e3/st[3][2];
	st[5][2] = total/cray;

	ntimes = NTIMES;
	tm2 = 0.0;
	UserTimerOn( &second_timer );

   for (i = 0; i < ntimes; i++) {
	TimerOn();
	matgen(a,lda,n,b,&norma);
	TimerOff();
	tm2 = tm2 + TimerElapsed();
	dgefa(a,lda,n,ipvt,&info);
	}

	UserTimerOff( &second_timer );
	st[0][3] = ( UserTimerElapsed( &second_timer ) - tm2)/ntimes;
	Report( "clinpack(dgefa#4)", st[0][3] );

	TimerOn();
   for (i = 0; i < ntimes; i++) {
		 dgesl(a,lda,n,ipvt,b,0);
	}
	TimerOff();

	st[1][3] = TimerElapsed()/ntimes;
	Report( "clinpack(dgesl#4)", st[1][3] );
	total = st[0][3] + st[1][3];
	st[2][3] = total;
	st[3][3] = ops/(1.0e3*total);
	st[4][3] = 2.0e3/st[3][3];
	st[5][3] = total/cray;

   print_time(1);
   print_time(2);
   print_time(3);

	matgen(aa,ldaa,n,b,&norma);
	TimerOn();
	dgefa(aa,ldaa,n,ipvt,&info);
	TimerOff();
	st[0][4] = TimerElapsed();
	Report( "clinpack(dgefa#5)", st[0][4] );
	
	TimerOn();
	dgesl(aa,ldaa,n,ipvt,b,0);
	TimerOff();
	st[1][4] = TimerElapsed();
	Report( "clinpack(dgesl#5)", st[1][4] );

	total = st[0][4] + st[1][4];
	st[2][4] = total;
	st[3][4] = ops/(1.0e3*total);
	st[4][4] = 2.0e3/st[3][4];
	st[5][4] = total/cray;

	matgen(aa,ldaa,n,b,&norma);
	TimerOn();
	dgefa(aa,ldaa,n,ipvt,&info);
	TimerOff();
	st[0][5] = TimerElapsed();
	Report( "clinpack(dgefa#6)", st[0][5] );

	TimerOn();
	dgesl(aa,ldaa,n,ipvt,b,0);
	TimerOff();
	st[1][5] = TimerElapsed();
	Report( "clinpack(dgesl#6)", st[1][5] );

	total = st[0][5] + st[1][5];
	st[2][5] = total;
	st[3][5] = ops/(1.0e3*total);
	st[4][5] = 2.0e3/st[3][5];
	st[5][5] = total/cray;

   matgen(aa,ldaa,n,b,&norma);
   TimerOn();
   dgefa(aa,ldaa,n,ipvt,&info);
   TimerOff();
   st[0][6] = TimerElapsed();
   Report( "clinpack(dgefa#7)", st[0][6] );

   TimerOn();
   dgesl(aa,ldaa,n,ipvt,b,0);
   TimerOff();
   st[1][6] = TimerElapsed();
   Report( "clinpack(dgesl#7)", st[1][6] );

   total = st[0][6] + st[1][6];
   st[2][6] = total;
   st[3][6] = ops/(1.0e3*total);
   st[4][6] = 2.0e3/st[3][6];
   st[5][6] = total/cray;

   ntimes = NTIMES;
   tm2 = 0;
   UserTimerOn( &second_timer );
   for (i = 0; i < ntimes; i++) {
	TimerOn();
	matgen(aa,ldaa,n,b,&norma);
	TimerOff();
	tm2 = tm2 + TimerElapsed();
	dgefa(aa,ldaa,n,ipvt,&info);
	}
   UserTimerOff( &second_timer );

   st[0][7] = ( UserTimerElapsed( &second_timer ) - tm2 ) / ntimes;
   Report( "clinpack(dgefa#8)", st[0][7] );
   
   TimerOn();
   for (i = 0; i < ntimes; i++) {
	dgesl(aa,ldaa,n,ipvt,b,0);
	}
   TimerOff();

   st[1][7] = TimerElapsed()/ntimes;
   Report( "clinpack(dgesl#8)", st[1][7] );

   total = st[0][7] + st[1][7];
   st[2][7] = total;
   st[3][7] = ops/(1.0e3*total);
   st[4][7] = 2.0e3/st[3][7];
   st[5][7] = total/cray;

   /* the following code sequence implements the semantics of
	the Fortran intrinsics "nint(min(st[3][3],st[3][7]))"   */
/*
   kf = (st[3][3] < st[3][7]) ? st[3][3] : st[3][7];
   kf = (kf > ZERO) ? (kf + .5) : (kf - .5);
   if (fabs((double)kf) < ONE) 
	kflops = 0;
   else {
	kflops = floor(fabs((double)kf));
	if (kf < ZERO) kflops = -kflops;
   }
*/
   if ( st[3][3] < ZERO ) st[3][3] = ZERO;
   if ( st[3][7] < ZERO ) st[3][7] = ZERO;
   kf = st[3][3];
   if ( st[3][7] < st[3][3] ) kf = st[3][7];
   kflops = (int)(kf + 0.5);

   printf(" times for array with leading dimension of%4d\n",ldaa);
   print_time(4);
   print_time(5);
   print_time(6);
   print_time(7);
   printf(ROLLING); printf(PREC);
   printf(" Precision %5d Kflops ; %d Reps \n",kflops,NTIMES);
   exit( EXIT_SUCCESS );
}
     
/*----------------------*/ 
void print_time (row)
int row;
{
printf("%11.2f%11.2f%11.2f%11.0f%11.2f%11.2f\n",
	 (double)st[0][row], (double)st[1][row], (double)st[2][row], 
	 (double)st[3][row], (double)st[4][row], (double)st[5][row]);
}
	
/*----------------------*/ 
void matgen(a,lda,n,b,norma)
REAL a[],b[],*norma;
int lda, n;

/* We would like to declare a[][lda], but c does not allow it.  In this
function, references to a[i][j] are written a[lda*i+j].  */

{
   int init, i, j;

   init = 1325;
   *norma = 0.0;
   for (j = 0; j < n; j++) {
	for (i = 0; i < n; i++) {
	 init = 3125*init % 65536;
	 a[lda*j+i] = (init - 32768.0)/16384.0;
	 *norma = (a[lda*j+i] > *norma) ? a[lda*j+i] : *norma;
	}
   }
   for (i = 0; i < n; i++) {
	  b[i] = 0.0;
   }
   for (j = 0; j < n; j++) {
	for (i = 0; i < n; i++) {
	 b[i] = b[i] + a[lda*j+i];
	}
   }
}

/*----------------------*/ 
void dgefa(a,lda,n,ipvt,info)
REAL a[];
int lda,n,ipvt[],*info;

/* We would like to declare a[][lda], but c does not allow it.  In this
function, references to a[i][j] are written a[lda*i+j].  
*/

/*
     dgefa factors a double precision matrix by gaussian elimination.

     dgefa is usually called by dgeco, but it can be called
     directly with a saving in time if  rcond  is not needed.
     (time for dgeco) = (1 + 9/n)*(time for dgefa) .

     on entry

	a       REAL precision[n][lda]
		the matrix to be factored.

	lda     integer
		the leading dimension of the array  a .

	n       integer
		the order of the matrix  a .

     on return

	a       an upper triangular matrix and the multipliers
		which were used to obtain it.
		the factorization can be written  a = l*u  where
		l  is a product of permutation and unit lower
		triangular matrices and  u  is upper triangular.

	ipvt    integer[n]
		an integer vector of pivot indices.

	info    integer
		= 0  normal value.
		= k  if  u[k][k] .eq. 0.0 .  this is not an error
		     condition for this subroutine, but it does
		     indicate that dgesl or dgedi will divide by zero
		     if called.  use  rcond  in dgeco for a reliable
		     indication of singularity.

     linpack. this version dated 08/14/78 .
     cleve moler, university of new mexico, argonne national lab.

     functions

     blas daxpy,dscal,idamax
*/


{
/*     internal variables   */

REAL t;
int idamax(),j,k,kp1,l,nm1;


/*     gaussian elimination with partial pivoting   */

   *info = 0;
   nm1 = n - 1;
   if (nm1 >=  0) {
	for (k = 0; k < nm1; k++) {
	 kp1 = k + 1;

		/* find l = pivot index   */

	 l = idamax(n-k,&a[lda*k+k],1) + k;
	 ipvt[k] = l;

	 /* zero pivot implies this column already 
	    triangularized */

	 if (a[lda*k+l] != ZERO) {

	    /* interchange if necessary */

	    if (l != k) {