home *** CD-ROM | disk | FTP | other *** search
/ rtsi.com / 2014.01.www.rtsi.com.tar / www.rtsi.com / OS9 / BENCHMARKS / linpack.c < prev    next >
C/C++ Source or Header  |  2009-11-06  |  21KB  |  907 lines
  1. /*
  2.  
  3. Translated to C by Bonnie Toy 5/88
  4.  
  5. To compile single precision version for Sun-4:
  6.  
  7.     cc -DSP -O4 -fsingle -fsingle2 clinpack.c -lm
  8.  
  9. To compile double precision version for Sun-4:
  10.  
  11.     cc -DDP -O4 clinpack.c -lm
  12.  
  13. To obtain rolled source BLAS, add -DROLL to the command lines.
  14. To obtain unrolled source BLAS, add -DUNROLL to the command lines.
  15.  
  16. You must specify one of -DSP or -DDP to compile correctly.
  17.  
  18. You must specify one of -DROLL or -DUNROLL to compile correctly.
  19.  
  20. */
  21.  
  22. #ifdef SP
  23. #define REAL float
  24. #define ZERO 0.0
  25. #define ONE 1.0
  26. #define PREC "Single "
  27. #endif
  28.  
  29. #ifdef DP
  30. #define REAL double
  31. #define ZERO 0.0e0
  32. #define ONE 1.0e0
  33. #define PREC "Double "
  34. #endif
  35.  
  36. #define NTIMES 10
  37.  
  38. #ifdef ROLL
  39. #define ROLLING "Rolled "
  40. #endif
  41. #ifdef UNROLL
  42. #define ROLLING "Unrolled "
  43. #endif
  44.  
  45. #include <stdio.h>
  46. #include <math.h>
  47.  
  48. static REAL time_array[8][6];
  49.  
  50. main ()
  51. {
  52.     static REAL aa[200][200],a[200][201],b[200],x[200];
  53.     REAL cray,ops,total,norma,normx;
  54.     REAL resid,residn,eps,t1,tm,tm2;
  55.     REAL epslon(),second(),kf;
  56.     static int ipvt[200],n,i,ntimes,info,lda,ldaa,kflops;
  57.  
  58.     lda = 201;
  59.     ldaa = 200;
  60.     cray = .056; 
  61.     n = 100;
  62.  
  63.     fprintf(stdout,ROLLING);fprintf(stdout,PREC);fprintf(stdout,"Precision Linpack\n\n");
  64.     fprintf(stderr,ROLLING);fprintf(stderr,PREC);fprintf(stderr,"Precision Linpack\n\n");
  65.  
  66.         ops = (2.0e0*(n*n*n))/3.0 + 2.0*(n*n);
  67.  
  68.         matgen(a,lda,n,b,&norma);
  69.         t1 = second();
  70.         dgefa(a,lda,n,ipvt,&info);
  71.         time_array[0][0] = second() - t1;
  72.         t1 = second();
  73.         dgesl(a,lda,n,ipvt,b,0);
  74.         time_array[1][0] = second() - t1;
  75.         total = time_array[0][0] + time_array[1][0];
  76.  
  77. /*     compute a residual to verify results.  */ 
  78.  
  79.         for (i = 0; i < n; i++) {
  80.                 x[i] = b[i];
  81.     }
  82.         matgen(a,lda,n,b,&norma);
  83.         for (i = 0; i < n; i++) {
  84.                 b[i] = -b[i];
  85.     }
  86.         dmxpy(n,b,n,lda,x,a);
  87.         resid = 0.0;
  88.         normx = 0.0;
  89.         for (i = 0; i < n; i++) {
  90.                 resid = (resid > fabs((double)b[i])) 
  91.             ? resid : fabs((double)b[i]);
  92.                 normx = (normx > fabs((double)x[i])) 
  93.             ? normx : fabs((double)x[i]);
  94.     }
  95.         eps = epslon((REAL)ONE);
  96.         residn = resid/( n*norma*normx*eps );
  97.     
  98.        printf("     norm. resid      resid           machep");
  99.         printf("         x[0]-1        x[n-1]-1\n");
  100.     printf("  %8.1f      %16.8e%16.8e%16.8e%16.8e\n",
  101.            (double)residn, (double)resid, (double)eps, 
  102.                (double)x[0]-1, (double)x[n-1]-1);
  103.  
  104.        fprintf(stderr,"    times are reported for matrices of order %5d\n",n);
  105.     fprintf(stderr,"      dgefa      dgesl      total       kflops     unit");
  106.     fprintf(stderr,"      ratio\n");
  107.  
  108.         time_array[2][0] = total;
  109.         time_array[3][0] = ops/(1.0e3*total);
  110.         time_array[4][0] = 2.0e3/time_array[3][0];
  111.         time_array[5][0] = total/cray;
  112.  
  113.        fprintf(stderr," times for array with leading dimension of%5d\n",lda);
  114.     print_time(0);
  115.  
  116.         matgen(a,lda,n,b,&norma);
  117.         t1 = second();
  118.         dgefa(a,lda,n,ipvt,&info);
  119.         time_array[0][1] = second() - t1;
  120.         t1 = second();
  121.         dgesl(a,lda,n,ipvt,b,0);
  122.         time_array[1][1] = second() - t1;
  123.         total = time_array[0][1] + time_array[1][1];
  124.         time_array[2][1] = total;
  125.         time_array[3][1] = ops/(1.0e3*total);
  126.         time_array[4][1] = 2.0e3/time_array[3][1];
  127.         time_array[5][1] = total/cray;
  128.  
  129.         matgen(a,lda,n,b,&norma);
  130.         t1 = second();
  131.         dgefa(a,lda,n,ipvt,&info);
  132.         time_array[0][2] = second() - t1;
  133.         t1 = second();
  134.         dgesl(a,lda,n,ipvt,b,0);
  135.         time_array[1][2] = second() - t1;
  136.         total = time_array[0][2] + time_array[1][2];
  137.         time_array[2][2] = total;
  138.         time_array[3][2] = ops/(1.0e3*total);
  139.         time_array[4][2] = 2.0e3/time_array[3][2];
  140.         time_array[5][2] = total/cray;
  141.  
  142.         ntimes = NTIMES;
  143.         tm2 = 0.0;
  144.         t1 = second();
  145.  
  146.     for (i = 0; i < ntimes; i++) {
  147.                 tm = second();
  148.         matgen(a,lda,n,b,&norma);
  149.         tm2 = tm2 + second() - tm;
  150.         dgefa(a,lda,n,ipvt,&info);
  151.     }
  152.  
  153.         time_array[0][3] = (second() - t1 - tm2)/ntimes;
  154.         t1 = second();
  155.  
  156.     for (i = 0; i < ntimes; i++) {
  157.                 dgesl(a,lda,n,ipvt,b,0);
  158.     }
  159.  
  160.         time_array[1][3] = (second() - t1)/ntimes;
  161.         total = time_array[0][3] + time_array[1][3];
  162.         time_array[2][3] = total;
  163.         time_array[3][3] = ops/(1.0e3*total);
  164.         time_array[4][3] = 2.0e3/time_array[3][3];
  165.         time_array[5][3] = total/cray;
  166.  
  167.     print_time(1);
  168.     print_time(2);
  169.     print_time(3);
  170.  
  171.         matgen(aa,ldaa,n,b,&norma);
  172.         t1 = second();
  173.         dgefa(aa,ldaa,n,ipvt,&info);
  174.         time_array[0][4] = second() - t1;
  175.         t1 = second();
  176.         dgesl(aa,ldaa,n,ipvt,b,0);
  177.         time_array[1][4] = second() - t1;
  178.         total = time_array[0][4] + time_array[1][4];
  179.         time_array[2][4] = total;
  180.         time_array[3][4] = ops/(1.0e3*total);
  181.         time_array[4][4] = 2.0e3/time_array[3][4];
  182.         time_array[5][4] = total/cray;
  183.  
  184.         matgen(aa,ldaa,n,b,&norma);
  185.         t1 = second();
  186.         dgefa(aa,ldaa,n,ipvt,&info);
  187.         time_array[0][5] = second() - t1;
  188.         t1 = second();
  189.         dgesl(aa,ldaa,n,ipvt,b,0);
  190.         time_array[1][5] = second() - t1;
  191.         total = time_array[0][5] + time_array[1][5];
  192.         time_array[2][5] = total;
  193.         time_array[3][5] = ops/(1.0e3*total);
  194.         time_array[4][5] = 2.0e3/time_array[3][5];
  195.         time_array[5][5] = total/cray;
  196.  
  197.     matgen(aa,ldaa,n,b,&norma);
  198.     t1 = second();
  199.     dgefa(aa,ldaa,n,ipvt,&info);
  200.     time_array[0][6] = second() - t1;
  201.     t1 = second();
  202.     dgesl(aa,ldaa,n,ipvt,b,0);
  203.     time_array[1][6] = second() - t1;
  204.     total = time_array[0][6] + time_array[1][6];
  205.     time_array[2][6] = total;
  206.     time_array[3][6] = ops/(1.0e3*total);
  207.     time_array[4][6] = 2.0e3/time_array[3][6];
  208.     time_array[5][6] = total/cray;
  209.  
  210.     ntimes = NTIMES;
  211.     tm2 = 0;
  212.     t1 = second();
  213.     for (i = 0; i < ntimes; i++) {
  214.         tm = second();
  215.         matgen(aa,ldaa,n,b,&norma);
  216.         tm2 = tm2 + second() - tm;
  217.         dgefa(aa,ldaa,n,ipvt,&info);
  218.     }
  219.     time_array[0][7] = (second() - t1 - tm2)/ntimes;
  220.     t1 = second();
  221.     for (i = 0; i < ntimes; i++) {
  222.         dgesl(aa,ldaa,n,ipvt,b,0);
  223.     }
  224.     time_array[1][7] = (second() - t1)/ntimes;
  225.     total = time_array[0][7] + time_array[1][7];
  226.     time_array[2][7] = total;
  227.     time_array[3][7] = ops/(1.0e3*total);
  228.     time_array[4][7] = 2.0e3/time_array[3][7];
  229.     time_array[5][7] = total/cray;
  230.  
  231.     /* the following code sequence implements the semantics of
  232.        the Fortran intrinsics "nint(min(time_array[3][3],time_array[3][7]))"    */
  233.  
  234.     kf = (time_array[3][3] < time_array[3][7]) ? time_array[3][3] : time_array[3][7];
  235.     kf = (kf > ZERO) ? (kf + .5) : (kf - .5);
  236.     if (fabs((double)kf) < ONE) 
  237.         kflops = 0;
  238.     else {
  239.         kflops = floor(fabs((double)kf));
  240.         if (kf < ZERO) kflops = -kflops;
  241.     }
  242.  
  243.     fprintf(stderr," times for array with leading dimension of%4d\n",ldaa);
  244.     print_time(4);
  245.     print_time(5);
  246.     print_time(6);
  247.     print_time(7);
  248.     fprintf(stderr,ROLLING);fprintf(stderr,PREC);
  249.     fprintf(stderr," Precision %5d Kflops ; %d Reps \n",kflops,NTIMES);
  250. }
  251.      
  252. /*----------------------*/ 
  253. print_time (row)
  254. int row;
  255. {
  256. fprintf(stderr,"%11.2f%11.2f%11.2f%11.0f%11.2f%11.2f\n",   (double)time_array[0][row],
  257.        (double)time_array[1][row], (double)time_array[2][row], (double)time_array[3][row], 
  258.        (double)time_array[4][row], (double)time_array[5][row]);
  259. }
  260.       
  261. /*----------------------*/ 
  262. matgen(a,lda,n,b,norma)
  263. REAL a[],b[],*norma;
  264. int lda, n;
  265.  
  266. /* We would like to declare a[][lda], but c does not allow it.  In this
  267. function, references to a[i][j] are written a[lda*i+j].  */
  268.  
  269. {
  270.     int init, i, j;
  271.  
  272.     init = 1325;
  273.     *norma = 0.0;
  274.     for (j = 0; j < n; j++) {
  275.         for (i = 0; i < n; i++) {
  276.             init = 3125*init % 65536;
  277.             a[lda*j+i] = (init - 32768.0)/16384.0;
  278.             *norma = (a[lda*j+i] > *norma) ? a[lda*j+i] : *norma;
  279.         }
  280.     }
  281.     for (i = 0; i < n; i++) {
  282.           b[i] = 0.0;
  283.     }
  284.     for (j = 0; j < n; j++) {
  285.         for (i = 0; i < n; i++) {
  286.             b[i] = b[i] + a[lda*j+i];
  287.         }
  288.     }
  289. }
  290.  
  291. /*----------------------*/ 
  292. dgefa(a,lda,n,ipvt,info)
  293. REAL a[];
  294. int lda,n,ipvt[],*info;
  295.  
  296. /* We would like to declare a[][lda], but c does not allow it.  In this
  297. function, references to a[i][j] are written a[lda*i+j].  */
  298. /*
  299.      dgefa factors a double precision matrix by gaussian elimination.
  300.  
  301.      dgefa is usually called by dgeco, but it can be called
  302.      directly with a saving in time if  rcond  is not needed.
  303.      (time for dgeco) = (1 + 9/n)*(time for dgefa) .
  304.  
  305.      on entry
  306.  
  307.         a       REAL precision[n][lda]
  308.                 the matrix to be factored.
  309.  
  310.         lda     integer
  311.                 the leading dimension of the array  a .
  312.  
  313.         n       integer
  314.                 the order of the matrix  a .
  315.  
  316.      on return
  317.  
  318.         a       an upper triangular matrix and the multipliers
  319.                 which were used to obtain it.
  320.                 the factorization can be written  a = l*u  where
  321.                 l  is a product of permutation and unit lower
  322.                 triangular matrices and  u  is upper triangular.
  323.  
  324.         ipvt    integer[n]
  325.                 an integer vector of pivot indices.
  326.  
  327.         info    integer
  328.                 = 0  normal value.
  329.                 = k  if  u[k][k] .eq. 0.0 .  this is not an error
  330.                      condition for this subroutine, but it does
  331.                      indicate that dgesl or dgedi will divide by zero
  332.                      if called.  use  rcond  in dgeco for a reliable
  333.                      indication of singularity.
  334.  
  335.      linpack. this version dated 08/14/78 .
  336.      cleve moler, university of new mexico, argonne national lab.
  337.  
  338.      functions
  339.  
  340.      blas daxpy,dscal,idamax
  341. */
  342.  
  343. {
  344. /*     internal variables    */
  345.  
  346. REAL t;
  347. int idamax(),j,k,kp1,l,nm1;
  348.  
  349.  
  350. /*     gaussian elimination with partial pivoting    */
  351.  
  352.     *info = 0;
  353.     nm1 = n - 1;
  354.     if (nm1 >=  0) {
  355.         for (k = 0; k < nm1; k++) {
  356.             kp1 = k + 1;
  357.  
  358.                   /* find l = pivot index    */
  359.  
  360.             l = idamax(n-k,&a[lda*k+k],1) + k;
  361.             ipvt[k] = l;
  362.  
  363.             /* zero pivot implies this column already 
  364.                triangularized */
  365.  
  366.             if (a[lda*k+l] != ZERO) {
  367.  
  368.                 /* interchange if necessary */
  369.  
  370.                 if (l != k) {
  371.                     t = a[lda*k+l];
  372.                     a[lda*k+l] = a[lda*k+k];
  373.                     a[lda*k+k] = t; 
  374.                 }
  375.  
  376.                 /* compute multipliers */
  377.  
  378.                 t = -ONE/a[lda*k+k];
  379.                 dscal(n-(k+1),t,&a[lda*k+k+1],1);
  380.  
  381.                 /* row elimination with column indexing */
  382.  
  383.                 for (j = kp1; j < n; j++) {
  384.                     t = a[lda*j+l];
  385.                     if (l != k) {
  386.                         a[lda*j+l] = a[lda*j+k];
  387.                         a[lda*j+k] = t;
  388.                     }
  389.                     daxpy(n-(k+1),t,&a[lda*k+k+1],1,
  390.                           &a[lda*j+k+1],1);
  391.                   } 
  392.               }
  393.             else { 
  394.                         *info = k;
  395.             }
  396.         } 
  397.     }
  398.     ipvt[n-1] = n-1;
  399.     if (a[lda*(n-1)+(n-1)] == ZERO) *info = n-1;
  400. }
  401.  
  402. /*----------------------*/ 
  403.  
  404. dgesl(a,lda,n,ipvt,b,job)
  405. int lda,n,ipvt[],job;
  406. REAL a[],b[];
  407.  
  408. /* We would like to declare a[][lda], but c does not allow it.  In this
  409. function, references to a[i][j] are written a[lda*i+j].  */
  410.  
  411. /*
  412.      dgesl solves the double precision system
  413.      a * x = b  or  trans(a) * x = b
  414.      using the factors computed by dgeco or dgefa.
  415.  
  416.      on entry
  417.  
  418.         a       double precision[n][lda]
  419.                 the output from dgeco or dgefa.
  420.  
  421.         lda     integer
  422.                 the leading dimension of the array  a .
  423.  
  424.         n       integer
  425.                 the order of the matrix  a .
  426.  
  427.         ipvt    integer[n]
  428.                 the pivot vector from dgeco or dgefa.
  429.  
  430.         b       double precision[n]
  431.                 the right hand side vector.
  432.  
  433.         job     integer
  434.                 = 0         to solve  a*x = b ,
  435.                 = nonzero   to solve  trans(a)*x = b  where
  436.                             trans(a)  is the transpose.
  437.  
  438.     on return
  439.  
  440.         b       the solution vector  x .
  441.  
  442.      error condition
  443.  
  444.         a division by zero will occur if the input factor contains a
  445.         zero on the diagonal.  technically this indicates singularity
  446.         but it is often caused by improper arguments or improper
  447.         setting of lda .  it will not occur if the subroutines are
  448.         called correctly and if dgeco has set rcond .gt. 0.0
  449.         or dgefa has set info .eq. 0 .
  450.  
  451.      to compute  inverse(a) * c  where  c  is a matrix
  452.      with  p  columns
  453.            dgeco(a,lda,n,ipvt,rcond,z)
  454.            if (!rcond is too small){
  455.                for (j=0,j<p,j++)
  456.                       dgesl(a,lda,n,ipvt,c[j][0],0);
  457.        }
  458.  
  459.      linpack. this version dated 08/14/78 .
  460.      cleve moler, university of new mexico, argonne national lab.
  461.  
  462.      functions
  463.  
  464.      blas daxpy,ddot
  465. */
  466. {
  467. /*     internal variables    */
  468.  
  469.     REAL ddot(),t;
  470.     int k,kb,l,nm1;
  471.  
  472.     nm1 = n - 1;
  473.     if (job == 0) {
  474.  
  475.         /* job = 0 , solve  a * x = b
  476.            first solve  l*y = b        */
  477.  
  478.         if (nm1 >= 1) {
  479.             for (k = 0; k < nm1; k++) {
  480.                 l = ipvt[k];
  481.                 t = b[l];
  482.                 if (l != k){ 
  483.                     b[l] = b[k];
  484.                     b[k] = t;
  485.                 }    
  486.                 daxpy(n-(k+1),t,&a[lda*k+k+1],1,&b[k+1],1);
  487.             }
  488.         } 
  489.  
  490.         /* now solve  u*x = y */
  491.  
  492.         for (kb = 0; kb < n; kb++) {
  493.             k = n - (kb + 1);
  494.             b[k] = b[k]/a[lda*k+k];
  495.             t = -b[k];
  496.             daxpy(k,t,&a[lda*k+0],1,&b[0],1);
  497.         }
  498.     }
  499.     else { 
  500.  
  501.         /* job = nonzero, solve  trans(a) * x = b
  502.            first solve  trans(u)*y = b             */
  503.  
  504.         for (k = 0; k < n; k++) {
  505.             t = ddot(k,&a[lda*k+0],1,&b[0],1);
  506.             b[k] = (b[k] - t)/a[lda*k+k];
  507.         }
  508.  
  509.         /* now solve trans(l)*x = y    */
  510.  
  511.         if (nm1 >= 1) {
  512.             for (kb = 1; kb < nm1; kb++) {
  513.                 k = n - (kb+1);
  514.                 b[k] = b[k] + ddot(n-(k+1),&a[lda*k+k+1],1,&b[k+1],1);
  515.                 l = ipvt[k];
  516.                 if (l != k) {
  517.                     t = b[l];
  518.                     b[l] = b[k];
  519.                     b[k] = t;
  520.                 }
  521.             }
  522.         }
  523.     }
  524. }
  525.  
  526. /*----------------------*/ 
  527.  
  528. daxpy(n,da,dx,incx,dy,incy)
  529. /*
  530.      constant times a vector plus a vector.
  531.      jack dongarra, linpack, 3/11/78.
  532. */
  533. REAL dx[],dy[],da;
  534. int incx,incy,n;
  535. {
  536.     int i,ix,iy,m,mp1;
  537.  
  538.     if(n <= 0) return;
  539.     if (da == ZERO) return;
  540.  
  541.     if(incx != 1 || incy != 1) {
  542.  
  543.         /* code for unequal increments or equal increments
  544.            not equal to 1                     */
  545.  
  546.         ix = 1;
  547.         iy = 1;
  548.         if(incx < 0) ix = (-n+1)*incx + 1;
  549.         if(incy < 0)iy = (-n+1)*incy + 1;
  550.         for (i = 0;i < n; i++) {
  551.             dy[iy] = dy[iy] + da*dx[ix];
  552.             ix = ix + incx;
  553.             iy = iy + incy;
  554.         }
  555.               return;
  556.     }
  557.  
  558.     /* code for both increments equal to 1 */
  559.  
  560. #ifdef ROLL
  561.     for (i = 0;i < n; i++) {
  562.         dy[i] = dy[i] + da*dx[i];
  563.     }
  564. #endif
  565. #ifdef UNROLL
  566.  
  567.     m = n % 4;
  568.     if ( m != 0) {
  569.         for (i = 0; i < m; i++) 
  570.             dy[i] = dy[i] + da*dx[i];
  571.         if (n < 4) return;
  572.     }
  573.     for (i = m; i < n; i = i + 4) {
  574.         dy[i] = dy[i] + da*dx[i];
  575.         dy[i+1] = dy[i+1] + da*dx[i+1];
  576.         dy[i+2] = dy[i+2] + da*dx[i+2];
  577.         dy[i+3] = dy[i+3] + da*dx[i+3];
  578.     }
  579. #endif
  580. }
  581.    
  582. /*----------------------*/ 
  583.  
  584. REAL ddot(n,dx,incx,dy,incy)
  585. /*
  586.      forms the dot product of two vectors.
  587.      jack dongarra, linpack, 3/11/78.
  588. */
  589. REAL dx[],dy[];
  590.  
  591. int incx,incy,n;
  592. {
  593.     REAL dtemp;
  594.     int i,ix,iy,m,mp1;
  595.  
  596.     dtemp = ZERO;
  597.  
  598.     if(n <= 0) return(ZERO);
  599.  
  600.     if(incx != 1 || incy != 1) {
  601.  
  602.         /* code for unequal increments or equal increments
  603.            not equal to 1                    */
  604.  
  605.         ix = 0;
  606.         iy = 0;
  607.         if (incx < 0) ix = (-n+1)*incx;
  608.         if (incy < 0) iy = (-n+1)*incy;
  609.         for (i = 0;i < n; i++) {
  610.             dtemp = dtemp + dx[ix]*dy[iy];
  611.             ix = ix + incx;
  612.             iy = iy + incy;
  613.         }
  614.         return(dtemp);
  615.     }
  616.  
  617.     /* code for both increments equal to 1 */
  618.  
  619. #ifdef ROLL
  620.     for (i=0;i < n; i++)
  621.         dtemp = dtemp + dx[i]*dy[i];
  622.     return(dtemp);
  623. #endif
  624. #ifdef UNROLL
  625.  
  626.     m = n % 5;
  627.     if (m != 0) {
  628.         for (i = 0; i < m; i++)
  629.             dtemp = dtemp + dx[i]*dy[i];
  630.         if (n < 5) return(dtemp);
  631.     }
  632.     for (i = m; i < n; i = i + 5) {
  633.         dtemp = dtemp + dx[i]*dy[i] +
  634.         dx[i+1]*dy[i+1] + dx[i+2]*dy[i+2] +
  635.         dx[i+3]*dy[i+3] + dx[i+4]*dy[i+4];
  636.     }
  637.     return(dtemp);
  638. #endif
  639. }
  640.  
  641. /*----------------------*/ 
  642. dscal(n,da,dx,incx)
  643.  
  644. /*     scales a vector by a constant.
  645.       jack dongarra, linpack, 3/11/78.
  646. */
  647. REAL da,dx[];
  648. int n, incx;
  649. {
  650.     int i,m,mp1,nincx;
  651.  
  652.     if(n <= 0)return;
  653.     if(incx != 1) {
  654.  
  655.         /* code for increment not equal to 1 */
  656.  
  657.         nincx = n*incx;
  658.         for (i = 0; i < nincx; i = i + incx)
  659.             dx[i] = da*dx[i];
  660.         return;
  661.     }
  662.  
  663.     /* code for increment equal to 1 */
  664.  
  665. #ifdef ROLL
  666.     for (i = 0; i < n; i++)
  667.         dx[i] = da*dx[i];
  668. #endif
  669. #ifdef UNROLL
  670.  
  671.     m = n % 5;
  672.     if (m != 0) {
  673.         for (i = 0; i < m; i++)
  674.             dx[i] = da*dx[i];
  675.         if (n < 5) return;
  676.     }
  677.     for (i = m; i < n; i = i + 5){
  678.         dx[i] = da*dx[i];
  679.         dx[i+1] = da*dx[i+1];
  680.         dx[i+2] = da*dx[i+2];
  681.         dx[i+3] = da*dx[i+3];
  682.         dx[i+4] = da*dx[i+4];
  683.     }
  684. #endif
  685.  
  686. }
  687.  
  688. /*----------------------*/ 
  689. int idamax(n,dx,incx)
  690.  
  691. /*
  692.      finds the index of element having max. absolute value.
  693.      jack dongarra, linpack, 3/11/78.
  694. */
  695.  
  696. REAL dx[];
  697. int incx,n;
  698. {
  699.     REAL dmax;
  700.     int i, ix, itemp;
  701.  
  702.     if( n < 1 ) return(-1);
  703.     if(n ==1 ) return(0);
  704.     if(incx != 1) {
  705.  
  706.         /* code for increment not equal to 1 */
  707.  
  708.         ix = 1;
  709.         dmax = fabs((double)dx[0]);
  710.         ix = ix + incx;
  711.         for (i = 1; i < n; i++) {
  712.             if(fabs((double)dx[ix]) > dmax)  {
  713.                 itemp = i;
  714.                 dmax = fabs((double)dx[ix]);
  715.             }
  716.             ix = ix + incx;
  717.         }
  718.     }
  719.     else {
  720.  
  721.         /* code for increment equal to 1 */
  722.  
  723.         itemp = 0;
  724.         dmax = fabs((double)dx[0]);
  725.         for (i = 1; i < n; i++) {
  726.             if(fabs((double)dx[i]) > dmax) {
  727.                 itemp = i;
  728.                 dmax = fabs((double)dx[i]);
  729.             }
  730.         }
  731.     }
  732.     return (itemp);
  733. }
  734.  
  735. /*----------------------*/ 
  736. REAL epslon (x)
  737. REAL x;
  738. /*
  739.      estimate unit roundoff in quantities of size x.
  740. */
  741.  
  742. {
  743.     REAL a,b,c,eps;
  744. /*
  745.      this program should function properly on all systems
  746.      satisfying the following two assumptions,
  747.         1.  the base used in representing dfloating point
  748.             numbers is not a power of three.
  749.         2.  the quantity  a  in statement 10 is represented to 
  750.             the accuracy used in dfloating point variables
  751.             that are stored in memory.
  752.      the statement number 10 and the go to 10 are intended to
  753.      force optimizing compilers to generate code satisfying 
  754.      assumption 2.
  755.      under these assumptions, it should be true that,
  756.             a  is not exactly equal to four-thirds,
  757.             b  has a zero for its last bit or digit,
  758.             c  is not exactly equal to one,
  759.             eps  measures the separation of 1.0 from
  760.                  the next larger dfloating point number.
  761.      the developers of eispack would appreciate being informed
  762.      about any systems where these assumptions do not hold.
  763.  
  764.      *****************************************************************
  765.      this routine is one of the auxiliary routines used by eispack iii
  766.      to avoid machine dependencies.
  767.      *****************************************************************
  768.  
  769.      this version dated 4/6/83.
  770. */
  771.  
  772.     a = 4.0e0/3.0e0;
  773.     eps = ZERO;
  774.     while (eps == ZERO) {
  775.         b = a - ONE;
  776.         c = b + b + b;
  777.         eps = fabs((double)(c-ONE));
  778.     }
  779.     return(eps*fabs((double)x));
  780. }
  781.  
  782. /*----------------------*/ 
  783. dmxpy (n1, y, n2, ldm, x, m)
  784. REAL y[], x[], m[];
  785. int n1, n2, ldm;
  786.  
  787. /* We would like to declare m[][ldm], but c does not allow it.  In this
  788. function, references to m[i][j] are written m[ldm*i+j].  */
  789.  
  790. /*
  791.    purpose:
  792.      multiply matrix m times vector x and add the result to vector y.
  793.  
  794.    parameters:
  795.  
  796.      n1 integer, number of elements in vector y, and number of rows in
  797.          matrix m
  798.  
  799.      y double [n1], vector of length n1 to which is added 
  800.          the product m*x
  801.  
  802.      n2 integer, number of elements in vector x, and number of columns
  803.          in matrix m
  804.  
  805.      ldm integer, leading dimension of array m
  806.  
  807.      x double [n2], vector of length n2
  808.  
  809.      m double [ldm][n2], matrix of n1 rows and n2 columns
  810.  
  811.  ----------------------------------------------------------------------
  812. */
  813. {
  814.     int j,i,jmin;
  815.     /* cleanup odd vector */
  816.  
  817.     j = n2 % 2;
  818.     if (j >= 1) {
  819.         j = j - 1;
  820.         for (i = 0; i < n1; i++) 
  821.                     y[i] = (y[i]) + x[j]*m[ldm*j+i];
  822.     } 
  823.  
  824.     /* cleanup odd group of two vectors */
  825.  
  826.     j = n2 % 4;
  827.     if (j >= 2) {
  828.         j = j - 1;
  829.         for (i = 0; i < n1; i++)
  830.                     y[i] = ( (y[i])
  831.                              + x[j-1]*m[ldm*(j-1)+i]) + x[j]*m[ldm*j+i];
  832.     } 
  833.  
  834.     /* cleanup odd group of four vectors */
  835.  
  836.     j = n2 % 8;
  837.     if (j >= 4) {
  838.         j = j - 1;
  839.         for (i = 0; i < n1; i++)
  840.             y[i] = ((( (y[i])
  841.                    + x[j-3]*m[ldm*(j-3)+i]) 
  842.                    + x[j-2]*m[ldm*(j-2)+i])
  843.                    + x[j-1]*m[ldm*(j-1)+i]) + x[j]*m[ldm*j+i];
  844.     } 
  845.  
  846.     /* cleanup odd group of eight vectors */
  847.  
  848.     j = n2 % 16;
  849.     if (j >= 8) {
  850.         j = j - 1;
  851.         for (i = 0; i < n1; i++)
  852.             y[i] = ((((((( (y[i])
  853.                    + x[j-7]*m[ldm*(j-7)+i]) + x[j-6]*m[ldm*(j-6)+i])
  854.                      + x[j-5]*m[ldm*(j-5)+i]) + x[j-4]*m[ldm*(j-4)+i])
  855.                    + x[j-3]*m[ldm*(j-3)+i]) + x[j-2]*m[ldm*(j-2)+i])
  856.                    + x[j-1]*m[ldm*(j-1)+i]) + x[j]  *m[ldm*j+i];
  857.     } 
  858.     
  859.     /* main loop - groups of sixteen vectors */
  860.  
  861.     jmin = (n2%16)+16;
  862.     for (j = jmin-1; j < n2; j = j + 16) {
  863.         for (i = 0; i < n1; i++) 
  864.             y[i] = ((((((((((((((( (y[i])
  865.                        + x[j-15]*m[ldm*(j-15)+i]) 
  866.                 + x[j-14]*m[ldm*(j-14)+i])
  867.                     + x[j-13]*m[ldm*(j-13)+i]) 
  868.                 + x[j-12]*m[ldm*(j-12)+i])
  869.                     + x[j-11]*m[ldm*(j-11)+i]) 
  870.                 + x[j-10]*m[ldm*(j-10)+i])
  871.                     + x[j- 9]*m[ldm*(j- 9)+i]) 
  872.                 + x[j- 8]*m[ldm*(j- 8)+i])
  873.                     + x[j- 7]*m[ldm*(j- 7)+i]) 
  874.                 + x[j- 6]*m[ldm*(j- 6)+i])
  875.                     + x[j- 5]*m[ldm*(j- 5)+i]) 
  876.                 + x[j- 4]*m[ldm*(j- 4)+i])
  877.                     + x[j- 3]*m[ldm*(j- 3)+i]) 
  878.                 + x[j- 2]*m[ldm*(j- 2)+i])
  879.                     + x[j- 1]*m[ldm*(j- 1)+i]) 
  880.                 + x[j]   *m[ldm*j+i];
  881.     }
  882. } 
  883.  
  884. /*----------------------*/ 
  885. REAL second()
  886. {
  887. #ifdef UNIX
  888. #include <sys/time.h>
  889. #include <sys/resource.h>
  890.  
  891. struct rusage ru;
  892. REAL t ;
  893.  
  894. void getrusage();
  895.  
  896. getrusage(RUSAGE_SELF,&ru) ;
  897.  
  898. t = (REAL) (ru.ru_utime.tv_sec+ru.ru_stime.tv_sec) + 
  899.     ((REAL) (ru.ru_utime.tv_usec+ru.ru_stime.tv_usec))/1.0e6 ;
  900. return t ;
  901. #endif
  902. #ifdef OSK
  903. #include <time.h>
  904. return((REAL)clock() / (REAL)CLK_TCK);
  905. #endif
  906. }
  907.