Running JAPHYC produces a failure to converge in EIGSYS error.
authorilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Fri, 9 Feb 2018 22:55:35 +0000 (22:55 +0000)
committerilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Fri, 9 Feb 2018 22:55:35 +0000 (22:55 +0000)
git-svn-id: https://citdcbmipav.cit.nih.gov/repos-pub/mipav/trunk@15371 ba61647d-9d00-f842-95cd-605cb4296b96

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

index 0a6da24..2427239 100644 (file)
@@ -27,12 +27,126 @@ public class SymmsIntegralMapping extends AlgorithmBase {
        // If SYMTY is true, what are the coordinates of the center of symmetry?\r
        private final int MNARC = 100;\r
        private double CENSY[] = new double[2];\r
-       // If SYMTY is true, number of arcs on the fundamental boundary section\r
-       // If SYMTY is false, number of arcs on the boundary\r
        // NARCS <= MNARC-1 = 99\r
+       // NARCS - INTEGER\r
+    // THE NUMBER OF ANALYTIC ARCS THAT MAKE UP THE\r
+       // W H O L E BOUNDARY OF THE PHYSICAL DOMAIN.\r
        private int NARCS;\r
+       // ISYGP - INTEGER\r
+       // THE MAGNITUDE OF ISYGP IS THE ORDER OF THE\r
+       // SYMMETRY GROUP OF THE PHYSICAL DOMAIN.\r
+    // ISYGP.EQ.1 -THE SYMMETRY GROUP HAS ONLY ONE ELE-\r
+       // MENT,THE IDENTITY TRANSFORMATION; IN\r
+       // OTHER WORDS, THE DOMAIN HAS 'NO\r
+       // SYMMETRY'.\r
+       // ISYGP.GT.1 -THE SYMMETRY GROUP CONTAINS ONLY\r
+       // PROPER (IN-PLANE) ROTATIONS; IN OTHER\r
+       // WORDS, THE DOMAIN HAS ONLY ROTATIONAL\r
+       // SYMMETRIES.\r
+       // ISYGP.LT.-1 -THE SYMMETRY GROUP CONTAINS IMPROPER\r
+       // (OUT-OF-PLANE) ROATIONS; IN OTHER\r
+       // WORDS, THE DOMAIN HAS REFLECTIONAL\r
+       // SYMMETRY AND MAY ALSO HAVE ROTATIONAL\r
+       // SYMMETRIES.\r
+       // AN INPUT VALUE OF -1 OR 0 IS TREATED AS IF IT WERE\r
+       // 1.\r
+       private int ISYGP;\r
+       // RFARC - INTEGER\r
+    // THE REFERENCE ARC USED TO DEFINE THE ORIENTATION\r
+    // THAT IS GIVEN TO THE MAP. THE CONVENTION IS THAT\r
+    // THE POINT AT THE START OF ANALYTIC ARC NUMBER\r
+    // RFARC IS MAPPED TO THE POINT WITH ARGUMENT\r
+    // RFARG*PI ON THE UNIT DISC.\r
+    private int RFARC;\r
+    // RFARG - DOUBLE\r
+    // THE REFERENCE ARGUMENT/PI USED TO DEFINE THE\r
+    // ORIENTATION THAT IS GIVEN TO THE MAP. SEE RFARC\r
+    // ABOVE.\r
+    private double RFARG[] = new double[1];\r
+    // INCST - BOOLEAN\r
+    // IF INCST IS TRUE THEN AN INCREMENTAL STRATEGY IS\r
+    // USED TO TRY TO ACHIEVE THE ACCURACY SPECIFIED BY\r
+    // MAXER; VERY ROUGHLY SPEAKING, THIS MEANS THAT THE\r
+    // METHOD SUCCESSIVELY ACHIEVES THE TARGET ACCURACIES\r
+    // 1E-1,1E-2,...UNTIL MAXER HAS BEEN ACHIEVED. IF\r
+    // THE PROBLEM IS THOUGHT TO BE EITHER PARTICULARLY\r
+    // DIFFICULT OR PARTICULARLY SIMPLE, THEN INCST\r
+    // SHOULD BE SET TO .TRUE. FOR PROBLEMS OF 'AVERAGE'\r
+    // DIFFICULTY, SETTING INCST TO .FALSE. IS USUALLY\r
+    // MORE EFFICIENT.\r
+    private boolean INCST;\r
+       // TSTNG - INTEGER\r
+       // EITHER 0 OR 1.\r
+    // ON SUCCESSFUL COMPLETION OF THE NUMERICAL SOLUTION\r
+    // OF SYMM'S EQUATION, A MODULE IS PROVIDED FOR\r
+       // TESTING THE ERROR IN THE MODULUS OF THE COMPUTED\r
+       // MAP ON THE BOUNDARY OF THE DOMAIN.\r
+       // TSTNG=0 - TEST ONLY AT SUB-ARC END POINTS\r
+    // TSTNG=1 - IN ADDITION TO TESTING AT SUB-ARC END\r
+    // POINTS TEST ALSO AT INTERIOR POINTS\r
+       // ON EACH SUB-ARC.\r
+       private int TSTNG;\r
+       // MAXER - DOUBLE\r
+    // RELATIVE ACCURACY REQUESTED FOR THE CONFORMAL MAP;\r
+    // THIS IS THE SAME AS THE ABSOLUTE ACCURACY ON THE\r
+    // BOUNDARY OF THE PHYSICAL DOMAIN.\r
+       private double MAXER;\r
+       // OULVL - INTEGER\r
+    // EITHER 0,1,2,3,4 OR 5.\r
+       // CONTROLS THE AMOUNT OF OUTPUT IN THE LISTING FILE\r
+       // <JBNM>pl.\r
+       // OULVL=0 - OUTPUT A SOLUTION SUMMARY AT EACH STAGE\r
+       // IN THE ADAPTIVE PROCESS AND A SHORT\r
+       // SUMMARY OF THE ERRORS IN MODULUS.\r
+       // OULVL=1 - AS 0, BUT ALSO OUTPUT A DETAILED LIST OF\r
+       // THE ERRORS IN MODULUS.\r
+       // OULVL=2 - AS 0, BUT ALSO OUTPUT FULL DETAILS OF\r
+       // THE FINAL COMPUTED JACOBI COEFFICIENTS\r
+    // ON SUCCESSFUL COMPLETION.\r
+       // OULVL=3 - AS 2, BUT ALSO OUTPUT A DETAILED LIST OF\r
+       // THE ERRORS IN MODULUS.\r
+    // OULVL=4 - OUTPUT FULL DETAILS OF THE OF THE COMPU-\r
+       // TED JACOBI COEFFICIENTS AT EVERY STAGE\r
+       // IN THE ADAPTIVE PROCESS AND A SHORT\r
+       // SUMMARY OF THE ERRORS IN MODULUS.\r
+       // OULVL=5 - AS 4, BUT ALSO OUTPUT A DETAILED LIST OF\r
+       // THE ERRORS IN MODULUS.\r
+       private int OULVL;\r
        // NUMDER is of length MNARC\r
        // NUMDER is initially false\r
+       // MNSUA= MAXIMUM NUMBER OF SUBARCS ALLOWED ON THE PHYSICAL BOUNDARY\r
+    private int MNSUA = 150; // 150 in DRIVE1\r
+    // NUMBER OF QUADRATURE POINTS\r
+    private int NQPTS = 8; // 8 in DRIVE1\r
+    // MNEQN= MAXIMUM NUMBER OF LINEAR EQUATIONS ALLOWED IN THE \r
+    //        COLLOCATION SOLUTION OF SYMM'S EQUATION\r
+    private int MNEQN = 150; // 150 in DRIVE1\r
+    // MNJXS= (MAXIMUM NUMBER OF CORNERS ALLOWED) + 1\r
+    private int MNJXS = 50; // 50 in DRIVE1\r
+    // MQIN1= (MAXIMUM NUMBER OF PANELS ALLOWED IN A SINGLE\r
+    //        COMPOSITE GAUSSIAN RULE) + 1\r
+    // From GQCANP:\r
+       // MQIN1 - INTEGER\r
+       // DEFINES THE NUMBER OF PANELS ALLOWED IN A\r
+       // COMPOSITE RULE. SPECIFICALLY, MQIN1 = 1 + (THE\r
+       // MAXIMUM NUMBER OF PANELS IN A COMPOSITE RULE FOR\r
+       // A SINGLE SUB-ARC ON THE BOUNDARY)\r
+    private int MQIN1 = 11; // 11 in DRIVE0\r
+    // MNQUA= MAXIMUM TOTAL NUMBER OF QUADRATURE POINTS ALLOWED OVER\r
+    //        ALL COMPOSITE GAUSSIAN RULES AT THE PROCESSING STAGE OF\r
+    // SOLVING SYMM'S EQUATION\r
+    private int MNQUA = 2000; // 2000 in DRIVE0\r
+    // SET UP THE ARRAY OF BOUNDS IBNDS FOR JAPHYC\r
+    private int IBNDS[] = new int[]{MNSUA,MNJXS,MQIN1,MNQUA,0};\r
+    // MATRX - DOUBLE ARRAY\r
+    // A 3-DIMENSIONAL MATRIX OF SIZE\r
+    // MNEQN X MNEQN X 2 .\r
+       // (IN THE ADAPTIVE PROCESS, MATRX(*,*,2) WILL STORE\r
+    // THE COEFFICIENT MATRIX OF THE CURRENT COLLOCATION\r
+    // SYSTEM AND MATRX(*,*,1) WILL STORE THE COEFFICIENT\r
+       // MATRIX OF THE PREVIOUS SYSTEM)\r
+    private double MATRX[][][] = new double[MNEQN][MNEQN][2];\r
+    private int IER[] = new int[1];\r
        private boolean NUMDER[] = new boolean[MNARC];\r
        // ARCTY is of length MNARC\r
        // Type of arc, 1 = LINE SEGMENT, 2 = CIRCULAR ARC SEGMENT, 3 = CARTESIAN\r
@@ -95,7 +209,6 @@ public class SymmsIntegralMapping extends AlgorithmBase {
        private boolean traditionalInput = true;\r
        Scanner input = new Scanner(System.in);\r
        private double zzset[][] = new double[400][2];\r
-       private int IBNDS[] = new int[5];\r
        private int IGEOM[] = new int[4]; // Written in original JAPHYC\r
        private int PARNT[] = new int[IBNDS[0]]; // Written in original JAPHYC\r
        private double RGEOM[] = new double[2]; // Written in original JAPHYC\r
@@ -107,10 +220,8 @@ public class SymmsIntegralMapping extends AlgorithmBase {
        private int JATYP[] = new int[IBNDS[0]]; // Written in original JAPHYC\r
        private int LOSUB[] = new int[IBNDS[0]]; // Written in original JAPHYC\r
        // Parts of RSNPH\r
-       private int NQPTS;\r
        private int NJIND = NARCS + 1;\r
        private int TNGQP = NQPTS * NJIND;\r
-       private int MNEQN;\r
        private double ACOEF[] = new double[TNGQP]; // Written in original JAPHYC\r
        private double BCOEF[] = new double[TNGQP]; // Written in original JAPHYC\r
        private double AICOF[] = new double[TNGQP]; // Written in original JAPHYC\r
@@ -258,13 +369,6 @@ public class SymmsIntegralMapping extends AlgorithmBase {
        private double ZSNCA[] = new double[2]; // At first location of ZSNCA\r
        private double BFSNC[][] = new double[IBNDS[1]][2];\r
        private double SOLNC[][] = new double[IBNDS[1]][2];\r
-       // From GQCANP:\r
-       // MQIN1 - INTEGER\r
-       // DEFINES THE NUMBER OF PANELS ALLOWED IN A\r
-       // COMPOSITE RULE. SPECIFICALLY, MQIN1 = 1 + (THE\r
-       // MAXIMUM NUMBER OF PANELS IN A COMPOSITE RULE FOR\r
-       // A SINGLE SUB-ARC ON THE BOUNDARY)\r
-       private int MQIN1;\r
        // MQUCA - INTEGER\r
        // THE MAXIMUM NUMBER OF QUADRATURE POINTS ALLOWED\r
        // IN THE FINAL GLOBAL RULE. THE VALUE OF THIS\r
@@ -299,7 +403,87 @@ public class SymmsIntegralMapping extends AlgorithmBase {
        private int example = 1;\r
 \r
        public SymmsIntegralMapping() {\r
-\r
+               if (example == 1) {\r
+                       INTER = true;\r
+               NARCS = 5;\r
+               ISYGP = 1;\r
+               RFARC = 2;\r
+               RFARG[0] = 0.5;\r
+               INCST = false;\r
+               TSTNG = 1;\r
+               CENTR[0] = 0.0;\r
+               CENTR[1] = 0.0;\r
+               MAXER = 1.0E-4;\r
+               OULVL = 2;\r
+        }\r
+        else if (example == 2) {\r
+               INTER = true;\r
+               NARCS = 4;\r
+               ISYGP = 1;\r
+               RFARC = 1;\r
+               RFARG[0] = 0.0;\r
+               INCST = true;\r
+               TSTNG = 1;\r
+               CENTR[0] = 0.6;\r
+               CENTR[1] = 0.0;\r
+               MAXER = 1.0E-4;\r
+               OULVL = 2;\r
+        }\r
+        else if (example == 3) {\r
+               INTER = false;\r
+               NARCS = 8;\r
+               ISYGP = 4;\r
+               RFARC = 1;\r
+               RFARG[0] = 0.0;\r
+               INCST = false;\r
+               TSTNG = 1;\r
+               CENTR[0] = 0.0;\r
+               CENTR[1] = 0.0;\r
+               MAXER = 1.0E-4;\r
+               OULVL = 2;\r
+        }\r
+        else if (example == 4) {\r
+               INTER = true;\r
+               NARCS = 48;\r
+               ISYGP = -12;\r
+               RFARC = 1;\r
+               RFARG[0] = 0.0;\r
+               INCST = true;\r
+               TSTNG = 1;\r
+               CENTR[0] = 0.0;\r
+               CENTR[1] = -5.196152422706632;\r
+               MAXER = 1.0E-4;\r
+               OULVL = 2;\r
+        }\r
+               NJIND = NARCS + 1;\r
+               TNGQP = NQPTS * NJIND;\r
+               // Parts of RSNPH\r
+               ACOEF = new double[TNGQP]; // Written in original JAPHYC\r
+               BCOEF = new double[TNGQP]; // Written in original JAPHYC\r
+               AICOF = new double[TNGQP]; // Written in original JAPHYC\r
+               BICOF = new double[TNGQP]; // Written in original JAPHYC\r
+               QUPTS = new double[TNGQP]; // Written in original JAPHYC\r
+               QUWTS = new double[TNGQP]; // Written in original JAPHYC\r
+               H0VAL = new double[NJIND]; // Written in original JAPHYC\r
+               HIVAL = new double[NJIND]; // Written in original JAPHYC\r
+               JACIN = new double[NJIND]; // Written in original JAPHYC\r
+               // Part of RWORK\r
+               AQCOF = new double[TNGQP];\r
+               BQCOF = new double[TNGQP];\r
+               CQCOF = new double[TNGQP];\r
+               COLSC = new double[TNGQP];\r
+               RIGLL = new double[TNGQP];\r
+               // Part of RSNCA\r
+               ACOFC = new double[TNGQP];\r
+               BCOFC = new double[TNGQP];\r
+               AICOC = new double[TNGQP];\r
+               BICOC = new double[TNGQP];\r
+               QUPTC = new double[TNGQP];\r
+               QUWTC = new double[TNGQP];\r
+               H0VLC = new double[NJIND];\r
+               HIVLC = new double[NJIND];\r
+               JAINC = new double[NJIND];\r
+               COARG = new double[NJIND];\r
        }\r
 \r
        public SymmsIntegralMapping(ModelImage destImg, ModelImage srcImg, String FORTFL, boolean SYMTY, boolean REFLN,\r
@@ -867,26 +1051,77 @@ public class SymmsIntegralMapping extends AlgorithmBase {
                // EXAMPLE 2: 4.78E-16\r
                // EXAMPLE 3: 2.449E-16;\r
                // EXAMPLE 4: 9.86E-16\r
-        int IER[] = new int[1];\r
         double MXMIS[] = new double[1];\r
         double MXDIF[] = new double[1];\r
         double PSD[] = new double[1];\r
         double MINPD = 0.0;\r
         double MAXPD = 0.0;\r
-        if (example == 1) {\r
-               NARCS = 5;\r
-        }\r
-        else if (example == 2) {\r
-               NARCS = 4;\r
-        }\r
-        else if (example == 3) {\r
-               NARCS = 8;\r
-        }\r
-        else if (example == 4) {\r
-               NARCS = 48;\r
-        }\r
-        TSTPLT(MXMIS,MXDIF,NARCS,PSD,MINPD,MAXPD,IER);\r
-       }\r
+        TSTPLT(MXMIS,MXDIF,PSD,MINPD,MAXPD);\r
+       } // public void DRIVE0()\r
+       \r
+\r
+    public void DRIVE1() {\r
+        \r
+        // .......................................................................\r
+        // EXAMPLE PROGRAM TO SHOW THE USE OF THE SYMM EQUATION SOLVING\r
+        // SUBROUTINE JAPHYC FROM THE\r
+    \r
+        //           C O N F P A C K    L I B R A R Y .\r
+        // .......................................................................\r
+\r
+        // SCALAR VARIABLES REQUIRED FOR MAIN SUBROUTINE\r
+\r
+        int IER[] = new int[1];\r
+\r
+        // MNEQN= MAXIMUM NUMBER OF LINEAR EQUATIONS ALLOWED IN THE \r
+        //        COLLOCATION SOLUTION OF SYMM'S EQUATION\r
+\r
+        // MNJXS= (MAXIMUM NUMBER OF CORNERS ALLOWED) + 1\r
+\r
+        // MQIN1= (MAXIMUM NUMBER OF PANELS ALLOWED IN A SINGLE\r
+        //        COMPOSITE GAUSSIAN RULE) + 1\r
+\r
+        // MNQUA= MAXIMUM TOTAL NUMBER OF QUADRATURE POINTS ALLOWED OVER\r
+        //        ALL COMPOSITE GAUSSIAN RULES AT THE PROCESSING STAGE OF\r
+        // SOLVING SYMM'S EQUATION\r
+\r
+        // MNSUA= MAXIMUM NUMBER OF SUBARCS ALLOWED ON THE PHYSICAL \r
+        // BOUNDARY\r
+\r
+        // THE ABOVE BOUNDS ARE RELATED TO THE NUMBER OF QUADRATURE POINTS\r
+        // NQPTS VIA\r
+        //       MNEQN <= MNSUA*NQPTS + 1\r
+        //       MNQUA <= (MQIN1-1)*MNJXS*NQPTS \r
+    \r
+        // BUT IN PRACTICE MNEQN AND MNQUA MAY BE SIGNIFICANTLY\r
+        // LESS THAN THE VALUES ON THE RIGHT OF THESE INEQUALITIES.\r
+\r
+        // SUGGESTED VALUES TO COVER A WIDE RANGE OF PROBLEMS ARE SET IN\r
+        // THE FOLLOWING PARAMETER STATEMENT.\r
+\r
+        //NQPTS=8;\r
+        //MNSUA=150;\r
+        //MNEQN=500;\r
+        //MNJXS=50;\r
+        //MQIN1=11;\r
+        //MNQUA=2000;\r
+\r
+        // INTEGER \r
+        // IBNDS(4)\r
+        // SET UP THE ARRAY OF BOUNDS IBNDS FOR JAPHYC\r
+        //IBNDS[0]=MNSUA;\r
+        //IBNDS[1]=MNJXS;\r
+        //IBNDS[2]=MQIN1;\r
+        //IBNDS[3]=MNQUA;\r
+   \r
+        // SOLVE SYMM'S INTEGRAL EQUATION, TO ESTIMATE THE JACOBI \r
+        // COEFFICIENTS FOR THE DENSITY ASSOCIATED WITH THE MAP : PHYSICAL\r
+        // CANONICAL, WITH AUTOMATIC CREATION OF OUTPUT FILES \r
+        // jbnm, <JBNM>pl, <JBNM>gm, <JBNM>ph.\r
+        JAPHYC();\r
+\r
+    } // public void DRIVE1()\r
+\r
 \r
 \r
        private void HEADER(String TXT, String REDD, RandomAccessFile raFile) {\r
@@ -2630,8 +2865,7 @@ public class SymmsIntegralMapping extends AlgorithmBase {
 \r
        } // private void WRSYM3\r
 \r
-       private void TSTPLT(double MXMIS[], double MXDIF[], int NARCS, double PSD[], double MINPD, double MAXPD,\r
-                       int IER[]) {\r
+       private void TSTPLT(double MXMIS[], double MXDIF[], double PSD[], double MINPD, double MAXPD) {\r
                // CHARACTER*4 JBNM\r
 \r
                // ......................................................................\r
@@ -3127,8 +3361,7 @@ public class SymmsIntegralMapping extends AlgorithmBase {
                return result;\r
        }\r
 \r
-       private void JAPHYC(double MAXER, int ISYGP, boolean INCST, int RFARC, double RFARG[], int TSTNG, int OULVL,\r
-                       double MATRX[][][], int IER[]) {\r
+       private void JAPHYC() {\r
 \r
                // INTEGER IBNDS(*),IGEOM(*),ISNPH(*),IWORK(*)\r
                // REAL RGEOM(*),MATRX(MNEQN,MNEQN,2),RSNPH(*),RWORK(*)\r
@@ -4212,7 +4445,7 @@ public class SymmsIntegralMapping extends AlgorithmBase {
                Preferences.debug("\n", Preferences.DEBUG_ALGORITHM);\r
                Preferences.debug("MINIMUM NUMBER OF QUADRATURE POINTS : " + NQPTS + "\n", Preferences.DEBUG_ALGORITHM);\r
                Preferences.debug("MAXIMUM DEGREE OF POLYNOMIAL        : " + (NQPTS - 1) + "\n", Preferences.DEBUG_ALGORITHM);\r
-               Preferences.debug("INITIAL DEGREE OF POLYNOMIAL        : " + INDEG + (NQPTS - 1) + "\n",\r
+               Preferences.debug("INITIAL DEGREE OF POLYNOMIAL        : " + INDEG + "\n",\r
                                Preferences.DEBUG_ALGORITHM);\r
                Preferences.debug("INCREMENTAL STRATEGY                : " + INCST + "\n", Preferences.DEBUG_ALGORITHM);\r
 \r