EPS now set in public SymmsIntegralMapping(). In CSCAL3 changed double PP[] = new...
authorilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Thu, 15 Feb 2018 19:17:07 +0000 (19:17 +0000)
committerilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Thu, 15 Feb 2018 19:17:07 +0000 (19:17 +0000)
git-svn-id: https://citdcbmipav.cit.nih.gov/repos-pub/mipav/trunk@15375 ba61647d-9d00-f842-95cd-605cb4296b96

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

index 2427239..0005bbc 100644 (file)
@@ -484,6 +484,26 @@ public class SymmsIntegralMapping extends AlgorithmBase {
                HIVLC = new double[NJIND];\r
                JAINC = new double[NJIND];\r
                COARG = new double[NJIND];\r
+               // eps returns the distance from 1.0 to the next larger double-precision\r
+               // number, that is, eps = 2^-52.\r
+               // epsilon = D1MACH(4)\r
+               // Machine epsilon is the smallest positive epsilon such that\r
+               // (1.0 + epsilon) != 1.0.\r
+               // epsilon = 2**(1 - doubleDigits) = 2**(1 - 53) = 2**(-52)\r
+               // epsilon = 2.2204460e-16\r
+               // epsilon is called the largest relative spacing\r
+               EPS = 1.0;\r
+               double neweps = 1.0;\r
+\r
+               while (true) {\r
+\r
+                       if (1.0 == (1.0 + neweps)) {\r
+                               break;\r
+                       } else {\r
+                               EPS = neweps;\r
+                               neweps = neweps / 2.0;\r
+                       }\r
+               } // while(true)\r
        }\r
 \r
        public SymmsIntegralMapping(ModelImage destImg, ModelImage srcImg, String FORTFL, boolean SYMTY, boolean REFLN,\r
@@ -3920,10 +3940,10 @@ public class SymmsIntegralMapping extends AlgorithmBase {
                        double anorm;\r
                        int iwork[] = new int[NROWS];\r
                        int info[] = new int[1];\r
-                       anorm = ge.dlange('1', MNEQN, MNEQN, A, MNEQN, WORK2);\r
-                       le2.dgetrf(MNEQN, MNEQN, A, MNEQN, IPIVT, info);\r
+                       anorm = ge.dlange('1', NROWS, NROWS, A, MNEQN, WORK2);\r
+                       le2.dgetrf(NROWS, NROWS, A, MNEQN, IPIVT, info);\r
                        if (info[0] > 0) {\r
-                               MipavUtil.displayError("In dgetrf factor U is exactly singular");\r
+                               System.err.println("In dgetrf factor U is exactly singular");\r
                        }\r
                        le2.dgecon('1', NROWS, A, MNEQN, anorm, RCOND, WORK2, iwork, info);\r
                        // Reciprocal condition number of A in 1-norm.\r
@@ -4901,16 +4921,18 @@ public class SymmsIntegralMapping extends AlgorithmBase {
                        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 PP[] = new double[NQPTS];\r
                                double AA[] = new double[NQPTS - 1];\r
                                double BB[] = new double[NQPTS - 1];\r
-                               for (N = 0; N < NQPTS - 1; N++) {\r
+                               for (N = 0; N < NQPTS; N++) {\r
                                        PP[N] = QQ[(J - 1) * NQPTS + N];\r
+                               }\r
+                               for (N = 0; N < NQPTS - 1; 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
+                               for (N = 0; N < NQPTS; N++) {\r
                                        QQ[(J - 1) * NQPTS + N] = PP[N];\r
                                }\r
                                TT[J - 1] = 1.0;\r