home *** CD-ROM | disk | FTP | other *** search
/ Inside Multimedia 1995 July / IMM0795.ISO / share / os2 / sysbench / src / pmb_linp.ckr < prev    next >
Encoding:
Text File  |  1994-11-05  |  26.2 KB  |  1,232 lines
  1. /*
  2. Translated to C by Bonnie Toy 5/88
  3.  
  4. You MUST specify one of -DSP   or -DDP     to compile correctly.
  5. You MUST specify one of -DROLL or -DUNROLL to compile correctly.
  6. You MUST specify a timer option(see below) to compile correctly.
  7.  
  8. To compile double precision version for Sun-4:
  9.    cc -DUNIX -DDP -DROLL -O4 clinpack.c
  10.  
  11. To compile single precision version for Sun-4:
  12.    cc -DUNIX -DSP -DROLL -O4 -fsingle -fsingle2 clinpack.c
  13.  
  14. To obtain   rolled source BLAS, add -DROLL   to the command lines.
  15. To obtain unrolled source BLAS, add -DUNROLL to the command lines.
  16.  
  17. PLEASE NOTE: You can also just 'uncomment' one of the options below.
  18. */
  19.  
  20. /* #define SP     */
  21. #define DP
  22. #define ROLL
  23. /* #define UNROLL */
  24.  
  25. /***************************************************************/
  26. /* Timer options. You MUST uncomment one of the options below  */
  27. /* or compile, for example, with the '-DUNIX' option.          */
  28. /***************************************************************/
  29. /* #define Amiga       */
  30. /* #define UNIX        */
  31. /* #define UNIX_Old    */
  32. /* #define VMS         */
  33. /* #define BORLAND_C   */
  34. /* #define MSC         */
  35. /* #define MAC         */
  36. /* #define IPSC        */
  37. /* #define FORTRAN_SEC */
  38. /* #define GTODay      */
  39. /* #define CTimer      */
  40. /* #define UXPM        */
  41.  
  42. #ifdef SP
  43. #define REAL float
  44. #define ZERO 0.0
  45. #define ONE  1.0
  46. #define PREC "Single "
  47. #endif
  48.  
  49. #ifdef DP
  50. #define REAL double
  51. #define ZERO 0.0e0
  52. #define ONE  1.0e0
  53. #define PREC "Double "
  54. #endif
  55.  
  56. #define NTIMES 100
  57.  
  58. #ifdef ROLL
  59. #define ROLLING "Rolled "
  60. #endif
  61.  
  62. #ifdef UNROLL
  63. #define ROLLING "Unrolled "
  64. #endif
  65.  
  66. //#include <stdio.h>
  67. #include <math.h>
  68.  
  69. static double st[8][6];
  70.  
  71. double pmb_linpack ()
  72. {
  73.    static REAL aa[200][200],a[200][201],b[200],x[200];
  74.    REAL cray,ops,total,norma,normx;
  75.    REAL resid,residn,eps;
  76.    REAL epslon(),kf;
  77.    double t1,tm,tm2,dtime();
  78.    static int ipvt[200],n,i,ntimes,info,lda,ldaa,kflops;
  79.  
  80.    lda = 201;
  81.    ldaa = 200;
  82.    cray = .056; 
  83.    n = 100;
  84.  
  85. /*
  86.    fprintf(stdout,ROLLING);fprintf(stdout,PREC);
  87.    fprintf(stdout,"Precision Linpack\n\n");
  88. */   
  89.    
  90. //   fprintf(stderr,ROLLING);fprintf(stderr,PREC);
  91. //   fprintf(stderr,"Precision Linpack\n\n");
  92.  
  93.     ops = (2.0e0*(n*n*n))/3.0 + 2.0*(n*n);
  94.  
  95.     matgen(a,lda,n,b,&norma);
  96.     t1 = dtime();
  97.     dgefa(a,lda,n,ipvt,&info);
  98.     st[0][0] = dtime() - t1;
  99.     
  100.     t1 = dtime();
  101.     dgesl(a,lda,n,ipvt,b,0);
  102.     st[1][0] = dtime() - t1;
  103.     total = st[0][0] + st[1][0];
  104.  
  105. /*     compute a residual to verify results.  */ 
  106.  
  107.     for (i = 0; i < n; i++) {
  108.            x[i] = b[i];
  109.    }
  110.     matgen(a,lda,n,b,&norma);
  111.     for (i = 0; i < n; i++) {
  112.            b[i] = -b[i];
  113.    }
  114.     dmxpy(n,b,n,lda,x,a);
  115.     resid = 0.0;
  116.     normx = 0.0;
  117.     for (i = 0; i < n; i++) {
  118.            resid = (resid > fabs((double)b[i])) 
  119.      ? resid : fabs((double)b[i]);
  120.            normx = (normx > fabs((double)x[i])) 
  121.      ? normx : fabs((double)x[i]);
  122.    }
  123.     eps = epslon((REAL)ONE);
  124.     residn = resid/( n*norma*normx*eps );
  125.    
  126. /*   printf("   norm. resid      resid           machep");
  127.    printf("         x[0]-1        x[n-1]-1\n");
  128.    printf("%8.1f      %16.8e%16.8e%16.8e%16.8e\n",
  129.       (double)residn, (double)resid, (double)eps, 
  130.            (double)x[0]-1, (double)x[n-1]-1);
  131.  
  132. fprintf(stderr," times are reported for matrices of order %5d\n",n);
  133. fprintf(stderr,"      dgefa      dgesl      total       kflops     unit");
  134. fprintf(stderr,"      ratio\n");
  135. */
  136.     st[2][0] = total;
  137.     st[3][0] = ops/(1.0e3*total);
  138.     st[4][0] = 2.0e3/st[3][0];
  139.     st[5][0] = total/cray;
  140.  
  141. //   fprintf(stderr," times for array with leading dimension of%5d\n",lda);
  142.    print_time(0);
  143.  
  144.     matgen(a,lda,n,b,&norma);
  145.     t1 = dtime();
  146.     dgefa(a,lda,n,ipvt,&info);
  147.     st[0][1] = dtime() - t1;
  148.     
  149.     t1 = dtime();
  150.     dgesl(a,lda,n,ipvt,b,0);
  151.     st[1][1] = dtime() - t1;
  152.     total = st[0][1] + st[1][1];
  153.     
  154.     st[2][1] = total;
  155.     st[3][1] = ops/(1.0e3*total);
  156.     st[4][1] = 2.0e3/st[3][1];
  157.     st[5][1] = total/cray;
  158.  
  159.     matgen(a,lda,n,b,&norma);
  160.     
  161.     t1 = dtime();
  162.     dgefa(a,lda,n,ipvt,&info);
  163.     st[0][2] = dtime() - t1;
  164.     
  165.     t1 = dtime();
  166.     dgesl(a,lda,n,ipvt,b,0);
  167.     st[1][2] = dtime() - t1;
  168.     
  169.     total = st[0][2] + st[1][2];
  170.     st[2][2] = total;
  171.     st[3][2] = ops/(1.0e3*total);
  172.     st[4][2] = 2.0e3/st[3][2];
  173.     st[5][2] = total/cray;
  174.  
  175.     ntimes = NTIMES;
  176.     tm2 = 0.0;
  177.     t1 = dtime();
  178.  
  179.    for (i = 0; i < ntimes; i++) {
  180.            tm = dtime();
  181.       matgen(a,lda,n,b,&norma);
  182.       tm2 = tm2 + dtime() - tm;
  183.       dgefa(a,lda,n,ipvt,&info);
  184.    }
  185.  
  186.     st[0][3] = (dtime() - t1 - tm2)/ntimes;
  187.     t1 = dtime();
  188.  
  189.    for (i = 0; i < ntimes; i++) {
  190.            dgesl(a,lda,n,ipvt,b,0);
  191.    }
  192.  
  193.     st[1][3] = (dtime() - t1)/ntimes;
  194.     total = st[0][3] + st[1][3];
  195.     st[2][3] = total;
  196.     st[3][3] = ops/(1.0e3*total);
  197.     st[4][3] = 2.0e3/st[3][3];
  198.     st[5][3] = total/cray;
  199.  
  200.    print_time(1);
  201.    print_time(2);
  202.    print_time(3);
  203.  
  204.     matgen(aa,ldaa,n,b,&norma);
  205.     t1 = dtime();
  206.     dgefa(aa,ldaa,n,ipvt,&info);
  207.     st[0][4] = dtime() - t1;
  208.     
  209.     t1 = dtime();
  210.     dgesl(aa,ldaa,n,ipvt,b,0);
  211.     st[1][4] = dtime() - t1;
  212.  
  213.     total = st[0][4] + st[1][4];
  214.     st[2][4] = total;
  215.     st[3][4] = ops/(1.0e3*total);
  216.     st[4][4] = 2.0e3/st[3][4];
  217.     st[5][4] = total/cray;
  218.  
  219.     matgen(aa,ldaa,n,b,&norma);
  220.     t1 = dtime();
  221.     dgefa(aa,ldaa,n,ipvt,&info);
  222.     st[0][5] = dtime() - t1;
  223.  
  224.     t1 = dtime();
  225.     dgesl(aa,ldaa,n,ipvt,b,0);
  226.     st[1][5] = dtime() - t1;
  227.  
  228.     total = st[0][5] + st[1][5];
  229.     st[2][5] = total;
  230.     st[3][5] = ops/(1.0e3*total);
  231.     st[4][5] = 2.0e3/st[3][5];
  232.     st[5][5] = total/cray;
  233.  
  234.    matgen(aa,ldaa,n,b,&norma);
  235.    t1 = dtime();
  236.    dgefa(aa,ldaa,n,ipvt,&info);
  237.    st[0][6] = dtime() - t1;
  238.  
  239.    t1 = dtime();
  240.    dgesl(aa,ldaa,n,ipvt,b,0);
  241.    st[1][6] = dtime() - t1;
  242.  
  243.    total = st[0][6] + st[1][6];
  244.    st[2][6] = total;
  245.    st[3][6] = ops/(1.0e3*total);
  246.    st[4][6] = 2.0e3/st[3][6];
  247.    st[5][6] = total/cray;
  248.  
  249.    ntimes = NTIMES;
  250.    tm2 = 0;
  251.    t1 = dtime();
  252.    for (i = 0; i < ntimes; i++) {
  253.       tm = dtime();
  254.       matgen(aa,ldaa,n,b,&norma);
  255.       tm2 = tm2 + dtime() - tm;
  256.       dgefa(aa,ldaa,n,ipvt,&info);
  257.    }
  258.    st[0][7] = (dtime() - t1 - tm2)/ntimes;
  259.    
  260.    t1 = dtime();
  261.    for (i = 0; i < ntimes; i++) {
  262.       dgesl(aa,ldaa,n,ipvt,b,0);
  263.    }
  264.    st[1][7] = (dtime() - t1)/ntimes;
  265.    total = st[0][7] + st[1][7];
  266.    st[2][7] = total;
  267.    st[3][7] = ops/(1.0e3*total);
  268.    st[4][7] = 2.0e3/st[3][7];
  269.    st[5][7] = total/cray;
  270.  
  271.    /* the following code sequence implements the semantics of
  272.       the Fortran intrinsics "nint(min(st[3][3],st[3][7]))"   */
  273. /*
  274.    kf = (st[3][3] < st[3][7]) ? st[3][3] : st[3][7];
  275.    kf = (kf > ZERO) ? (kf + .5) : (kf - .5);
  276.    if (fabs((double)kf) < ONE) 
  277.       kflops = 0;
  278.    else {
  279.       kflops = floor(fabs((double)kf));
  280.       if (kf < ZERO) kflops = -kflops;
  281.    }
  282. */
  283.    if ( st[3][3] < ZERO ) st[3][3] = ZERO;
  284.    if ( st[3][7] < ZERO ) st[3][7] = ZERO;
  285.    kf = st[3][3];
  286.    if ( st[3][7] < st[3][3] ) kf = st[3][7];
  287.    kflops = (int)(kf + 0.5);
  288.  
  289. //   fprintf(stderr," times for array with leading dimension of%4d\n",ldaa);
  290.    print_time(4);
  291.    print_time(5);
  292.    print_time(6);
  293.    print_time(7);
  294. //   fprintf(stderr,ROLLING);fprintf(stderr,PREC);
  295. //   fprintf(stderr," Precision %5d Kflops ; %d Reps \n",kflops,NTIMES);
  296.    return kflops;
  297. }
  298.      
  299. /*----------------------*/ 
  300. static print_time (row)
  301. int row;
  302. {
  303. /*fprintf(stderr,"%11.2f%11.2f%11.2f%11.0f%11.2f%11.2f\n",
  304.        (double)st[0][row], (double)st[1][row], (double)st[2][row], 
  305.        (double)st[3][row], (double)st[4][row], (double)st[5][row]);
  306. */
  307. }
  308.       
  309. /*----------------------*/ 
  310. static matgen(a,lda,n,b,norma)
  311. REAL a[],b[],*norma;
  312. int lda, n;
  313.  
  314. /* We would like to declare a[][lda], but c does not allow it.  In this
  315. function, references to a[i][j] are written a[lda*i+j].  */
  316.  
  317. {
  318.    int init, i, j;
  319.  
  320.    init = 1325;
  321.    *norma = 0.0;
  322.    for (j = 0; j < n; j++) {
  323.       for (i = 0; i < n; i++) {
  324.      init = 3125*init % 65536;
  325.      a[lda*j+i] = (init - 32768.0)/16384.0;
  326.      *norma = (a[lda*j+i] > *norma) ? a[lda*j+i] : *norma;
  327.       }
  328.    }
  329.    for (i = 0; i < n; i++) {
  330.       b[i] = 0.0;
  331.    }
  332.    for (j = 0; j < n; j++) {
  333.       for (i = 0; i < n; i++) {
  334.      b[i] = b[i] + a[lda*j+i];
  335.       }
  336.    }
  337. }
  338.  
  339. /*----------------------*/ 
  340. static dgefa(a,lda,n,ipvt,info)
  341. REAL a[];
  342. int lda,n,ipvt[],*info;
  343.  
  344. /* We would like to declare a[][lda], but c does not allow it.  In this
  345. function, references to a[i][j] are written a[lda*i+j].  
  346. */
  347.  
  348. /*
  349.      dgefa factors a double precision matrix by gaussian elimination.
  350.  
  351.      dgefa is usually called by dgeco, but it can be called
  352.      directly with a saving in time if  rcond  is not needed.
  353.      (time for dgeco) = (1 + 9/n)*(time for dgefa) .
  354.  
  355.      on entry
  356.  
  357.     a       REAL precision[n][lda]
  358.         the matrix to be factored.
  359.  
  360.     lda     integer
  361.         the leading dimension of the array  a .
  362.  
  363.     n       integer
  364.         the order of the matrix  a .
  365.  
  366.      on return
  367.  
  368.     a       an upper triangular matrix and the multipliers
  369.         which were used to obtain it.
  370.         the factorization can be written  a = l*u  where
  371.         l  is a product of permutation and unit lower
  372.         triangular matrices and  u  is upper triangular.
  373.  
  374.     ipvt    integer[n]
  375.         an integer vector of pivot indices.
  376.  
  377.     info    integer
  378.         = 0  normal value.
  379.         = k  if  u[k][k] .eq. 0.0 .  this is not an error
  380.              condition for this subroutine, but it does
  381.              indicate that dgesl or dgedi will divide by zero
  382.              if called.  use  rcond  in dgeco for a reliable
  383.              indication of singularity.
  384.  
  385.      linpack. this version dated 08/14/78 .
  386.      cleve moler, university of new mexico, argonne national lab.
  387.  
  388.      functions
  389.  
  390.      blas daxpy,dscal,idamax
  391. */
  392.  
  393.  
  394. {
  395. /*     internal variables   */
  396.  
  397. REAL t;
  398. int idamax(),j,k,kp1,l,nm1;
  399.  
  400.  
  401. /*     gaussian elimination with partial pivoting   */
  402.  
  403.    *info = 0;
  404.    nm1 = n - 1;
  405.    if (nm1 >=  0) {
  406.       for (k = 0; k < nm1; k++) {
  407.      kp1 = k + 1;
  408.  
  409.         /* find l = pivot index   */
  410.  
  411.      l = idamax(n-k,&a[lda*k+k],1) + k;
  412.      ipvt[k] = l;
  413.  
  414.      /* zero pivot implies this column already 
  415.         triangularized */
  416.  
  417.      if (a[lda*k+l] != ZERO) {
  418.  
  419.         /* interchange if necessary */
  420.  
  421.         if (l != k) {
  422.            t = a[lda*k+l];
  423.            a[lda*k+l] = a[lda*k+k];
  424.            a[lda*k+k] = t; 
  425.         }
  426.  
  427.         /* compute multipliers */
  428.  
  429.         t = -ONE/a[lda*k+k];
  430.         dscal(n-(k+1),t,&a[lda*k+k+1],1);
  431.  
  432.         /* row elimination with column indexing */
  433.  
  434.         for (j = kp1; j < n; j++) {
  435.            t = a[lda*j+l];
  436.            if (l != k) {
  437.           a[lda*j+l] = a[lda*j+k];
  438.           a[lda*j+k] = t;
  439.            }
  440.            daxpy(n-(k+1),t,&a[lda*k+k+1],1,
  441.              &a[lda*j+k+1],1);
  442.           } 
  443.        }
  444.      else { 
  445.              *info = k;
  446.      }
  447.       } 
  448.    }
  449.    ipvt[n-1] = n-1;
  450.    if (a[lda*(n-1)+(n-1)] == ZERO) *info = n-1;
  451. }
  452.  
  453.  
  454. /*----------------------*/ 
  455.  
  456. static dgesl(a,lda,n,ipvt,b,job)
  457. int lda,n,ipvt[],job;
  458. REAL a[],b[];
  459.  
  460. /* We would like to declare a[][lda], but c does not allow it.  In this
  461. function, references to a[i][j] are written a[lda*i+j].  */
  462.  
  463. /*
  464.      dgesl solves the double precision system
  465.      a * x = b  or  trans(a) * x = b
  466.      using the factors computed by dgeco or dgefa.
  467.  
  468.      on entry
  469.  
  470.     a       double precision[n][lda]
  471.         the output from dgeco or dgefa.
  472.  
  473.     lda     integer
  474.         the leading dimension of the array  a .
  475.  
  476.     n       integer
  477.         the order of the matrix  a .
  478.  
  479.     ipvt    integer[n]
  480.         the pivot vector from dgeco or dgefa.
  481.  
  482.     b       double precision[n]
  483.         the right hand side vector.
  484.  
  485.     job     integer
  486.         = 0         to solve  a*x = b ,
  487.         = nonzero   to solve  trans(a)*x = b  where
  488.                 trans(a)  is the transpose.
  489.  
  490.     on return
  491.  
  492.     b       the solution vector  x .
  493.  
  494.      error condition
  495.  
  496.     a division by zero will occur if the input factor contains a
  497.     zero on the diagonal.  technically this indicates singularity
  498.     but it is often caused by improper arguments or improper
  499.     setting of lda .  it will not occur if the subroutines are
  500.     called correctly and if dgeco has set rcond .gt. 0.0
  501.     or dgefa has set info .eq. 0 .
  502.  
  503.      to compute  inverse(a) * c  where  c  is a matrix
  504.      with  p  columns
  505.        dgeco(a,lda,n,ipvt,rcond,z)
  506.        if (!rcond is too small){
  507.           for (j=0,j<p,j++)
  508.             dgesl(a,lda,n,ipvt,c[j][0],0);
  509.       }
  510.  
  511.      linpack. this version dated 08/14/78 .
  512.      cleve moler, university of new mexico, argonne national lab.
  513.  
  514.      functions
  515.  
  516.      blas daxpy,ddot
  517. */
  518.  
  519.  
  520. {
  521. /*     internal variables   */
  522.  
  523.    REAL ddot(),t;
  524.    int k,kb,l,nm1;
  525.  
  526.    nm1 = n - 1;
  527.    if (job == 0) {
  528.  
  529.       /* job = 0 , solve  a * x = b
  530.      first solve  l*y = b       */
  531.  
  532.       if (nm1 >= 1) {
  533.      for (k = 0; k < nm1; k++) {
  534.         l = ipvt[k];
  535.         t = b[l];
  536.         if (l != k){ 
  537.            b[l] = b[k];
  538.            b[k] = t;
  539.         }   
  540.         daxpy(n-(k+1),t,&a[lda*k+k+1],1,&b[k+1],1);
  541.      }
  542.       } 
  543.  
  544.       /* now solve  u*x = y */
  545.  
  546.       for (kb = 0; kb < n; kb++) {
  547.       k = n - (kb + 1);
  548.       b[k] = b[k]/a[lda*k+k];
  549.       t = -b[k];
  550.       daxpy(k,t,&a[lda*k+0],1,&b[0],1);
  551.       }
  552.    }
  553.    else { 
  554.  
  555.       /* job = nonzero, solve  trans(a) * x = b
  556.      first solve  trans(u)*y = b          */
  557.  
  558.       for (k = 0; k < n; k++) {
  559.      t = ddot(k,&a[lda*k+0],1,&b[0],1);
  560.      b[k] = (b[k] - t)/a[lda*k+k];
  561.       }
  562.  
  563.       /* now solve trans(l)*x = y   */
  564.  
  565.       if (nm1 >= 1) {
  566.      for (kb = 1; kb < nm1; kb++) {
  567.         k = n - (kb+1);
  568.         b[k] = b[k] + ddot(n-(k+1),&a[lda*k+k+1],1,&b[k+1],1);
  569.         l = ipvt[k];
  570.         if (l != k) {
  571.            t = b[l];
  572.            b[l] = b[k];
  573.            b[k] = t;
  574.         }
  575.      }
  576.       }
  577.    }
  578. }
  579.  
  580. /*----------------------*/ 
  581.  
  582. static daxpy(n,da,dx,incx,dy,incy)
  583. /*
  584.      constant times a vector plus a vector.
  585.      jack dongarra, linpack, 3/11/78.
  586. */
  587. REAL dx[],dy[],da;
  588. int incx,incy,n;
  589. {
  590.    int i,ix,iy,m,mp1;
  591.  
  592.    if(n <= 0) return;
  593.    if (da == ZERO) return;
  594.  
  595.    if(incx != 1 || incy != 1) {
  596.  
  597.       /* code for unequal increments or equal increments
  598.      not equal to 1                */
  599.  
  600.       ix = 1;
  601.       iy = 1;
  602.       if(incx < 0) ix = (-n+1)*incx + 1;
  603.       if(incy < 0) iy = (-n+1)*incy + 1;
  604.       for (i = 0;i < n; i++) {
  605.      dy[iy] = dy[iy] + da*dx[ix];
  606.      ix = ix + incx;
  607.      iy = iy + incy;
  608.       }
  609.         return;
  610.    }
  611.  
  612.    /* code for both increments equal to 1 */
  613.  
  614. #ifdef ROLL
  615.    for (i = 0;i < n; i++) {
  616.       dy[i] = dy[i] + da*dx[i];
  617.    }
  618. #endif
  619. #ifdef UNROLL
  620.  
  621.    m = n % 4;
  622.    if ( m != 0) {
  623.       for (i = 0; i < m; i++) 
  624.      dy[i] = dy[i] + da*dx[i];
  625.       if (n < 4) return;
  626.    }
  627.    for (i = m; i < n; i = i + 4) {
  628.       dy[i]   = dy[i]   + da*dx[i];
  629.       dy[i+1] = dy[i+1] + da*dx[i+1];
  630.       dy[i+2] = dy[i+2] + da*dx[i+2];
  631.       dy[i+3] = dy[i+3] + da*dx[i+3];
  632.    }
  633. #endif
  634. }
  635.    
  636. /*----------------------*/ 
  637.  
  638. static REAL ddot(n,dx,incx,dy,incy)
  639. /*
  640.      forms the dot product of two vectors.
  641.      jack dongarra, linpack, 3/11/78.
  642. */
  643. REAL dx[],dy[];
  644.  
  645. int incx,incy,n;
  646. {
  647.    REAL dtemp;
  648.    int i,ix,iy,m,mp1;
  649.  
  650.    dtemp = ZERO;
  651.  
  652.    if(n <= 0) return(ZERO);
  653.  
  654.    if(incx != 1 || incy != 1) {
  655.  
  656.       /* code for unequal increments or equal increments
  657.      not equal to 1               */
  658.  
  659.       ix = 0;
  660.       iy = 0;
  661.       if (incx < 0) ix = (-n+1)*incx;
  662.       if (incy < 0) iy = (-n+1)*incy;
  663.       for (i = 0;i < n; i++) {
  664.      dtemp = dtemp + dx[ix]*dy[iy];
  665.      ix = ix + incx;
  666.      iy = iy + incy;
  667.       }
  668.       return(dtemp);
  669.    }
  670.  
  671.    /* code for both increments equal to 1 */
  672.  
  673. #ifdef ROLL
  674.    for (i=0;i < n; i++)
  675.       dtemp = dtemp + dx[i]*dy[i];
  676.    return(dtemp);
  677. #endif
  678. #ifdef UNROLL
  679.  
  680.    m = n % 5;
  681.    if (m != 0) {
  682.       for (i = 0; i < m; i++)
  683.      dtemp = dtemp + dx[i]*dy[i];
  684.       if (n < 5) return(dtemp);
  685.    }
  686.    for (i = m; i < n; i = i + 5) {
  687.       dtemp = dtemp + dx[i]*dy[i] +
  688.       dx[i+1]*dy[i+1] + dx[i+2]*dy[i+2] +
  689.       dx[i+3]*dy[i+3] + dx[i+4]*dy[i+4];
  690.    }
  691.    return(dtemp);
  692. #endif
  693. }
  694.  
  695. /*----------------------*/ 
  696. static dscal(n,da,dx,incx)
  697.  
  698. /*     scales a vector by a constant.
  699.       jack dongarra, linpack, 3/11/78.
  700. */
  701. REAL da,dx[];
  702. int n, incx;
  703. {
  704.    int i,m,mp1,nincx;
  705.  
  706.    if(n <= 0)return;
  707.    if(incx != 1) {
  708.  
  709.       /* code for increment not equal to 1 */
  710.  
  711.       nincx = n*incx;
  712.       for (i = 0; i < nincx; i = i + incx)
  713.      dx[i] = da*dx[i];
  714.       return;
  715.    }
  716.  
  717.    /* code for increment equal to 1 */
  718.  
  719. #ifdef ROLL
  720.    for (i = 0; i < n; i++)
  721.       dx[i] = da*dx[i];
  722. #endif
  723. #ifdef UNROLL
  724.  
  725.    m = n % 5;
  726.    if (m != 0) {
  727.       for (i = 0; i < m; i++)
  728.      dx[i] = da*dx[i];
  729.       if (n < 5) return;
  730.    }
  731.    for (i = m; i < n; i = i + 5){
  732.       dx[i] = da*dx[i];
  733.       dx[i+1] = da*dx[i+1];
  734.       dx[i+2] = da*dx[i+2];
  735.       dx[i+3] = da*dx[i+3];
  736.       dx[i+4] = da*dx[i+4];
  737.    }
  738. #endif
  739.  
  740. }
  741.  
  742. /*----------------------*/ 
  743. static int idamax(n,dx,incx)
  744.  
  745. /*
  746.      finds the index of element having max. absolute value.
  747.      jack dongarra, linpack, 3/11/78.
  748. */
  749.  
  750. REAL dx[];
  751. int incx,n;
  752. {
  753.    REAL dmax;
  754.    int i, ix, itemp;
  755.  
  756.    if( n < 1 ) return(-1);
  757.    if(n ==1 ) return(0);
  758.    if(incx != 1) {
  759.  
  760.       /* code for increment not equal to 1 */
  761.  
  762.       ix = 1;
  763.       dmax = fabs((double)dx[0]);
  764.       ix = ix + incx;
  765.       for (i = 1; i < n; i++) {
  766.      if(fabs((double)dx[ix]) > dmax)  {
  767.         itemp = i;
  768.         dmax = fabs((double)dx[ix]);
  769.      }
  770.      ix = ix + incx;
  771.       }
  772.    }
  773.    else {
  774.  
  775.       /* code for increment equal to 1 */
  776.  
  777.       itemp = 0;
  778.       dmax = fabs((double)dx[0]);
  779.       for (i = 1; i < n; i++) {
  780.      if(fabs((double)dx[i]) > dmax) {
  781.         itemp = i;
  782.         dmax = fabs((double)dx[i]);
  783.      }
  784.       }
  785.    }
  786.    return (itemp);
  787. }
  788.  
  789. /*----------------------*/ 
  790. static REAL epslon (x)
  791. REAL x;
  792. /*
  793.      estimate unit roundoff in quantities of size x.
  794. */
  795.  
  796. {
  797.    REAL a,b,c,eps;
  798. /*
  799.      this program should function properly on all systems
  800.      satisfying the following two assumptions,
  801.     1.  the base used in representing dfloating point
  802.         numbers is not a power of three.
  803.     2.  the quantity  a  in statement 10 is represented to 
  804.         the accuracy used in dfloating point variables
  805.         that are stored in memory.
  806.      the statement number 10 and the go to 10 are intended to
  807.      force optimizing compilers to generate code satisfying 
  808.      assumption 2.
  809.      under these assumptions, it should be true that,
  810.         a  is not exactly equal to four-thirds,
  811.         b  has a zero for its last bit or digit,
  812.         c  is not exactly equal to one,
  813.         eps  measures the separation of 1.0 from
  814.          the next larger dfloating point number.
  815.      the developers of eispack would appreciate being informed
  816.      about any systems where these assumptions do not hold.
  817.  
  818.      *****************************************************************
  819.      this routine is one of the auxiliary routines used by eispack iii
  820.      to avoid machine dependencies.
  821.      *****************************************************************
  822.  
  823.      this version dated 4/6/83.
  824. */
  825.  
  826.    a = 4.0e0/3.0e0;
  827.    eps = ZERO;
  828.    while (eps == ZERO) {
  829.       b = a - ONE;
  830.       c = b + b + b;
  831.       eps = fabs((double)(c-ONE));
  832.    }
  833.    return(eps*fabs((double)x));
  834. }
  835.  
  836. /*----------------------*/ 
  837. static dmxpy (n1, y, n2, ldm, x, m)
  838. REAL y[], x[], m[];
  839. int n1, n2, ldm;
  840.  
  841. /* We would like to declare m[][ldm], but c does not allow it.  In this
  842. function, references to m[i][j] are written m[ldm*i+j].  */
  843.  
  844. /*
  845.    purpose:
  846.      multiply matrix m times vector x and add the result to vector y.
  847.  
  848.    parameters:
  849.  
  850.      n1 integer, number of elements in vector y, and number of rows in
  851.      matrix m
  852.  
  853.      y double [n1], vector of length n1 to which is added 
  854.      the product m*x
  855.  
  856.      n2 integer, number of elements in vector x, and number of columns
  857.      in matrix m
  858.  
  859.      ldm integer, leading dimension of array m
  860.  
  861.      x double [n2], vector of length n2
  862.  
  863.      m double [ldm][n2], matrix of n1 rows and n2 columns
  864.  
  865.  ----------------------------------------------------------------------
  866. */
  867. {
  868.    int j,i,jmin;
  869.    /* cleanup odd vector */
  870.  
  871.    j = n2 % 2;
  872.    if (j >= 1) {
  873.       j = j - 1;
  874.       for (i = 0; i < n1; i++) 
  875.           y[i] = (y[i]) + x[j]*m[ldm*j+i];
  876.    } 
  877.  
  878.    /* cleanup odd group of two vectors */
  879.  
  880.    j = n2 % 4;
  881.    if (j >= 2) {
  882.       j = j - 1;
  883.       for (i = 0; i < n1; i++)
  884.           y[i] = ( (y[i])
  885.                 + x[j-1]*m[ldm*(j-1)+i]) + x[j]*m[ldm*j+i];
  886.    } 
  887.  
  888.    /* cleanup odd group of four vectors */
  889.  
  890.    j = n2 % 8;
  891.    if (j >= 4) {
  892.       j = j - 1;
  893.       for (i = 0; i < n1; i++)
  894.      y[i] = ((( (y[i])
  895.         + x[j-3]*m[ldm*(j-3)+i]) 
  896.         + x[j-2]*m[ldm*(j-2)+i])
  897.         + x[j-1]*m[ldm*(j-1)+i]) + x[j]*m[ldm*j+i];
  898.    } 
  899.  
  900.    /* cleanup odd group of eight vectors */
  901.  
  902.    j = n2 % 16;
  903.    if (j >= 8) {
  904.       j = j - 1;
  905.       for (i = 0; i < n1; i++)
  906.      y[i] = ((((((( (y[i])
  907.         + x[j-7]*m[ldm*(j-7)+i]) + x[j-6]*m[ldm*(j-6)+i])
  908.           + x[j-5]*m[ldm*(j-5)+i]) + x[j-4]*m[ldm*(j-4)+i])
  909.         + x[j-3]*m[ldm*(j-3)+i]) + x[j-2]*m[ldm*(j-2)+i])
  910.         + x[j-1]*m[ldm*(j-1)+i]) + x[j]  *m[ldm*j+i];
  911.    } 
  912.    
  913.    /* main loop - groups of sixteen vectors */
  914.  
  915.    jmin = (n2%16)+16;
  916.    for (j = jmin-1; j < n2; j = j + 16) {
  917.       for (i = 0; i < n1; i++) 
  918.      y[i] = ((((((((((((((( (y[i])
  919.            + x[j-15]*m[ldm*(j-15)+i]) 
  920.         + x[j-14]*m[ldm*(j-14)+i])
  921.          + x[j-13]*m[ldm*(j-13)+i]) 
  922.         + x[j-12]*m[ldm*(j-12)+i])
  923.          + x[j-11]*m[ldm*(j-11)+i]) 
  924.         + x[j-10]*m[ldm*(j-10)+i])
  925.          + x[j- 9]*m[ldm*(j- 9)+i]) 
  926.         + x[j- 8]*m[ldm*(j- 8)+i])
  927.          + x[j- 7]*m[ldm*(j- 7)+i]) 
  928.         + x[j- 6]*m[ldm*(j- 6)+i])
  929.          + x[j- 5]*m[ldm*(j- 5)+i]) 
  930.         + x[j- 4]*m[ldm*(j- 4)+i])
  931.          + x[j- 3]*m[ldm*(j- 3)+i]) 
  932.         + x[j- 2]*m[ldm*(j- 2)+i])
  933.          + x[j- 1]*m[ldm*(j- 1)+i]) 
  934.         + x[j]   *m[ldm*j+i];
  935.    }
  936. } 
  937.  
  938. /*****************************************************/
  939. /* Various timer routines.                           */
  940. /* Al Aburto, aburto@marlin.nosc.mil, 26 Sep 1992    */
  941. /*                                                   */
  942. /* t = dtime() outputs the current time in seconds.  */
  943. /* Use CAUTION as some of these routines will mess   */
  944. /* up when timing across the hour mark!!!            */
  945. /*                                                   */
  946. /* For timing I use the 'user' time whenever         */
  947. /* possible. Using 'user+sys' time is a separate     */
  948. /* issue.                                            */
  949. /*                                                   */
  950. /*****************************************************/
  951.  
  952. /*********************************/
  953. /* Timer code.                   */
  954. /*********************************/
  955. /*******************/
  956. /*  Amiga dtime()  */
  957. /*******************/
  958. #ifdef Amiga
  959. #include <ctype.h>
  960. #define HZ 50
  961.  
  962. double dtime()
  963. {
  964.    double q;
  965.  
  966.    struct   tt {
  967.       long  days;
  968.       long  minutes;
  969.       long  ticks;
  970.    } tt;
  971.  
  972.    DateStamp(&tt);
  973.  
  974.    q = ((double)(tt.ticks + (tt.minutes * 60L * 50L))) / (double)HZ;
  975.  
  976.    return q;
  977. }
  978. #endif
  979.  
  980. /*****************************************************/
  981. /*  UNIX dtime(). This is the preferred UNIX timer.  */
  982. /*  Provided by: Markku Kolkka, mk59200@cc.tut.fi    */
  983. /*  HP-UX Addition by: Bo Thide', bt@irfu.se         */
  984. /*****************************************************/
  985. #ifdef UNIX
  986. #include <sys/time.h>
  987. #include <sys/resource.h>
  988.  
  989. #ifdef __hpux
  990. #include <sys/syscall.h>
  991. #define getrusage(a,b) syscall(SYS_getrusage,a,b)
  992. #endif
  993.  
  994. struct rusage rusage;
  995.  
  996. double dtime()
  997. {
  998.    double q;
  999.  
  1000.    getrusage(RUSAGE_SELF,&rusage);
  1001.  
  1002.    q = (double)(rusage.ru_utime.tv_sec);
  1003.    q = q + (double)(rusage.ru_utime.tv_usec) * 1.0e-06;
  1004.    
  1005.    return q;
  1006. }
  1007. #endif
  1008.  
  1009. /***************************************************/
  1010. /*  UNIX_Old dtime(). This is the old UNIX timer.  */
  1011. /*  Use only if absolutely necessary as HZ may be  */
  1012. /*  ill defined on your system.                    */
  1013. /***************************************************/
  1014. #ifdef UNIX_Old
  1015. #include <sys/types.h>
  1016. #include <sys/times.h>
  1017. #include <sys/param.h>
  1018.  
  1019. #ifndef HZ
  1020. #define HZ 60
  1021. #endif
  1022.  
  1023. struct tms tms;
  1024.  
  1025. double dtime()
  1026. {
  1027.    double q;
  1028.  
  1029.    times(&tms);
  1030.  
  1031.    q = (double)(tms.tms_utime) / (double)HZ;
  1032.    
  1033.    return q;
  1034. }
  1035. #endif
  1036.  
  1037. /*********************************************************/
  1038. /*  VMS dtime() for VMS systems.                         */
  1039. /*  Provided by: RAMO@uvphys.phys.UVic.CA                */
  1040. /*  Some people have run into problems with this timer.  */
  1041. /*********************************************************/
  1042. #ifdef VMS
  1043. #include time
  1044.  
  1045. #ifndef HZ
  1046. #define HZ 100
  1047. #endif
  1048.  
  1049. struct tbuffer_t
  1050.        {
  1051.     int proc_user_time;
  1052.     int proc_system_time;
  1053.     int child_user_time;
  1054.     int child_system_time;
  1055.        };
  1056. struct tbuffer_t tms;
  1057.  
  1058. double dtime()
  1059. {
  1060.    double q;
  1061.  
  1062.    times(&tms);
  1063.  
  1064.    q = (double)(tms.proc_user_time) / (double)HZ;
  1065.    
  1066.    return q;
  1067. }
  1068. #endif
  1069.  
  1070. /******************************/
  1071. /*  BORLAND C dtime() for DOS */
  1072. /******************************/
  1073. #ifdef BORLAND_C
  1074. #include <ctype.h>
  1075. #include <dos.h>
  1076. #include <time.h>
  1077.  
  1078. #define HZ 100
  1079. struct time tnow;
  1080.  
  1081. double dtime()
  1082. {
  1083.    double q;
  1084.  
  1085.    gettime(&tnow);
  1086.  
  1087.    q = 60.0 * (double)(tnow.ti_min);
  1088.    q = q + (double)(tnow.ti_sec);
  1089.    q = q + (double)(tnow.ti_hund)/(double)HZ;
  1090.    
  1091.    return q;
  1092. }
  1093. #endif
  1094.  
  1095. /**************************************/
  1096. /*  Microsoft C (MSC) dtime() for DOS */
  1097. /**************************************/
  1098. #ifdef MSC
  1099. #include <time.h>
  1100. #include <ctype.h>
  1101.  
  1102. #define HZ CLK_TCK
  1103. clock_t tnow;
  1104.  
  1105. double dtime()
  1106. {
  1107.    double q;
  1108.  
  1109.    tnow = clock();
  1110.  
  1111.    q = (double)tnow / (double)HZ;
  1112.    
  1113.    return q;
  1114. }
  1115. #endif
  1116.  
  1117. /*************************************/
  1118. /*  Macintosh (MAC) Think C dtime()  */
  1119. /*************************************/
  1120. #ifdef MAC
  1121. #include <time.h>
  1122.  
  1123. #define HZ 60
  1124.  
  1125. double dtime()
  1126. {
  1127.    double q;
  1128.  
  1129.    q = (double)clock() / (double)HZ;
  1130.    
  1131.    return q;
  1132. }
  1133. #endif
  1134.  
  1135. /************************************************************/
  1136. /*  iPSC/860 (IPSC) dtime() for i860.                       */
  1137. /*  Provided by: Dan Yergeau, yergeau@gloworm.Stanford.EDU  */
  1138. /************************************************************/
  1139. #ifdef IPSC
  1140. extern double dclock();
  1141.  
  1142. double dtime()
  1143. {
  1144.    double q;
  1145.  
  1146.    q = dclock();
  1147.    
  1148.    return q;
  1149. }
  1150. #endif
  1151.  
  1152. /**************************************************/
  1153. /*  FORTRAN dtime() for Cray type systems.        */
  1154. /*  This is the preferred timer for Cray systems. */
  1155. /**************************************************/
  1156. #ifdef FORTRAN_SEC
  1157.  
  1158. fortran double second();
  1159.  
  1160. double dtime()
  1161. {
  1162.    double q;
  1163.  
  1164.    second(&q);
  1165.    
  1166.    return q;
  1167. }
  1168. #endif
  1169.  
  1170. /***********************************************************/
  1171. /*  UNICOS C dtime() for Cray UNICOS systems.  Don't use   */
  1172. /*  unless absolutely necessary as returned time includes  */
  1173. /*  'user+system' time.  Provided by: R. Mike Dority,      */
  1174. /*  dority@craysea.cray.com                                */
  1175. /***********************************************************/
  1176. #ifdef CTimer
  1177. #include <time.h>
  1178.  
  1179. double dtime()
  1180. {
  1181.    double    q;
  1182.    clock_t   t;
  1183.  
  1184.        t = clock();
  1185.  
  1186.        q = (double)t / (double)CLOCKS_PER_SEC;
  1187.  
  1188.        return q;
  1189. }
  1190. #endif
  1191.  
  1192. /********************************************/
  1193. /* Another UNIX timer using gettimeofday(). */
  1194. /* However, getrusage() is preferred.       */
  1195. /********************************************/
  1196. #ifdef GTODay
  1197. #include <sys/time.h>
  1198.  
  1199. struct timeval tnow;
  1200.  
  1201. double dtime()
  1202. {
  1203.    double q;
  1204.  
  1205.    gettimeofday(&tnow,NULL);
  1206.    q = (double)tnow.tv_sec + (double)tnow.tv_usec * 1.0e-6;
  1207.  
  1208.    return q;
  1209. }
  1210. #endif
  1211.  
  1212. /*****************************************************/
  1213. /*  Fujitsu UXP/M timer.                             */
  1214. /*  Provided by: Mathew Lim, ANUSF, M.Lim@anu.edu.au */
  1215. /*****************************************************/
  1216. #ifdef UXPM
  1217. #include <sys/types.h>
  1218. #include <sys/timesu.h>
  1219. struct tmsu rusage;
  1220.  
  1221. double dtime()
  1222. {
  1223.    double q;
  1224.  
  1225.    timesu(&rusage);
  1226.  
  1227.    q = (double)(rusage.tms_utime) * 1.0e-06;
  1228.    
  1229.    return q;
  1230. }
  1231. #endif
  1232.