Fix issue #133

tzaeschke · Apr 21, 2024 · de839f1 · de839f1
1 parent ec56f9e
commit de839f1
Show file tree

Hide file tree

Showing 2 changed files with 39 additions and 34 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -25,6 +25,9 @@
 
 --> See TODO.txt
 
+## 0.5.3-SNAPSHOT
+- Fix issue in DCLP.solve1() causing ODE INTERNAL ERROR. [#133](https://github.com/tzaeschke/ode4j/issues/133)
+
 ## 0.5.2 - 2023-10-07
 - Fix DVector3.cross() returning always 0. [#128](https://github.com/tzaeschke/ode4j/issues/128)
 - Fix JUnit test lint. [#130](https://github.com/tzaeschke/ode4j/pull/130)

diff --git a/core/src/main/java/org/ode4j/ode/internal/DLCP.java b/core/src/main/java/org/ode4j/ode/internal/DLCP.java
@@ -916,41 +916,43 @@ void solve1 (double[] a, int i, boolean dir_positive, boolean only_transfer)
 
 		final int nC = m_nC;
 		if (nC > 0) {
-		    double[] Dell = m_Dell;
-		    int[] C = m_C;
-			//double[] aptr = AROW(i);
-			int aPos = AROWp(i);
-			if (NUB_OPTIMIZATIONS) {//#   ifdef NUB_OPTIMIZATIONS
-				// if nub>0, initial part of aptr[] is guaranteed unpermuted
-			    final int nub = m_nub;
-			    int j=0;
-				for ( ; j<nub; ++j) Dell[j] = m_A[aPos+j];//aptr[j];
-				for ( ; j<nC; ++j) Dell[j] = m_A[aPos+C[j]];//aptr[C[j]];
-			} else {//#   else
-        		for (int j=0; j<nC; j++) Dell[j] = m_A[aPos+C[j]];//aptr[C[j]];
-			}//#   endif
-		}
-		solveL1Straight(m_L, m_Dell, 0, nC,  m_nskip, 1);
-		{
-		    double[] ell = m_ell, Dell = m_Dell, d = m_d;
-		    for (int j=0; j<nC; j++) ell[j] = Dell[j] * d[j];
-		}
+			{
+				double[] Dell = m_Dell;
+				int[] C = m_C;
+				//double[] aptr = AROW(i);
+				int aPos = AROWp(i);
+				if (NUB_OPTIMIZATIONS) {//#   ifdef NUB_OPTIMIZATIONS
+					// if nub>0, initial part of aptr[] is guaranteed unpermuted
+					final int nub = m_nub;
+					int j = 0;
+					for (; j < nub; ++j) Dell[j] = m_A[aPos + j];//aptr[j];
+					for (; j < nC; ++j) Dell[j] = m_A[aPos + C[j]];//aptr[C[j]];
+				} else {//#   else
+					for (int j = 0; j < nC; j++) Dell[j] = m_A[aPos + C[j]];//aptr[C[j]];
+				}//#   endif
+			}
+			solveL1Straight(m_L, m_Dell, 0, nC, m_nskip, 1);
+			{
+				double[] ell = m_ell, Dell = m_Dell, d = m_d;
+				for (int j = 0; j < nC; j++) ell[j] = Dell[j] * d[j];
+			}
 
-		if (!only_transfer) {
-		    double[] tmp = m_tmp, ell = m_ell;
-		    {
-		        for (int j=0; j<nC; ++j) tmp[j] = ell[j];
-		    }
-			solveL1Transposed(m_L,tmp,0, nC,m_nskip, 1);
-		    if (dir_positive) {
-		        int[] C = m_C;
-		        //double[] tmp = m_tmp;
-		        for (int j=0; j<nC; ++j) a[C[j]] = -tmp[j];
-		    } else {
-		        int[] C = m_C;
-		        //double[] tmp = m_tmp;
-		        for (int j=0; j<nC; ++j) a[C[j]] = tmp[j];
-		    }
+			if (!only_transfer) {
+				double[] tmp = m_tmp, ell = m_ell;
+				{
+					for (int j = 0; j < nC; ++j) tmp[j] = ell[j];
+				}
+				solveL1Transposed(m_L, tmp, 0, nC, m_nskip, 1);
+				if (dir_positive) {
+					int[] C = m_C;
+					//double[] tmp = m_tmp;
+					for (int j = 0; j < nC; ++j) a[C[j]] = -tmp[j];
+				} else {
+					int[] C = m_C;
+					//double[] tmp = m_tmp;
+					for (int j = 0; j < nC; ++j) a[C[j]] = tmp[j];
+				}
+			}
 		}
 	}