author ilb@NIH.GOV Thu, 15 Feb 2018 19:17:07 +0000 (19:17 +0000) committer ilb@NIH.GOV 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

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;\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) {\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