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

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

index 8edddf6..78cef3b 100644 (file)
@@ -1822,8 +1822,9 @@ double neweps;
                //      EXTERNAL FCN\r
                //     ..\r
                //     .. Local Scalars ..\r
-                     double ACTRED,DELTA,EPSMCH,FNORM,FNORM1,ONE,P0001,P001,\r
+                     double ACTRED,EPSMCH,FNORM,FNORM1,ONE,P0001,P001,\r
                                      P1,P5,PNORM,PRERED,RATIO,SUM,TEMP,XNORM,ZERO;\r
+                     double DELTA = 0.0;\r
                      int IFLAG[] = new int[1];\r
                      int I,ITER,J,JM1,L,MSUM,NCFAIL,NCSUC,NSLOW1,NSLOW2;\r
                      boolean JEVAL,SING;\r
@@ -1843,7 +1844,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
@@ -1921,129 +1922,150 @@ double neweps;
                \r
                // BEGINNING OF THE OUTER LOOP.\r
                \r
-               while (true) {\r
+               outer: 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
-               C        COMPUTE THE QR FACTORIZATION OF THE JACOBIAN.\r
-               C\r
-                     CALL QRFAC(N,N,FJAC,LDFJAC,.FALSE.,IWA,1,WA1,WA2,WA3)\r
-               C\r
-               C        ON THE FIRST ITERATION AND IF MODE IS 1, SCALE ACCORDING\r
-               C        TO THE NORMS OF THE COLUMNS OF THE INITIAL JACOBIAN.\r
-               C\r
-                     IF (ITER.NE.1) GO TO 70\r
-                     IF (MODE.EQ.2) GO TO 50\r
-                     DO 40 J = 1,N\r
-                         DIAG(J) = WA2(J)\r
-                         IF (WA2(J).EQ.ZERO) DIAG(J) = ONE\r
-                  40 CONTINUE\r
-                  50 CONTINUE\r
-               C\r
-               C        ON THE FIRST ITERATION, CALCULATE THE NORM OF THE SCALED X\r
-               C        AND INITIALIZE THE STEP BOUND DELTA.\r
-               C\r
-                     DO 60 J = 1,N\r
-                         WA3(J) = DIAG(J)*X(J)\r
-                  60 CONTINUE\r
-                     XNORM = ENORM(N,WA3)\r
-                     DELTA = FACTOR*XNORM\r
-                     IF (DELTA.EQ.ZERO) DELTA = FACTOR\r
-                  70 CONTINUE\r
-               C\r
-               C        FORM (Q TRANSPOSE)*FVEC AND STORE IN QTF.\r
-               C\r
-                     DO 80 I = 1,N\r
-                         QTF(I) = FVEC(I)\r
-                  80 CONTINUE\r
-                     DO 120 J = 1,N\r
-                         IF (FJAC(J,J).EQ.ZERO) GO TO 110\r
-                         SUM = ZERO\r
-                         DO 90 I = J,N\r
-                             SUM = SUM + FJAC(I,J)*QTF(I)\r
-                  90     CONTINUE\r
-                         TEMP = -SUM/FJAC(J,J)\r
-                         DO 100 I = J,N\r
-                             QTF(I) = QTF(I) + FJAC(I,J)*TEMP\r
-                 100     CONTINUE\r
-                 110     CONTINUE\r
-                 120 CONTINUE\r
-               C\r
-               C        COPY THE TRIANGULAR FACTOR OF THE QR FACTORIZATION INTO R.\r
-               C\r
-                     SING = .FALSE.\r
-                     DO 150 J = 1,N\r
-                         L = J\r
-                         JM1 = J - 1\r
-                         IF (JM1.LT.1) GO TO 140\r
-                         DO 130 I = 1,JM1\r
-                             R(L) = FJAC(I,J)\r
-                             L = L + N - I\r
-                 130     CONTINUE\r
-                 140     CONTINUE\r
-                         R(L) = WA1(J)\r
-                         IF (WA1(J).EQ.ZERO) SING = .TRUE.\r
-                 150 CONTINUE\r
-               C\r
-               C        ACCUMULATE THE ORTHOGONAL FACTOR IN FJAC.\r
-               C\r
-                     CALL QFORM(N,N,FJAC,LDFJAC,WA1)\r
-               C\r
-               C        RESCALE IF NECESSARY.\r
-               C\r
-                     IF (MODE.EQ.2) GO TO 170\r
-                     DO 160 J = 1,N\r
-                         DIAG(J) = DMAX1(DIAG(J),WA2(J))\r
-                 160 CONTINUE\r
-                 170 CONTINUE\r
-               C\r
-               C        BEGINNING OF THE INNER LOOP.\r
-               C\r
-                 180 CONTINUE\r
-               C\r
-               C           IF REQUESTED, CALL FCN TO ENABLE PRINTING OF ITERATES.\r
-               C\r
-                     IF (NPRINT.LE.0) GO TO 190\r
-                     IFLAG = 0\r
-                     IF (MOD(ITER-1,NPRINT).EQ.0) CALL FCN(N,X,FVEC,IFLAG)\r
-                     IF (IFLAG.LT.0) GO TO 300\r
-                 190 CONTINUE\r
-               C\r
-               C           DETERMINE THE DIRECTION P.\r
-               C\r
-                     CALL DOGLEG(N,R,LR,DIAG,QTF,DELTA,WA1,WA2,WA3)\r
-               C\r
-               C           STORE THE DIRECTION P AND X + P. CALCULATE THE NORM OF P.\r
-               C\r
-                     DO 200 J = 1,N\r
-                         WA1(J) = -WA1(J)\r
-                         WA2(J) = X(J) + WA1(J)\r
-                         WA3(J) = DIAG(J)*WA1(J)\r
-                 200 CONTINUE\r
-                     PNORM = ENORM(N,WA3)\r
-               C\r
-               C           ON THE FIRST ITERATION, ADJUST THE INITIAL STEP BOUND.\r
-               C\r
-                     IF (ITER.EQ.1) DELTA = DMIN1(DELTA,PNORM)\r
-               C\r
-               C           EVALUATE THE FUNCTION AT X + P AND CALCULATE ITS NORM.\r
-               C\r
-                     IFLAG = 1\r
-                     CALL FCN(N,WA2,WA4,IFLAG)\r
-                     NFEV = NFEV + 1\r
-                     IF (IFLAG.LT.0) GO TO 300\r
-                     FNORM1 = ENORM(N,WA4)\r
-               C\r
-               C           COMPUTE THE SCALED ACTUAL REDUCTION.\r
-               C\r
-                     ACTRED = -ONE\r
-                     IF (FNORM1.LT.FNORM) ACTRED = ONE - (FNORM1/FNORM)**2\r
+                   NFEV[0] = NFEV[0] + MSUM;\r
+                   if (IFLAG[0] < 0) {\r
+                       break;\r
+                   }\r
+               \r
+                   // COMPUTE THE QR FACTORIZATION OF THE JACOBIAN.\r
+               \r
+                   QRFAC(N,N,FJAC,LDFJAC,false,IWA,1,WA1,WA2,WA3);\r
+               \r
+                   // ON THE FIRST ITERATION AND IF MODE IS 1, SCALE ACCORDING\r
+                   // TO THE NORMS OF THE COLUMNS OF THE INITIAL JACOBIAN.\r
+               \r
+                   if (ITER == 1) {\r
+                       if (MODE != 2) {\r
+                           for (J = 0; J < N; J++) {\r
+                               DIAG[J] = WA2[J];\r
+                               if (WA2[J] == ZERO){ \r
+                                       DIAG[J] = ONE;  \r
+                               }\r
+                           } // for (J = 0; J < N; J++)\r
+                       } // if (MODE != 2)\r
+               \r
+                       // ON THE FIRST ITERATION, CALCULATE THE NORM OF THE SCALED X\r
+                       // AND INITIALIZE THE STEP BOUND DELTA.\r
+               \r
+                       for (J = 0; J < N; J++) {\r
+                           WA3[J] = DIAG[J]*X[J];\r
+                       } // for (J = 0; J < N; J++)\r
+                       XNORM = ENORM(N,WA3);\r
+                       DELTA = FACTOR*XNORM;\r
+                       if (DELTA == ZERO) {\r
+                               DELTA = FACTOR;\r
+                       }\r
+                   } // if (ITER == 1)\r
+               \r
+                   // FORM (Q TRANSPOSE)*FVEC AND STORE IN QTF.\r
+               \r
+                   for (I = 0; I < N; I++) {\r
+                         QTF[I] = FVEC[I];\r
+                   } // for (I = 0; I < N; I++)\r
+                   for (J = 0; J < N; J++) {\r
+                       if (FJAC[J][J] == ZERO) {\r
+                               continue;\r
+                       }\r
+                       SUM = ZERO;\r
+                       for (I = J; I < N; I++) {\r
+                           SUM = SUM + FJAC[I][J]*QTF[I];\r
+                       } // for (I = J; I < N; I++)\r
+                       TEMP = -SUM/FJAC[J][J];\r
+                       for (I = J; I < N; I++) {\r
+                           QTF[I] = QTF[I] + FJAC[I][J]*TEMP;\r
+                       } // for (I = J; I < N; I++)\r
+                   } // for (J = 0; J < N; J++)\r
+               \r
+                   // COPY THE TRIANGULAR FACTOR OF THE QR FACTORIZATION INTO R.\r
+               \r
+                     SING = false;\r
+                     for (J = 1; J <= N; J++) {\r
+                         L = J;\r
+                         JM1 = J - 1;\r
+                         if (JM1 >= 1) {\r
+                             for (I = 1; I <= JM1; I++) {\r
+                                 R[L-1] = FJAC[I-1][J-1];\r
+                                 L = L + N - I;\r
+                             } // for (I = 1; I <= JM1; I++)\r
+                         } // if (JM1 >= 1) \r
+                         R[L-1] = WA1[J-1];\r
+                         if (WA1[J-1] == ZERO) {\r
+                                 SING = true;\r
+                         }\r
+                     } // for (J = 1; J <= N; J++)\r
+               \r
+                 // ACCUMULATE THE ORTHOGONAL FACTOR IN FJAC.\r
+               \r
+                     QFORM(N,N,FJAC,LDFJAC,WA1);\r
+               \r
+                     // RESCALE IF NECESSARY.\r
+               \r
+                    if (MODE != 2) {\r
+                        for (J = 0; J < N; J++) {\r
+                             DIAG[J] = Math.max(DIAG[J],WA2[J]);\r
+                        } // for (J = 0; J < N; J++)\r
+                    } // if (MODE != 2)\r
+               \r
+                    // BEGINNING OF THE INNER LOOP.\r
+               \r
+                    while (true) {\r
+               \r
+                        // IF REQUESTED, CALL FCN TO ENABLE PRINTING OF ITERATES.\r
+               \r
+                        if  (NPRINT > 0) {\r
+                            IFLAG[0] = 0;\r
+                            if ((ITER-1%NPRINT) == 0) {\r
+                                if (FCN == dscfun) {\r
+                                    DSCFUN(N,X,FVEC,IFLAG);\r
+                                }\r
+                            }\r
+                            if (IFLAG[0] < 0) {\r
+                                break outer;\r
+                            }\r
+                        } // if (NPRINT > 0)\r
+               \r
+                        // DETERMINE THE DIRECTION P.\r
+               \r
+                        DOGLEG(N,R,LR,DIAG,QTF,DELTA,WA1,WA2,WA3);\r
+               \r
+                        // STORE THE DIRECTION P AND X + P. CALCULATE THE NORM OF P.\r
+               \r
+                        for (J = 0; J < N; J++) {\r
+                            WA1[J] = -WA1[J];\r
+                            WA2[J] = X[J] + WA1[J];\r
+                            WA3[J] = DIAG[J]*WA1[J];\r
+                        } // for (J = 0; J < N; J++)\r
+                        PNORM = ENORM(N,WA3);\r
+               \r
+                        // ON THE FIRST ITERATION, ADJUST THE INITIAL STEP BOUND.\r
+               \r
+                        if (ITER == 1) {\r
+                                DELTA = Math.min(DELTA,PNORM);\r
+                        }\r
+               \r
+                        // EVALUATE THE FUNCTION AT X + P AND CALCULATE ITS NORM.\r
+               \r
+                        IFLAG[0] = 1;\r
+                        if (FCN == dscfun) {\r
+                            DSCFUN(N,WA2,WA4,IFLAG);\r
+                        }\r
+                        NFEV[0] = NFEV[0] + 1;\r
+                        if (IFLAG[0] < 0) {\r
+                                break outer;\r
+                        }\r
+                        FNORM1 = ENORM(N,WA4);\r
+               \r
+                        // 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
@@ -2139,28 +2161,221 @@ double neweps;
                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
-               C\r
-               C           END OF THE INNER LOOP.\r
-               C\r
-                     JEVAL = .FALSE.\r
-                     GO TO 180\r
-\r
-                 290 CONTINUE\r
-               C\r
-               C        END OF THE OUTER LOOP.\r
-               C\r
+                     CALL R1MPYQ(1,N,QTF,1,WA2,WA3)*/\r
+               \r
+                        // END OF THE INNER LOOP.\r
+               \r
+                         JEVAL = false;\r
+                    } // while (true)\r
+               \r
+                   // END OF THE OUTER LOOP.\r
+               \r
                } // while (true)\r
 \r
-                 300 CONTINUE\r
-               C\r
-               C     TERMINATION, EITHER NORMAL OR USER IMPOSED.\r
-               C\r
-                     IF (IFLAG.LT.0) INFO = IFLAG\r
-                     IFLAG = 0\r
-                     IF (NPRINT.GT.0) CALL FCN(N,X,FVEC,IFLAG)\r
-               return;*/\r
+               //     TERMINATION, EITHER NORMAL OR USER IMPOSED.\r
+               \r
+               if (IFLAG[0] < 0) {\r
+                       INFO[0] = IFLAG[0];\r
+               }\r
+               IFLAG[0] = 0;\r
+               if (NPRINT > 0) {\r
+                       if (FCN == dscfun) {\r
+                           DSCFUN(N,X,FVEC,IFLAG);\r
+                       }\r
+               }\r
+               return;\r
       }\r
+       \r
+       private void DOGLEG(int N, double R[], int LR, double DIAG[], double QTB[],\r
+                       double DELTA, double X[], double WA1[], double WA2[]) {\r
+           //     **********\r
+       \r
+           // SUBROUTINE DOGLEG\r
+       \r
+           // GIVEN AN M BY N MATRIX A, AN N BY N NONSINGULAR DIAGONAL\r
+           // MATRIX D, AN M-VECTOR B, AND A POSITIVE NUMBER DELTA, THE\r
+           // PROBLEM IS TO DETERMINE THE CONVEX COMBINATION X OF THE\r
+           // GAUSS-NEWTON AND SCALED GRADIENT DIRECTIONS THAT MINIMIZES\r
+           // (A*X - B) IN THE LEAST SQUARES SENSE, SUBJECT TO THE\r
+        // RESTRICTION THAT THE EUCLIDEAN NORM OF D*X BE AT MOST DELTA.\r
+       \r
+        // THIS SUBROUTINE COMPLETES THE SOLUTION OF THE PROBLEM\r
+           // IF IT IS PROVIDED WITH THE NECESSARY INFORMATION FROM THE\r
+           // QR FACTORIZATION OF A. THAT IS, IF A = Q*R, WHERE Q HAS\r
+           // ORTHOGONAL COLUMNS AND R IS AN UPPER TRIANGULAR MATRIX,\r
+           // THEN DOGLEG EXPECTS THE FULL UPPER TRIANGLE OF R AND\r
+           // THE FIRST N COMPONENTS OF (Q TRANSPOSE)*B.\r
+       \r
+           // THE SUBROUTINE STATEMENT IS\r
+       \r
+           // SUBROUTINE DOGLEG(N,R,LR,DIAG,QTB,DELTA,X,WA1,WA2)\r
+       \r
+           // WHERE\r
+       \r
+           // N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE ORDER OF R.\r
+       \r
+           // R IS AN INPUT ARRAY OF LENGTH LR WHICH MUST CONTAIN THE UPPER\r
+           //   TRIANGULAR MATRIX R STORED BY ROWS.\r
+       \r
+           // LR IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN\r
+           //    (N*(N+1))/2.\r
+       \r
+           // DIAG IS AN INPUT ARRAY OF LENGTH N WHICH MUST CONTAIN THE\r
+           //   DIAGONAL ELEMENTS OF THE MATRIX D.\r
+       \r
+           // QTB IS AN INPUT ARRAY OF LENGTH N WHICH MUST CONTAIN THE FIRST\r
+           //   N ELEMENTS OF THE VECTOR (Q TRANSPOSE)*B.\r
+       \r
+           // DELTA IS A POSITIVE INPUT VARIABLE WHICH SPECIFIES AN UPPER\r
+           //   BOUND ON THE EUCLIDEAN NORM OF D*X.\r
+       \r
+           // X IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS THE DESIRED\r
+           //   CONVEX COMBINATION OF THE GAUSS-NEWTON DIRECTION AND THE\r
+           //   SCALED GRADIENT DIRECTION.\r
+       \r
+           // WA1 AND WA2 ARE WORK ARRAYS OF LENGTH N.\r
+       \r
+           // SUBPROGRAMS CALLED\r
+       \r
+           // MINPACK-SUPPLIED ... DPMPAR,ENORM\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 DELTA\r
+           // INTEGER LR,N\r
+       \r
+           //.. Array Arguments ..\r
+           // DOUBLE PRECISION DIAG(N),QTB(N),R(LR),WA1(N),WA2(N),X(N)\r
+       \r
+           // .. Local Scalars ..\r
+           double ALPHA,BNORM,EPSMCH,GNORM,ONE,QNORM,SGNORM,SUM,TEMP,ZERO;\r
+           double ratio, ratio2, ratio3;\r
+           int I,J,JJ,JP1,K,L;\r
+\r
+           // .. External Functions ..\r
+           // DOUBLE PRECISION DPMPAR,ENORM\r
+           // EXTERNAL DPMPAR,ENORM\r
+       \r
+           // .. Data statements ..\r
+           ONE = 1.0;\r
+           ZERO = 0.0;\r
+       \r
+           // EPSMCH IS THE MACHINE PRECISION.\r
+       \r
+             EPSMCH = MACHEP;\r
+       \r
+             // FIRST, CALCULATE THE GAUSS-NEWTON DIRECTION.\r
+       \r
+             JJ = (N* (N+1))/2 + 1;\r
+             for (K = 1; K <= N; K++) {\r
+                 J = N - K + 1;\r
+                 JP1 = J + 1;\r
+                 JJ = JJ - K;\r
+                 L = JJ + 1;\r
+                 SUM = ZERO;\r
+                 if (N >= JP1) {\r
+                     for (I = JP1; I <= N; I++) {\r
+                         SUM = SUM + R[L-1]*X[I-1];\r
+                         L = L + 1;\r
+                     } // for (I = JP1; I <= N; I++)\r
+                 } // if (N >= JP1)\r
+                 TEMP = R[JJ-1];\r
+                 if (TEMP == ZERO) {\r
+                     L = J;\r
+                     for (I = 1; I <= J; I++) {\r
+                         TEMP = Math.max(TEMP,Math.abs(R[L-1]));\r
+                         L = L + N - I;\r
+                     } // for (I = 1; I <= J; I++)\r
+                     TEMP = EPSMCH*TEMP;\r
+                     if (TEMP == ZERO) {\r
+                         TEMP = EPSMCH;\r
+                     }\r
+                 } // if (TEMP == ZERO)\r
+                 X[J-1] = (QTB[J-1]-SUM)/TEMP;\r
+             } // for (K = 1; K <= N; K++)\r
+       \r
+             // TEST WHETHER THE GAUSS-NEWTON DIRECTION IS ACCEPTABLE.\r
+       \r
+             for (J = 0; J < N; J++) {\r
+                 WA1[J] = ZERO;\r
+                 WA2[J] = DIAG[J]*X[J];\r
+             } // for (J = 0; J < N; J++)\r
+             QNORM = ENORM(N,WA2);\r
+             if (QNORM <= DELTA) {\r
+                 return;\r
+             }\r
+       \r
+             // THE GAUSS-NEWTON DIRECTION IS NOT ACCEPTABLE.\r
+             // NEXT, CALCULATE THE SCALED GRADIENT DIRECTION.\r
+       \r
+             L = 0;\r
+             for (J = 0; J < N; J++) {\r
+                 TEMP = QTB[J];\r
+                 for (I = J; I < N; I++) {\r
+                     WA1[I] = WA1[I] + R[L]*TEMP;\r
+                     L = L + 1;\r
+                 } // for (I = J; I < N; I++)\r
+                 WA1[J] = WA1[J]/DIAG[J];\r
+             } // for (J = 0; J < N; J++)\r
+       \r
+             // CALCULATE THE NORM OF THE SCALED GRADIENT AND TEST FOR\r
+             // THE SPECIAL CASE IN WHICH THE SCALED GRADIENT IS ZERO.\r
+       \r
+             GNORM = ENORM(N,WA1);\r
+             SGNORM = ZERO;\r
+             ALPHA = DELTA/QNORM;\r
+             if (GNORM != ZERO) {\r
+       \r
+                 // CALCULATE THE POINT ALONG THE SCALED GRADIENT\r
+                 // AT WHICH THE QUADRATIC IS MINIMIZED.\r
+       \r
+                 for (J = 0; J < N; J++) {\r
+                     WA1[J] = (WA1[J]/GNORM)/DIAG[J];\r
+                 } // for (J = 0; J < N; J++)\r
+                 L = 0;\r
+                 for (J = 0; J < N; J++) {\r
+                     SUM = ZERO;\r
+                     for (I = J; I < N; I++) {\r
+                         SUM = SUM + R[L]*WA1[I];\r
+                         L = L + 1;\r
+                     } // for (I = J; I < N; I++)\r
+                     WA2[J] = SUM;\r
+                 } // for (J = 0; J < N; J++)\r
+                 TEMP = ENORM(N,WA2);\r
+                 SGNORM = (GNORM/TEMP)/TEMP;\r
+       \r
+                 // TEST WHETHER THE SCALED GRADIENT DIRECTION IS ACCEPTABLE.\r
+       \r
+                 ALPHA = ZERO;\r
+                 if  (SGNORM < DELTA) {\r
+       \r
+                     // THE SCALED GRADIENT DIRECTION IS NOT ACCEPTABLE.\r
+                     // FINALLY, CALCULATE THE POINT ALONG THE DOGLEG\r
+                     // AT WHICH THE QUADRATIC IS MINIMIZED.\r
+       \r
+                     BNORM = ENORM(N,QTB);\r
+                     TEMP = (BNORM/GNORM)* (BNORM/QNORM)* (SGNORM/DELTA);\r
+                     ratio = SGNORM/DELTA;\r
+                     ratio2 = TEMP - (DELTA/QNORM);\r
+                     ratio3 = DELTA/QNORM;\r
+                     TEMP = TEMP - (DELTA/QNORM)* ratio*ratio +\r
+                         Math.sqrt(ratio2*ratio2+(ONE- ratio3*ratio3)* (ONE- ratio*ratio));\r
+                     ALPHA = ((DELTA/QNORM)* (ONE- ratio*ratio))/TEMP;\r
+                 } // if  (SGNORM < DELTA)\r
+             } // if (GNORM != ZERO)\r
+       \r
+             // FORM APPROPRIATE CONVEX COMBINATION OF THE GAUSS-NEWTON\r
+             // DIRECTION AND THE SCALED GRADIENT DIRECTION.\r
+       \r
+             TEMP = (ONE-ALPHA)*Math.min(SGNORM,DELTA);\r
+             for (J = 0; J < N; J++) {\r
+                 X[J] = TEMP*WA1[J] + ALPHA*X[J];\r
+             } // for (J = 0; J < N; J++)\r
+             return;\r
+    }\r
 \r
        \r
        private void QINIT(int M, int N, double ALFA0[], double ALFA1[], int NPTQ, double QWORK[]) {\r
@@ -2755,32 +2970,336 @@ double neweps;
                \r
                // COMPUTATION OF BANDED APPROXIMATE JACOBIAN.\r
                \r
-               /*for (K = 1; K <= MSUM; K++) {\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
+                       if (H == ZERO) {\r
+                               H = EPS;\r
+                       }\r
+                       X[J-1] = WA2[J-1] + H;\r
+                   } // for (J = K; J <= N; J += MSUM)\r
+                   if (FCN == dscfun) {\r
+                       DSCFUN(N,X,WA1,IFLAG);\r
+                   }\r
+                   if (IFLAG[0] < 0) {\r
+                       return;\r
+                   }\r
+                   for (J = K; J <= N; J += MSUM) {\r
+                       X[J-1] = WA2[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 (I = 1; I <= N; I++) {\r
+                           FJAC[I-1][J-1] = ZERO;\r
+                           if (I >= J-MU && I <= J+ML) {\r
+                               FJAC[I-1][J-1] = (WA1[I-1] - FVEC[I-1])/H;\r
+                           }\r
+                       } // for (I = 1; I <= N; I++)\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
+               return;\r
        }\r
+       \r
+       private void QFORM(int M,int N, double Q[][], int LDQ, double WA[]) {\r
+           //     **********\r
+       \r
+           // SUBROUTINE QFORM\r
+       \r
+           // THIS SUBROUTINE PROCEEDS FROM THE COMPUTED QR FACTORIZATION OF\r
+           // AN M BY N MATRIX A TO ACCUMULATE THE M BY M ORTHOGONAL MATRIX\r
+           // Q FROM ITS FACTORED FORM.\r
+       \r
+           // THE SUBROUTINE STATEMENT IS\r
+       \r
+           // SUBROUTINE QFORM(M,N,Q,LDQ,WA)\r
+       \r
+           // WHERE\r
+       \r
+           // M IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER\r
+           //   OF ROWS OF A AND THE ORDER OF Q.\r
+       \r
+           // N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER\r
+           //   OF COLUMNS OF A.\r
+       \r
+           // Q IS AN M BY M ARRAY. ON INPUT THE FULL LOWER TRAPEZOID IN\r
+        //   THE FIRST MIN(M,N) COLUMNS OF Q CONTAINS THE FACTORED FORM.\r
+           //   ON OUTPUT Q HAS BEEN ACCUMULATED INTO A SQUARE MATRIX.\r
+       \r
+           // LDQ IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN M\r
+           //   WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY Q.\r
+       \r
+           // WA IS A WORK ARRAY OF LENGTH M.\r
+       \r
+           // SUBPROGRAMS CALLED\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 LDQ,M,N\r
+       \r
+           // .. Array Arguments ..\r
+           // DOUBLE PRECISION Q(LDQ,M),WA(M)\r
+\r
+           //      .. Local Scalars ..\r
+           double ONE,SUM,TEMP,ZERO;\r
+           int I,J,JM1,K,L,MINMN,NP1;\r
+       \r
+           //     .. Data statements ..\r
+           ONE = 1.0;\r
+           ZERO = 0.0;\r
+       \r
+           // ZERO OUT UPPER TRIANGLE OF Q IN THE FIRST MIN(M,N) COLUMNS.\r
+       \r
+           MINMN = Math.min(M,N);\r
+           if (MINMN >= 2) {\r
+               for (J = 2; J <= MINMN; J++) {\r
+                 JM1 = J - 1;\r
+                 for (I = 1; I <= JM1; I++) {\r
+                     Q[I-1][J-1] = ZERO;\r
+                 } // for (I = 1; I <= JM1; I++)\r
+               } // for (J = 2; J <= MINMN; J++)\r
+           } // if (MINMN >= 2) \r
+       \r
+        // INITIALIZE REMAINING COLUMNS TO THOSE OF THE IDENTITY MATRIX.\r
+       \r
+           NP1 = N + 1;\r
+           if (M >= NP1) {\r
+               for (J = NP1; J <= M; J++) {\r
+                   for (I = 1; I <= M; I++) {\r
+                     Q[I-1][J-1] = ZERO;\r
+                   } // for (I = 1; I <= M; I++)\r
+                   Q[J-1][J-1] = ONE;\r
+               } // for (J = NP1; J <= M; J++)\r
+           } // if (M >= NP1)\r
+       \r
+           // ACCUMULATE Q FROM ITS FACTORED FORM.\r
+       \r
+           for (L = 1; L <= MINMN; L++) {\r
+               K = MINMN - L + 1;\r
+               for (I = K; I <= M; I++) {\r
+                   WA[I-1] = Q[I-1][K-1];\r
+                   Q[I-1][K-1] = ZERO;\r
+               } // for (I = K; I <= M; I++)\r
+               Q[K-1][K-1] = ONE;\r
+               if (WA[K-1] == ZERO) {\r
+                       continue;\r
+               }\r
+               for (J = K; J <= M; J++) {\r
+                   SUM = ZERO;\r
+                   for (I = K; I <= M; I++) {\r
+                       SUM = SUM + Q[I-1][J-1]*WA[I-1];\r
+                   } // for (I = K; I <= M; I++)\r
+                   TEMP = SUM/WA[K-1];\r
+                   for (I = K; I <= M; I++) {\r
+                       Q[I-1][J-1] = Q[I-1][J-1] - TEMP*WA[I-1];\r
+                   } // for (I = K; I <= M; I++)\r
+               } // for (J = K; J <= M; J++)\r
+           } // for (L = 1; L <= MINMN; L++)\r
+             return;\r
+       }\r
+\r
+       \r
+       private void QRFAC(int M,int N, double A[][], int LDA, boolean PIVOT,\r
+                       int IPVT[], int LIPVT, double RDIAG[], double ACNORM[], double WA[]) {\r
+           //     **********\r
+       \r
+           // SUBROUTINE QRFAC\r
+       \r
+           // THIS SUBROUTINE USES HOUSEHOLDER TRANSFORMATIONS WITH COLUMN\r
+           // PIVOTING (OPTIONAL) TO COMPUTE A QR FACTORIZATION OF THE\r
+           // M BY N MATRIX A. THAT IS, QRFAC DETERMINES AN ORTHOGONAL\r
+           // MATRIX Q, A PERMUTATION MATRIX P, AND AN UPPER TRAPEZOIDAL\r
+           // MATRIX R WITH DIAGONAL ELEMENTS OF NONINCREASING MAGNITUDE,\r
+           // SUCH THAT A*P = Q*R. THE HOUSEHOLDER TRANSFORMATION FOR\r
+           // COLUMN K, K = 1,2,...,MIN(M,N), IS OF THE FORM\r
+       \r
+           //                           T\r
+           //           I - (1/U(K))*U*U\r
+       \r
+        // WHERE U HAS ZEROS IN THE FIRST K-1 POSITIONS. THE FORM OF\r
+           // THIS TRANSFORMATION AND THE METHOD OF PIVOTING FIRST\r
+           // APPEARED IN THE CORRESPONDING LINPACK SUBROUTINE.\r
+       \r
+           // THE SUBROUTINE STATEMENT IS\r
+       \r
+           // SUBROUTINE QRFAC(M,N,A,LDA,PIVOT,IPVT,LIPVT,RDIAG,ACNORM,WA)\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 CONTAINS THE MATRIX FOR\r
+           //   WHICH THE QR FACTORIZATION IS TO BE COMPUTED. ON OUTPUT\r
+           // THE STRICT UPPER TRAPEZOIDAL PART OF A CONTAINS THE STRICT\r
+           // UPPER TRAPEZOIDAL PART OF R, AND THE LOWER TRAPEZOIDAL\r
+           // PART OF A CONTAINS A FACTORED FORM OF Q (THE NON-TRIVIAL\r
+           // ELEMENTS OF THE U VECTORS DESCRIBED ABOVE).\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
+           // PIVOT IS A LOGICAL INPUT VARIABLE. IF PIVOT IS SET TRUE,\r
+           //   THEN COLUMN PIVOTING IS ENFORCED. IF PIVOT IS SET FALSE,\r
+           //   THEN NO COLUMN PIVOTING IS DONE.\r
+       \r
+           // IPVT IS AN INTEGER OUTPUT ARRAY OF LENGTH LIPVT. IPVT\r
+           //   DEFINES THE PERMUTATION MATRIX P SUCH THAT A*P = Q*R.\r
+           //   COLUMN J OF P IS COLUMN IPVT(J) OF THE IDENTITY MATRIX.\r
+           //   IF PIVOT IS FALSE, IPVT IS NOT REFERENCED.\r
+       \r
+           // LIPVT IS A POSITIVE INTEGER INPUT VARIABLE. IF PIVOT IS FALSE,\r
+           //   THEN LIPVT MAY BE AS SMALL AS 1. IF PIVOT IS TRUE, THEN\r
+           //   LIPVT MUST BE AT LEAST N.\r
+       \r
+           // RDIAG IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS THE\r
+           //   DIAGONAL ELEMENTS OF R.\r
+       \r
+        // ACNORM IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS THE\r
+           //   NORMS OF THE CORRESPONDING COLUMNS OF THE INPUT MATRIX A.\r
+           //   IF THIS INFORMATION IS NOT NEEDED, THEN ACNORM CAN COINCIDE\r
+           //   WITH RDIAG.\r
+       \r
+           // WA IS A WORK ARRAY OF LENGTH N. IF PIVOT IS FALSE, THEN WA\r
+           //   CAN COINCIDE WITH RDIAG.\r
+       \r
+           // SUBPROGRAMS CALLED\r
+       \r
+           // MINPACK-SUPPLIED ... DPMPAR,ENORM\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,LIPVT,M,N\r
+           // LOGICAL PIVOT\r
+       \r
+           // .. Array Arguments ..\r
+           // DOUBLE PRECISION A(LDA,N),ACNORM(N),RDIAG(N),WA(N)\r
+           // INTEGER IPVT(LIPVT)\r
+       \r
+           //      .. Local Scalars ..\r
+               double ratio;\r
+           double AJNORM,EPSMCH,ONE,P05,SUM,TEMP,ZERO;\r
+           double vec[];\r
+           int I,J,JP1,K,KMAX,MINMN;\r
+       \r
+           //     .. External Functions ..\r
+           //      DOUBLE PRECISION DPMPAR,ENORM\r
+           //      EXTERNAL DPMPAR,ENORM\r
+       \r
+           //     .. Data statements ..\r
+           ONE = 1.0;\r
+           P05 = 5.0E-2;\r
+           ZERO = 0.0;\r
+       \r
+        // EPSMCH IS THE MACHINE PRECISION.\r
+       \r
+           EPSMCH = MACHEP;\r
+       \r
+           // COMPUTE THE INITIAL COLUMN NORMS AND INITIALIZE SEVERAL ARRAYS.\r
+       \r
+           for (J = 1; J <= N; J++) {\r
+               vec = new double[M];\r
+               for (I = 0; I < M; I++) {\r
+                       vec[I] = A[I][J-1];\r
+               }\r
+               ACNORM[J-1] = ENORM(M,vec);\r
+               RDIAG[J-1] = ACNORM[J-1];\r
+               WA[J-1] = RDIAG[J-1];\r
+               if (PIVOT) {\r
+                       IPVT[J-1] = J;\r
+               }\r
+           } // for (J = 1; J <= N; J++)\r
+       \r
+           // REDUCE A TO R WITH HOUSEHOLDER TRANSFORMATIONS.\r
+       \r
+           MINMN = Math.min(M,N);\r
+           for (J = 1; J <= MINMN; J++) {\r
+               if (PIVOT) {\r
+       \r
+                   // BRING THE COLUMN OF LARGEST NORM INTO THE PIVOT POSITION.\r
+       \r
+                   KMAX = J;\r
+                   for (K = J; K <= N; K++) {\r
+                       if (RDIAG[K-1] > RDIAG[KMAX-1]) {\r
+                               KMAX = K;\r
+                       }\r
+                   } // for (K = J; K <= N; K++)\r
+                   if  (KMAX != J) {\r
+                       for (I = 1; I <= M; I++) {\r
+                           TEMP = A[I-1][J-1];\r
+                           A[I-1][J-1] = A[I-1][KMAX-1];\r
+                           A[I-1][KMAX-1] = TEMP;\r
+                       } // for (I = 1; I <= M; I++)\r
+                       RDIAG[KMAX-1] = RDIAG[J-1];\r
+                       WA[KMAX-1] = WA[J-1];\r
+                       K = IPVT[J-1];\r
+                       IPVT[J-1] = IPVT[KMAX-1];\r
+                       IPVT[KMAX-1] = K;\r
+                   } // if (KMAX != J)\r
+               } // if (PIVOT)\r
+       \r
+               // COMPUTE THE HOUSEHOLDER TRANSFORMATION TO REDUCE THE\r
+               // J-TH COLUMN OF A TO A MULTIPLE OF THE J-TH UNIT VECTOR.\r
+       \r
+               vec = new double[M-J+1];\r
+               for (I = 0; I < M-J+1; I++) {\r
+                       vec[I] = A[J-1+I][J-1];\r
+               }\r
+               AJNORM = ENORM(M-J+1,vec);\r
+               if (AJNORM != 0) {\r
+                   if (A[J-1][J-1] < ZERO) {\r
+                       AJNORM = -AJNORM;\r
+                   }\r
+                   for (I = J; I <= M; I++) {\r
+                     A[I-1][J-1] = A[I-1][J-1]/AJNORM;\r
+                   } // for (I = J; I <= M; I++)\r
+                   A[J-1][J-1] = A[J-1][J-1] + ONE;\r
+       \r
+                   // APPLY THE TRANSFORMATION TO THE REMAINING COLUMNS\r
+                   // AND UPDATE THE NORMS.\r
+       \r
+                   JP1 = J + 1;\r
+                   if (N >= JP1) {\r
+                       for (K = JP1; K <= N; K++) {\r
+                           SUM = ZERO;\r
+                           for (I = J; I <= M; I++) {\r
+                               SUM = SUM + A[I-1][J-1]*A[I-1][K-1];\r
+                           } // for (I = J; I <= M; I++)\r
+                           TEMP = SUM/A[J-1][J-1];\r
+                           for (I = J; I <= M; I++) {\r
+                               A[I-1][K-1] = A[I-1][K-1] - TEMP*A[I-1][J-1];\r
+                           } // for (I = J; I <= M; I++)\r
+                           if (PIVOT && RDIAG[K-1] != ZERO) {\r
+                               TEMP = A[J-1][K-1]/RDIAG[K-1];\r
+                               RDIAG[K-1] = RDIAG[K-1]*Math.sqrt(Math.max(ZERO,ONE-TEMP*TEMP));\r
+                               ratio = RDIAG[K-1]/WA[K-1];\r
+                               if (P05*ratio*ratio <= EPSMCH) {\r
+                                       vec = new double[M-J];\r
+                                       for (I = 0; I < M-J; I++) {\r
+                                               vec[I] = A[JP1-1+I][K-1];\r
+                                       }\r
+                                   RDIAG[K-1] = ENORM(M-J,vec);\r
+                                   WA[K-1] = RDIAG[K-1];\r
+                               } // if (P05*ratio*ratio <= EPSMCH)\r
+                           } // if (PIVOT && RDIAG[K-1] != ZERO)\r
+                       } // for (K = JP1; K <= N; K++)\r
+                   } // if (N >= JP1)\r
+               } // if (AJNORM != 0)\r
+               RDIAG[J-1] = -AJNORM;\r
+           } // for (J = 1; J <= MINMN; J++)\r
+             return;\r
+       }\r
+\r
 \r
 \r
 }
\ No newline at end of file