home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / kaffe-0.5p4-src.tgz / tar.out / contrib / kaffe / test / floats / linpack / LinpackApp.java
Text File  |  1996-09-28  |  14KB  |  627 lines
  1. /*
  2.  
  3. Optimized by Jonathan Hardwick (jch@cs.cmu.edu), 3/28/96
  4. Compare to Linkpack.java.
  5. Optimizations performed:
  6.  - added "final" modifier to performance-critical methods.
  7.  - changed lines of the form "a[i] = a[i] + x" to "a[i] += x".
  8.  - minimized array references using common subexpression elimination.
  9.  - eliminated unused variables.
  10.  - undid an unrolled loop.
  11.  - added temporary 1D arrays to hold frequently-used columns of 2D arrays.
  12.  - wrote my own abs() method
  13. See http://www.cs.cmu.edu/~jch/java/linpack.html for more details.
  14.  
  15.  
  16. Ported to Java by Reed Wade  (wade@cs.utk.edu) 2/96
  17. built using JDK 1.0 on solaris
  18. using "javac -O Linpack.java"
  19.  
  20.  
  21. Translated to C by Bonnie Toy 5/88
  22.   (modified on 2/25/94  to fix a problem with daxpy  for
  23.    unequal increments or equal increments not equal to 1.
  24.      Jack Dongarra)
  25.  
  26. */
  27.  
  28.  
  29. import java.io.*;
  30. import java.util.*;
  31.  
  32.  
  33. class LinpackClass {
  34.  
  35.     String applet_version = "LinpackJavaV2.0";
  36.  
  37.     double mflops_result = 0.0;
  38.     double residn_result = 0.0;
  39.     double time_result = 0.0;
  40.     double eps_result = 0.0;
  41.  
  42.  
  43.  
  44.  
  45.   final double abs (double d) {
  46.     return (d >= 0) ? d : -d;
  47.   }
  48.  
  49.   
  50.  
  51.   double second_orig = -1;
  52.   double second()
  53.   {
  54.     if (second_orig==-1) {
  55.       second_orig = System.currentTimeMillis();
  56.     }
  57.     return (System.currentTimeMillis() - second_orig)/1000;
  58.   }
  59.  
  60.   
  61.   
  62.   public void run_benchmark ()
  63.   {
  64.     double a[][] = new double[200][201];
  65.     double b[] = new double[200];
  66.     double x[] = new double[200];
  67.     double cray,ops,total,norma,normx;
  68.     double resid,time;
  69.     double kf;
  70.     int n,i,ntimes,info,lda,ldaa,kflops;
  71.     int ipvt[] = new int[200];
  72.     
  73.     //double mflops_result;
  74.     //double residn_result;
  75.     //double time_result;
  76.     //double eps_result;
  77.  
  78.     lda = 201;
  79.     ldaa = 200;
  80.     cray = .056;
  81.     n = 100;
  82.     
  83.     ops = (2.0e0*(n*n*n))/3.0 + 2.0*(n*n);
  84.     
  85.     norma = matgen(a,lda,n,b);
  86.     time = second();
  87.     info = dgefa(a,lda,n,ipvt);
  88.     dgesl(a,lda,n,ipvt,b,0);
  89.     total = second() - time;
  90.     
  91.     for (i = 0; i < n; i++) {
  92.       x[i] = b[i];
  93.     }
  94.     norma = matgen(a,lda,n,b);
  95.     for (i = 0; i < n; i++) {
  96.       b[i] = -b[i];
  97.     }
  98.     dmxpy(n,b,n,lda,x,a);
  99.     resid = 0.0;
  100.     normx = 0.0;
  101.     for (i = 0; i < n; i++) {
  102.       resid = (resid > abs(b[i])) ? resid : abs(b[i]);
  103.       normx = (normx > abs(x[i])) ? normx : abs(x[i]);
  104.     }
  105.     
  106.     eps_result = epslon((double)1.0);
  107. /*
  108.  
  109.     residn_result = resid/( n*norma*normx*eps_result );
  110.     time_result = total;
  111.     mflops_result = ops/(1.0e6*total);
  112.  
  113.     return ("Mflops/s: " + mflops_result +
  114.         "  Time: " + time_result + " secs" +
  115.         "  Norm Res: " + residn_result +
  116.         "  Precision: " + eps_result);
  117. */
  118.     residn_result = resid/( n*norma*normx*eps_result );
  119.     residn_result += 0.005; // for rounding
  120.     residn_result = (int)(residn_result*100);
  121.     residn_result /= 100;
  122.  
  123.     time_result = total;
  124.     time_result += 0.005; // for rounding
  125.     time_result = (int)(time_result*100);
  126.     time_result /= 100;
  127.  
  128.     mflops_result = ops/(1.0e6*total);
  129.     mflops_result += 0.0005; // for rounding
  130.     mflops_result = (int)(mflops_result*1000);
  131.     mflops_result /= 1000;
  132.  
  133.   }
  134.   
  135.  
  136.   
  137.   final double matgen (double a[][], int lda, int n, double b[])
  138.   {
  139.     double norma;
  140.     int init, i, j;
  141.     
  142.     init = 1325;
  143.     norma = 0.0;
  144.     for (j = 0; j < n; j++) {
  145.       for (i = 0; i < n; i++) {
  146.     init = 3125*init % 65536;
  147.     a[j][i] = (init - 32768.0)/16384.0;
  148.     norma = (a[j][i] > norma) ? a[j][i] : norma;
  149.       }
  150.     }
  151.     for (i = 0; i < n; i++) {
  152.       b[i] = 0.0;
  153.     }
  154.     for (j = 0; j < n; j++) {
  155.       for (i = 0; i < n; i++) {
  156.     b[i] += a[j][i];
  157.       }
  158.     }
  159.     
  160.     return norma;
  161.   }
  162.   
  163.  
  164.   
  165.   /*
  166.     dgefa factors a double precision matrix by gaussian elimination.
  167.     
  168.     dgefa is usually called by dgeco, but it can be called
  169.     directly with a saving in time if  rcond  is not needed.
  170.     (time for dgeco) = (1 + 9/n)*(time for dgefa) .
  171.     
  172.     on entry
  173.     
  174.     a       double precision[n][lda]
  175.     the matrix to be factored.
  176.     
  177.     lda     integer
  178.     the leading dimension of the array  a .
  179.     
  180.     n       integer
  181.     the order of the matrix  a .
  182.     
  183.     on return
  184.     
  185.     a       an upper triangular matrix and the multipliers
  186.     which were used to obtain it.
  187.     the factorization can be written  a = l*u  where
  188.     l  is a product of permutation and unit lower
  189.     triangular matrices and  u  is upper triangular.
  190.     
  191.     ipvt    integer[n]
  192.     an integer vector of pivot indices.
  193.     
  194.     info    integer
  195.     = 0  normal value.
  196.     = k  if  u[k][k] .eq. 0.0 .  this is not an error
  197.     condition for this subroutine, but it does
  198.     indicate that dgesl or dgedi will divide by zero
  199.     if called.  use  rcond  in dgeco for a reliable
  200.     indication of singularity.
  201.     
  202.     linpack. this version dated 08/14/78.
  203.     cleve moler, university of new mexico, argonne national lab.
  204.     
  205.     functions
  206.     
  207.     blas daxpy,dscal,idamax
  208.   */
  209.   final int dgefa( double a[][], int lda, int n, int ipvt[])
  210.   {
  211.     double[] col_k, col_j;
  212.     double t;
  213.     int j,k,kp1,l,nm1;
  214.     int info;
  215.     
  216.     // gaussian elimination with partial pivoting
  217.     
  218.     info = 0;
  219.     nm1 = n - 1;
  220.     if (nm1 >=  0) {
  221.       for (k = 0; k < nm1; k++) {
  222.     col_k = a[k];
  223.     kp1 = k + 1;
  224.     
  225.     // find l = pivot index
  226.     
  227.     l = idamax(n-k,col_k,k,1) + k;
  228.     ipvt[k] = l;
  229.     
  230.     // zero pivot implies this column already triangularized
  231.     
  232.     if (col_k[l] != 0) {
  233.       
  234.       // interchange if necessary
  235.       
  236.       if (l != k) {
  237.         t = col_k[l];
  238.         col_k[l] = col_k[k];
  239.         col_k[k] = t;
  240.       }
  241.       
  242.       // compute multipliers
  243.       
  244.       t = -1.0/col_k[k];
  245.       dscal(n-(kp1),t,col_k,kp1,1);
  246.       
  247.       // row elimination with column indexing
  248.       
  249.       for (j = kp1; j < n; j++) {
  250.         col_j = a[j];
  251.         t = col_j[l];
  252.         if (l != k) {
  253.           col_j[l] = col_j[k];
  254.           col_j[k] = t;
  255.         }
  256.         daxpy(n-(kp1),t,col_k,kp1,1,
  257.           col_j,kp1,1);
  258.       }
  259.     }
  260.     else {
  261.       info = k;
  262.     }
  263.       }
  264.     }
  265.     ipvt[n-1] = n-1;
  266.     if (a[(n-1)][(n-1)] == 0) info = n-1;
  267.     
  268.     return info;
  269.   }
  270.  
  271.   
  272.   
  273.   /*
  274.     dgesl solves the double precision system
  275.     a * x = b  or  trans(a) * x = b
  276.     using the factors computed by dgeco or dgefa.
  277.   
  278.     on entry
  279.   
  280.     a       double precision[n][lda]
  281.     the output from dgeco or dgefa.
  282.   
  283.     lda     integer
  284.     the leading dimension of the array  a .
  285.     
  286.     n       integer
  287.     the order of the matrix  a .
  288.   
  289.     ipvt    integer[n]
  290.     the pivot vector from dgeco or dgefa.
  291.  
  292.     b       double precision[n]
  293.     the right hand side vector.
  294.     
  295.     job     integer
  296.     = 0         to solve  a*x = b ,
  297.     = nonzero   to solve  trans(a)*x = b  where
  298.     trans(a)  is the transpose.
  299.     
  300.     on return
  301.     
  302.     b       the solution vector  x .
  303.     
  304.     error condition
  305.     
  306.     a division by zero will occur if the input factor contains a
  307.     zero on the diagonal.  technically this indicates singularity
  308.     but it is often caused by improper arguments or improper
  309.     setting of lda .  it will not occur if the subroutines are
  310.     called correctly and if dgeco has set rcond .gt. 0.0
  311.     or dgefa has set info .eq. 0 .
  312.     
  313.     to compute  inverse(a) * c  where  c  is a matrix
  314.     with  p  columns
  315.     dgeco(a,lda,n,ipvt,rcond,z)
  316.     if (!rcond is too small){
  317.     for (j=0,j<p,j++)
  318.     dgesl(a,lda,n,ipvt,c[j][0],0);
  319.     }
  320.     
  321.     linpack. this version dated 08/14/78 .
  322.     cleve moler, university of new mexico, argonne national lab.
  323.     
  324.     functions
  325.     
  326.     blas daxpy,ddot
  327.   */
  328.   final void dgesl( double a[][], int lda, int n, int ipvt[], double b[], int job)
  329.   {
  330.     double t;
  331.     int k,kb,l,nm1,kp1;
  332.  
  333.     nm1 = n - 1;
  334.     if (job == 0) {
  335.  
  336.       // job = 0 , solve  a * x = b.  first solve  l*y = b
  337.  
  338.       if (nm1 >= 1) {
  339.     for (k = 0; k < nm1; k++) {
  340.       l = ipvt[k];
  341.       t = b[l];
  342.       if (l != k){
  343.         b[l] = b[k];
  344.         b[k] = t;
  345.       }
  346.       kp1 = k + 1;
  347.       daxpy(n-(kp1),t,a[k],kp1,1,b,kp1,1);
  348.     }
  349.       }
  350.  
  351.       // now solve  u*x = y
  352.  
  353.       for (kb = 0; kb < n; kb++) {
  354.     k = n - (kb + 1);
  355.     b[k] /= a[k][k];
  356.     t = -b[k];
  357.     daxpy(k,t,a[k],0,1,b,0,1);
  358.       }
  359.     }
  360.     else {
  361.  
  362.       // job = nonzero, solve  trans(a) * x = b.  first solve  trans(u)*y = b
  363.  
  364.       for (k = 0; k < n; k++) {
  365.     t = ddot(k,a[k],0,1,b,0,1);
  366.     b[k] = (b[k] - t)/a[k][k];
  367.       }
  368.  
  369.       // now solve trans(l)*x = y 
  370.  
  371.       if (nm1 >= 1) {
  372.     for (kb = 1; kb < nm1; kb++) {
  373.       k = n - (kb+1);
  374.       kp1 = k + 1;
  375.       b[k] += ddot(n-(kp1),a[k],kp1,1,b,kp1,1);
  376.       l = ipvt[k];
  377.       if (l != k) {
  378.         t = b[l];
  379.         b[l] = b[k];
  380.         b[k] = t;
  381.       }
  382.     }
  383.       }
  384.     }
  385.   }
  386.  
  387.  
  388.  
  389.   /*
  390.     constant times a vector plus a vector.
  391.     jack dongarra, linpack, 3/11/78.
  392.   */
  393.   final void daxpy( int n, double da, double dx[], int dx_off, int incx,
  394.           double dy[], int dy_off, int incy)
  395.   {
  396.     int i,ix,iy;
  397.  
  398.     if ((n > 0) && (da != 0)) {
  399.       if (incx != 1 || incy != 1) {
  400.  
  401.     // code for unequal increments or equal increments not equal to 1
  402.  
  403.     ix = 0;
  404.     iy = 0;
  405.     if (incx < 0) ix = (-n+1)*incx;
  406.     if (incy < 0) iy = (-n+1)*incy;
  407.     for (i = 0;i < n; i++) {
  408.       dy[iy +dy_off] += da*dx[ix +dx_off];
  409.       ix += incx;
  410.       iy += incy;
  411.     }
  412.     return;
  413.       } else {
  414.  
  415.     // code for both increments equal to 1
  416.  
  417.     for (i=0; i < n; i++)
  418.       dy[i +dy_off] += da*dx[i +dx_off];
  419.       }
  420.     }
  421.   }
  422.  
  423.  
  424.  
  425.   /*
  426.     forms the dot product of two vectors.
  427.     jack dongarra, linpack, 3/11/78.
  428.   */
  429.   final double ddot( int n, double dx[], int dx_off, int incx, double dy[],
  430.            int dy_off, int incy)
  431.   {
  432.     double dtemp;
  433.     int i,ix,iy;
  434.  
  435.     dtemp = 0;
  436.  
  437.     if (n > 0) {
  438.       
  439.       if (incx != 1 || incy != 1) {
  440.  
  441.     // code for unequal increments or equal increments not equal to 1
  442.  
  443.     ix = 0;
  444.     iy = 0;
  445.     if (incx < 0) ix = (-n+1)*incx;
  446.     if (incy < 0) iy = (-n+1)*incy;
  447.     for (i = 0;i < n; i++) {
  448.       dtemp += dx[ix +dx_off]*dy[iy +dy_off];
  449.       ix += incx;
  450.       iy += incy;
  451.     }
  452.       } else {
  453.  
  454.     // code for both increments equal to 1
  455.     
  456.     for (i=0;i < n; i++)
  457.       dtemp += dx[i +dx_off]*dy[i +dy_off];
  458.       }
  459.     }
  460.     return(dtemp);
  461.   }
  462.  
  463.   
  464.   
  465.   /*
  466.     scales a vector by a constant.
  467.     jack dongarra, linpack, 3/11/78.
  468.   */
  469.   final void dscal( int n, double da, double dx[], int dx_off, int incx)
  470.   {
  471.     int i,nincx;
  472.  
  473.     if (n > 0) {
  474.       if (incx != 1) {
  475.  
  476.     // code for increment not equal to 1
  477.  
  478.     nincx = n*incx;
  479.     for (i = 0; i < nincx; i += incx)
  480.       dx[i +dx_off] *= da;
  481.       } else {
  482.  
  483.     // code for increment equal to 1
  484.  
  485.     for (i = 0; i < n; i++)
  486.       dx[i +dx_off] *= da;
  487.       }
  488.     }
  489.   }
  490.  
  491.   
  492.   
  493.   /*
  494.     finds the index of element having max. absolute value.
  495.     jack dongarra, linpack, 3/11/78.
  496.   */
  497.   final int idamax( int n, double dx[], int dx_off, int incx)
  498.   {
  499.     double dmax, dtemp;
  500.     int i, ix, itemp=0;
  501.  
  502.     if (n < 1) {
  503.       itemp = -1;
  504.     } else if (n ==1) {
  505.       itemp = 0;
  506.     } else if (incx != 1) {
  507.  
  508.       // code for increment not equal to 1
  509.  
  510.       dmax = abs(dx[0 +dx_off]);
  511.       ix = 1 + incx;
  512.       for (i = 1; i < n; i++) {
  513.     dtemp = abs(dx[ix + dx_off]);
  514.     if (dtemp > dmax)  {
  515.       itemp = i;
  516.       dmax = dtemp;
  517.     }
  518.     ix += incx;
  519.       }
  520.     } else {
  521.  
  522.       // code for increment equal to 1
  523.  
  524.       itemp = 0;
  525.       dmax = abs(dx[0 +dx_off]);
  526.       for (i = 1; i < n; i++) {
  527.     dtemp = abs(dx[i + dx_off]);
  528.     if (dtemp > dmax) {
  529.       itemp = i;
  530.       dmax = dtemp;
  531.     }
  532.       }
  533.     }
  534.     return (itemp);
  535.   }
  536.  
  537.  
  538.   
  539.   /*
  540.     estimate unit roundoff in quantities of size x.
  541.     
  542.     this program should function properly on all systems
  543.     satisfying the following two assumptions,
  544.     1.  the base used in representing dfloating point
  545.     numbers is not a power of three.
  546.     2.  the quantity  a  in statement 10 is represented to
  547.     the accuracy used in dfloating point variables
  548.     that are stored in memory.
  549.     the statement number 10 and the go to 10 are intended to
  550.     force optimizing compilers to generate code satisfying
  551.     assumption 2.
  552.     under these assumptions, it should be true that,
  553.     a  is not exactly equal to four-thirds,
  554.     b  has a zero for its last bit or digit,
  555.     c  is not exactly equal to one,
  556.     eps  measures the separation of 1.0 from
  557.     the next larger dfloating point number.
  558.     the developers of eispack would appreciate being informed
  559.     about any systems where these assumptions do not hold.
  560.     
  561.     *****************************************************************
  562.     this routine is one of the auxiliary routines used by eispack iii
  563.     to avoid machine dependencies.
  564.     *****************************************************************
  565.   
  566.     this version dated 4/6/83.
  567.   */
  568.   final double epslon (double x)
  569.   {
  570.     double a,b,c,eps;
  571.  
  572.     a = 4.0e0/3.0e0;
  573.     eps = 0;
  574.     while (eps == 0) {
  575.       b = a - 1.0;
  576.       c = b + b + b;
  577.       eps = abs(c-1.0);
  578.     }
  579.     return(eps*abs(x));
  580.   }
  581.  
  582.   
  583.  
  584.   /*
  585.     purpose:
  586.     multiply matrix m times vector x and add the result to vector y.
  587.     
  588.     parameters:
  589.     
  590.     n1 integer, number of elements in vector y, and number of rows in
  591.     matrix m
  592.     
  593.     y double [n1], vector of length n1 to which is added
  594.     the product m*x
  595.     
  596.     n2 integer, number of elements in vector x, and number of columns
  597.     in matrix m
  598.     
  599.     ldm integer, leading dimension of array m
  600.     
  601.     x double [n2], vector of length n2
  602.     
  603.     m double [ldm][n2], matrix of n1 rows and n2 columns
  604.   */
  605.   final void dmxpy ( int n1, double y[], int n2, int ldm, double x[], double m[][])
  606.   {
  607.     int j,i;
  608.  
  609.     // cleanup odd vector
  610.     for (j = 0; j < n2; j++) {
  611.       for (i = 0; i < n1; i++) {
  612.     y[i] += x[j]*m[j][i];
  613.       }
  614.     }
  615.   }
  616. }
  617.  
  618. public class LinpackApp {
  619.  
  620.     public static void main(String args[]) {
  621.     LinpackClass l = new LinpackClass();
  622.  
  623.         l.run_benchmark();
  624.     System.out.println("Mflop/s: " + l.mflops_result + "  Time: " + l.time_result + " secs" + "  Norm Res: " + l.residn_result + "  Precision: " + l.eps_result);
  625.     }
  626. }
  627.