Porting continues.
authorilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Tue, 21 Nov 2017 20:34:43 +0000 (20:34 +0000)
committerilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Tue, 21 Nov 2017 20:34:43 +0000 (20:34 +0000)
git-svn-id: https://citdcbmipav.cit.nih.gov/repos-pub/mipav/trunk@15268 ba61647d-9d00-f842-95cd-605cb4296b96

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

index 116aebb..8edddf6 100644 (file)
@@ -1021,6 +1021,29 @@ double neweps;
                  return;\r
        }\r
        \r
+       private double FMAX(int MN, double FVAL[]) {\r
+           //   -------------------------\r
+           //    DETERMINES THE MAXIMUM-NORM OF THE VECTOR FVAL.\r
+       \r
+           //     .. Array Arguments ..\r
+           // DOUBLE PRECISION FVAL(MN)\r
+               \r
+           //      .. Local Scalars ..\r
+               double maxVal;\r
+           double FV;\r
+           int K;\r
+       \r
+           maxVal = 0.0;\r
+           for (K = 0; K < MN; K++) {\r
+               FV = Math.abs(FVAL[K]);\r
+               if (FV > maxVal) {\r
+                       maxVal = FV;\r
+               }\r
+           } // for (K = 0; K < MN; K++)\r
+           return maxVal;\r
+       }\r
+\r
+       \r
        private void DSCFUN(int NDIM, double X[], double FVAL[],int IFLAG[]) {\r
            //    --------------------------------------\r
            //  FORMS THE NONLINEAR SYSTEM SATISFIED BY D-SC PARAMETERS.THE\r
@@ -1048,6 +1071,10 @@ double neweps;
             //.. Local Scalars ..\r
                double cr[] = new double[1];\r
                double ci[] = new double[1];\r
+           double wout[];\r
+           double base;\r
+           double diffR;\r
+           double diffI;\r
                double WA[] = new double[2];\r
                double WAI[] = new double[2];\r
                double WARC[] = new double[2];\r
@@ -1093,7 +1120,7 @@ double neweps;
              ICOUNT = ICOUNT + 1;\r
        \r
              //  TWO EQUATIONS TO ELIMINATE POSSIBLE ROTATION OF THE INNER POLYGON:\r
-             double wout[] = WQUAD(W12[N2-1],PHI12[N2-1],N2,1,W12[0],PHI12[0],1,\r
+             wout = WQUAD(W12[N2-1],PHI12[N2-1],N2,1,W12[0],PHI12[0],1,\r
                             1,U2[0],M2,N2,U2[0],W02,W12,ALFA02,ALFA12,NPTQ2,QWORK2,LINEARC2,2);\r
              scm.zmlt(C2[0], C2[1], wout[0], wout[1], cr, ci);\r
              WIN1[0] = Z12[0][0] - Z12[N2-1][0] - cr[0];\r
@@ -1112,110 +1139,137 @@ double neweps;
        \r
              //  TWO EQUATIONS TO FIX THE RELATIVE POSITION OF THE INNER POLYGON:\r
              TEST1 = Math.cos(PHI12[N2-1]);\r
-            /* if (TEST1 < U2[0]) {\r
+             if (TEST1 < U2[0]) {\r
        \r
                  // IF THE LINE PATH FROM W02[M2-1] TO W12[N2-1] IS OUT OF DOMAIN,THE\r
               //  COMBINATION OF TWO DIFFERENT PATHS WILL BE USED INSTEAD:\r
-             WX = DCMPLX(U2,0.D0)\r
-             WLINE = WQUAD(W02(M2),0.D0,M2,0,WX,0.D0,0,2,0.D0,M2,N2,U2,W02,W12,ALFA02,\r
-            +        ALFA12,NPTQ2,QWORK2,0,2)\r
-             IF (PHI12(N2).LE.0.D0) THEN\r
-                 WARC = WQUAD(W12(N2),PHI12(N2),N2,1,WX,0.D0,0,2,U2,M2,N2,U2,W02,W12,\r
-            +           ALFA02,ALFA12,NPTQ2,QWORK2,1,2)\r
-                 WIN2 = WLINE - WARC\r
-\r
-             ELSE\r
-                 WARC = WQUAD(WX,0.D0,0,2,W12(N2),PHI12(N2),N2,1,U2,M2,N2,U2,W02,W12,\r
-            +           ALFA02,ALFA12,NPTQ2,QWORK2,1,2)\r
-                 WIN2 = WLINE + WARC\r
-             END IF\r
-\r
-             GO TO 30\r
+                 WX[0] = U2[0];\r
+                 WX[1] = 0.0;\r
+                 WLINE = WQUAD(W02[M2-1],0.0,M2,0,WX,0.0,0,2,0.0,M2,N2,U2[0],W02,W12,ALFA02,\r
+                    ALFA12,NPTQ2,QWORK2,0,2);\r
+                 if (PHI12[N2-1] <= 0.0) {\r
+                     WARC = WQUAD(W12[N2-1],PHI12[N2-1],N2,1,WX,0.0,0,2,U2[0],M2,N2,U2[0],W02,W12,\r
+                       ALFA02,ALFA12,NPTQ2,QWORK2,1,2);\r
+                     WIN2[0] = WLINE[0] - WARC[0];\r
+                     WIN2[1] = WLINE[1] - WARC[1];\r
+                 }\r
+                 else {\r
+                     WARC = WQUAD(WX,0.0,0,2,W12[N2-1],PHI12[N2-1],N2,1,U2[0],M2,N2,U2[0],W02,W12,\r
+                       ALFA02,ALFA12,NPTQ2,QWORK2,1,2);\r
+                     WIN2[0] = WLINE[0] + WARC[0];\r
+                     WIN2[1] = WLINE[1] + WARC[1];\r
+                 }\r
              } // if (TEST1 < U2[0])\r
              else {\r
 \r
-                 WIN2 = WQUAD(W02(M2),0.D0,M2,0,W12(N2),0.D0,N2,1,0.D0,M2,N2,U2,W02,W12,ALFA02,\r
-            +       ALFA12,NPTQ2,QWORK2,0,2)\r
+                 WIN2 = WQUAD(W02[M2-1],0.0,M2,0,W12[N2-1],0.0,N2,1,0.0,M2,N2,U2[0],W02,W12,ALFA02,\r
+                   ALFA12,NPTQ2,QWORK2,0,2);\r
              }\r
-             FVAL(N2+2) = DIMAG(ZI* (Z12(N2)-Z02(M2)-C2*WIN2))\r
-             FVAL(N2+3) = DIMAG(Z12(N2)-Z02(M2)-C2*WIN2)\r
-       C\r
-       C  TWO EQUATIONS TO ELIMINATE POSSIBLE ROTATION OF THE OUTER POLYGON:\r
-             WIN3 = Z02(1) - Z02(M2) - C2*WQUAD(W02(M2),PHI02(M2),M2,0,W02(1),PHI02(1),1,\r
-            +       0,1.D0,M2,N2,U2,W02,W12,ALFA02,ALFA12,NPTQ2,QWORK2,LINEARC2,2)\r
-             FVAL(N2+4) = DIMAG(ZI*WIN3)\r
-             FVAL(N2+5) = DIMAG(WIN3)\r
-       C\r
-             IF (M2.EQ.3) THEN\r
-       C\r
-       C  CALCULATE THE MAXIMUM-NORM OF THE FUNCTION FVAL:\r
-                 IF (ISPRT.NE.1) GO TO 40\r
-                 FMAXN = FMAX(NDIM,FVAL)\r
-                 WRITE (6,FMT=9000) ICOUNT,FMAXN\r
-          40     CONTINUE\r
-                 RETURN\r
-\r
-             END IF\r
+             scm.zmlt(C2[0], C2[1], WIN2[0], WIN2[1], cr, ci);\r
+             diffR = Z12[N2-1][0] - Z02[M2-1][0] - cr[0];\r
+             diffI = Z12[N2-1][1] - Z02[M2-1][1] - ci[0];\r
+             scm.zmlt(ZI[0], ZI[1], diffR, diffI, cr, ci);\r
+             FVAL[N2+1] = ci[0];\r
+             FVAL[N2+2] = diffI;\r
+       \r
+             //  TWO EQUATIONS TO ELIMINATE POSSIBLE ROTATION OF THE OUTER POLYGON:\r
+             wout = WQUAD(W02[M2-1],PHI02[M2-1],M2,0,W02[0],PHI02[0],1,\r
+                                   0,1.0,M2,N2,U2[0],W02,W12,ALFA02,ALFA12,NPTQ2,QWORK2,LINEARC2,2);\r
+             scm.zmlt(C2[0], C2[1], wout[0], wout[1], cr, ci);\r
+             WIN3[0] = Z02[0][0] - Z02[M2-1][0] - cr[0];\r
+             WIN3[1] = Z02[0][1] - Z02[M2-1][1] - ci[0];\r
+             scm.zmlt(ZI[0],ZI[1], WIN3[0], WIN3[1], cr, ci);\r
+             FVAL[N2+3] = ci[0];\r
+             FVAL[N2+4] = WIN3[1];\r
+       \r
+             if (M2 == 3) {\r
+       \r
+                 // CALCULATE THE MAXIMUM-NORM OF THE FUNCTION FVAL:\r
+                 if (ISPRT == 1) {\r
+                     FMAXN = FMAX(NDIM,FVAL);\r
+                     System.out.println("# of iterations = " + ICOUNT);\r
+                     System.out.println("Maximum norm of the residuals = " + FMAXN);\r
+                 } // if (ISPRT == 1)\r
+                 return;\r
 \r
-             IF (ISHAPE2.EQ.1) GO TO 70\r
-       C\r
-       C  M-3 SIDE LENGTH CONDITIONS OF THE OUTER POLYGON:\r
-             DO 50 J = 1,M2 - 3\r
-                 WINT2 = WQUAD(W02(J),PHI02(J),J,0,W02(J+1),PHI02(J+1),J+1,0,1.D0,\r
-            +            M2,N2,U2,W02,W12,ALFA02,ALFA12,NPTQ2,QWORK2,LINEARC2,2)\r
-                 FVAL(N2+5+J) = ABS(Z02(J+1)-Z02(J)) - ABS(C2*WINT2)\r
-          50 CONTINUE\r
-       C\r
-       C  CALCULATE THE MAXIMUM-NORM OF THE FUNCTION FVAL:\r
-             IF (ISPRT.NE.1) GO TO 60\r
-             FMAXN = FMAX(NDIM,FVAL)\r
-             WRITE (6,FMT=9010) ICOUNT,FMAXN\r
-          60 CONTINUE\r
-             RETURN\r
-       C\r
-       C  OUTER POLYGON CONTAINS SOME INFINITE VERTICES & FOR EACH OF THEM\r
-       C  TWO LENGTH CONDITIONS WILL BE REPLACED BY A COMPLEX INTEGRAL:\r
-          70 DO 100 K = 1,NSHAPE - 1\r
-                 IF (IND(K+1).EQ.2 .OR. IND(K).GE.IND(K+1)-2) GO TO 90\r
-                 DO 80 J = IND(K) + 1,IND(K+1) - 2\r
-                     WINT3 = WQUAD(W02(J),PHI02(J),J,0,W02(J+1),PHI02(J+1),J+1,0,\r
-            +                1.D0,M2,N2,U2,W02,W12,ALFA02,ALFA12,NPTQ2,QWORK2,LINEARC2,2)\r
-                     FVAL(N+5+J) = ABS(Z02(J+1)-Z02(J)) - ABS(C2*WINT3)\r
-          80     CONTINUE\r
-          90     IF (K.EQ.NSHAPE-1 .OR. IND(K+1).EQ.M2-1) GO TO 100\r
-       C\r
-       C  THE COMBINATION  OF THREE DIFFERENT PATHS  IS USED TO INTEGRATE\r
-       C  FROM WA TO WB TO AVOID DOMAIN PROBLEM.THE BY-PRODUCT OF THIS IS\r
-       C  THAT IT IS NUMERICALLY  MORE EFFICIENT THAN A SINGLE LINE PATH:\r
-                 WA = W02(IND(K+1)-1)\r
-                 WB = W02(IND(K+1)+1)\r
-                 PHIA = PHI02(IND(K+1)-1)\r
-                 PHIB = PHI02(IND(K+1)+1)\r
-                 RADIUS = (1.D0+U2)/2.D0\r
-                 WAI = RADIUS*EXP(ZI*PHIA)\r
-                 WBI = RADIUS*EXP(ZI*PHIB)\r
-                 WLINE1 = WQUAD(WA,0.D0,IND(K+1)-1,0,WAI,0.D0,0,2,0.D0,M2,N2,U2,\r
-            +             W02,W12,ALFA02,ALFA12,NPTQ2,QWORK2,0,2)\r
-                 WLINE2 = WQUAD(WB,0.D0,IND(K+1)+1,0,WBI,0.D0,0,2,0.D0,M2,N2,U2,\r
-            +             W02,W12,ALFA02,ALFA12,NPTQ2,QWORK2,0,2)\r
-                 WCIRCLE = WQUAD(WAI,PHIA,0,2,WBI,PHIB,0,2,RADIUS,M2,N2,U2,W02,W12,\r
-            +              ALFA02,ALFA12,NPTQ2,QWORK2,1,2)\r
-                 WIN4 = C2* (WLINE1+WCIRCLE-WLINE2)\r
-                 FVAL(N2+5+IND(K+1)-1) = DIMAG(ZI*\r
-            +                           (Z02(IND(K+1)+1)-Z02(IND(K+1)-1)-WIN4))\r
-                 FVAL(N2+5+IND(K+1)) = DIMAG(Z02(IND(K+1)+1)-Z02(IND(K+1)-1)-WIN4)\r
-         100 CONTINUE\r
-       C\r
-       C  CALCULATE THE MAXIMUM-NORM OF THE FUNCTION FVAL:\r
-             IF (ISPRT.NE.1) GO TO 110\r
-             FMAXN = FMAX(NDIM,FVAL)\r
-             WRITE (6,FMT=9020) ICOUNT,FMAXN\r
-         110 CONTINUE\r
-             RETURN\r
+             } // if (M2 == 3)\r
 \r
-        9000 FORMAT (2X,I5,6X,E10.4)\r
-        9010 FORMAT (2X,I5,6X,E10.4)\r
-        9020 FORMAT (2X,I5,6X,E10.4)*/\r
+             if (ISHAPE2 != 1) {\r
+       \r
+                 //  M2-3 SIDE LENGTH CONDITIONS OF THE OUTER POLYGON:\r
+                 for (J = 1; J <= M2 - 3; J++) {\r
+                     WINT2 = WQUAD(W02[J-1],PHI02[J-1],J,0,W02[J],PHI02[J],J+1,0,1.0,\r
+                         M2,N2,U2[0],W02,W12,ALFA02,ALFA12,NPTQ2,QWORK2,LINEARC2,2);\r
+                     scm.zmlt(C2[0], C2[1], WINT2[0], WINT2[1], cr, ci);\r
+                     FVAL[N2+4+J] = scm.zabs(Z02[J][0]-Z02[J-1][0], Z02[J][1]-Z02[J-1][1]) - scm.zabs(cr[0], ci[0]);\r
+                 } // for (J = 1; J <= M2 - 3; J++)\r
+       \r
+                 //  CALCULATE THE MAXIMUM-NORM OF THE FUNCTION FVAL:\r
+                 if (ISPRT == 1) {\r
+                     FMAXN = FMAX(NDIM,FVAL);\r
+                     System.out.println("# of iterations = " + ICOUNT);\r
+                     System.out.println("Maximum norm of the residuals = " + FMAXN);\r
+                 } // if (ISPRT == 1)\r
+                 return;\r
+             } // if (ISHAPE2 != 1)\r
+       \r
+             //  OUTER POLYGON CONTAINS SOME INFINITE VERTICES & FOR EACH OF THEM\r
+             //  TWO LENGTH CONDITIONS WILL BE REPLACED BY A COMPLEX INTEGRAL:\r
+             for (K = 1; K <= NSHAPE - 1; K++) {\r
+                 if (IND[K] != 1 && IND[K-1] < IND[K]-2) {\r
+                     for (J = IND[K-1] + 2; J <= IND[K] - 1; J++) {\r
+                         WINT3 = WQUAD(W02[J-1],PHI02[J-1],J,0,W02[J],PHI02[J],J+1,0,\r
+                            1.0,M2,N2,U2[0],W02,W12,ALFA02,ALFA12,NPTQ2,QWORK2,LINEARC2,2);\r
+                         scm.zmlt(C2[0], C2[1], WINT3[0], WINT3[1], cr, ci);\r
+                         FVAL[N2+4+J] = scm.zabs(Z02[J][0]-Z02[J-1][0], Z02[J][1]-Z02[J-1][1]) - scm.zabs(cr[0], ci[0]);\r
+                     } // for (J = IND[K-1] + 2; J <= IND[K] - 1; J++) \r
+                 } // if (IND[K] != 1 && IND[K-1] < IND[K]-2)\r
+                 if (K == NSHAPE-1 || IND[K] == M2-2) {\r
+                         continue;\r
+                 }\r
+       \r
+                 //  THE COMBINATION  OF THREE DIFFERENT PATHS  IS USED TO INTEGRATE\r
+                 //  FROM WA TO WB TO AVOID DOMAIN PROBLEM.THE BY-PRODUCT OF THIS IS\r
+                 //  THAT IT IS NUMERICALLY  MORE EFFICIENT THAN A SINGLE LINE PATH:\r
+                 WA[0] = W02[IND[K]-1][0];\r
+                 WA[1] = W02[IND[K]-1][1];\r
+                 WB[0] = W02[IND[K]+1][0];\r
+                 WB[1] = W02[IND[K]+1][1];\r
+                 PHIA = PHI02[IND[K]-1];\r
+                 PHIB = PHI02[IND[K]+1];\r
+                 RADIUS = (1.0+U2[0])/2.0;\r
+                 scm.zmlt(ZI[0], ZI[1], PHIA, 0.0, cr, ci);\r
+                 base = RADIUS * Math.exp(cr[0]);\r
+                 WAI[0] = base * Math.cos(ci[0]);\r
+                 WAI[1] = base * Math.sin(ci[0]);\r
+                 scm.zmlt(ZI[0], ZI[1], PHIB, 0.0, cr, ci);\r
+                 base = RADIUS * Math.exp(cr[0]);\r
+                 WBI[0] = base * Math.cos(ci[0]);\r
+                 WBI[1] = base * Math.sin(ci[0]);\r
+                 WLINE1 = WQUAD(WA,0.0,IND[K]-1,0,WAI,0.0,0,2,0.0,M2,N2,U2[0],\r
+                            W02,W12,ALFA02,ALFA12,NPTQ2,QWORK2,0,2);\r
+                 WLINE2 = WQUAD(WB,0.0,IND[K]+1,0,WBI,0.0,0,2,0.0,M2,N2,U2[0],\r
+                            W02,W12,ALFA02,ALFA12,NPTQ2,QWORK2,0,2);\r
+                 WCIRCLE = WQUAD(WAI,PHIA,0,2,WBI,PHIB,0,2,RADIUS,M2,N2,U2[0],W02,W12,\r
+                          ALFA02,ALFA12,NPTQ2,QWORK2,1,2);\r
+                 scm.zmlt(C2[0], C2[1], (WLINE1[0]+WCIRCLE[0]-WLINE2[0]),\r
+                                 (WLINE1[1]+WCIRCLE[1]-WLINE2[1]), cr, ci);\r
+                 WIN4[0] = cr[0];\r
+                 WIN4[1] = ci[0];\r
+                 diffR = Z02[IND[K]+1][0]-Z02[IND[K]-1][0]-WIN4[0];\r
+                 diffI = Z02[IND[K]+1][1]-Z02[IND[K]-1][1]-WIN4[1];\r
+                 scm.zmlt(ZI[0], ZI[1], diffR, diffI, cr, ci);\r
+                 FVAL[N2+5+IND[K]-1] = ci[0];\r
+                 FVAL[N2+5+IND[K]] = diffI;\r
+             } // for (K = 1; K <= NSHAPE - 1; K++)\r
+       \r
+             //  CALCULATE THE MAXIMUM-NORM OF THE FUNCTION FVAL:\r
+             if (ISPRT == 1) {\r
+              FMAXN = FMAX(NDIM,FVAL);\r
+              System.out.println("# of iterations = " + ICOUNT);\r
+              System.out.println("Maximum norm of the residuals = " + FMAXN);\r
+          } // if (ISPRT == 1)\r
+          return;\r
        }\r
 \r
 \r
@@ -1789,7 +1843,7 @@ double neweps;
                //     ..\r
            //     .. Data statements ..\r
                //\r
-               ONE = 1.0;\r
+               /*ONE = 1.0;\r
                P1 = 1.0E-1;\r
                P5 = 5.0E-1;\r
                P001 = 1.0E-3;\r
@@ -1850,30 +1904,30 @@ double neweps;
                        } // if (NPRINT > 0)\r
                return;\r
                } // if (IFLAG[0] < 0)\r
-               /*      FNORM = ENORM(N,FVEC)\r
-               C\r
-               C     DETERMINE THE NUMBER OF CALLS TO FCN NEEDED TO COMPUTE\r
-               C     THE JACOBIAN MATRIX.\r
-               C\r
-                     MSUM = MIN0(ML+MU+1,N)\r
-               C\r
-               C     INITIALIZE ITERATION COUNTER AND MONITORS.\r
-               C\r
-                     ITER = 1\r
-                     NCSUC = 0\r
-                     NCFAIL = 0\r
-                     NSLOW1 = 0\r
-                     NSLOW2 = 0\r
-               C\r
-               C     BEGINNING OF THE OUTER LOOP.\r
-               C\r
-                  30 CONTINUE\r
-                     JEVAL = .TRUE.\r
-               C\r
-               C        CALCULATE THE JACOBIAN MATRIX.\r
-               C\r
-                     IFLAG = 2\r
-                     CALL FDJAC1(FCN,N,X,FVEC,FJAC,LDFJAC,IFLAG,ML,MU,EPSFCN,WA1,WA2)\r
+               FNORM = ENORM(N,FVEC);\r
+               \r
+               // DETERMINE THE NUMBER OF CALLS TO FCN NEEDED TO COMPUTE\r
+               // THE JACOBIAN MATRIX.\r
+               \r
+               MSUM = Math.min(ML+MU+1,N);\r
+               \r
+           // INITIALIZE ITERATION COUNTER AND MONITORS.\r
+               \r
+               ITER = 1;\r
+               NCSUC = 0;\r
+               NCFAIL = 0;\r
+               NSLOW1 = 0;\r
+               NSLOW2 = 0;\r
+               \r
+               // BEGINNING OF THE OUTER LOOP.\r
+               \r
+               while (true) {\r
+                   JEVAL = true;\r
+               \r
+            // CALCULATE THE JACOBIAN MATRIX.\r
+               \r
+                   IFLAG[0] = 2;\r
+                   FDJAC1(FCN,N,X,FVEC,FJAC,LDFJAC,IFLAG,ML,MU,EPSFCN,WA1,WA2);\r
                      NFEV = NFEV + MSUM\r
                      IF (IFLAG.LT.0) GO TO 300\r
                C\r
@@ -2096,7 +2150,7 @@ double neweps;
                C\r
                C        END OF THE OUTER LOOP.\r
                C\r
-                     GO TO 30\r
+               } // while (true)\r
 \r
                  300 CONTINUE\r
                C\r
@@ -2445,4 +2499,288 @@ double neweps;
          result = FAC/G;\r
          return result;\r
        }\r
+       \r
+       private double ENORM(int N, double X[]) {\r
+           //     **********\r
+           //\r
+           //     FUNCTION ENORM\r
+       \r
+           // GIVEN AN N-VECTOR X, THIS FUNCTION CALCULATES THE\r
+           // EUCLIDEAN NORM OF X.\r
+       \r
+           // THE EUCLIDEAN NORM IS COMPUTED BY ACCUMULATING THE SUM OF\r
+           // SQUARES IN THREE DIFFERENT SUMS. THE SUMS OF SQUARES FOR THE\r
+           // SMALL AND LARGE COMPONENTS ARE SCALED SO THAT NO OVERFLOWS\r
+           // OCCUR. NON-DESTRUCTIVE UNDERFLOWS ARE PERMITTED. UNDERFLOWS\r
+           // AND OVERFLOWS DO NOT OCCUR IN THE COMPUTATION OF THE UNSCALED\r
+           // SUM OF SQUARES FOR THE INTERMEDIATE COMPONENTS.\r
+           // THE DEFINITIONS OF SMALL, INTERMEDIATE AND LARGE COMPONENTS\r
+           // DEPEND ON TWO CONSTANTS, RDWARF AND RGIANT. THE MAIN\r
+           // RESTRICTIONS ON THESE CONSTANTS ARE THAT RDWARF**2 NOT\r
+           // UNDERFLOW AND RGIANT**2 NOT OVERFLOW. THE CONSTANTS\r
+           // GIVEN HERE ARE SUITABLE FOR EVERY KNOWN COMPUTER.\r
+       \r
+           // THE FUNCTION STATEMENT IS\r
+       \r
+           // DOUBLE PRECISION FUNCTION ENORM(N,X)\r
+       \r
+        // WHERE\r
+       \r
+           // N IS A POSITIVE INTEGER INPUT VARIABLE.\r
+       \r
+           // X IS AN INPUT ARRAY OF LENGTH N.\r
+       \r
+           // SUBPROGRAMS CALLED\r
+       \r
+           // FORTRAN-SUPPLIED ... DABS,DSQRT\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 N\r
+\r
+           // .. Array Arguments ..\r
+           // DOUBLE PRECISION X(N)\r
+\r
+           // .. Local Scalars ..\r
+               double result = 0.0;\r
+           double AGIANT,FLOATN,S1,S2,S3,X1MAX,\r
+                             X3MAX,XABS,ratio;\r
+           int I;\r
+       \r
+             // .. Data statements ..\r
+           double ONE = 1.0;\r
+           double ZERO = 0.0;\r
+           double RDWARF = 3.834E-20;\r
+           double RGIANT = 1.304E19;\r
+\r
+           S1 = ZERO;\r
+           S2 = ZERO;\r
+           S3 = ZERO;\r
+           X1MAX = ZERO;\r
+           X3MAX = ZERO;\r
+           FLOATN = N;\r
+           AGIANT = RGIANT/FLOATN;\r
+           for (I = 0; I < N; I++) {\r
+               XABS = Math.abs(X[I]);\r
+                 if (XABS <= RDWARF || XABS >= AGIANT) {\r
+                     if (XABS > RDWARF) {\r
+       \r
+                         // SUM FOR LARGE COMPONENTS.\r
+       \r
+                         if  (XABS > X1MAX) {\r
+                                 ratio = X1MAX/XABS;\r
+                             S1 = ONE + S1* ratio * ratio;\r
+                             X1MAX = XABS;\r
+                         }\r
+                         else {\r
+                                 ratio = XABS/X1MAX;\r
+                             S1 = S1 + ratio * ratio;\r
+                         }\r
+                         continue;\r
+\r
+                     } //  if (XABS > RDWARF) \r
+       \r
+                     // SUM FOR SMALL COMPONENTS.\r
+       \r
+                     if (XABS > X3MAX) {\r
+                         ratio = X3MAX/XABS;\r
+                         S3 = ONE + S3* ratio * ratio;\r
+                         X3MAX = XABS;\r
+                 }\r
+                 else {\r
+                     if (XABS != ZERO) {\r
+                         ratio = XABS/X3MAX;\r
+                         S3 = S3 + ratio * ratio;\r
+                     }\r
+                 }\r
+                 continue;\r
+\r
+                 } // if (XABS <= RDWARF || XABS >= AGIANT)\r
+       \r
+                 // SUM FOR INTERMEDIATE COMPONENTS.\r
+       \r
+                 S2 = S2 + XABS * XABS;\r
+           } // for (I = 0; I < N; I++)\r
+       \r
+        // CALCULATION OF NORM.\r
+       \r
+             if (S1 != ZERO) {\r
+                 result = X1MAX*Math.sqrt(S1+ (S2/X1MAX)/X1MAX);\r
+             }\r
+             else {\r
+                 if (S2 != ZERO) {\r
+                     if (S2 >= X3MAX) result = Math.sqrt(S2* (ONE+ (X3MAX/S2)* (X3MAX*S3)));\r
+                     if (S2 < X3MAX) result = Math.sqrt(X3MAX* ((S2/X3MAX)+ (X3MAX*S3)));\r
+                 }\r
+                 else {\r
+                     result = X3MAX*Math.sqrt(S3);\r
+                 } // else\r
+             } // else\r
+             return result;\r
+       }\r
+       \r
+       private void FDJAC1(int FCN,int N, double X[], double FVEC[], double FJAC[][], int LDFJAC,\r
+                       int IFLAG[] ,int ML,int MU,double EPSFCN, double WA1[], double WA2[]) {\r
+        // **********\r
+               \r
+               //     SUBROUTINE FDJAC1\r
+               \r
+               //     THIS SUBROUTINE COMPUTES A FORWARD-DIFFERENCE APPROXIMATION\r
+               //     TO THE N BY N JACOBIAN MATRIX ASSOCIATED WITH A SPECIFIED\r
+               //     PROBLEM OF N FUNCTIONS IN N VARIABLES. IF THE JACOBIAN HAS\r
+               //     A BANDED FORM, THEN FUNCTION EVALUATIONS ARE SAVED BY ONLY\r
+               //     APPROXIMATING THE NONZERO TERMS.\r
+               \r
+               //     THE SUBROUTINE STATEMENT IS\r
+               \r
+               //       SUBROUTINE FDJAC1(FCN,N,X,FVEC,FJAC,LDFJAC,IFLAG,ML,MU,EPSFCN,\r
+               //                         WA1,WA2)\r
+               \r
+               //     WHERE\r
+        //         dscfun = 1 for DSCFUN\r
+               //         FCN IS THE NAME OF THE USER-SUPPLIED SUBROUTINE WHICH\r
+               //         CALCULATES THE FUNCTIONS. FCN MUST BE DECLARED\r
+               //         IN AN EXTERNAL STATEMENT IN THE USER CALLING\r
+               //         PROGRAM, AND SHOULD BE WRITTEN AS FOLLOWS.\r
+               \r
+               //         SUBROUTINE FCN(N,X,FVEC,IFLAG)\r
+               //         INTEGER N,IFLAG\r
+               //         DOUBLE PRECISION X(N),FVEC(N)\r
+               //         ----------\r
+               //         CALCULATE THE FUNCTIONS AT X AND\r
+               //         RETURN THIS VECTOR IN FVEC.\r
+               //         ----------\r
+               //         RETURN\r
+               //         END\r
+               \r
+               //         THE VALUE OF IFLAG SHOULD NOT BE CHANGED BY FCN UNLESS\r
+               //         THE USER WANTS TO TERMINATE EXECUTION OF FDJAC1.\r
+               //         IN THIS CASE SET IFLAG TO A NEGATIVE INTEGER.\r
+               \r
+               //       N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER\r
+               //         OF FUNCTIONS AND VARIABLES.\r
+               \r
+               //       X IS AN INPUT ARRAY OF LENGTH N.\r
+               \r
+               //       FVEC IS AN INPUT ARRAY OF LENGTH N WHICH MUST CONTAIN THE\r
+               //         FUNCTIONS EVALUATED AT X.\r
+               \r
+               //       FJAC IS AN OUTPUT N BY N ARRAY WHICH CONTAINS THE\r
+               //         APPROXIMATION TO THE JACOBIAN MATRIX EVALUATED AT X.\r
+               \r
+               //       LDFJAC IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN N\r
+               //         WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY FJAC.\r
+               \r
+               //       IFLAG IS AN INTEGER VARIABLE WHICH CAN BE USED TO TERMINATE\r
+               //         THE EXECUTION OF FDJAC1. SEE DESCRIPTION OF FCN.\r
+               \r
+               //       ML IS A NONNEGATIVE INTEGER INPUT VARIABLE WHICH SPECIFIES\r
+               //         THE NUMBER OF SUBDIAGONALS WITHIN THE BAND OF THE\r
+               //         JACOBIAN MATRIX. IF THE JACOBIAN IS NOT BANDED, SET\r
+               //         ML TO AT LEAST N - 1.\r
+               \r
+               //       EPSFCN IS AN INPUT VARIABLE USED IN DETERMINING A SUITABLE\r
+               //         STEP LENGTH FOR THE FORWARD-DIFFERENCE APPROXIMATION. THIS\r
+               //         APPROXIMATION ASSUMES THAT THE RELATIVE ERRORS IN THE\r
+               //         FUNCTIONS ARE OF THE ORDER OF EPSFCN. IF EPSFCN IS LESS\r
+               //         THAN THE MACHINE PRECISION, IT IS ASSUMED THAT THE RELATIVE\r
+               //         ERRORS IN THE FUNCTIONS ARE OF THE ORDER OF THE MACHINE\r
+               //         PRECISION.\r
+               \r
+               //       MU IS A NONNEGATIVE INTEGER INPUT VARIABLE WHICH SPECIFIES\r
+               //         THE NUMBER OF SUPERDIAGONALS WITHIN THE BAND OF THE\r
+               //         JACOBIAN MATRIX. IF THE JACOBIAN IS NOT BANDED, SET\r
+               //         MU TO AT LEAST N - 1.\r
+               \r
+               //       WA1 AND WA2 ARE WORK ARRAYS OF LENGTH N. IF ML + MU + 1 IS AT\r
+               //         LEAST N, THEN THE JACOBIAN IS CONSIDERED DENSE, AND WA2 IS\r
+               //         NOT REFERENCED.\r
+               \r
+               //     SUBPROGRAMS CALLED\r
+               //\r
+               //       MINPACK-SUPPLIED ... DPMPAR\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
+               //      DOUBLE PRECISION EPSFCN\r
+               //      INTEGER IFLAG,LDFJAC,ML,MU,N\r
+               \r
+               //     .. Array Arguments ..\r
+               //      DOUBLE PRECISION FJAC(LDFJAC,N),FVEC(N),WA1(N),WA2(N),X(N)\r
+               \r
+               //     .. Subroutine Arguments ..\r
+               //      EXTERNAL FCN\r
+               \r
+               //     .. Local Scalars ..\r
+               double EPS,EPSMCH,H,TEMP,ZERO;\r
+               int I,J,K,MSUM;\r
+               \r
+               ZERO = 0.0;\r
+               \r
+               //     EPSMCH IS THE MACHINE PRECISION.\r
+               \r
+               EPSMCH = MACHEP;\r
+               \r
+               EPS = Math.sqrt(Math.max(EPSFCN,EPSMCH));\r
+               MSUM = ML + MU + 1;\r
+               if (MSUM >= N) {\r
+               \r
+                   // COMPUTATION OF DENSE APPROXIMATE JACOBIAN.\r
+               \r
+                   for (J = 0; J < N; J++) {\r
+                       TEMP = X[J];\r
+                       H = EPS*Math.abs(TEMP);\r
+                       if (H == ZERO) H = EPS;\r
+                       X[J] = TEMP + H;\r
+                       if (FCN == dscfun) {\r
+                           DSCFUN(N,X,WA1,IFLAG);\r
+                       }\r
+                       if (IFLAG[0] < 0) {\r
+                               return;\r
+                       }\r
+                       X[J] = TEMP;\r
+                       for (I = 0; I < N; I++) {\r
+                             FJAC[I][J] = (WA1[I]-FVEC[I])/H;\r
+                       } // for (I = 0; I < N; I++)\r
+                   } // for (J = 0; J < N; J++)\r
+                   return;\r
+\r
+               } // if (MSUM >= N)\r
+               \r
+               // COMPUTATION OF BANDED APPROXIMATE JACOBIAN.\r
+               \r
+               /*for (K = 1; K <= MSUM; K++) {\r
+                   for (J = K; J <= N; J += MSUM) {\r
+                       WA2[J-1] = X[J-1];\r
+                       H = EPS*Math.abs(WA2[J-1]);\r
+                       IF (H == ZERO) {\r
+                               H = EPS;\r
+                       }\r
+                       X[J-1] = WA2[J-1] + H\r
+                   } // for (J = K; J <= N; J += MSUM)\r
+                         CALL FCN(N,X,WA1,IFLAG)\r
+                         IF (IFLAG.LT.0) GO TO 90\r
+                         DO 70 J = K,N,MSUM\r
+                             X(J) = WA2(J)\r
+                             H = EPS*DABS(WA2(J))\r
+                             IF (H.EQ.ZERO) H = EPS\r
+                             DO 60 I = 1,N\r
+                                 FJAC(I,J) = ZERO\r
+                                 IF (I.GE.J-MU .AND. I.LE.J+ML) FJAC(I,\r
+                    +                J) = (WA1(I)-FVEC(I))/H\r
+                  60         CONTINUE\r
+                  70     CONTINUE\r
+               } // for (K = 1; K <= MSUM; K++)\r
+                  90 CONTINUE\r
+                 100 CONTINUE*/\r
+                     return;\r
+       }\r
+\r
+\r
 }
\ No newline at end of file