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
        @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
 
 
 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 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() ;
                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() ;
                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()
         * 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.
        
        /**
         * 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 ) 
        {
         */
        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.
        
        /**
         * Add a task into the graph.
-        * @param t Task to be add
+        * @param _t Task to be add
         */
        public void addGTask( GTask _t )
        {
         */
        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;
 
     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);
            "  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
 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 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) .
     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
     on entry
-    
+
     a       double precision[n][lda]
     the matrix to be factored.
     a       double precision[n][lda]
     the matrix to be factored.
-    
+
     lda     integer
     the leading dimension of the array  a .
     lda     integer
     the leading dimension of the array  a .
-    
+
     n       integer
     the order of the matrix  a .
     n       integer
     the order of the matrix  a .
-    
+
     on return
     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.
     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.
     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
     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.
     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.
     linpack. this version dated 08/14/78.
     cleve moler, university of new mexico, argonne national lab.
-    
+
     functions
     functions
-    
+
     blas daxpy,dscal,idamax
     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.
     dgesl solves the double precision system
     a * x = b  or  trans(a) * x = b
     using the factors computed by dgeco or dgefa.
-  
+
     on entry
     on entry
-  
+
     a       double precision[n][lda]
     the output from dgeco or dgefa.
     a       double precision[n][lda]
     the output from dgeco or dgefa.
-  
+
     lda     integer
     the leading dimension of the array  a .
     lda     integer
     the leading dimension of the array  a .
-    
+
     n       integer
     the order of the matrix  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.
     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.
     job     integer
     = 0         to solve  a*x = b ,
     = nonzero   to solve  trans(a)*x = b  where
     trans(a)  is the transpose.
-    
+
     on return
     on return
-    
+
     b       the solution vector  x .
     b       the solution vector  x .
-    
+
     error condition
     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 .
     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)
     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);
     }
     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.
     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.
     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.
     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.
     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.
     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.
     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
     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.
     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 routine is one of the auxiliary routines used by eispack iii
     to avoid machine dependencies.
-    *****************************************************************
-  
+        *****************************************************************
+
     this version dated 4/6/83.
     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.
     purpose:
     multiply matrix m times vector x and add the result to vector y.
-    
+
     parameters:
     parameters:
-    
+
     n1 integer, number of elements in vector y, and number of rows in
     matrix m
     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
     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
     n2 integer, number of elements in vector x, and number of columns
     in matrix m
-    
+
     ldm integer, leading dimension of array m
     ldm integer, leading dimension of array m
-    
+
     x double [n2], vector of length n2
     x double [n2], vector of length n2
-    
+
     m double [ldm][n2], matrix of n1 rows and n2 columns
     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];
+                       }
+               }
+       }
 
 }
 
 }