Logo AND Algorithmique Numérique Distribuée

Public GIT Repository
Correction of some bugs.
authorSébastien Miquée <sebastien.miquee@univ-fcomte.fr>
Fri, 26 Feb 2010 07:09:36 +0000 (08:09 +0100)
committerSébastien Miquée <sebastien.miquee@univ-fcomte.fr>
Fri, 26 Feb 2010 07:09:36 +0000 (08:09 +0100)
Makefile
src/and/Mapping/Algo.java
src/and/Mapping/GNode.java
src/and/Mapping/Graph.java
src/and/Mapping/Linpack.java

index 292a7c7..5e7081a 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -23,7 +23,7 @@ jar:clean compile
        @echo
        @echo "## Creation of Mapping jar ..."
        @echo
-       jar -cvfm ./$(JAR) Manifest ./$(EXT) ./$(JAVADOC) -C ./$(BIN) ./$(PACK)/
+       jar -cvfm ./$(JAR) Manifest ./$(EXT) -C ./$(BIN) ./$(PACK)/
 
 
 javadoc:cleanDoc
index 608c667..142b8a6 100644 (file)
@@ -16,6 +16,7 @@ public abstract class Algo implements Serializable
        protected Graph gr ;
        protected Grid gl ;
        protected Mapping mp ;
+       protected String ids ;
        
        
        /**
@@ -26,6 +27,7 @@ public abstract class Algo implements Serializable
                gr = new Graph() ;
                gl = new Grid() ;
                mp = new Mapping() ;
+               ids = "" ;
        }
        
        
@@ -39,6 +41,7 @@ public abstract class Algo implements Serializable
                gr =  _gr ;
                gl = _gl ;
                mp = new Mapping() ;
+               ids = "" ;
        }
        
        
@@ -97,6 +100,26 @@ public abstract class Algo implements Serializable
        
        
        /**
+        * Set the string identifier for the algorithm.
+        * @param _s The algorithm's identifier
+        */
+       public void setIdS( String _s )
+       {
+               ids = _s ;
+       }
+       
+       
+       /**
+        * Return the string identifier of the algorithm.
+        * @return The algorithm's identifier
+        */
+       public String getIdS()
+       {
+               return ids ;
+       }
+       
+       
+       /**
         * Update the grid status after having done the mapping.
         */
        public void updateGrid()
index c1015db..b76d0c1 100644 (file)
@@ -210,7 +210,7 @@ public class GNode implements Serializable
        
        /**
         * Set the MFlops of each computing core of the computing node.
-        * @param _freq The MFlops of cores
+        * @param _mflops The MFlops of cores
         */
        public void setMFlops( int _mflops ) 
        {
index 3387f60..79efb2f 100644 (file)
@@ -45,7 +45,7 @@ public class Graph implements Serializable
        
        /**
         * Add a task into the graph.
-        * @param t Task to be add
+        * @param _t Task to be add
         */
        public void addGTask( GTask _t )
        {
index 08978b4..b92d7d3 100644 (file)
@@ -32,71 +32,88 @@ Translated to C by Bonnie Toy 5/88
 
 */
 
+/**
+ * Class allowing to retrieve the computing power of a 
+ * computer; this power is represented by MFlops
+ * 
+ * @author Bonnie Toy, Reed Wade
+ */
+
+public class Linpack 
+{
+
+       /**
+        * Empty default constructor.
+        */
+       public Linpack(){}
+
+       final double abs (double d) {
+               return (d >= 0) ? d : -d;
+       }
 
-public class Linpack {
+       double second_orig = -1;
 
-  final double abs (double d) {
-    return (d >= 0) ? d : -d;
-  }
+       double second()
+       {
+               if (second_orig==-1) {
+                       second_orig = System.currentTimeMillis();
+               }
+               return (System.currentTimeMillis() - second_orig)/1000;
+       }
 
-  double second_orig = -1;
-
-  double second()
-  {
-    if (second_orig==-1) {
-      second_orig = System.currentTimeMillis();
-    }
-    return (System.currentTimeMillis() - second_orig)/1000;
-  }
-
-  public double getMFlops() 
-  {
-    double mflops_result = 0.0;
-    double residn_result = 0.0;
-    double time_result = 0.0;
-    double eps_result = 0.0;
-
-    double a[][] = new double[1000][1001];
-    double b[] = new double[1000];
-    double x[] = new double[1000];
-    double ops,total,norma,normx;
-    double resid,time;
-    int n,i,lda;
-    int ipvt[] = new int[1000];
-    
-    //double mflops_result;
-    //double residn_result;
-    //double time_result;
-    //double eps_result;
-
-    lda = 1001;
-    n = 500;
-    
-    ops = (2.0e0*(n*n*n))/3.0 + 2.0*(n*n);
-    
-    norma = matgen(a,lda,n,b);
-    time = second();
-    dgefa(a,lda,n,ipvt);
-    dgesl(a,lda,n,ipvt,b,0);
-    total = second() - time;
-    
-    for (i = 0; i < n; i++) {
-      x[i] = b[i];
-    }
-    norma = matgen(a,lda,n,b);
-    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 > abs(b[i])) ? resid : abs(b[i]);
-      normx = (normx > abs(x[i])) ? normx : abs(x[i]);
-    }
-    
-    eps_result = epslon((double)1.0);
-/*
+       /**
+        * Estimate the computing power of the local computer
+        * and return as result the obtained MFlops.
+        * @return The MFlops of the computer
+        */
+       public double getMFlops() 
+       {
+               double mflops_result = 0.0;
+               double residn_result = 0.0;
+               double time_result = 0.0;
+               double eps_result = 0.0;
+
+               double a[][] = new double[1000][1001];
+               double b[] = new double[1000];
+               double x[] = new double[1000];
+               double ops,total,norma,normx;
+               double resid,time;
+               int n,i,lda;
+               int ipvt[] = new int[1000];
+
+               //double mflops_result;
+               //double residn_result;
+               //double time_result;
+               //double eps_result;
+
+               lda = 1001;
+               n = 500;
+
+               ops = (2.0e0*(n*n*n))/3.0 + 2.0*(n*n);
+
+               norma = matgen(a,lda,n,b);
+               time = second();
+               dgefa(a,lda,n,ipvt);
+               dgesl(a,lda,n,ipvt,b,0);
+               total = second() - time;
+
+               for (i = 0; i < n; i++) {
+                       x[i] = b[i];
+               }
+               norma = matgen(a,lda,n,b);
+               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 > abs(b[i])) ? resid : abs(b[i]);
+                       normx = (normx > abs(x[i])) ? normx : abs(x[i]);
+               }
+
+               eps_result = epslon((double)1.0);
+               /*
 
     residn_result = resid/( n*norma*normx*eps_result );
     time_result = total;
@@ -106,92 +123,92 @@ public class Linpack {
            "  Time: " + time_result + " secs" +
            "  Norm Res: " + residn_result +
            "  Precision: " + eps_result);
-*/
-       residn_result = resid/( n*norma*normx*eps_result );
-       residn_result += 0.005; // for rounding
-       residn_result = (int)(residn_result*100);
-       residn_result /= 100;
-
-       time_result = total;
-       time_result += 0.005; // for rounding
-       time_result = (int)(time_result*100);
-       time_result /= 100;
-
-       mflops_result = ops/(1.0e6*total);
-       mflops_result += 0.0005; // for rounding
-       mflops_result = (int)(mflops_result*1000);
-       mflops_result /= 1000;
-
-//    System.out.println("Mflops/s: " + mflops_result +
-//         "  Time: " + time_result + " secs" +
-//         "  Norm Res: " + residn_result +
-//         "  Precision: " + eps_result);
-       
-       return mflops_result ;
-  }
-  
-
-  
-  final double matgen (double a[][], int lda, int n, double b[])
-  {
-    double norma;
-    int init, i, j;
-    
-    init = 1325;
-    norma = 0.0;
-/*  Next two for() statements switched.  Solver wants
+                */
+               residn_result = resid/( n*norma*normx*eps_result );
+               residn_result += 0.005; // for rounding
+               residn_result = (int)(residn_result*100);
+               residn_result /= 100;
+
+               time_result = total;
+               time_result += 0.005; // for rounding
+               time_result = (int)(time_result*100);
+               time_result /= 100;
+
+               mflops_result = ops/(1.0e6*total);
+               mflops_result += 0.0005; // for rounding
+               mflops_result = (int)(mflops_result*1000);
+               mflops_result /= 1000;
+
+               //    System.out.println("Mflops/s: " + mflops_result +
+               //          "  Time: " + time_result + " secs" +
+               //          "  Norm Res: " + residn_result +
+               //          "  Precision: " + eps_result);
+
+               return mflops_result ;
+       }
+
+
+
+       final double matgen (double a[][], int lda, int n, double b[])
+       {
+               double norma;
+               int init, i, j;
+
+               init = 1325;
+               norma = 0.0;
+               /*  Next two for() statements switched.  Solver wants
 matrix in column order. --dmd 3/3/97
-*/
-      for (i = 0; i < n; i++) {
-    for (j = 0; j < n; j++) {
-       init = 3125*init % 65536;
-       a[j][i] = (init - 32768.0)/16384.0;
-       norma = (a[j][i] > norma) ? a[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] += a[j][i];
-      }
-    }
-    
-    return norma;
-  }
-  
+                */
+               for (i = 0; i < n; i++) {
+                       for (j = 0; j < n; j++) {
+                               init = 3125*init % 65536;
+                               a[j][i] = (init - 32768.0)/16384.0;
+                               norma = (a[j][i] > norma) ? a[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] += a[j][i];
+                       }
+               }
+
+               return norma;
+       }
+
+
 
-  
-  /*
+       /*
     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       double 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
@@ -199,118 +216,118 @@ matrix in column order. --dmd 3/3/97
     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
-  */
-  final int dgefa( double a[][], int lda, int n, int ipvt[])
-  {
-    double[] col_k, col_j;
-    double t;
-    int j,k,kp1,l,nm1;
-    int info;
-    
-    // gaussian elimination with partial pivoting
-    
-    info = 0;
-    nm1 = n - 1;
-    if (nm1 >=  0) {
-      for (k = 0; k < nm1; k++) {
-       col_k = a[k];
-       kp1 = k + 1;
-       
-       // find l = pivot index
-       
-       l = idamax(n-k,col_k,k,1) + k;
-       ipvt[k] = l;
-       
-       // zero pivot implies this column already triangularized
-       
-       if (col_k[l] != 0) {
-         
-         // interchange if necessary
-         
-         if (l != k) {
-           t = col_k[l];
-           col_k[l] = col_k[k];
-           col_k[k] = t;
-         }
-         
-         // compute multipliers
-         
-         t = -1.0/col_k[k];
-         dscal(n-(kp1),t,col_k,kp1,1);
-         
-         // row elimination with column indexing
-         
-         for (j = kp1; j < n; j++) {
-           col_j = a[j];
-           t = col_j[l];
-           if (l != k) {
-             col_j[l] = col_j[k];
-             col_j[k] = t;
-           }
-           daxpy(n-(kp1),t,col_k,kp1,1,
-                 col_j,kp1,1);
-         }
-       }
-       else {
-         info = k;
+        */
+       final int dgefa( double a[][], int lda, int n, int ipvt[])
+       {
+               double[] col_k, col_j;
+               double t;
+               int j,k,kp1,l,nm1;
+               int info;
+
+               // gaussian elimination with partial pivoting
+
+               info = 0;
+               nm1 = n - 1;
+               if (nm1 >=  0) {
+                       for (k = 0; k < nm1; k++) {
+                               col_k = a[k];
+                               kp1 = k + 1;
+
+                               // find l = pivot index
+
+                               l = idamax(n-k,col_k,k,1) + k;
+                               ipvt[k] = l;
+
+                               // zero pivot implies this column already triangularized
+
+                               if (col_k[l] != 0) {
+
+                                       // interchange if necessary
+
+                                       if (l != k) {
+                                               t = col_k[l];
+                                               col_k[l] = col_k[k];
+                                               col_k[k] = t;
+                                       }
+
+                                       // compute multipliers
+
+                                       t = -1.0/col_k[k];
+                                       dscal(n-(kp1),t,col_k,kp1,1);
+
+                                       // row elimination with column indexing
+
+                                       for (j = kp1; j < n; j++) {
+                                               col_j = a[j];
+                                               t = col_j[l];
+                                               if (l != k) {
+                                                       col_j[l] = col_j[k];
+                                                       col_j[k] = t;
+                                               }
+                                               daxpy(n-(kp1),t,col_k,kp1,1,
+                                                               col_j,kp1,1);
+                                       }
+                               }
+                               else {
+                                       info = k;
+                               }
+                       }
+               }
+               ipvt[n-1] = n-1;
+               if (a[(n-1)][(n-1)] == 0) info = n-1;
+
+               return info;
        }
-      }
-    }
-    ipvt[n-1] = n-1;
-    if (a[(n-1)][(n-1)] == 0) info = n-1;
-    
-    return info;
-  }
-
-  
-  
-  /*
+
+
+
+       /*
     dgesl solves the double precision system
     a * x = b  or  trans(a) * x = b
     using the factors computed by dgeco or dgefa.
-  
+
     on entry
-  
+
     a       double precision[n][lda]
     the output from dgeco or dgefa.
-  
+
     lda     integer
     the leading dimension of the array  a .
-    
+
     n       integer
     the order of the matrix  a .
-  
+
     ipvt    integer[n]
     the pivot vector from dgeco or dgefa.
 
     b       double precision[n]
     the right hand side vector.
-    
+
     job     integer
     = 0         to solve  a*x = b ,
     = nonzero   to solve  trans(a)*x = b  where
     trans(a)  is the transpose.
-    
+
     on return
-    
+
     b       the solution vector  x .
-    
+
     error condition
-    
+
     a division by zero will occur if the input factor contains a
     zero on the diagonal.  technically this indicates singularity
     but it is often caused by improper arguments or improper
     setting of lda .  it will not occur if the subroutines are
     called correctly and if dgeco has set rcond .gt. 0.0
     or dgefa has set info .eq. 0 .
-    
+
     to compute  inverse(a) * c  where  c  is a matrix
     with  p  columns
     dgeco(a,lda,n,ipvt,rcond,z)
@@ -318,228 +335,228 @@ matrix in column order. --dmd 3/3/97
     for (j=0,j<p,j++)
     dgesl(a,lda,n,ipvt,c[j][0],0);
     }
-    
+
     linpack. this version dated 08/14/78 .
     cleve moler, university of new mexico, argonne national lab.
-    
-    functions
-    
-    blas daxpy,ddot
-  */
-  final void dgesl( double a[][], int lda, int n, int ipvt[], double b[], int job)
-  {
-    double t;
-    int k,kb,l,nm1,kp1;
-
-    nm1 = n - 1;
-    if (job == 0) {
-
-      // job = 0 , solve  a * x = b.  first solve  l*y = b
-
-      if (nm1 >= 1) {
-       for (k = 0; k < nm1; k++) {
-         l = ipvt[k];
-         t = b[l];
-         if (l != k){
-           b[l] = b[k];
-           b[k] = t;
-         }
-         kp1 = k + 1;
-         daxpy(n-(kp1),t,a[k],kp1,1,b,kp1,1);
-       }
-      }
 
-      // now solve  u*x = y
+    functions
 
-      for (kb = 0; kb < n; kb++) {
-       k = n - (kb + 1);
-       b[k] /= a[k][k];
-       t = -b[k];
-       daxpy(k,t,a[k],0,1,b,0,1);
-      }
-    }
-    else {
-
-      // job = nonzero, solve  trans(a) * x = b.  first solve  trans(u)*y = b
-
-      for (k = 0; k < n; k++) {
-       t = ddot(k,a[k],0,1,b,0,1);
-       b[k] = (b[k] - t)/a[k][k];
-      }
-
-      // now solve trans(l)*x = y 
-
-      if (nm1 >= 1) {
-       for (kb = 1; kb < nm1; kb++) {
-         k = n - (kb+1);
-         kp1 = k + 1;
-         b[k] += ddot(n-(kp1),a[k],kp1,1,b,kp1,1);
-         l = ipvt[k];
-         if (l != k) {
-           t = b[l];
-           b[l] = b[k];
-           b[k] = t;
-         }
+    blas daxpy,ddot
+        */
+       final void dgesl( double a[][], int lda, int n, int ipvt[], double b[], int job)
+       {
+               double t;
+               int k,kb,l,nm1,kp1;
+
+               nm1 = n - 1;
+               if (job == 0) {
+
+                       // job = 0 , solve  a * x = b.  first solve  l*y = b
+
+                       if (nm1 >= 1) {
+                               for (k = 0; k < nm1; k++) {
+                                       l = ipvt[k];
+                                       t = b[l];
+                                       if (l != k){
+                                               b[l] = b[k];
+                                               b[k] = t;
+                                       }
+                                       kp1 = k + 1;
+                                       daxpy(n-(kp1),t,a[k],kp1,1,b,kp1,1);
+                               }
+                       }
+
+                       // now solve  u*x = y
+
+                       for (kb = 0; kb < n; kb++) {
+                               k = n - (kb + 1);
+                               b[k] /= a[k][k];
+                               t = -b[k];
+                               daxpy(k,t,a[k],0,1,b,0,1);
+                       }
+               }
+               else {
+
+                       // job = nonzero, solve  trans(a) * x = b.  first solve  trans(u)*y = b
+
+                       for (k = 0; k < n; k++) {
+                               t = ddot(k,a[k],0,1,b,0,1);
+                               b[k] = (b[k] - t)/a[k][k];
+                       }
+
+                       // now solve trans(l)*x = y 
+
+                       if (nm1 >= 1) {
+                               for (kb = 1; kb < nm1; kb++) {
+                                       k = n - (kb+1);
+                                       kp1 = k + 1;
+                                       b[k] += ddot(n-(kp1),a[k],kp1,1,b,kp1,1);
+                                       l = ipvt[k];
+                                       if (l != k) {
+                                               t = b[l];
+                                               b[l] = b[k];
+                                               b[k] = t;
+                                       }
+                               }
+                       }
+               }
        }
-      }
-    }
-  }
 
 
 
-  /*
+       /*
     constant times a vector plus a vector.
     jack dongarra, linpack, 3/11/78.
-  */
-  final void daxpy( int n, double da, double dx[], int dx_off, int incx,
-             double dy[], int dy_off, int incy)
-  {
-    int i,ix,iy;
-
-    if ((n > 0) && (da != 0)) {
-      if (incx != 1 || incy != 1) {
-
-       // code for unequal increments or equal increments not equal to 1
-
-       ix = 0;
-       iy = 0;
-       if (incx < 0) ix = (-n+1)*incx;
-       if (incy < 0) iy = (-n+1)*incy;
-       for (i = 0;i < n; i++) {
-         dy[iy +dy_off] += da*dx[ix +dx_off];
-         ix += incx;
-         iy += incy;
+        */
+       final void daxpy( int n, double da, double dx[], int dx_off, int incx,
+                       double dy[], int dy_off, int incy)
+       {
+               int i,ix,iy;
+
+               if ((n > 0) && (da != 0)) {
+                       if (incx != 1 || incy != 1) {
+
+                               // code for unequal increments or equal increments not equal to 1
+
+                               ix = 0;
+                               iy = 0;
+                               if (incx < 0) ix = (-n+1)*incx;
+                               if (incy < 0) iy = (-n+1)*incy;
+                               for (i = 0;i < n; i++) {
+                                       dy[iy +dy_off] += da*dx[ix +dx_off];
+                                       ix += incx;
+                                       iy += incy;
+                               }
+                               return;
+                       } else {
+
+                               // code for both increments equal to 1
+
+                               for (i=0; i < n; i++)
+                                       dy[i +dy_off] += da*dx[i +dx_off];
+                       }
+               }
        }
-       return;
-      } else {
-
-       // code for both increments equal to 1
-
-       for (i=0; i < n; i++)
-         dy[i +dy_off] += da*dx[i +dx_off];
-      }
-    }
-  }
 
 
 
-  /*
+       /*
     forms the dot product of two vectors.
     jack dongarra, linpack, 3/11/78.
-  */
-  final double ddot( int n, double dx[], int dx_off, int incx, double dy[],
-              int dy_off, int incy)
-  {
-    double dtemp;
-    int i,ix,iy;
-
-    dtemp = 0;
-
-    if (n > 0) {
-      
-      if (incx != 1 || incy != 1) {
-
-       // code for unequal increments or equal increments not equal to 1
-
-       ix = 0;
-       iy = 0;
-       if (incx < 0) ix = (-n+1)*incx;
-       if (incy < 0) iy = (-n+1)*incy;
-       for (i = 0;i < n; i++) {
-         dtemp += dx[ix +dx_off]*dy[iy +dy_off];
-         ix += incx;
-         iy += incy;
+        */
+       final double ddot( int n, double dx[], int dx_off, int incx, double dy[],
+                       int dy_off, int incy)
+       {
+               double dtemp;
+               int i,ix,iy;
+
+               dtemp = 0;
+
+               if (n > 0) {
+
+                       if (incx != 1 || incy != 1) {
+
+                               // code for unequal increments or equal increments not equal to 1
+
+                               ix = 0;
+                               iy = 0;
+                               if (incx < 0) ix = (-n+1)*incx;
+                               if (incy < 0) iy = (-n+1)*incy;
+                               for (i = 0;i < n; i++) {
+                                       dtemp += dx[ix +dx_off]*dy[iy +dy_off];
+                                       ix += incx;
+                                       iy += incy;
+                               }
+                       } else {
+
+                               // code for both increments equal to 1
+
+                               for (i=0;i < n; i++)
+                                       dtemp += dx[i +dx_off]*dy[i +dy_off];
+                       }
+               }
+               return(dtemp);
        }
-      } else {
 
-       // code for both increments equal to 1
-       
-       for (i=0;i < n; i++)
-         dtemp += dx[i +dx_off]*dy[i +dy_off];
-      }
-    }
-    return(dtemp);
-  }
 
-  
-  
-  /*
+
+       /*
     scales a vector by a constant.
     jack dongarra, linpack, 3/11/78.
-  */
-  final void dscal( int n, double da, double dx[], int dx_off, int incx)
-  {
-    int i,nincx;
+        */
+       final void dscal( int n, double da, double dx[], int dx_off, int incx)
+       {
+               int i,nincx;
 
-    if (n > 0) {
-      if (incx != 1) {
+               if (n > 0) {
+                       if (incx != 1) {
 
-       // code for increment not equal to 1
+                               // code for increment not equal to 1
 
-       nincx = n*incx;
-       for (i = 0; i < nincx; i += incx)
-         dx[i +dx_off] *= da;
-      } else {
+                               nincx = n*incx;
+                               for (i = 0; i < nincx; i += incx)
+                                       dx[i +dx_off] *= da;
+                       } else {
+
+                               // code for increment equal to 1
+
+                               for (i = 0; i < n; i++)
+                                       dx[i +dx_off] *= da;
+                       }
+               }
+       }
 
-       // code for increment equal to 1
 
-       for (i = 0; i < n; i++)
-         dx[i +dx_off] *= da;
-      }
-    }
-  }
 
-  
-  
-  /*
+       /*
     finds the index of element having max. absolute value.
     jack dongarra, linpack, 3/11/78.
-  */
-  final int idamax( int n, double dx[], int dx_off, int incx)
-  {
-    double dmax, dtemp;
-    int i, ix, itemp=0;
-
-    if (n < 1) {
-      itemp = -1;
-    } else if (n ==1) {
-      itemp = 0;
-    } else if (incx != 1) {
-
-      // code for increment not equal to 1
-
-      dmax = abs(dx[0 +dx_off]);
-      ix = 1 + incx;
-      for (i = 1; i < n; i++) {
-       dtemp = abs(dx[ix + dx_off]);
-       if (dtemp > dmax)  {
-         itemp = i;
-         dmax = dtemp;
+        */
+       final int idamax( int n, double dx[], int dx_off, int incx)
+       {
+               double dmax, dtemp;
+               int i, ix, itemp=0;
+
+               if (n < 1) {
+                       itemp = -1;
+               } else if (n ==1) {
+                       itemp = 0;
+               } else if (incx != 1) {
+
+                       // code for increment not equal to 1
+
+                       dmax = abs(dx[0 +dx_off]);
+                       ix = 1 + incx;
+                       for (i = 1; i < n; i++) {
+                               dtemp = abs(dx[ix + dx_off]);
+                               if (dtemp > dmax)  {
+                                       itemp = i;
+                                       dmax = dtemp;
+                               }
+                               ix += incx;
+                       }
+               } else {
+
+                       // code for increment equal to 1
+
+                       itemp = 0;
+                       dmax = abs(dx[0 +dx_off]);
+                       for (i = 1; i < n; i++) {
+                               dtemp = abs(dx[i + dx_off]);
+                               if (dtemp > dmax) {
+                                       itemp = i;
+                                       dmax = dtemp;
+                               }
+                       }
+               }
+               return (itemp);
        }
-       ix += incx;
-      }
-    } else {
-
-      // code for increment equal to 1
-
-      itemp = 0;
-      dmax = abs(dx[0 +dx_off]);
-      for (i = 1; i < n; i++) {
-       dtemp = abs(dx[i + dx_off]);
-       if (dtemp > dmax) {
-         itemp = i;
-         dmax = dtemp;
-       }
-      }
-    }
-    return (itemp);
-  }
 
 
-  
-  /*
+
+       /*
     estimate unit roundoff in quantities of size x.
-    
+
     this program should function properly on all systems
     satisfying the following two assumptions,
     1.  the base used in representing dfloating point
@@ -558,61 +575,61 @@ matrix in column order. --dmd 3/3/97
     the next larger dfloating point number.
     the developers of eispack would appreciate being informed
     about any systems where these assumptions do not hold.
-    
-    *****************************************************************
+
+        *****************************************************************
     this routine is one of the auxiliary routines used by eispack iii
     to avoid machine dependencies.
-    *****************************************************************
-  
+        *****************************************************************
+
     this version dated 4/6/83.
-  */
-  final double epslon (double x)
-  {
-    double a,b,c,eps;
-
-    a = 4.0e0/3.0e0;
-    eps = 0;
-    while (eps == 0) {
-      b = a - 1.0;
-      c = b + b + b;
-      eps = abs(c-1.0);
-    }
-    return(eps*abs(x));
-  }
+        */
+       final double epslon (double x)
+       {
+               double a,b,c,eps;
+
+               a = 4.0e0/3.0e0;
+               eps = 0;
+               while (eps == 0) {
+                       b = a - 1.0;
+                       c = b + b + b;
+                       eps = abs(c-1.0);
+               }
+               return(eps*abs(x));
+       }
 
-  
 
-  /*
+
+       /*
     purpose:
     multiply matrix m times vector x and add the result to vector y.
-    
+
     parameters:
-    
+
     n1 integer, number of elements in vector y, and number of rows in
     matrix m
-    
+
     y double [n1], vector of length n1 to which is added
     the product m*x
-    
+
     n2 integer, number of elements in vector x, and number of columns
     in matrix m
-    
+
     ldm integer, leading dimension of array m
-    
+
     x double [n2], vector of length n2
-    
+
     m double [ldm][n2], matrix of n1 rows and n2 columns
-  */
-  final void dmxpy ( int n1, double y[], int n2, int ldm, double x[], double m[][])
-  {
-    int j,i;
-
-    // cleanup odd vector
-    for (j = 0; j < n2; j++) {
-      for (i = 0; i < n1; i++) {
-       y[i] += x[j]*m[j][i];
-      }
-    }
-  }
+        */
+       final void dmxpy ( int n1, double y[], int n2, int ldm, double x[], double m[][])
+       {
+               int j,i;
+
+               // cleanup odd vector
+               for (j = 0; j < n2; j++) {
+                       for (i = 0; i < n1; i++) {
+                               y[i] += x[j]*m[j][i];
+                       }
+               }
+       }
 
 }