Porting continues.
authorilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Fri, 24 Nov 2017 23:14:23 +0000 (23:14 +0000)
committerilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Fri, 24 Nov 2017 23:14:23 +0000 (23:14 +0000)
git-svn-id: https://citdcbmipav.cit.nih.gov/repos-pub/mipav/trunk@15274 ba61647d-9d00-f842-95cd-605cb4296b96

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

index 78cef3b..58d8ef4 100644 (file)
@@ -84,6 +84,8 @@ public class DoublyConnectedSC extends AlgorithmBase {
        private int NPTQ;\r
        // 1 for solving the nonlinear system, 2 for not solving the nonlinear system\r
        private int ISOLV;\r
+       // Inverse points\r
+       private double INVERSE_POINTS[][] = null; ;\r
        private boolean testRoutine = false;\r
        private double MACHEP = 2.2204460E-16;\r
        \r
@@ -91,15 +93,17 @@ public class DoublyConnectedSC extends AlgorithmBase {
                \r
        }\r
        \r
-       public DoublyConnectedSC(int IPOLY, int NPTQ, int ISOLV, int ISPRT) {\r
+       public DoublyConnectedSC(int IPOLY, int NPTQ, int ISOLV, int ISPRT, double INVERSE_POINTS[][]) {\r
                this.IPOLY = IPOLY;\r
                this.NPTQ = NPTQ;\r
                this.ISOLV = ISOLV;\r
                this.ISPRT = ISPRT;\r
+               this.INVERSE_POINTS = INVERSE_POINTS;\r
                testRoutine = true;\r
        }\r
        \r
        public void runAlgorithm() {\r
+               int I;\r
                int M[] = new int[1];\r
                int N[] = new int[1];\r
                double Z0[][] = new double[30][2];\r
@@ -117,6 +121,9 @@ public class DoublyConnectedSC extends AlgorithmBase {
                double W1[][] = new double[30][2];\r
                double PHI0[] = new double[30];\r
                double PHI1[] = new double[30];\r
+               double ZZ[] = new double[2];\r
+               double EPS;\r
+               double WW[];\r
                \r
 double neweps;\r
                \r
@@ -166,6 +173,24 @@ double neweps;
                                DSCSOLV(TOL, IGUESS, M[0], N[0], U, C, W0, W1, PHI0, PHI1,\r
                                                Z0, Z1, ALFA0, ALFA1, NPTQ, QWORK, ISHAPE, LINEARC);\r
                        }\r
+                       \r
+                       // Output will be arranged in DSCSOLV\r
+                       \r
+                       // Compute data for theta-function and test the accuracy\r
+                       THDATA(U);\r
+                   DSCTEST(M[0],N[0],U[0],C,W0,W1,Z0,Z1,ALFA0,ALFA1,NPTQ,QWORK);\r
+                   \r
+            // Test inverse evaluations\r
+                   if ((INVERSE_POINTS != null) && (INVERSE_POINTS.length != 0)) {\r
+                       for (I = 0; I < INVERSE_POINTS.length; I++) {\r
+                           ZZ[0] = INVERSE_POINTS[I][0];       \r
+                           ZZ[1] = INVERSE_POINTS[I][1];\r
+                           EPS = 1.0E-6;\r
+                           WW = WDSC(ZZ,M[0],N[0],U[0],C,W0,W1,Z0,Z1,ALFA0,ALFA1,PHI0,PHI1,\r
+                                                 NPTQ,QWORK,EPS,1);\r
+\r
+                       } // for (I = 0; I < INVERSE_POINTS.length; I++)\r
+                   } // if ((INVERSE_POINTS != null) && (INVERSE_POINTS.length != 0))\r
                } // if (testRoutine)\r
                \r
                \r
@@ -247,6 +272,11 @@ double neweps;
             ALFA0[9] = 1.39199980651013222;\r
             ALFA0[10] = 0.819604487273064009;\r
             ALFA0[11] = 0.295854728586457991;\r
+            INVERSE_POINTS = new double[1][2];\r
+            INVERSE_POINTS[0][0] = 0.5;\r
+            INVERSE_POINTS[0][1] = 0.5;\r
+            // Gives:\r
+            // THE PREIMAGE OF ZZ=          (-0.623005768209842,0.782172159414472)\r
                } // else if (IPOLY == 2)\r
                else if (IPOLY == 3) {\r
                        M[0] = 11;\r
@@ -1495,8 +1525,190 @@ double neweps;
         for (I = 0; I < NM; I++) {\r
                QW[I+903] = QTF[I];\r
         }\r
+        \r
+        //  CHECK ERROR INFORMATION. THE DESCRIPTION OF INFO CAN BE\r
+        //  FOUND IN THE DOCUMENTATION OF HYBRD.(INFO=1: SUCCESSFUL EXIT)\r
+        if (INFO[0] == 0) {\r
+               System.out.println("INFO[0] == 0, IMPROPER INPUT PARAMETERS.");\r
+        }\r
+        else if (INFO[0] == 1) {\r
+               System.out.println("INFO[0] == 1, SUCESSFUL EXIT");\r
+               System.out.println("RELATIVE ERROR BETWEEN TWO CONSECUTIVE ITERATES IS AT MOST XTOL.");\r
+        }\r
+        else if (INFO[0] == 2) {\r
+               System.out.println("INFO[0] == 2, UNSUCESSFUL EXIT");\r
+               System.out.println("NUMBER OF CALLS TO FCN HAS REACHED OR EXCEEDED MAXFEV");\r
+        }\r
+        else if (INFO[0] == 3) {\r
+               System.out.println("INFO[0] == 3, UNSUCESSFUL EXIT");\r
+               System.out.println("XTOL IS TOO SMALL. NO FURTHER IMPROVEMENT IN THE APPROXIMATE SOLUTION X IS POSSIBLE.");\r
+        }\r
+        else if (INFO[0] == 4) {\r
+               System.out.println("INFO[0] == 4, UNSUCESSFUL EXIT");\r
+               System.out.println("ITERATION IS NOT MAKING GOOD PROGRESS, ");\r
+               System.out.println("AS MEASURED BY THE IMPROVEMENT FROM THE LAST FIVE JACOBIAN EVALUATIONS.");\r
+        }\r
+        else if (INFO[0] == 5) {\r
+               System.out.println("INFO[0] == 5, UNSUCESSFUL EXIT");\r
+               System.out.println("ITERATION IS NOT MAKING GOOD PROGRESS, ");\r
+            System.out.println("AS MEASURED BY THE IMPROVEMENT FROM THE LAST TEN ITERATIONS.");\r
+        }\r
+               \r
+        // COPY OUTPUT DATA FROM COMMON BLOCK AND PRINT THE RESULTS:\r
+        XWTRAN(M,N,X,U,C,W0,W1,PHI0,PHI1);\r
+        DSCPRINT(M,N,C,U[0],W0,W1,PHI0,PHI1,TOL,NPTQ);\r
+        return;\r
+\r
+       }\r
+       \r
+       private double[] WDSC(double ZZ[], int M,int N, double U, double C[], double W0[][], double W1[][],\r
+                       double Z0[][], double Z1[][], double ALFA0[], double ALFA1[], double PHI0[], double PHI1[],\r
+                       int NPTQ, double QWORK[], double EPS,int IOPT) {\r
+        //    ----------------------------------------------------\r
+               // COMPUTES THE INVERSE MAP AFTER  ALL ACCESSARY PARAMETERS HAVE\r
+               // BEEN DETERMINED.EULER'S SCHEME FOR SOLVING ODE IS USED TO GIVE\r
+           // THE INITIAL GUESS WHICH IS THEN REFINED BY NEWTON'S ITERATION.\r
+               // ZZ IS THE POINT AT WHICH INVERSE MAP IS EVALUATED. EPS IS THE\r
+               // REQUIRED ACCURACY SUPPLIED BY THE USER.\r
+               // NOTE: ZZ IS NOT ALLOWED TO BE A VERTEX!\r
+               \r
+               \r
+               // IF THE INVERSE EVALUATION IS NOT SUCCESSFUL, THE NAME OF THE\r
+               // FUNCTION, WDSC, WILL BE ASSIGNED WITH:\r
+               //      .. Scalar Arguments ..\r
+               // DOUBLE COMPLEX C,ZZ\r
+               // DOUBLE PRECISION EPS,U\r
+               // INTEGER IOPT,M,N,NPTQ\r
+       \r
+               // .. Array Arguments ..\r
+               // DOUBLE COMPLEX W0(M),W1(N),Z0(M),Z1(N)\r
+               // DOUBLE PRECISION ALFA0(M),ALFA1(N),PHI0(M),PHI1(N),\r
+               //                  QWORK(NPTQ* (2* (M+N)+3))\r
+               \r
+               // .. Scalars in Common ..\r
+               // DOUBLE PRECISION DLAM\r
+               // INTEGER IU\r
+               \r
+               // .. Arrays in Common ..\r
+               // DOUBLE PRECISION UARY(8),VARY(3)\r
+       \r
+               // .. Local Scalars ..\r
+               double result[] = new double[2];\r
+               double WFN[] = new double[2];\r
+               double WI[] = new double[2];\r
+               double ZS[] = new double[2];\r
+               double ZZ1[] = new double[2];\r
+               // DOUBLE COMPLEX WFN,WI,ZS,ZZ1\r
+               int INZ[] = new int[1];\r
+               int KNZ[] = new int[1];\r
+               int I,IT,K;\r
+       \r
+               // .. Local Arrays ..\r
+               double ZS0[][] = new double[50][2];\r
+               double ZS1[][] = new double[50][2];\r
+               // DOUBLE COMPLEX ZS0(50),ZS1(50)\r
+               \r
+               //     .. External Functions ..\r
+               //      DOUBLE COMPLEX WPROD,ZDSC\r
+               //      EXTERNAL WPROD,ZDSC\r
+       \r
+               //     .. External Subroutines ..\r
+               //      EXTERNAL NEARZ\r
+               \r
+               //     .. Common blocks ..\r
+               //      COMMON /PARAM4/UARY,VARY,DLAM,IU\r
+               \r
+               result[0] = 0.0;\r
+               result[1] = 0.0;\r
+               \r
+               //  1.FORM WORK ARRAYS:\r
+               /*for (I = 0; I < M; I++) {\r
+                   ZS0[I][0] = Z0[I][0];\r
+                   ZS0[I][1] = Z0[I][1];\r
+               } // for (I = 0; I < M; I++)\r
+               for (I = 0; I < N; I++) {\r
+                   ZS1[I][0] = Z1[I][0];\r
+                   ZS1[I][1] = Z1[I][1];\r
+               } // for (I = 0; I < N; I++)\r
+               \r
+               // 2.GENERATE THE INITIAL VALUE FOR SOLVING ODE:\r
+               while (true) {\r
+                   NEARZ(M,N,ZS0,ZS1,ALFA0,ZZ,KNZ,INZ);\r
+                   if (INZ[0] == 2) {\r
+                       System.err.println("THE INVERSE EVALUATION FAILED");    \r
+                       System.err.println("AT POINT Z = (" + ZZ[0] + ", " + ZZ[1] + ")");\r
+                       System.err.println("THE DESIGNED INVERSION PROCEDURE EXPERIENCED");\r
+                       System.err.println("SOME ESSENTIAL DIFFICULATIES");\r
+                       return result;\r
+                   } // if (INZ[0] == 2)\r
+                   if (INZ[0] == 0) {\r
+                       if (KNZ[0] >= 2) {\r
+                           WI[0] = (W0[KNZ[0]-1][0]+W0[KNZ[0]-2][0])/2.0;\r
+                           WI[1] = (W0[KNZ[0]-1][1]+W0[KNZ[0]-2][1])/2.0;\r
+                       }\r
+                       else {\r
+                             WI[0] = (W0[KNZ[0]-1][0]+W0[M-1][0]))/2.0\r
+                      WI[1] = (W0[KNZ[0]-1][1]+W0[M-1][1]))/2.0\r
+                       }\r
+\r
+                   } // if (INZ[0] == 0)\r
+                   else {\r
+                         IF (KNZ.GE.2) THEN\r
+                             WI = (W1(KNZ)+W1(KNZ-1))/2.D0\r
+\r
+                         ELSE\r
+                             WI = (W0(KNZ)+W1(N))/2.D0\r
+                         END IF\r
+\r
+                         WI = U*WI/ABS(WI)\r
+                   } //else\r
+\r
+                     ZS = ZDSC(WI,0,2,M,N,U,C,W0,W1,Z0,Z1,ALFA0,ALFA1,PHI0,PHI1,NPTQ,\r
+                    +     QWORK,1)\r
+               C\r
+               C  3.SOLVE ODE INITIAL VALUE PROBLEM (ALONG LINE SEGMENT FROM\r
+               C    ZS TO ZZ) BY EULER'S SCHEME TO GENERATE THE INITIAL GUESS:\r
+                     DO 40 K = 1,20\r
+                         WI = WI + (ZZ-ZS)/ (20.D0*C*WPROD(WI,M,N,U,W0,W1,ALFA0,ALFA1))\r
+                  40 CONTINUE\r
+                     IF (ABS(WI).GT.1.D0) WI = WI/ (ABS(WI)+ABS(1.D0-ABS(WI)))\r
+                     IF (ABS(WI).LT.U) WI = U*WI* (ABS(WI)+ABS(U-ABS(WI)))/ABS(WI)\r
+               C\r
+               C  4.REFINE THE SOLUTION BY NEWTON'S ITERATION:\r
+                     DO 50 IT = 1,15\r
+                         WFN = ZZ - ZDSC(WI,0,2,M,N,U,C,W0,W1,Z0,Z1,ALFA0,ALFA1,PHI0,\r
+                    +          PHI1,NPTQ,QWORK,IOPT)\r
+                         WI = WI + WFN/ (C*WPROD(WI,M,N,U,W0,W1,ALFA0,ALFA1))\r
+                         IF (ABS(WI).GT.1.D0) WI = WI/ (ABS(WI)+ABS(1.D0-ABS(WI)))\r
+                         IF (ABS(WI).LT.U) WI = U*WI* (ABS(WI)+ABS(U-ABS(WI)))/ABS(WI)\r
+                         IF (ABS(WFN).LT.EPS) GO TO 60\r
+                  50 CONTINUE\r
+               C\r
+               C  THE ITERATION FAILED TO MEET THE TOLERANCE IN 15 ITERATIONS.\r
+               C  TRY A DIFFERENT VERTEX AS A REFERENCE POINT:\r
+                     ZZ1 = (1.D0,1.D0) + ZZ\r
+                     GO TO 70\r
+\r
+                  60 WDSC = WI\r
+               C\r
+               C  5.VERIFICATION:\r
+                     ZZ1 = ZDSC(WI,0,2,M,N,U,C,W0,W1,Z0,Z1,ALFA0,ALFA1,PHI0,PHI1,NPTQ,\r
+                    +      QWORK,1)\r
+                  70 IF (ABS(WI).GE.U .AND. ABS(WI).LE.1.D0 .AND.\r
+                    +    ABS(ZZ-ZZ1).LE.1.D-3) GO TO 90\r
+                     IF (INZ.EQ.0) THEN\r
+                         ZS0(KNZ) = ZZ\r
+\r
+                     ELSE\r
+                         ZS1(KNZ) = ZZ\r
+                     END IF\r
+\r
+               } // while (true)*/\r
+\r
+               return result;\r
        }\r
 \r
+\r
        \r
        private void ANGLES(int MN, double Z01[][], double ALFA01[], int I01) {\r
            // Computes the interior angles of a doubly\r
@@ -1530,6 +1742,57 @@ double neweps;
                return;\r
        }\r
        \r
+       private void NEARZ(int M,int N, double Z0[][], double Z1[][],\r
+                       double ALFA0[], double Z[], int KNZ[],int INZ[]) {\r
+           //   --------------------------------------------------\r
+           // GIVEN PT Z, THIS ROUTINE DETERMINES THE NEAREST VERTEX TO\r
+           // THE PT Z.ON RETURN INTEGER'INZ'INDICATES THAT THE NEAREST\r
+           // PT IS FOUND  EITHER IN Z0 (INZ=0 ) OR IN Z1 ( INZ=1). KNZ\r
+           // CONTAINS THE CORRESPONDING INDEX IN Z0 OR 1N Z1.\r
+           // (IF ON RETURN INZ=2, THAT MEANS NO APPROPRIATE VERTEX IS\r
+           // FOUND, WHICH IS A DEVICE USED FOR INVERSE MAPPING ROUTINE.\r
+           // NOTE: THE VERTICES CORRESPONDING TO INFINITE VERTICES\r
+        //       WILL BE SKIPPED (OUTER POLYGON ONLY).\r
+           // SEE ROUTINE DSCSOLV FOR THE REST OF CALLING SEQUENCE.\r
+       \r
+       \r
+           // .. Scalar Arguments ..\r
+           // DOUBLE COMPLEX Z\r
+           // INTEGER INZ,KNZ,M,N\r
+       \r
+           // .. Array Arguments ..\r
+           // DOUBLE COMPLEX Z0(M),Z1(N)\r
+           // DOUBLE PRECISION ALFA0(M)\r
+       \r
+           // .. Local Scalars ..\r
+           double D,DIST;\r
+           int I;\r
+       \r
+           INZ[0] = 2;\r
+           DIST = 99.0;\r
+           for (I = 1; I <= M; I++) {\r
+               D = scm.zabs(Z[0]-Z0[I-1][0],Z[1]-Z0[I-1][1]);\r
+               if (D <= 1.0E-11 || D >= DIST ||\r
+                    ALFA0[I-1] <= 0.0) {\r
+                       continue;\r
+               }\r
+               KNZ[0] = I;\r
+               DIST = D;\r
+               INZ[0] = 0;\r
+           } // for (I = 1; I <= M; I++)\r
+           for (I = 1; I <= N; I++) {\r
+                 D = scm.zabs(Z[0]-Z1[I-1][0],Z[1]-Z1[I-1][1]);\r
+                 if (D <= 1.0E-11 || D >= DIST) {\r
+                         continue;\r
+                 }\r
+                 DIST = D;\r
+                 KNZ[0] = I;\r
+                 INZ[0] = 1;\r
+           } // for (I = 1; I <= N; I++)\r
+           return;\r
+       }\r
+\r
+       \r
        private double DIST(int M,int N, double W0[][], double W1[][],\r
                        double W[], int KWA, int IC) {\r
        //   -----------------------------------\r
@@ -1649,6 +1912,131 @@ double neweps;
                return;\r
        }\r
        \r
+       private void DSCPRINT(int M,int N, double C[], double U, double W0[][], double W1[][],\r
+                       double PHI0[], double PHI1[], double TOL, int NPTQ) {\r
+           // ---------------------------------------------------------\r
+           // PRINTS THE COMPUTED SC-PARAMETERS AND SOME\r
+        // OTHER CONTROLING PARAMETERS FOR REFERENCE:\r
+       \r
+       \r
+         // OPEN FILE DSCPACK FOR OUTPUT:\r
+         //     .. Scalar Arguments ..\r
+         // DOUBLE COMPLEX C\r
+         // DOUBLE PRECISION TOL,U\r
+         // INTEGER M,N,NPTQ\r
+       \r
+         // .. Array Arguments ..\r
+         // DOUBLE COMPLEX W0(M),W1(N)\r
+         // DOUBLE PRECISION PHI0(M),PHI1(N)\r
+       \r
+         //.. Local Scalars ..\r
+         int K;\r
+       \r
+        System.out.println(" PARAMETERS DEFINING MAP:   (M = "+M+")   (N = "+N+")   (NPTQ = "+NPTQ+")   (TOL = "+TOL+")");\r
+        System.out.println(" U = " + U);\r
+        System.out.println(" C = (" + C[0] + ", " + C[1] + ")");\r
+     for (K = 0; K < M; K++) {\r
+        System.out.println("W0["+K+"] = (" + W0[K][0] + ", " + W0[K][1] + "   PHI0["+K+"] = " + PHI0[K]);\r
+     }\r
+     for (K = 0; K < N; K++) {\r
+        System.out.println("W1["+K+"] = (" + W1[K][0] + ", " + W1[K][1] + "   PHI1["+K+"] = " + PHI1[K]);\r
+     }\r
+            \r
+        return;\r
+       \r
+       }\r
+       \r
+       private void DSCTEST(int M,int N, double U, double C[], double W0[][], double W1[][],\r
+                       double Z0[][], double Z1[][], double ALFA0[], double ALFA1[], int NPTQ, double QWORK[]) {\r
+           //   ----------------------------------------------------------------\r
+           // TESTS THE COMPUTED PARAMETERS FOR ACCURACY BY COMPUTING VERTICES\r
+           // Z0(K) NUMERICALLY AND COMPARING WITH THE EXACT ONES. ON OUTPUT,\r
+           // THE MAXIMUM AND MINIMUM ERRORS ACHIEVED WILL BE GIVEN TOGETHER\r
+           // WITH THE CORRESPONDING INDICES KMAX AND KMIN RESPECTIVELY. SEE\r
+           // ROUTINE DSCSOLV FOR THE CALLING SEQUENCE.\r
+       \r
+       \r
+           // .. Scalar Arguments ..\r
+           // DOUBLE COMPLEX C\r
+           // DOUBLE PRECISION U\r
+           // INTEGER M,N,NPTQ\r
+       \r
+           // .. Array Arguments ..\r
+           // DOUBLE COMPLEX W0(M),W1(N),Z0(M),Z1(N)\r
+           // DOUBLE PRECISION ALFA0(M),ALFA1(N),QWORK(NPTQ* (2* (M+N)+3))\r
+       \r
+           // .. Scalars in Common ..\r
+           // DOUBLE PRECISION DLAM\r
+           // INTEGER IU\r
+       \r
+           // .. Arrays in Common ..\r
+           // DOUBLE PRECISION UARY(8),VARY(3)\r
+       \r
+           // .. Local Scalars ..\r
+               double WA[] = new double[2];\r
+               double ZC[] = new double[2];\r
+               double ZTEST[] = new double[2];\r
+               double wout[];\r
+               double cr[] = new double[1];\r
+               double ci[] = new double[1];\r
+           // DOUBLE COMPLEX WA,ZC,ZTEST\r
+           double D,D1,DIST,ERRMAX,ERRMIN;\r
+           int I,K;\r
+           int KMAX = -1;\r
+           int KMIN = -1;\r
+           int KWA = 0;\r
+       \r
+           // .. External Functions ..\r
+           // DOUBLE COMPLEX WQUAD\r
+           // EXTERNAL WQUAD\r
+       \r
+           // .. Common blocks ..\r
+           // COMMON /PARAM4/UARY,VARY,DLAM,IU\r
+       \r
+           ERRMAX = 0.0;\r
+           ERRMIN = 99.0;\r
+           for (K = 0; K < M - 1; K++) {\r
+               if (ALFA0[K] <= 0.0) {\r
+                       continue;\r
+               }\r
+               DIST = 2.0;\r
+               for (I = 0; I < N; I++) {\r
+                   D = scm.zabs(W0[K][0]-W1[I][0], W0[K][1]-W1[I][1]);\r
+                   if (D >= DIST) {\r
+                       continue;\r
+                   }\r
+                   DIST = D;\r
+                   WA[0] = W1[I][0];\r
+                   WA[1] = W1[I][1];\r
+                   KWA = I;\r
+                   ZC[0] = Z1[I][0];\r
+                   ZC[1] = Z1[I][1];\r
+               } // for (I = 0; I < N; I++)\r
+               wout = WQUAD(WA,0.0,KWA,1,W0[K],0.0,K,0,0.0,M,N,U,\r
+                                W0,W1,ALFA0,ALFA1,NPTQ,QWORK,0,1);\r
+               scm.zmlt(C[0], C[1], wout[0], wout[1], cr, ci);\r
+               ZTEST[0] = ZC[0] + cr[0];\r
+               ZTEST[1] = ZC[1] + ci[0];\r
+               D1 = scm.zabs(Z0[K][0]-ZTEST[0],Z0[K][1]-ZTEST[1]);\r
+               if (D1 > ERRMAX) {\r
+                   ERRMAX = D1;\r
+                   KMAX = K;\r
+               } // if (D1 > ERRMAX)\r
+\r
+               if (D1 < ERRMIN) {\r
+                   ERRMIN = D1;\r
+                   KMIN = K;\r
+               } // if (D1 < ERRMIN)\r
+\r
+           } // for (K = 0; K < M - 1; K++)\r
+           System.out.println("ACCURACY TEST: ");\r
+           System.out.println("MAXIMUM ERROR = " + ERRMAX + " ACHIEVED AT KMAX = " + KMAX);\r
+           System.out.println("MINIMUM ERROR = " + ERRMIN + " ACHIEVED AT KMIN = " + KMIN);\r
+           return;\r
+    }\r
+\r
+\r
+       \r
        private void HYBRD(int FCN,int N, double X[], double FVEC[], double XTOL, int MAXFEV,\r
                        int ML, int MU, double EPSFCN, double DIAG[], int MODE,\r
                    double FACTOR, int NPRINT,int INFO[], int NFEV[], double FJAC[][],\r
@@ -1823,11 +2211,15 @@ double neweps;
                //     ..\r
                //     .. Local Scalars ..\r
                      double ACTRED,EPSMCH,FNORM,FNORM1,ONE,P0001,P001,\r
-                                     P1,P5,PNORM,PRERED,RATIO,SUM,TEMP,XNORM,ZERO;\r
+                                     P1,P5,PNORM,PRERED,RATIO,SUM,TEMP,ZERO;\r
                      double DELTA = 0.0;\r
+                     double XNORM = 0.0;\r
+                     double ratio;\r
+                     double A[][];\r
                      int IFLAG[] = new int[1];\r
                      int I,ITER,J,JM1,L,MSUM,NCFAIL,NCSUC,NSLOW1,NSLOW2;\r
-                     boolean JEVAL,SING;\r
+                     boolean JEVAL;\r
+                     boolean SING[] = new boolean[1];\r
                //     ..\r
            //     .. Local Arrays ..\r
                      int IWA[] = new int[1];\r
@@ -1985,7 +2377,7 @@ double neweps;
                \r
                    // COPY THE TRIANGULAR FACTOR OF THE QR FACTORIZATION INTO R.\r
                \r
-                     SING = false;\r
+                     SING[0] = false;\r
                      for (J = 1; J <= N; J++) {\r
                          L = J;\r
                          JM1 = J - 1;\r
@@ -1997,7 +2389,7 @@ double neweps;
                          } // if (JM1 >= 1) \r
                          R[L-1] = WA1[J-1];\r
                          if (WA1[J-1] == ZERO) {\r
-                                 SING = true;\r
+                                 SING[0] = true;\r
                          }\r
                      } // for (J = 1; J <= N; J++)\r
                \r
@@ -2065,107 +2457,149 @@ double neweps;
                         // COMPUTE THE SCALED ACTUAL REDUCTION.\r
                \r
                         ACTRED = -ONE;\r
-                     /*IF (FNORM1.LT.FNORM) ACTRED = ONE - (FNORM1/FNORM)**2\r
-               C\r
-               C           COMPUTE THE SCALED PREDICTED REDUCTION.\r
-               C\r
-                     L = 1\r
-                     DO 220 I = 1,N\r
-                         SUM = ZERO\r
-                         DO 210 J = I,N\r
-                             SUM = SUM + R(L)*WA1(J)\r
-                             L = L + 1\r
-                 210     CONTINUE\r
-                         WA3(I) = QTF(I) + SUM\r
-                 220 CONTINUE\r
-                     TEMP = ENORM(N,WA3)\r
-                     PRERED = ZERO\r
-                     IF (TEMP.LT.FNORM) PRERED = ONE - (TEMP/FNORM)**2\r
-               C\r
-               C           COMPUTE THE RATIO OF THE ACTUAL TO THE PREDICTED\r
-               C           REDUCTION.\r
-               C\r
-                     RATIO = ZERO\r
-                     IF (PRERED.GT.ZERO) RATIO = ACTRED/PRERED\r
-               C\r
-               C           UPDATE THE STEP BOUND.\r
-               C\r
-                     IF (RATIO.GE.P1) GO TO 230\r
-                     NCSUC = 0\r
-                     NCFAIL = NCFAIL + 1\r
-                     DELTA = P5*DELTA\r
-                     GO TO 240\r
-\r
-                 230 CONTINUE\r
-                     NCFAIL = 0\r
-                     NCSUC = NCSUC + 1\r
-                     IF (RATIO.GE.P5 .OR. NCSUC.GT.1) DELTA = DMAX1(DELTA,PNORM/P5)\r
-                     IF (DABS(RATIO-ONE).LE.P1) DELTA = PNORM/P5\r
-                 240 CONTINUE\r
-               C\r
-               C           TEST FOR SUCCESSFUL ITERATION.\r
-               C\r
-                     IF (RATIO.LT.P0001) GO TO 260\r
-               C\r
-               C           SUCCESSFUL ITERATION. UPDATE X, FVEC, AND THEIR NORMS.\r
-               C\r
-                     DO 250 J = 1,N\r
-                         X(J) = WA2(J)\r
-                         WA2(J) = DIAG(J)*X(J)\r
-                         FVEC(J) = WA4(J)\r
-                 250 CONTINUE\r
-                     XNORM = ENORM(N,WA2)\r
-                     FNORM = FNORM1\r
-                     ITER = ITER + 1\r
-                 260 CONTINUE\r
-               C\r
-               C           DETERMINE THE PROGRESS OF THE ITERATION.\r
-               C\r
-                     NSLOW1 = NSLOW1 + 1\r
-                     IF (ACTRED.GE.P001) NSLOW1 = 0\r
-                     IF (JEVAL) NSLOW2 = NSLOW2 + 1\r
-                     IF (ACTRED.GE.P1) NSLOW2 = 0\r
-               C\r
-               C           TEST FOR CONVERGENCE.\r
-               C\r
-                     IF (DELTA.LE.XTOL*XNORM .OR. FNORM.EQ.ZERO) INFO = 1\r
-                     IF (INFO.NE.0) GO TO 300\r
-               C\r
-               C           TESTS FOR TERMINATION AND STRINGENT TOLERANCES.\r
-               C\r
-                     IF (NFEV.GE.MAXFEV) INFO = 2\r
-                     IF (P1*DMAX1(P1*DELTA,PNORM).LE.EPSMCH*XNORM) INFO = 3\r
-                     IF (NSLOW2.EQ.5) INFO = 4\r
-                     IF (NSLOW1.EQ.10) INFO = 5\r
-                     IF (INFO.NE.0) GO TO 300\r
-               C\r
-               C           CRITERION FOR RECALCULATING JACOBIAN APPROXIMATION\r
-               C           BY FORWARD DIFFERENCES.\r
-               C\r
-                     IF (NCFAIL.EQ.2) GO TO 290\r
-               C\r
-               C           CALCULATE THE RANK ONE MODIFICATION TO THE JACOBIAN\r
-               C           AND UPDATE QTF IF NECESSARY.\r
-               C\r
-                     DO 280 J = 1,N\r
-                         SUM = ZERO\r
-                         DO 270 I = 1,N\r
-                             SUM = SUM + FJAC(I,J)*WA4(I)\r
-                 270     CONTINUE\r
-                         WA2(J) = (SUM-WA3(J))/PNORM\r
-                         WA1(J) = DIAG(J)* ((DIAG(J)*WA1(J))/PNORM)\r
-                         IF (RATIO.GE.P0001) QTF(J) = SUM\r
-                 280 CONTINUE\r
-               C\r
-               C           COMPUTE THE QR FACTORIZATION OF THE UPDATED JACOBIAN.\r
-               C\r
-                     CALL R1UPDT(N,N,R,LR,WA1,WA2,WA3,SING)\r
-                     CALL R1MPYQ(N,N,FJAC,LDFJAC,WA2,WA3)\r
-                     CALL R1MPYQ(1,N,QTF,1,WA2,WA3)*/\r
+                       if (FNORM1 < FNORM) {\r
+                               ratio = FNORM1/FNORM;\r
+                               ACTRED = ONE - ratio * ratio;\r
+                       }\r
+               \r
+                // COMPUTE THE SCALED PREDICTED REDUCTION.\r
+\r
+                       L = 0;\r
+                       for (I = 0; I < N; I++) {\r
+                           SUM = ZERO;\r
+                           for (J = I; J < N; J++) {\r
+                               SUM = SUM + R[L]*WA1[J];\r
+                               L = L + 1;\r
+                           } // for (J = I; J < N; J++)\r
+                           WA3[I] = QTF[I] + SUM;\r
+                       } // for (I = 0; I < N; I++)\r
+                       TEMP = ENORM(N,WA3);\r
+                       PRERED = ZERO;\r
+                       if (TEMP< FNORM) {\r
+                               ratio = TEMP/FNORM;\r
+                               PRERED = ONE - ratio * ratio;\r
+                       }\r
+               \r
+                       // COMPUTE THE RATIO OF THE ACTUAL TO THE PREDICTED\r
+                       // REDUCTION.\r
+               \r
+                       RATIO = ZERO;\r
+                       if (PRERED > ZERO) {\r
+                               RATIO = ACTRED/PRERED;\r
+                       }\r
+               \r
+                       // UPDATE THE STEP BOUND.\r
+               \r
+                       if (RATIO < P1) {\r
+                           NCSUC = 0;\r
+                           NCFAIL = NCFAIL + 1;\r
+                           DELTA = P5*DELTA;\r
+                     }\r
+                     else {\r
+                           NCFAIL = 0;\r
+                           NCSUC = NCSUC + 1;\r
+                           if (RATIO >= P5 || NCSUC > 1) {\r
+                               DELTA = Math.max(DELTA,PNORM/P5);\r
+                           }\r
+                           if (Math.abs(RATIO-ONE) <= P1) {\r
+                               DELTA = PNORM/P5;\r
+                           }\r
+                     }\r
+               \r
+                     // TEST FOR SUCCESSFUL ITERATION.\r
+               \r
+                     if (RATIO >= P0001) {\r
+               \r
+                         // SUCCESSFUL ITERATION. UPDATE X, FVEC, AND THEIR NORMS.\r
+               \r
+                         for (J = 0; J < N; J++) {\r
+                             X[J] = WA2[J];\r
+                             WA2[J] = DIAG[J]*X[J];\r
+                             FVEC[J] = WA4[J];\r
+                         } // for (J = 0; J < N; J++)\r
+                         XNORM = ENORM(N,WA2);\r
+                         FNORM = FNORM1;\r
+                         ITER = ITER + 1;\r
+                     } // if (RATIO >= P0001)\r
+               \r
+              // DETERMINE THE PROGRESS OF THE ITERATION.\r
+               \r
+                     NSLOW1 = NSLOW1 + 1;\r
+                     if (ACTRED >= P001) {\r
+                         NSLOW1 = 0;\r
+                     }\r
+                     if (JEVAL) {\r
+                         NSLOW2 = NSLOW2 + 1;\r
+                     }\r
+                     if (ACTRED >= P1) {\r
+                         NSLOW2 = 0;\r
+                     }\r
+               \r
+                     // TEST FOR CONVERGENCE.\r
+               \r
+                     if (DELTA <= XTOL*XNORM || FNORM == ZERO) {\r
+                         INFO[0] = 1;\r
+                     }\r
+                     if (INFO[0] != 0) {\r
+                         break outer;\r
+                     }\r
+               \r
+                     // TESTS FOR TERMINATION AND STRINGENT TOLERANCES.\r
+               \r
+                     if (NFEV[0] >= MAXFEV) {\r
+                         INFO[0] = 2;\r
+                     }\r
+                     if (P1*Math.max(P1*DELTA,PNORM) <= EPSMCH*XNORM) {\r
+                         INFO[0] = 3;\r
+                     }\r
+                     if (NSLOW2 == 5) {\r
+                         INFO[0] = 4;\r
+                     }\r
+                     if (NSLOW1 == 10) {\r
+                         INFO[0] = 5;\r
+                     }\r
+                     if (INFO[0] != 0) {\r
+                         break outer;\r
+                     }\r
+               \r
+                     // CRITERION FOR RECALCULATING JACOBIAN APPROXIMATION\r
+                     // BY FORWARD DIFFERENCES.\r
+               \r
+                     if (NCFAIL == 2) {\r
+                         continue outer;\r
+                     }\r
+               \r
+                     // CALCULATE THE RANK ONE MODIFICATION TO THE JACOBIAN\r
+                     // AND UPDATE QTF IF NECESSARY.\r
+               \r
+                     for (J = 0; J < N; J++) {\r
+                         SUM = ZERO;\r
+                         for (I = 0; I < N; I++) {\r
+                             SUM = SUM + FJAC[I][J]*WA4[I];\r
+                         } // for (I = 0; I < N; I++)\r
+                         WA2[J] = (SUM-WA3[J])/PNORM;\r
+                         WA1[J] = DIAG[J]* ((DIAG[J]*WA1[J])/PNORM);\r
+                         if (RATIO >= P0001) {\r
+                                 QTF[J] = SUM;\r
+                         }\r
+                     } // for (J = 0; J < N; J++)\r
+               \r
+                     // COMPUTE THE QR FACTORIZATION OF THE UPDATED JACOBIAN.\r
+               \r
+                     R1UPDT(N,N,R,LR,WA1,WA2,WA3,SING);\r
+                     R1MPYQ(N,N,FJAC,LDFJAC,WA2,WA3);\r
+                     A = new double[1][N];\r
+                     for (I = 0; I < N; I++) {\r
+                         A[0][I] = QTF[I];\r
+                     }\r
+                     R1MPYQ(1,N,A,1,WA2,WA3);\r
+                     for (I = 0; I < N; I++) {\r
+                         QTF[I] = A[0][I];\r
+                     }\r
                \r
                         // END OF THE INNER LOOP.\r
                \r
-                         JEVAL = false;\r
+                       JEVAL = false;\r
                     } // while (true)\r
                \r
                    // END OF THE OUTER LOOP.\r
@@ -2445,6 +2879,338 @@ double neweps;
                return;\r
        }\r
        \r
+       private void R1MPYQ(int M,int N, double A[][], int LDA,\r
+                       double V[], double W[]) {\r
+           //     **********\r
+       \r
+           // SUBROUTINE R1MPYQ\r
+       \r
+           // GIVEN AN M BY N MATRIX A, THIS SUBROUTINE COMPUTES A*Q WHERE\r
+           // Q IS THE PRODUCT OF 2*(N - 1) TRANSFORMATIONS\r
+       \r
+        // GV(N-1)*...*GV(1)*GW(1)*...*GW(N-1)\r
+       \r
+           // AND GV(I), GW(I) ARE GIVENS ROTATIONS IN THE (I,N) PLANE WHICH\r
+           // ELIMINATE ELEMENTS IN THE I-TH AND N-TH PLANES, RESPECTIVELY.\r
+           // Q ITSELF IS NOT GIVEN, RATHER THE INFORMATION TO RECOVER THE\r
+           // GV, GW ROTATIONS IS SUPPLIED.\r
+       \r
+           // THE SUBROUTINE STATEMENT IS\r
+       \r
+           // SUBROUTINE R1MPYQ(M,N,A,LDA,V,W)\r
+       \r
+           // WHERE\r
+       \r
+           //   M IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER\r
+           //     OF ROWS OF A.\r
+       \r
+           //   N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER\r
+           //     OF COLUMNS OF A.\r
+       \r
+           //   A IS AN M BY N ARRAY. ON INPUT A MUST CONTAIN THE MATRIX\r
+           //     TO BE POSTMULTIPLIED BY THE ORTHOGONAL MATRIX Q\r
+           //     DESCRIBED ABOVE. ON OUTPUT A*Q HAS REPLACED A.\r
+       \r
+           //   LDA IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN M\r
+           //     WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY A.\r
+       \r
+           //   V IS AN INPUT ARRAY OF LENGTH N. V(I) MUST CONTAIN THE\r
+           //     INFORMATION NECESSARY TO RECOVER THE GIVENS ROTATION GV(I)\r
+           //     DESCRIBED ABOVE.\r
+       \r
+           //  W IS AN INPUT ARRAY OF LENGTH N. W(I) MUST CONTAIN THE\r
+           //    INFORMATION NECESSARY TO RECOVER THE GIVENS ROTATION GW(I)\r
+           //    DESCRIBED ABOVE.\r
+       \r
+           // ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980.\r
+           // BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE\r
+       \r
+           // **********\r
+        // .. Scalar Arguments ..\r
+           // INTEGER LDA,M,N\r
+       \r
+           // .. Array Arguments ..\r
+           // DOUBLE PRECISION A(LDA,N),V(N),W(N)\r
+       \r
+           // .. Local Scalars ..\r
+           double ONE,TEMP;\r
+           double COS = 0.0;\r
+           double SIN = 0.0;\r
+           int I,J,NM1,NMJ;\r
+       \r
+           // .. Data statements ..\r
+           ONE = 1.0;\r
+       \r
+           // APPLY THE FIRST SET OF GIVENS ROTATIONS TO A.\r
+       \r
+           NM1 = N - 1;\r
+           if (NM1 < 1) {\r
+               return;\r
+           }\r
+           for (NMJ = 1; NMJ <= NM1; NMJ++) {\r
+               J = N - NMJ;\r
+               if (Math.abs(V[J-1]) > ONE) {\r
+                       COS = ONE/V[J-1];\r
+               }\r
+               if (Math.abs(V[J-1]) > ONE) {\r
+                       SIN = Math.sqrt(ONE-COS*COS);\r
+               }\r
+               if (Math.abs(V[J-1]) <= ONE) {\r
+                       SIN = V[J-1];\r
+               }\r
+               if (Math.abs(V[J-1]) <= ONE) {\r
+                       COS = Math.sqrt(ONE-SIN*SIN);\r
+               }\r
+               for (I = 1; I <= M; I++) {\r
+                   TEMP = COS*A[I-1][J-1] - SIN*A[I-1][N-1];\r
+                   A[I-1][N-1] = SIN*A[I-1][J-1] + COS*A[I-1][N-1];\r
+                   A[I-1][J-1] = TEMP;\r
+               } // for (I = 1; I <= M; I++)\r
+           } // for (NMJ = 1; NMJ <= NM1; NMJ++) \r
+       \r
+           // APPLY THE SECOND SET OF GIVENS ROTATIONS TO A.\r
+       \r
+           for (J = 1; J <= NM1; J++) {\r
+               if (Math.abs(W[J-1]) > ONE) {\r
+                       COS = ONE/W[J-1];\r
+               }\r
+               if (Math.abs(W[J-1]) > ONE) {\r
+                       SIN = Math.sqrt(ONE-COS*COS);\r
+               }\r
+               if (Math.abs(W[J-1]) <= ONE) {\r
+                       SIN = W[J-1];\r
+               }\r
+               if (Math.abs(W[J-1]) <= ONE) {\r
+                       COS = Math.sqrt(ONE-SIN*SIN);\r
+               }\r
+               for (I = 1; I <= M; I++) {\r
+                   TEMP = COS*A[I-1][J-1] + SIN*A[I-1][N-1];\r
+                   A[I-1][N-1] = -SIN*A[I-1][J-1] + COS*A[I-1][N-1];\r
+                   A[I-1][J-1] = TEMP;\r
+               } // for (I = 1; I <= M; I++)\r
+           } // for (J = 1; J <= NM1; J++)\r
+           return;\r
+       }\r
+\r
+       \r
+       private void R1UPDT(int M,int N, double S[], int LS, double U[],\r
+                       double V[], double W[], boolean SING[]) {\r
+           //     **********\r
+       \r
+           // SUBROUTINE R1UPDT\r
+       \r
+           // GIVEN AN M BY N LOWER TRAPEZOIDAL MATRIX S, AN M-VECTOR U,\r
+           // AND AN N-VECTOR V, THE PROBLEM IS TO DETERMINE AN\r
+           // ORTHOGONAL MATRIX Q SUCH THAT\r
+       \r
+           //                   T\r
+           //           (S + U*V )*Q\r
+       \r
+           //     IS AGAIN LOWER TRAPEZOIDAL.\r
+       \r
+           // THIS SUBROUTINE DETERMINES Q AS THE PRODUCT OF 2*(N - 1)\r
+           // TRANSFORMATIONS\r
+       \r
+           // GV(N-1)*...*GV(1)*GW(1)*...*GW(N-1)\r
+       \r
+           // WHERE GV(I), GW(I) ARE GIVENS ROTATIONS IN THE (I,N) PLANE\r
+           // WHICH ELIMINATE ELEMENTS IN THE I-TH AND N-TH PLANES,\r
+           // RESPECTIVELY. Q ITSELF IS NOT ACCUMULATED, RATHER THE\r
+           // INFORMATION TO RECOVER THE GV, GW ROTATIONS IS RETURNED.\r
+       \r
+           // THE SUBROUTINE STATEMENT IS\r
+       \r
+           // SUBROUTINE R1UPDT(M,N,S,LS,U,V,W,SING)\r
+       \r
+           // WHERE\r
+       \r
+           //   M IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER\r
+           //     OF ROWS OF S.\r
+       \r
+           //   N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER\r
+           //     OF COLUMNS OF S. N MUST NOT EXCEED M.\r
+       \r
+           //   S IS AN ARRAY OF LENGTH LS. ON INPUT S MUST CONTAIN THE LOWER\r
+           //     TRAPEZOIDAL MATRIX S STORED BY COLUMNS. ON OUTPUT S CONTAINS\r
+           //     THE LOWER TRAPEZOIDAL MATRIX PRODUCED AS DESCRIBED ABOVE.\r
+       \r
+           //   LS IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN\r
+           //     (N*(2*M-N+1))/2.\r
+       \r
+           //   U IS AN INPUT ARRAY OF LENGTH M WHICH MUST CONTAIN THE\r
+           //     VECTOR U.\r
+       \r
+           //   V IS AN ARRAY OF LENGTH N. ON INPUT V MUST CONTAIN THE VECTOR\r
+           //     V. ON OUTPUT V(I) CONTAINS THE INFORMATION NECESSARY TO\r
+           //     RECOVER THE GIVENS ROTATION GV(I) DESCRIBED ABOVE.\r
+       \r
+           //   W IS AN OUTPUT ARRAY OF LENGTH M. W(I) CONTAINS INFORMATION\r
+           //     NECESSARY TO RECOVER THE GIVENS ROTATION GW(I) DESCRIBED\r
+           //     ABOVE.\r
+       \r
+           //   SING IS A LOGICAL OUTPUT VARIABLE. SING IS SET TRUE IF ANY\r
+           //     OF THE DIAGONAL ELEMENTS OF THE OUTPUT S ARE ZERO. OTHERWISE\r
+           //     SING IS SET FALSE.\r
+       \r
+           // ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980.\r
+           // BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE,\r
+           // JOHN L. NAZARETH\r
+       \r
+           //     **********\r
+           // .. Scalar Arguments ..\r
+           // INTEGER LS,M,N\r
+       \r
+           // .. Array Arguments ..\r
+           // DOUBLE PRECISION S(LS),U(M),V(N),W(M)\r
+           // LOGICAL SING[]\r
+       \r
+           // .. Local Scalars ..\r
+           double COS,COTAN,GIANT,ONE,P25,P5,SIN,TAN,TAU,TEMP,ZERO;\r
+           int I,J,JJ,L,NM1,NMJ;\r
+       \r
+           // Data statements ..\r
+           ONE = 1.0;\r
+           P5 = 5.0E-1;\r
+           P25 = 2.5E-1;\r
+           ZERO = 0.0;\r
+       \r
+           // GIANT IS THE LARGEST MAGNITUDE.\r
+       \r
+           GIANT = Double.MAX_VALUE;\r
+       \r
+           // INITIALIZE THE DIAGONAL ELEMENT POINTER.\r
+       \r
+           JJ = (N* (2*M-N+1))/2 - (M-N);\r
+       \r
+           // MOVE THE NONTRIVIAL PART OF THE LAST COLUMN OF S INTO W.\r
+       \r
+           L = JJ;\r
+           for (I = N; I <= M; I++) {\r
+               W[I-1] = S[L-1];\r
+               L = L + 1;\r
+           } // for (I = N; I <= M; I++)\r
+       \r
+           // ROTATE THE VECTOR V INTO A MULTIPLE OF THE N-TH UNIT VECTOR\r
+           // IN SUCH A WAY THAT A SPIKE IS INTRODUCED INTO W.\r
+       \r
+           NM1 = N - 1;\r
+           if (NM1 >= 1) {\r
+               for (NMJ = 1; NMJ <= NM1; NMJ++) {\r
+                   J = N - NMJ;\r
+                   JJ = JJ - (M-J+1);\r
+                   W[J-1] = ZERO;\r
+                   if (V[J-1] == ZERO) {\r
+                       continue;\r
+                   }\r
+       \r
+                   // DETERMINE A GIVENS ROTATION WHICH ELIMINATES THE\r
+                   // J-TH ELEMENT OF V.\r
+       \r
+                 if (Math.abs(V[N-1]) < Math.abs(V[J-1])) {\r
+                     COTAN = V[N-1]/V[J-1];\r
+                     SIN = P5/Math.sqrt(P25+P25*COTAN*COTAN);\r
+                     COS = SIN*COTAN;\r
+                     TAU = ONE;\r
+                     if (Math.abs(COS)*GIANT > ONE) {\r
+                         TAU = ONE/COS;\r
+                     }\r
+                 } // if (Math.abs(V[N-1]) < Math.abs(V[J-1]))\r
+                 else {\r
+                     TAN = V[J-1]/V[N-1];\r
+                     COS = P5/Math.sqrt(P25+P25*TAN*TAN);\r
+                     SIN = COS*TAN;\r
+                     TAU = SIN;\r
+                 } // else\r
+       \r
+                 // APPLY THE TRANSFORMATION TO V AND STORE THE INFORMATION\r
+                 // NECESSARY TO RECOVER THE GIVENS ROTATION.\r
+       \r
+                 V[N-1] = SIN*V[J-1] + COS*V[N-1];\r
+                 V[J-1] = TAU;\r
+       \r
+                 // APPLY THE TRANSFORMATION TO S AND EXTEND THE SPIKE IN W.\r
+       \r
+                 L = JJ;\r
+                 for (I = J; I <= M; I++) {\r
+                     TEMP = COS*S[L-1] - SIN*W[I-1];\r
+                     W[I-1] = SIN*S[L-1] + COS*W[I-1];\r
+                     S[L-1] = TEMP;\r
+                     L = L + 1;\r
+                 } // for (I = J; I <= M; I++)\r
+               } // for (NMJ = 1; NMJ <= NM1; NMJ++)\r
+           } // if (NM1 >= 1)\r
+       \r
+           // ADD THE SPIKE FROM THE RANK 1 UPDATE TO W.\r
+       \r
+           for (I = 1; I <= M; I++) {\r
+               W[I-1] = W[I-1] + V[N-1]*U[I-1];\r
+           } // for (I = 1; I <= M; I++)\r
+       \r
+           // ELIMINATE THE SPIKE.\r
+       \r
+           SING[0] = false;\r
+           if (NM1 >= 1) {\r
+               for (J = 1; J <= NM1; J++) {\r
+                   if (W[J-1] != ZERO) {\r
+       \r
+                       // DETERMINE A GIVENS ROTATION WHICH ELIMINATES THE\r
+                       // J-TH ELEMENT OF THE SPIKE.\r
+       \r
+                       if (Math.abs(S[JJ-1]) < Math.abs(W[J-1])) {\r
+                           COTAN = S[JJ-1]/W[J-1];\r
+                           SIN = P5/Math.sqrt(P25+P25*COTAN*COTAN);\r
+                           COS = SIN*COTAN;\r
+                           TAU = ONE;\r
+                           if (Math.abs(COS)*GIANT > ONE) {\r
+                               TAU = ONE/COS;\r
+                           }\r
+                       } // if (Math.abs(S[JJ-1]) < Math.abs(W[J-1]))\r
+                       else {\r
+                           TAN = W[J-1]/S[JJ-1];\r
+                           COS = P5/Math.sqrt(P25+P25*TAN*TAN);\r
+                           SIN = COS*TAN;\r
+                           TAU = SIN;\r
+                       } // else\r
+       \r
+                    // APPLY THE TRANSFORMATION TO S AND REDUCE THE SPIKE IN W.\r
+       \r
+                       L = JJ;\r
+                       for (I = J; I <= M; I++) {\r
+                           TEMP = COS*S[L-1] + SIN*W[I-1];\r
+                           W[I-1] = -SIN*S[L-1] + COS*W[I-1];\r
+                           S[L-1] = TEMP;\r
+                           L = L + 1;\r
+                       } // for (I = J; I <= M; I++)\r
+       \r
+                       // STORE THE INFORMATION NECESSARY TO RECOVER THE\r
+                    // GIVENS ROTATION.\r
+       \r
+                       W[J-1] = TAU;\r
+                   } // if (W[J-1] != ZERO)\r
+       \r
+                   // TEST FOR ZERO DIAGONAL ELEMENTS IN THE OUTPUT S.\r
+       \r
+                   if (S[JJ-1] == ZERO) {\r
+                       SING[0] = true;\r
+                   }\r
+                   JJ = JJ + (M-J+1);\r
+               } // for (J = 1; J <= NM1; J++)\r
+           } // if (NM1 >= 1)\r
+       \r
+           // MOVE W BACK INTO THE LAST COLUMN OF THE OUTPUT S.\r
+       \r
+           L = JJ;\r
+           for (I = N; I <= M; I++) {\r
+               S[L-1] = W[I-1];\r
+               L = L + 1;\r
+           } // for (I = N; I <= M; I++) \r
+           if (S[JJ-1] == ZERO) {\r
+               SING[0] = true;\r
+           }\r
+           return;\r
+       } \r
+\r
+       \r
        private void GAUSSJ(int N, double ALPHA, double BETA, double B[], double T[], double W[]) {\r
            // B(N), T(N), W(N)\r
            //        THIS ROUTINE COMPUTES THE NODES T(J) AND WEIGHTS\r