Corrected looping in IMTQL2. 6 of the 7 examples give good results. Example
authorilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Thu, 7 Dec 2017 22:13:56 +0000 (22:13 +0000)
committerilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Thu, 7 Dec 2017 22:13:56 +0000 (22:13 +0000)
#3, A Heavily Crowded Region gives INFO[0] = 4, unsuccessful exit.

git-svn-id: https://citdcbmipav.cit.nih.gov/repos-pub/mipav/trunk@15297 ba61647d-9d00-f842-95cd-605cb4296b96

mipav/src/gov/nih/mipav/model/algorithms/DoublyConnectedSC.java

index d47c525..db7304b 100644 (file)
@@ -3551,7 +3551,6 @@ double neweps;
                double MUZERO[] = new double[1];\r
                int I;\r
                int IERR[] = new int[1];\r
-               \r
                CLASS(N, ALPHA, BETA, B, T, MUZERO);\r
         W[0] = 1.0;\r
         for (I = 1; I < N; I++) {\r
@@ -3620,87 +3619,82 @@ double neweps;
                double B, C, F, G, P, R, S;\r
                int I, II, K, L, M, MML;\r
                int J = 0;\r
-               boolean zeroJ = true;\r
                IERR[0] = 0;\r
                if (N == 1) {\r
                        return;\r
                }\r
                \r
                E[N-1] = 0.0;\r
-               for (L = 1; L <= N; L++) {\r
-                       if (zeroJ) {\r
-                       J = 0;\r
-                       }\r
-                       else {\r
-                               zeroJ = true;\r
-                       }\r
+               outer: for (L = 1; L <= N; L++) {\r
+                   J = 0;\r
                    // Look for small sub-diagonal element\r
-                   for (M = L; M <= N; M++) {\r
-                       if (M == N) {\r
-                               break;\r
-                       }\r
-                       if (Math.abs(E[M-1]) <= MACHEP * (Math.abs(D[M-1] + Math.abs(D[M])))) {\r
-                               break;\r
-                       }\r
-                   } // for (M = L; M <= N; M++)\r
-                   \r
-                   P = D[L-1];\r
-                   if (M == L) {\r
-                       continue;\r
-                   }\r
-                   if (J == 30) {\r
-                       // Set error -- no convergence to an eigenvalue after 30 iterations.\r
-                       IERR[0] = L;\r
-                       return;\r
-                   }\r
-                   J = J+1;\r
-                   // Form shift\r
-                   G = (D[L] - P)/(2.0*E[L-1]);\r
-                   R = Math.sqrt(G*G + 1.0);\r
-                   if (G >= 0) {\r
-                       G = D[M-1] - P + E[L-1]/(G + Math.abs(R));\r
-                   }\r
-                   else {\r
-                       G = D[M-1] - P + E[L-1]/(G - Math.abs(R));      \r
-                   }\r
-                   S = 1.0;\r
-                   C = 1.0;\r
-                   P = 0.0;\r
-                   MML = M - L;\r
-                   // For I=M-1 step -1 until L do\r
-                   for (II = 1; II <= MML; II++) {\r
-                       I = M - II;\r
-                       F = S * E[I-1];\r
-                       B = C * E[I-1];\r
-                       if (Math.abs(F) >= Math.abs(G)) {\r
-                           C = G/F;\r
-                           R = Math.sqrt(C*C + 1.0);\r
-                           E[I] = F*R;\r
-                           S = 1.0/R;\r
-                           C = C*S;\r
-                       }\r
-                       else {\r
-                           S= F/G;\r
-                           R = Math.sqrt(S*S + 1.0);\r
-                           E[I] = G*R;\r
-                           C = 1.0/R;\r
-                           S = S * C;\r
-                       }\r
-                       G = D[I] - P;\r
-                       R = (D[I-1] - G)*S + 2.0*C*B;\r
-                       P = S*R;\r
-                       D[I] = G + P;\r
-                       G = C*R - B;\r
-                       // Form first component of vector\r
-                       F = Z[I];\r
-                       Z[I] = S * Z[I-1] + C * F;\r
-                       Z[I-1] = C * Z[I-1] - S * F;\r
-                   } // for (II = 1; II <= MML; II++)\r
-                   \r
-                   D[L-1] = D[L-1] - P;\r
-                   E[L-1] = G;\r
-                   E[M-1] = 0.0;\r
-                   zeroJ = false;\r
+                   while (true) {\r
+                           for (M = L; M <= N; M++) {\r
+                               if (M == N) {\r
+                                       break;\r
+                               }\r
+                               if (Math.abs(E[M-1]) <= MACHEP * (Math.abs(D[M-1] + Math.abs(D[M])))) {\r
+                                       break;\r
+                               }\r
+                           } // for (M = L; M <= N; M++)\r
+                           \r
+                           P = D[L-1];\r
+                           if (M == L) {\r
+                               continue outer;\r
+                           }\r
+                           if (J == 30) {\r
+                               // Set error -- no convergence to an eigenvalue after 30 iterations.\r
+                               IERR[0] = L;\r
+                               return;\r
+                           }\r
+                           J = J+1;\r
+                           // Form shift\r
+                           G = (D[L] - P)/(2.0*E[L-1]);\r
+                           R = Math.sqrt(G*G + 1.0);\r
+                           if (G >= 0) {\r
+                               G = D[M-1] - P + E[L-1]/(G + Math.abs(R));\r
+                           }\r
+                           else {\r
+                               G = D[M-1] - P + E[L-1]/(G - Math.abs(R));      \r
+                           }\r
+                           S = 1.0;\r
+                           C = 1.0;\r
+                           P = 0.0;\r
+                           MML = M - L;\r
+                           // For I=M-1 step -1 until L do\r
+                           for (II = 1; II <= MML; II++) {\r
+                               I = M - II;\r
+                               F = S * E[I-1];\r
+                               B = C * E[I-1];\r
+                               if (Math.abs(F) >= Math.abs(G)) {\r
+                                   C = G/F;\r
+                                   R = Math.sqrt(C*C + 1.0);\r
+                                   E[I] = F*R;\r
+                                   S = 1.0/R;\r
+                                   C = C*S;\r
+                               }\r
+                               else {\r
+                                   S= F/G;\r
+                                   R = Math.sqrt(S*S + 1.0);\r
+                                   E[I] = G*R;\r
+                                   C = 1.0/R;\r
+                                   S = S * C;\r
+                               }\r
+                               G = D[I] - P;\r
+                               R = (D[I-1] - G)*S + 2.0*C*B;\r
+                               P = S*R;\r
+                               D[I] = G + P;\r
+                               G = C*R - B;\r
+                               // Form first component of vector\r
+                               F = Z[I];\r
+                               Z[I] = S * Z[I-1] + C * F;\r
+                               Z[I-1] = C * Z[I-1] - S * F;\r
+                           } // for (II = 1; II <= MML; II++)\r
+                           \r
+                           D[L-1] = D[L-1] - P;\r
+                           E[L-1] = G;\r
+                           E[M-1] = 0.0;\r
+                   } // while (true)\r
                } // for (L = 1; L <= N; L++)\r
                \r
                // Order eigenvalues and eigenvectors\r