--- trunk/OOPSE-3.0/src/math/SquareMatrix.hpp	2004/10/21 21:31:39	1630
+++ trunk/OOPSE-3.0/src/math/SquareMatrix.hpp	2004/10/22 23:09:57	1639
@@ -199,192 +199,155 @@ namespace oopse {
     // normalized.
     template<typename Real, int Dim>
     int SquareMatrix<Real, Dim>::jacobi(SquareMatrix<Real, Dim>& a, Vector<Real, Dim>& w, 
-                                  SquareMatrix<Real, Dim>& v) {
-      const int n = Dim;  
-      int i, j, k, iq, ip, numPos;
-      Real tresh, theta, tau, t, sm, s, h, g, c, tmp;
-      Real bspace[4], zspace[4];
-      Real *b = bspace;
-      Real *z = zspace;
+                                        SquareMatrix<Real, Dim>& v) {
+        const int n = Dim;  
+        int i, j, k, iq, ip, numPos;
+        Real tresh, theta, tau, t, sm, s, h, g, c, tmp;
+        Real bspace[4], zspace[4];
+        Real *b = bspace;
+        Real *z = zspace;
 
-      // only allocate memory if the matrix is large
-      if (n > 4)
-        {
-        b = new Real[n];
-        z = new Real[n]; 
+        // only allocate memory if the matrix is large
+        if (n > 4) {
+            b = new Real[n];
+            z = new Real[n]; 
         }
 
-      // initialize
-      for (ip=0; ip<n; ip++) 
-        {
-        for (iq=0; iq<n; iq++)
-          {
-          v(ip, iq) = 0.0;
-          }
-        v(ip, ip) = 1.0;
+        // initialize
+        for (ip=0; ip<n; ip++) {
+            for (iq=0; iq<n; iq++) {
+                v(ip, iq) = 0.0;
+            }
+            v(ip, ip) = 1.0;
         }
-      for (ip=0; ip<n; ip++) 
-        {
-        b[ip] = w[ip] = a(ip, ip);
-        z[ip] = 0.0;
+        for (ip=0; ip<n; ip++) {
+            b[ip] = w[ip] = a(ip, ip);
+            z[ip] = 0.0;
         }
 
-      // begin rotation sequence
-      for (i=0; i<VTK_MAX_ROTATIONS; i++) 
-        {
-        sm = 0.0;
-        for (ip=0; ip<n-1; ip++) 
-          {
-          for (iq=ip+1; iq<n; iq++)
-            {
-            sm += fabs(a(ip, iq));
+        // begin rotation sequence
+        for (i=0; i<VTK_MAX_ROTATIONS; i++) {
+            sm = 0.0;
+            for (ip=0; ip<n-1; ip++) {
+                for (iq=ip+1; iq<n; iq++) {
+                    sm += fabs(a(ip, iq));
+                }
             }
-          }
-        if (sm == 0.0)
-          {
-          break;
-          }
-
-        if (i < 3)                                // first 3 sweeps
-          {
-          tresh = 0.2*sm/(n*n);
-          }
-        else
-          {
-          tresh = 0.0;
-          }
+            if (sm == 0.0) {
+                break;
+            }
 
-        for (ip=0; ip<n-1; ip++) 
-          {
-          for (iq=ip+1; iq<n; iq++) 
-            {
-            g = 100.0*fabs(a(ip, iq));
+            if (i < 3) {                                // first 3 sweeps
+                tresh = 0.2*sm/(n*n);
+            } else {
+                tresh = 0.0;
+            }
 
-            // after 4 sweeps
-            if (i > 3 && (fabs(w[ip])+g) == fabs(w[ip])
-            && (fabs(w[iq])+g) == fabs(w[iq]))
-              {
-              a(ip, iq) = 0.0;
-              }
-            else if (fabs(a(ip, iq)) > tresh) 
-              {
-              h = w[iq] - w[ip];
-              if ( (fabs(h)+g) == fabs(h))
-                {
-                t = (a(ip, iq)) / h;
-                }
-              else 
-                {
-                theta = 0.5*h / (a(ip, iq));
-                t = 1.0 / (fabs(theta)+sqrt(1.0+theta*theta));
-                if (theta < 0.0)
-                  {
-                  t = -t;
-                  }
-                }
-              c = 1.0 / sqrt(1+t*t);
-              s = t*c;
-              tau = s/(1.0+c);
-              h = t*a(ip, iq);
-              z[ip] -= h;
-              z[iq] += h;
-              w[ip] -= h;
-              w[iq] += h;
-              a(ip, iq)=0.0;
+            for (ip=0; ip<n-1; ip++) {
+                for (iq=ip+1; iq<n; iq++) {
+                    g = 100.0*fabs(a(ip, iq));
 
-              // ip already shifted left by 1 unit
-              for (j = 0;j <= ip-1;j++) 
-                {
-                VTK_ROTATE(a,j,ip,j,iq);
+                    // after 4 sweeps
+                    if (i > 3 && (fabs(w[ip])+g) == fabs(w[ip])
+                        && (fabs(w[iq])+g) == fabs(w[iq])) {
+                        a(ip, iq) = 0.0;
+                    } else if (fabs(a(ip, iq)) > tresh) {
+                        h = w[iq] - w[ip];
+                        if ( (fabs(h)+g) == fabs(h)) {
+                            t = (a(ip, iq)) / h;
+                        } else {
+                            theta = 0.5*h / (a(ip, iq));
+                            t = 1.0 / (fabs(theta)+sqrt(1.0+theta*theta));
+                            if (theta < 0.0) {
+                                t = -t;
+                            }
+                        }
+                        c = 1.0 / sqrt(1+t*t);
+                        s = t*c;
+                        tau = s/(1.0+c);
+                        h = t*a(ip, iq);
+                        z[ip] -= h;
+                        z[iq] += h;
+                        w[ip] -= h;
+                        w[iq] += h;
+                        a(ip, iq)=0.0;
+
+                        // ip already shifted left by 1 unit
+                        for (j = 0;j <= ip-1;j++) {
+                            VTK_ROTATE(a,j,ip,j,iq);
+                        }
+                        // ip and iq already shifted left by 1 unit
+                        for (j = ip+1;j <= iq-1;j++) {
+                            VTK_ROTATE(a,ip,j,j,iq);
+                        }
+                        // iq already shifted left by 1 unit
+                        for (j=iq+1; j<n; j++) {
+                            VTK_ROTATE(a,ip,j,iq,j);
+                        }
+                        for (j=0; j<n; j++) {
+                            VTK_ROTATE(v,j,ip,j,iq);
+                        }
+                    }
                 }
-              // ip and iq already shifted left by 1 unit
-              for (j = ip+1;j <= iq-1;j++) 
-                {
-                VTK_ROTATE(a,ip,j,j,iq);
-                }
-              // iq already shifted left by 1 unit
-              for (j=iq+1; j<n; j++) 
-                {
-                VTK_ROTATE(a,ip,j,iq,j);
-                }
-              for (j=0; j<n; j++) 
-                {
-                VTK_ROTATE(v,j,ip,j,iq);
-                }
-              }
             }
-          }
 
-        for (ip=0; ip<n; ip++) 
-          {
-          b[ip] += z[ip];
-          w[ip] = b[ip];
-          z[ip] = 0.0;
-          }
+            for (ip=0; ip<n; ip++) {
+                b[ip] += z[ip];
+                w[ip] = b[ip];
+                z[ip] = 0.0;
+            }
         }
 
-      //// this is NEVER called
-      if ( i >= VTK_MAX_ROTATIONS )
-        {
-           std::cout << "vtkMath::Jacobi: Error extracting eigenfunctions" << std::endl;
-           return 0;
+        //// this is NEVER called
+        if ( i >= VTK_MAX_ROTATIONS ) {
+            std::cout << "vtkMath::Jacobi: Error extracting eigenfunctions" << std::endl;
+            return 0;
         }
 
-      // sort eigenfunctions                 these changes do not affect accuracy 
-      for (j=0; j<n-1; j++)                  // boundary incorrect
-        {
-        k = j;
-        tmp = w[k];
-        for (i=j+1; i<n; i++)                // boundary incorrect, shifted already
-          {
-          if (w[i] >= tmp)                   // why exchage if same?
-            {
-            k = i;
+        // sort eigenfunctions                 these changes do not affect accuracy 
+        for (j=0; j<n-1; j++) {                  // boundary incorrect
+            k = j;
             tmp = w[k];
+            for (i=j+1; i<n; i++) {                // boundary incorrect, shifted already
+                if (w[i] >= tmp) {                   // why exchage if same?
+                    k = i;
+                    tmp = w[k];
+                }
             }
-          }
-        if (k != j) 
-          {
-          w[k] = w[j];
-          w[j] = tmp;
-          for (i=0; i<n; i++) 
-            {
-            tmp = v(i, j);
-            v(i, j) = v(i, k);
-            v(i, k) = tmp;
+            if (k != j) {
+                w[k] = w[j];
+                w[j] = tmp;
+                for (i=0; i<n; i++) {
+                    tmp = v(i, j);
+                    v(i, j) = v(i, k);
+                    v(i, k) = tmp;
+                }
             }
-          }
         }
-      // insure eigenvector consistency (i.e., Jacobi can compute vectors that
-      // are negative of one another (.707,.707,0) and (-.707,-.707,0). This can
-      // reek havoc in hyperstreamline/other stuff. We will select the most
-      // positive eigenvector.
-      int ceil_half_n = (n >> 1) + (n & 1);
-      for (j=0; j<n; j++)
-        {
-        for (numPos=0, i=0; i<n; i++)
-          {
-          if ( v(i, j) >= 0.0 )
-            {
-            numPos++;
+        // insure eigenvector consistency (i.e., Jacobi can compute vectors that
+        // are negative of one another (.707,.707,0) and (-.707,-.707,0). This can
+        // reek havoc in hyperstreamline/other stuff. We will select the most
+        // positive eigenvector.
+        int ceil_half_n = (n >> 1) + (n & 1);
+        for (j=0; j<n; j++) {
+            for (numPos=0, i=0; i<n; i++) {
+                if ( v(i, j) >= 0.0 ) {
+                    numPos++;
+                }
             }
-          }
-    //    if ( numPos < ceil(double(n)/double(2.0)) )
-        if ( numPos < ceil_half_n)
-          {
-          for(i=0; i<n; i++)
-            {
-            v(i, j) *= -1.0;
+            //    if ( numPos < ceil(double(n)/double(2.0)) )
+            if ( numPos < ceil_half_n) {
+                for (i=0; i<n; i++) {
+                    v(i, j) *= -1.0;
+                }
             }
-          }
         }
 
-      if (n > 4)
-        {
-        delete [] b;
-        delete [] z;
+        if (n > 4) {
+            delete [] b;
+            delete [] z;
         }
-      return 1;
+        return 1;
     }