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

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

index fd90e54..de0966f 100644 (file)
@@ -88,6 +88,76 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
     private boolean traditionalInput = false;\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];\r
+    private int PARNT[] = new int[IBNDS[0]];\r
+    private double RGEOM[] = new double[2];\r
+    private double HALEN[]= new double[IBNDS[0]];\r
+       private double MIDPT[] = new double[IBNDS[0]];\r
+       private double VTARG[] = new double[IBNDS[0]]; \r
+    private int ISNPH[] = new int[6];\r
+    private int DGPOL[] = new int[IBNDS[0]];\r
+       private int JATYP[] = new int[IBNDS[0]];\r
+       private int LOSUB[] = new int[IBNDS[0]];\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];\r
+       private double BCOEF[] = new double[TNGQP];\r
+       private double AICOF[] = new double[TNGQP];\r
+       private double BICOF[] = new double[TNGQP];\r
+       private double QUPTS[] = new double[TNGQP];\r
+       private double QUWTS[] = new double[TNGQP];\r
+       private double H0VAL[] = new double[NJIND];\r
+       private double HIVAL[] = new double[NJIND];\r
+       private double JACIN[] = new double[NJIND];\r
+       private double ERARC[] = new double[IBNDS[0]];\r
+       private double BCFSN[] = new double[MNEQN];\r
+       private double SOLUN[] = new double[MNEQN];\r
+       // Parts of IWORK\r
+       private int IPIVT[] = new int[MNEQN];\r
+       private int LOQSB[] = new int[IBNDS[1]];\r
+       private int NQUAD[] = new int[IBNDS[1]];\r
+       private int HISUB[] = new int[IBNDS[0]];\r
+       private int LOTES[] = new int[IBNDS[0]];\r
+       private int HITES[] = new int[IBNDS[0]];\r
+       private int AXION[] = new int[IBNDS[0]];\r
+       private int NEWDG[] = new int[IBNDS[0]];\r
+       private int ICOPY[] = new int[IBNDS[0]];\r
+       private int LOOLD[] = new int[IBNDS[0]];\r
+       private int HIOLD[] = new int[IBNDS[0]];\r
+       // Parts of RWORK\r
+       private double WORK2[] = new double[MNEQN];\r
+       private double COLPR[] = new double[MNEQN];\r
+       private double A1COF[] = new double[IBNDS[1]];\r
+       private double B1COF[] = new double[IBNDS[1]];\r
+       private double TOLOU[] = new double[IBNDS[1]];\r
+       private double XIDST[] = new double[2*IBNDS[1]];\r
+       private double XENPT[] = new double[IBNDS[2]];\r
+       private double QCOMX[] = new double[IBNDS[3]];\r
+       private double QCOMW[] = new double[IBNDS[3]];\r
+       private double RCOPY[] = new double[IBNDS[0]];\r
+       private double NEWHL[] = new double[IBNDS[0]];\r
+       private double AQCOF[] = new double[TNGQP];\r
+       private double BQCOF[] = new double[TNGQP];\r
+       private double CQCOF[] = new double[TNGQP];\r
+       private double COLSC[] = new double[TNGQP];\r
+       private double RIGLL[] = new double[TNGQP];\r
+       private double WORK[] = new double[2*NQPTS];\r
+       private double DIAG[] = new double[NQPTS];\r
+       private double SDIAG[] = new double[NQPTS];\r
+       private double WORKT[] = new double[2*NQPTS*NQPTS];\r
+       private double WORKQ[] = new double[NQPTS*NQPTS];\r
+       // Parts of ZWORK\r
+       private double ZCOLL[][] = new double[MNEQN][2];\r
+       private double XIVAL[][] = new double[2*IBNDS[1]][2];\r
+       // Parts of LWORK\r
+       private boolean NEWQU[] = new boolean[IBNDS[1]];\r
+       private boolean LCOPY[] = new boolean[IBNDS[0]];\r
+       private boolean PNEWQ[] = new boolean[IBNDS[0]];\r
+       private boolean LNSEG[] = new boolean[IBNDS[0]];\r
     \r
        public SymmsIntegralMapping() {\r
                \r
@@ -2254,12 +2324,12 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        return result;\r
    }\r
 \r
-   private void JAPHYC(String JBNM, String HEAD, double MAXER,boolean INTER, int NARCS,\r
-                     int ISYGP, int NQPTS, boolean INCST,\r
-                     int RFARC, double RFARG[], double CENTR[], int TSTNG,int OULVL, int IBNDS[],\r
-                     int MNEQN, double MATRX[][][], int IWORK[], double RWORK[],\r
-                     double ZWORK[][], boolean LWORK[], int OCH, int IGEOM[], double RGEOM[],\r
-                     int ISNPH[], double RSNPH[], int IER[]) {\r
+   private void JAPHYC(String JBNM, String HEAD, double MAXER,boolean INTER,\r
+                     int ISYGP, boolean INCST,\r
+                     int RFARC, double RFARG[], double CENTR[], int TSTNG,int OULVL,\r
+                     double MATRX[][][],\r
+                     int OCH,\r
+                     int IER[]) {\r
                \r
                      //INTEGER IBNDS(*),IGEOM(*),ISNPH(*),IWORK(*)\r
                      //REAL RGEOM(*),MATRX(MNEQN,MNEQN,2),RSNPH(*),RWORK(*)\r
@@ -2567,23 +2637,9 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                     \r
                //     LOCAL VARAIBLES\r
                \r
-               //**** POINTERS FOR IGEOM,RGEOM,ISNPH,RSNPH\r
-               \r
-               int ACOEF,AICOF,BCFSN,BCOEF,BICOF,DGPOL,ERARC,H0VAL,HALEN,\r
-                    HIVAL,JACIN,JATYP,LOSUB,MIDPT,PARNT,QUPTS,QUWTS,SOLUN,VTARG;\r
-               \r
-               //**** POINTERS FOR IWORK,RWORK,ZWORK,LWORK\r
-               \r
-               int A1COF,AQCOF,AXION,B1COF,BQCOF,COLPR,COLSC,CQCOF,DIAG,\r
-                    HIOLD,HISUB,HITES,ICOPY,IPIVT,LCOPY,LNSEG,LOOLD,LOQSB,LOTES,NEWDG,\r
-                    NEWHL,NEWQU,NQUAD,PNEWQ,QCOMW,QCOMX,RCOPY,RIGLL,SDIAG,TOLOU,WORK2,\r
-                    WORKQ,WORK,WORKT,XENPT,XIDST,XIVAL,ZCOLL;\r
-               \r
-               //**** OTHER LOCAL VARIABLES\r
-               \r
                int I,IMXER,INDEG,J,L,MDGPO,MNJXS,MNQUA,MNSUA,MQIN1,NCOLL,\r
-                    NEFF,NEQNS,NJIND,NROWS,NTEST,\r
-                    TNSUA,TNGQP,ORDSG;\r
+                    NEFF,NEQNS,NROWS,NTEST,\r
+                    TNSUA,ORDSG;\r
                int SOLCO = 0;\r
                int QIERC[] = new int[7];\r
                int QIERR[] = new int[7];\r
@@ -2698,83 +2754,10 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                RGEOM[0]=GSUPE;\r
                RGEOM[1]=GLGTL;\r
                \r
-               //**** SET UP THE POINTERS TO ELEMENTS IN ARRAYS IGEOM AND RGEOM\r
-               \r
-               PARNT=5;\r
-               HALEN=3;\r
-               MIDPT=MNSUA+3;\r
-               VTARG=2*MNSUA+3; \r
-               \r
-               //**** SET UP THE POINTERS TO ELEMENTS IN ARRAYS ISNPH AND RSNPH\r
-               \r
-               DGPOL=7;\r
-               JATYP=MNSUA+7;\r
-               LOSUB=2*MNSUA+7;\r
-               ACOEF=1;\r
-               BCOEF=TNGQP+1;\r
-               AICOF=2*TNGQP+1;\r
-               BICOF=3*TNGQP+1;\r
-               QUPTS=4*TNGQP+1;\r
-               QUWTS=5*TNGQP+1;\r
-               H0VAL=6*TNGQP+1;\r
-               HIVAL=NJIND+6*TNGQP+1;\r
-               JACIN=2*NJIND+6*TNGQP+1;\r
-               ERARC=3*NJIND+6*TNGQP+1;\r
-               BCFSN=MNSUA+3*NJIND+6*TNGQP+1;\r
-               SOLUN=MNEQN+MNSUA+3*NJIND+6*TNGQP+1;\r
-               \r
-               //**** SET UP THE POINTERS TO ELEMENTS IN ARRAYS IWORK,RWORK,ZWORK AND \r
-               //**** LWORK\r
-               \r
-               IPIVT=1;\r
-               LOQSB=MNEQN+1;\r
-               NQUAD=MNJXS+MNEQN+1;\r
-               HISUB=2*MNJXS+MNEQN+1;\r
-               LOTES=MNSUA+2*MNJXS+MNEQN+1;\r
-               HITES=2*MNSUA+2*MNJXS+MNEQN+1;\r
-               AXION=3*MNSUA+2*MNJXS+MNEQN+1;\r
-               NEWDG=4*MNSUA+2*MNJXS+MNEQN+1;\r
-               ICOPY=5*MNSUA+2*MNJXS+MNEQN+1;\r
-               LOOLD=6*MNSUA+2*MNJXS+MNEQN+1;\r
-               HIOLD=7*MNSUA+2*MNJXS+MNEQN+1;\r
-               WORK2=1;\r
-               COLPR=MNEQN+1;\r
-               A1COF=2*MNEQN+1;\r
-               B1COF=MNJXS+2*MNEQN+1;\r
-               TOLOU=2*MNJXS+2*MNEQN+1;\r
-               XIDST=3*MNJXS+2*MNEQN+1;\r
-               XENPT=5*MNJXS+2*MNEQN+1;\r
-               QCOMX=MQIN1+5*MNJXS+2*MNEQN+1;\r
-               QCOMW=MNQUA+MQIN1+5*MNJXS+2*MNEQN+1;\r
-               RCOPY=2*MNQUA+MQIN1+5*MNJXS+2*MNEQN+1;\r
-               NEWHL=MNSUA+2*MNQUA+MQIN1+5*MNJXS+2*MNEQN+1;\r
-               AQCOF=2*MNSUA+2*MNQUA+MQIN1+5*MNJXS+2*MNEQN+1;\r
-               BQCOF=TNGQP+2*MNSUA+2*MNQUA+MQIN1+5*MNJXS+2*MNEQN+1;\r
-               CQCOF=2*TNGQP+2*MNSUA+2*MNQUA+MQIN1+5*MNJXS+2*MNEQN+1;\r
-               COLSC=3*TNGQP+2*MNSUA+2*MNQUA+MQIN1+5*MNJXS+2*MNEQN+1;\r
-               RIGLL=4*TNGQP+2*MNSUA+2*MNQUA+MQIN1+5*MNJXS+2*MNEQN+1;\r
-               WORK=5*TNGQP+2*MNSUA+2*MNQUA+MQIN1+5*MNJXS+2*MNEQN+1;\r
-               DIAG=2*NQPTS+5*TNGQP+2*MNSUA+2*MNQUA+MQIN1+5*MNJXS+2*MNEQN+1;\r
-               SDIAG=3*NQPTS+5*TNGQP+2*MNSUA+2*MNQUA+MQIN1+5*MNJXS+2*MNEQN+1;\r
-               WORKT=4*NQPTS+5*TNGQP+2*MNSUA+2*MNQUA+MQIN1+5*MNJXS+2*MNEQN+1;\r
-               WORKQ=2*NQPTS*NQPTS+4*NQPTS+5*TNGQP+2*MNSUA+2*MNQUA+MQIN1+5*MNJXS+\r
-                     2*MNEQN+1;\r
-               ZCOLL=1;\r
-               XIVAL=MNEQN+1;\r
-               NEWQU=1;\r
-               LCOPY=MNJXS+1;\r
-               PNEWQ=MNSUA+MNJXS+1;\r
-               LNSEG=2*MNSUA+MNJXS+1;\r
-               \r
                //**** ASSIGN THE JACOBI INDECES FOR EACH ARC.\r
                \r
-               //ANGLE7(RSNPH(JACIN),NARCS,INTER);\r
-               double BE[] = new double[NARCS];\r
-               ANGLE7(BE, NARCS, INTER);\r
-               for (I = 0; I < NARCS; I++) {\r
-                       RSNPH[JACIN + I - 1] = BE[I];\r
-               }\r
-               RSNPH[JACIN+NJIND-2]=0.0;\r
+               ANGLE7(JACIN,NARCS,INTER);\r
+               JACIN[NJIND-1]=0.0;\r
                \r
                //**** SET SUB-TOLERANCES AND INDEG\r
                \r
@@ -2805,23 +2788,12 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                \r
                //**** ASSIGN THE LOGICAL LINE SEGMENT TYPE TO EACH ARC.\r
                \r
-               boolean linout[] = new boolean[NARCS];\r
-               LINSEG(linout, NARCS);\r
-               for (I = 0; I < NARCS; I++) {\r
-                       LWORK[LNSEG-1+I]= linout[I];\r
-               }\r
-               //LINSEG(LWORK(LNSEG),NARCS);\r
+               LINSEG(LNSEG,NARCS);\r
                \r
                //**** LIST THE INPUT ARGUMENTS AND ASSOCIATED QUANTITIES\r
                \r
-               double beta[] = new double[NARCS];\r
-               boolean linear[] = new boolean[NARCS];\r
-               for (I = 0; I < NARCS; I++) {\r
-                       beta[I] = RSNPH[JACIN-1+I];\r
-                       linear[I] = LWORK[LNSEG-1+I];\r
-               }\r
                RSLT80(JBNM,HEAD,GSUPE,MAXER,GAQTL,INTER,NARCS,ORDSG,NQPTS,\r
-                    INCST,INDEG,RFARC,RFARG[0],CENTR,beta,linear,\r
+                    INCST,INDEG,RFARC,RFARG[0],CENTR,JACIN,LNSEG,\r
                     TSTNG,OULVL,IBNDS,MNEQN,OCH);\r
                PI=Math.PI;\r
                RFARG[0]=RFARG[0]*PI;\r
@@ -2834,43 +2806,9 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                //**** FOR THE INTEGRATED POLYNOMIALS, STORING RESULTS IN AICOF,BICOF\r
                //**** AND HIVAL.\r
                      \r
-               double JAC[] = new double[NJIND];\r
-               for (I = 0; I < NJIND; I++) {\r
-                       JAC[I] = RSNPH[JACIN+I-1];\r
-               }\r
-               double ACO[] = new double[NJIND*NQPTS];\r
-               double BCO[] = new double[NJIND*NQPTS];\r
-               double H0V[] = new double[NJIND];\r
-               double AIC[] = new double[NJIND*NQPTS+1];\r
-               double BIC[] = new double[NJIND*NQPTS+1];\r
-               double HIV[] = new double[NJIND];\r
-               double QUP[] = new double[NJIND*NQPTS];\r
-               double QUW[] = new double[NJIND*NQPTS];\r
-               double WOR[] = new double[NQPTS];\r
-               OPQUD1(NJIND,NQPTS,JAC,ACO,BCO,\r
-                    H0V,AIC,BIC,HIV,QUP,\r
-                    QUW,WOR,IER);\r
-               for (I = 0; I < NJIND*NQPTS; I++) {\r
-                       RSNPH[ACOEF+I-1] = ACO[I];\r
-                       RSNPH[BCOEF+I-1] = BCO[I];\r
-               }\r
-               for (I = 0; I < NJIND; I++) {\r
-                       RSNPH[H0VAL+I-1] = H0V[I];\r
-               }\r
-               for (I = 0; I < NJIND*NQPTS+1; I++) {\r
-                       RSNPH[AICOF+I-1] = AIC[I];\r
-                       RSNPH[BICOF+I-1] = BIC[I];\r
-               }\r
-               for (I = 0; I < NJIND; I++) {\r
-                       RSNPH[HIVAL+I-1] = HIV[I];\r
-               }\r
-               for (I = 0; I < NJIND*NQPTS; I++) {\r
-                       RSNPH[QUPTS+I-1] = QUP[I];\r
-                       RSNPH[QUWTS+I-1] = QUW[I];\r
-               }\r
-               for (I = 0; I < NQPTS; I++) {\r
-                       RWORK[WORK+I-1] = WOR[I];\r
-               }\r
+               OPQUD1(NJIND,NQPTS,JACIN,ACOEF,BCOEF,\r
+                    H0VAL,AICOF,BICOF,HIVAL,QUPTS,\r
+                    QUWTS,WORK,IER);\r
                if (IER[0] > 0) {\r
                        WRTAIL(1,0,IER[0],null);\r
                    return;            \r
@@ -2878,8 +2816,8 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                J=1-NQPTS;\r
                for (I=1; I <= NJIND; I++) {\r
                    J=J+NQPTS;\r
-                   RWORK[A1COF+I-2]=RSNPH[ACOEF+J-2];\r
-                   RWORK[B1COF+I-2]=RSNPH[BCOEF+J-2];\r
+                   A1COF[I-1]=ACOEF[J-1];\r
+                   B1COF[I-1]=BCOEF[J-1];\r
                } // for (I=1; I <= NJIND; I++)\r
                System.out.println("BASIC GAUSS QUADRATURE DATA DONE:");\r
                \r
@@ -2888,21 +2826,19 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                //**** JACOBI WEIGHT FUNCTIONS AND THEIR ORTHONORMAL POLYNOMIALS; STORE\r
                //**** THESE COEFFICIENTS IN AQCOF, BQCOF AND CQCOF\r
                \r
-/*                   CALL ASQUC7(RWORK(AQCOF),RWORK(BQCOF),RWORK(CQCOF),RSNPH(JACIN),\r
-                    +NJIND,NQPTS)\r
-                     WRITE(*,1) 'DATA FOR SINGULAR INTEGRALS DONE:'\r
-               C\r
-               C**** SET UP THE A PRIORI COLUMN SCALE FACTORS, STORED IN COLSC.\r
-               C\r
-                     CALL CSCAL3(RWORK(COLSC),NQPTS,NJIND,RSNPH(ACOEF),RSNPH(BCOEF),\r
-                    +RSNPH(H0VAL),RSNPH(QUPTS),RSNPH(QUWTS),RSNPH(JACIN),RWORK(WORK),\r
-                    +RWORK(WORKT),RWORK(WORKQ))\r
-               C\r
-               C**** SET UP THE ARRAY RIGLL OF REFERENCE IGNORE LEVELS.\r
-               C\r
-                     CALL IGNLVL(RWORK(RIGLL),RWORK(COLSC),RSNPH(ACOEF),RSNPH(BCOEF),\r
-                    +RSNPH(H0VAL),RSNPH(JACIN),NJIND,NQPTS,IER)\r
-                     IF (IER .GT. 0) THEN\r
+               ASQUC7(AQCOF,BQCOF,CQCOF,JACIN,NJIND,NQPTS);\r
+               System.out.println("DATA FOR SINGULAR INTEGRALS DONE:");\r
+               \r
+               //**** SET UP THE A PRIORI COLUMN SCALE FACTORS, STORED IN COLSC.\r
+               \r
+               CSCAL3(COLSC,NQPTS,NJIND,ACOEF,BCOEF,\r
+                    H0VAL,QUPTS,QUWTS,JACIN,WORK,\r
+                    WORKT,WORKQ);\r
+               \r
+               // **** SET UP THE ARRAY RIGLL OF REFERENCE IGNORE LEVELS.\r
+               \r
+               IGNLVL(RIGLL,COLSC,ACOEF,BCOEF,H0VAL,JACIN,NJIND,NQPTS,IER);\r
+               /*      IF (IER .GT. 0) THEN\r
                        GOTO 999\r
                      ENDIF\r
                C\r
@@ -3736,6 +3672,463 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
            result=5E-1*(B0+A[0]-B2)*FACTOR;\r
            return result;\r
     } // private double GAMMA\r
+    \r
+    private void ASQUC7(double AQCOF[], double BQCOF[], double CQCOF[],\r
+               double JACIN[], int NJIND, int NQPTS) {\r
+        //REAL AQCOF(*),BQCOF(*),CQCOF(*),JACIN(*)\r
+\r
+        // ..TO ASSIGN THE COEFFICIENTS A(J), B(J) AND C(J) ,J=1,MN IN THE\r
+        // ..3-TERM RECURRENCE FORMULA \r
+        // ..\r
+        // ..  Q   (Z) = (A(J)Z - B(J))Q (Z) - C(J)Q   (Z) , J=2,...,MN\r
+        // ..   J+1                     J           J-1\r
+        // ..\r
+        // ..    Q (Z) = (A(1)Z - B(1))Q (Z) - C(1)\r
+        // ..     2                     1 \r
+        // ..\r
+        // ..WHERE Q (Z):=<P ,LOG(Z-.)>, P  IS THE ORTHONORMAL JACOBI POLY\r
+        // ..       J       J             J\r
+        // ..OF DEGREE J AND THE INNER PRODUCT IS WITH RESPECT TO THE JACOBI\r
+        // ..DISTRIBUTION OVER (-1,1).  HERE A(J,I) STORES "A(J)" FOR THE ITH\r
+        // ..ARC ON THE BOUNDARY (WITH SIMILAR ROLE FOR ARRAYS B AND C) AND\r
+        // ..THE JACOBI WEIGHT FUNCTION FOR THE ITH ARC IS \r
+        // ..(1-X)**AB(I,1)*(1+X)**AB(I,2), J=1,MN, I=1,NA.  \r
+\r
+        // **** AUTHOR: DAVID HOUGH\r
+        // **** LAST UPDATE: 15.09.89\r
+\r
+        // **** LOCAL VARIABLES        \r
+        int I,J,K,LOSUB;\r
+        double BE,D,N,N1,N2,F;\r
+        double H[] = new double[1];\r
+        // EXTERNAL ASONJ7\r
+\r
+      for (I=1; I <= NJIND; I++) {\r
+          LOSUB=(I-1)*NQPTS+1;\r
+          BE=JACIN[I-1];\r
+          double A[] = new double[NQPTS];\r
+          double B[] = new double[NQPTS];\r
+          ASONJ7(1.0,BE+1.0,A,B,H,NQPTS);\r
+          for (K = 0; K < NQPTS; K++) {\r
+                 AQCOF[LOSUB+K-1] = A[K];\r
+                 BQCOF[LOSUB+K-1] = B[K];\r
+          }\r
+          for (K=1; K <= NQPTS; K++) {\r
+              J=LOSUB+K-1;\r
+              N=(double)(K);\r
+              D=(N+1.0)*(N+BE+2.0);\r
+              N1=N*(N+BE+1.0);\r
+              N2=(N-1.0)*(N+BE);\r
+              F=Math.sqrt(N1/D);\r
+              AQCOF[J-1]=F/AQCOF[J-1];\r
+              BQCOF[J-1]=BQCOF[J-1]*AQCOF[J-1];\r
+              if (K > 1) {\r
+                  CQCOF[J-1]=AQCOF[J-1]*N2/AQCOF[J-2]/N1;\r
+              }\r
+              else {\r
+                  CQCOF[J-1]=-AQCOF[J-1]*Math.sqrt(H[0]/N1);\r
+              }\r
+          } // for (K=1; K <= NQPTS; K++)\r
+      } // for (I=1; I <= NJIND; I++)\r
+\r
+    } // private void ASQUC7\r
+    \r
+    private void CSCAL3(double COLSC[], int NQPTS,int NJIND, double ACOEF[], double BCOEF[],\r
+               double H0VAL[], double QUPTS[], double QUWTS[], double JACIN[], double MU[],\r
+               double TT[], double QQ[]) {\r
+        // REAL COLSC(*),ACOEF(*),BCOEF(*),H0VAL(*),QUPTS(*),QUWTS(*),\r
+       // +JACIN(*),MU(*),TT(*),QQ(*)\r
+       \r
+       // TO SET UP THE A PRIORI COLUMN SCALE FACTORS USING TRUNCATED\r
+       // CHEBYSHEV EXAPANSIONS FOR THE LOGARITHMIC KERNEL, GAUSS-JACOBI\r
+       // QUADRATURE AND GAUSS-JACOBI TEST POINTS.\r
+       \r
+       // LOCAL VARIABLES\r
+       \r
+        int I,J,J1,JI,K,K1,KMAX,LO,LO1,M,N;\r
+       double BETA,ROOTH,P0SCL,X,MAXVL,SUM1,SUM2;\r
+       // EXTERNAL JAPAR7\r
+       \r
+       K1=0;\r
+       MU[0]=-Math.log(2.0);\r
+       for (I=2; I <= 2*NQPTS; I++) {\r
+           MU[I-1]=-2.0/(double)(I-1);\r
+       }\r
+       \r
+       for (JI=1; JI <= NJIND; JI++) {\r
+           BETA=JACIN[JI-1];\r
+           ROOTH=Math.sqrt(H0VAL[JI-1]);\r
+           P0SCL=1.0/ROOTH;\r
+           LO=(JI-1)*NQPTS;\r
+           LO1=LO+1;\r
+       \r
+           for (J=1; J <= NQPTS; J++) {\r
+               X=QUPTS[LO+J-1];\r
+               QQ[(J-1)*NQPTS]=P0SCL;\r
+               double PP[] = new double[NQPTS-1];\r
+               double AA[] = new double[NQPTS-1];\r
+               double BB[] = new double[NQPTS-1];\r
+               for (N = 0; N < NQPTS-1; N++) {\r
+                       PP[N] = QQ[(J-1)*NQPTS+N];\r
+                       AA[N] = ACOEF[LO1+N-1];\r
+                       BB[N] = BCOEF[LO1+N-1];\r
+               }\r
+               JAPAR7(PP,X,AA,BB,NQPTS-1);\r
+               for (N = 0; N < NQPTS-1; N++) {\r
+                       QQ[(J-1)*NQPTS+N] = PP[N];\r
+               }\r
+               TT[J-1]=1.0;\r
+               TT[J+NQPTS-1]=X;\r
+               for (K=3; K <= 2*NQPTS; K++) {\r
+                   TT[J+(K-1)*NQPTS-1]=2.0*X*TT[J+(K-2)*NQPTS-1]-TT[J+(K-3)*NQPTS-1];\r
+               }\r
+           } // for (J=1; J <= NQPTS; J++)\r
+       \r
+           for (K=1; K <= NQPTS; K++) {\r
+               MAXVL=0.0;\r
+               for (I=1; I <= NQPTS; I++) {\r
+                   SUM2=0.0;\r
+                   KMAX=2*NQPTS+1-K;\r
+                   for (M=K; M <= KMAX; M++) {\r
+                       SUM1=0.0;\r
+                       for (J=1; J <= NQPTS; J++) {\r
+                           J1=LO+J;\r
+                           SUM1=SUM1+QUWTS[J1-1]*QQ[K+(J-1)*NQPTS-1]*TT[J+(M-1)*NQPTS-1];\r
+                       } // for (J=1; J <= NQPTS; J++)\r
+                       SUM2=SUM2+MU[M-1]*TT[I+(M-1)*NQPTS-1]*SUM1;\r
+                   } // for (M=K; M <= KMAX; M++)\r
+                   MAXVL=Math.max(MAXVL,Math.abs(SUM2));\r
+               } // for (I=1; I <= NQPTS; I++)\r
+               K1=K1+1;\r
+               COLSC[K1-1]=1.0/MAXVL;\r
+           } // for (K=1; K <= NQPTS; K++)\r
+       \r
+       } // for (JI=1; JI <= NJIND; JI++)\r
+       \r
+    }\r
+\r
+    private void JAPAR7(double PP[] ,double X, double AA[], double BB[], int N) {\r
+        // REAL X,PP(*),AA(*),BB(*)\r
+\r
+        // ..................................................................\r
+\r
+        // 1.  JAPAR7\r
+        //     EVALUATES ORTHONORMAL JACOBI POLYNOMIALS.\r
+\r
+\r
+       // 2.  PURPOSE\r
+       //     TO COMPUTE PP(I), I=1,..,N+1, GIVEN PP(1) ON ENTRY, WHERE PP(I)\r
+       //     STORES THE VALUE OF THE ORTHONORMAL JACOBI POLYNOMIAL OF DEGREE\r
+       //     I-1 AT THE GIVEN PARAMETER VALUE X. \r
+\r
+\r
+       // 3.  CALLING SEQUENCE\r
+       //     CALL JAPAR7(PP,X,AA,BB,N)\r
+\r
+       // PARAMETERS (SEE DECLARATIONS ABOVE FOR TYPES)\r
+       //  ON ENTRY:\r
+       //  PP(1)  - THE VALUE OF THE POLYNOMIAL OF DEGREE ZERO.\r
+\r
+       //  X      - THE REAL NUMBER AT WHICH THE POLYNOMIALS ARE TO BE\r
+       //           EVALUATED.\r
+\r
+       //  AA     - ARRAY OF COEFFICIENTS IN THE 3-TERM RECURRENCE\r
+       //           AA(I)*PP(I+1)=(X-BB(I))*PP(I)-AA(I-1)*PP(I-1),\r
+       //           I=1,..,N, WITH "PP(0)" = 0.\r
+\r
+       //  BB     - ARRAY OF COEFFICIENTS  IN THE 3-TERM RECURRENCE ABOVE.\r
+\r
+       //  N      - THE LARGEST DEGREE OF POLYNOMIAL REQUIRING EVALUATION.\r
+\r
+       //  ON EXIT:\r
+       //  PP     - PP(I) IS THE VALUE OF THE POLYNOMIAL OF DEGREE I-1 AT\r
+       //           X, I=1,..,N+1.\r
+\r
+\r
+       // 4.  NO SUBROUTINES OR FUNCTIONS NEEDED\r
+        \r
+       // ..................................................................\r
+   \r
+       int I;\r
+\r
+       if (N > 0) {\r
+           PP[1]=(X-BB[0])*PP[0]/AA[0];\r
+           for (I=2; I <= N; I++) {\r
+               PP[I]=((X-BB[I-1])*PP[I-1]-AA[I-2]*PP[I-2])/AA[I-1];\r
+           }\r
+       }\r
+\r
+    } // private void JAPAR7\r
+\r
+    private void IGNLVL(double RIGLL[], double COLSC[], double ACOEF[], double BCOEF[],\r
+               double H0VAL[], double JACIN[], int NJIND, int NQPTS, int IER[]) {\r
+        // REAL ACOEF(*),BCOEF(*),H0VAL(*),JACIN(*),RIGLL(*),COLSC(*)\r
+       \r
+       // TO SET UP THE A PRIORI COLUMN SCALE FACTORS USING THE COEFFICIENT\r
+       // IGNORE LEVELS OBTAINED FROM BOUNDARY CORRESPONDENCE FUNCTION \r
+       // REPRESENTATION FOR THE BOUNDARY MAP; SEE RB #50 P141 ET SEQ.\r
+       // COMBINE WITH A PRIORI COLUMN SCALE FACTORS COLSC IN THE CASE OF\r
+       // SOLVING SYMM'S EQUATION.\r
+       \r
+       // IER=0 - NORMAL EXIT\r
+       // IER=5 - LOCAL PARAMETER MNQPT MUST BE INCREASED TO AT LEAST THE \r
+       //         VALUE OF NQPTS\r
+       // IER=6 - FAILURE OF IMTQLH EIGENVALUE ROUTINE\r
+       // IER=53- A CORNER ANGLE IS SO SMALL THAT IT MAY CAUSE OVERFLOW\r
+       \r
+       // LOCAL VARIABLES\r
+       \r
+       final int MNQPT = 20;\r
+       int JT,K,K1,LO,N;\r
+       double B1,B2,BETA,EXPON,H1VAL,LGTWO,OFLOW,PI,RH,RH1,SUP,T;\r
+       boolean SYMM;\r
+       double A1COF[] = new double[MNQPT];\r
+       double B1COF[] = new double[MNQPT];\r
+       double DIAG[] = new double[MNQPT];\r
+       double PP[] = new double[MNQPT];\r
+       double SDIAG[] = new double[MNQPT];\r
+       // EXTERNAL ASONJ7,IMTQLH,JAPAR7\r
+       \r
+       if (NQPTS > MNQPT) {\r
+           IER[0]=5;\r
+           return;\r
+       }\r
+       \r
+       if (COLSC[0] <= 0.0) {\r
+           SYMM= false;\r
+       }\r
+       else {\r
+           SYMM= true;\r
+       }\r
+       \r
+       K1=0;\r
+       PI= Math.PI;\r
+       LGTWO=Math.log(2.0);\r
+       OFLOW=Math.log(Double.MAX_VALUE);\r
+       /*for (JT=1; JT <= NJIND; JT++) {\r
+           BETA=JACIN[JT-1];\r
+           B1=BETA+1.0;\r
+           B2=BETA+2.0;\r
+           RH=Math.sqrt(H0VAL[JT-1]);\r
+           ASONJ7(1.0,B1,A1COF,B1COF,H1VAL,NQPTS);\r
+           LO=(JT-1)*NQPTS;\r
+           K1=K1+1;\r
+       \r
+           // ****   COMPUTE THE QUANTITY\r
+           // ****       PI*2E+0**B2/B1/RH\r
+           // ****   BUT CHECK FOR POSSIBLE OVERFLOW\r
+       \r
+           EXPON=B2*LGTWO+Math.log(PI/B1/RH);\r
+           if (EXPON >= OFLOW) {\r
+               IER[0]=53;\r
+               return;\r
+           }\r
+           else {\r
+               if (SYMM) {\r
+                   RIGLL[K1-1]=Math.exp(EXPON)*COLSC[K1-1];\r
+               }\r
+               else {\r
+                   RIGLL[K1-1]=Math.exp(EXPON);\r
+               }\r
+           }\r
+       \r
+           if (NQPTS == 1) {\r
+               continue;\r
+           }\r
+       \r
+           RH1=Math.sqrt(H1VAL);\r
+           K1=K1+1;\r
+       \r
+           // ****   COMPUTE THE QUANTITY\r
+           // ****       PI*2E+0**(B2+1E+0)*B1**B1/RH1/B2**(B2+5E-1)\r
+           // ****   BUT CHECK FOR POSSIBLE OVERFLOW\r
+       \r
+           EXPON=(B2+1.0)*LGTWO+B1*Math.log(B1)-(B2+0.5)*Math.log(B2)+Math.log(PI/RH1);\r
+           if (EXPON >= OFLOW) {\r
+               IER[0]=53;\r
+               return;\r
+           }\r
+           else {\r
+               if (SYMM) {\r
+                   RIGLL[K1-1]=Math.exp(EXPON)*COLSC[K1-1];\r
+               }\r
+               else {\r
+                   RIGLL[K1-1]=Math.exp(EXPON);\r
+               }\r
+           }\r
+       \r
+           for (N=2; N <= NQPTS-1; N++) {\r
+       \r
+               // FIND THE ZEROES OF THE JACOBI POLYNOMIAL OF DEGREE N FOR\r
+               // WEIGHT (1+X)**BETA\r
+       \r
+               for (K=1; K <= N; K++) {\r
+                   DIAG[K-1]=BCOEF[LO+K-1];\r
+                   if (K == 1) {\r
+                       SDIAG[K-1]=0.0;\r
+                   }\r
+                   else {\r
+                       SDIAG[K-1]=ACOEF[LO+K-2];\r
+                   }\r
+               } // for (K=1; K <= N; K++)\r
+       \r
+               IMTQLH(N,DIAG,SDIAG,IER);\r
+                 IF (IER.GT.0) THEN\r
+                   IER=6\r
+                   RETURN\r
+                 ENDIF\r
+       C\r
+                 SUP=0E+0\r
+                 DO 20 K=1,N\r
+                   T=DIAG(K)\r
+                   PP(1)=1E+0/RH1\r
+                   CALL JAPAR7(PP,T,A1COF,B1COF,N)\r
+                   T=(1E+0-T)*(1E+0+T)**B1*PP(N)\r
+                   SUP=MAX(SUP,ABS(T))\r
+       20        CONTINUE\r
+       C\r
+                 K1=K1+1\r
+                 IF (SYMM) THEN\r
+                   RIGLL(K1)=2E+0*PI*SUP*COLSC(K1)/SQRT(N*(N+B1))\r
+                 ELSE\r
+                   RIGLL(K1)=2E+0*PI*SUP/SQRT(N*(N+B1))\r
+                 ENDIF\r
+           } // for (N=2; N <= NQPTS-1; N++)\r
+       } // for (JT=1; JT <= NJIND; JT++)\r
+       \r
+       // NORMAL EXIT\r
+       \r
+       IER[0]=0;*/\r
+       \r
+    } // private void IGNLVL\r
+    \r
+    private void IMTQLH(int N, double D[], double E[], int IERR[]) {\r
+    \r
+        int I,J,L,M,II,MML;\r
+        //REAL D(N),E(N)\r
+        double B,C,F,G,P,R,S,TST1,TST2,PYTHAG;\r
+    \r
+        // THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE IMTQL1,\r
+        // NUM. MATH. 12, 377-383(1968) BY MARTIN AND WILKINSON,\r
+        // AS MODIFIED IN NUM. MATH. 15, 450(1970) BY DUBRULLE.\r
+        // HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 241-248(1971).\r
+    \r
+        // THIS SUBROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC\r
+        // TRIDIAGONAL MATRIX BY THE IMPLICIT QL METHOD.\r
+    \r
+        // ON INPUT\r
+    \r
+        //     N IS THE ORDER OF THE MATRIX.\r
+    \r
+        //     D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX.\r
+    \r
+        //     E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX\r
+        //       IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY.\r
+    \r
+        // ON OUTPUT\r
+    \r
+        //     D CONTAINS THE EIGENVALUES IN ASCENDING ORDER.  IF AN\r
+        //       ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND\r
+        //       ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE\r
+        //       THE SMALLEST EIGENVALUES.\r
+    \r
+        //    E HAS BEEN DESTROYED.\r
+    \r
+        //    IERR IS SET TO\r
+        //           ZERO       FOR NORMAL RETURN,\r
+        //           J          IF THE J-TH EIGENVALUE HAS NOT BEEN\r
+        //                      DETERMINED AFTER 30 ITERATIONS.\r
+    \r
+        // CALLS PYTHAG FOR  SQRT(A*A + B*B) .\r
+    \r
+        // QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,\r
+        // MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY\r
+    \r
+        // THIS VERSION ORIGINALLY DATED AUGUST 1983; RENDERED INACCURATE \r
+        // AND TRANSLATED INTO SINGLE PRECISION BY DAVID HOUGH, ETH, ZUERICH\r
+        // OCTOBER, 1989.\r
+    \r
+        // ------------------------------------------------------------------\r
+    \r
+        IERR[0] = 0;\r
+        if (N == 1) {\r
+            return;    \r
+        }\r
+    \r
+        for (I = 1; I <= N-1; I++) {\r
+            E[I-1] = E[I];\r
+        }\r
+        E[N-1] = 0.0;\r
+    \r
+        /*for (L = 1; L <= N; L++) {\r
+            J = 0;\r
+            // .......... LOOK FOR SMALL SUB-DIAGONAL ELEMENT ..........\r
+            while (true) {\r
+                for (M = L; M <= N; M++) {\r
+                IF (M .EQ. N) GO TO 120\r
+                TST1 = ABS(D(M)) + ABS(D(M+1))\r
+                TST2 = TST1 + ABS(E(M))\r
+                IF (TST2 .EQ. TST1) GO TO 120\r
+                } // for (M = L; M <= N; M++)\r
+    C\r
+      120    P = D(L)\r
+             IF (M .EQ. L) GO TO 215\r
+             IF (J .EQ. 30) GO TO 1000\r
+             J = J + 1\r
+    C     .......... FORM SHIFT ..........\r
+             G = (D(L+1) - P) / (2.0E0 * E(L))\r
+             R = PYTHAG(G,1.0E0)\r
+             G = D(M) - P + E(L) / (G + SIGN(R,G))\r
+             S = 1.0E0\r
+             C = 1.0E0\r
+             P = 0.0E0\r
+             MML = M - L\r
+    C     .......... FOR I=M-1 STEP -1 UNTIL L DO -- ..........\r
+             DO 200 II = 1, MML\r
+                I = M - II\r
+                F = S * E(I)\r
+                B = C * E(I)\r
+                R = PYTHAG(F,G)\r
+                E(I+1) = R\r
+                IF (R .EQ. 0.0E0) GO TO 210\r
+                S = F / R\r
+                C = G / R\r
+                G = D(I+1) - P\r
+                R = (D(I) - G) * S + 2.0E0 * C * B\r
+                P = S * R\r
+                D(I+1) = G + P\r
+                G = C * R - B\r
+      200    CONTINUE\r
+    C\r
+             D(L) = D(L) - P\r
+             E(L) = G\r
+             E(M) = 0.0E0\r
+             GO TO 105\r
+    C     .......... RECOVER FROM UNDERFLOW ..........\r
+      210    D(I+1) = D(I+1) - P\r
+             E(M) = 0.0E0\r
+             GO TO 105\r
+    C     .......... ORDER EIGENVALUES ..........\r
+      215    IF (L .EQ. 1) GO TO 250\r
+    C     .......... FOR I=L STEP -1 UNTIL 2 DO -- ..........\r
+             DO 230 II = 2, L\r
+                I = L + 2 - II\r
+                IF (P .GE. D(I-1)) GO TO 270\r
+                D(I) = D(I-1)\r
+      230    CONTINUE\r
+    C\r
+      250    I = 1\r
+      270    D(I) = P\r
+            } // while (true)\r
+        } // for (L = 1; L <= N; L++)\r
+    C\r
+          GO TO 1001\r
+    C     .......... SET ERROR -- NO CONVERGENCE TO AN\r
+    C                EIGENVALUE AFTER 30 ITERATIONS ..........\r
+     1000 IERR = L\r
+     1001 RETURN*/\r
+    }\r
+\r
 \r
 \r
             /* SUBROUTINE RSLT71(QIERC,RCOND,SOLUN,NEQNS,LOSUB,HISUB,COLSC,\r