From 6177c3efeb9d0fe9c41c983bdb9864573baa00b2 Mon Sep 17 00:00:00 2001 From: rwl Date: Wed, 2 May 2012 02:40:40 +0100 Subject: [PATCH] Fixing pointers up to Dklu_solve. --- .../edu/ufl/cise/klu/tdouble/Dklu_diagnostics.java | 97 ++++++++++--------- .../java/edu/ufl/cise/klu/tdouble/Dklu_dump.java | 17 ++-- .../edu/ufl/cise/klu/tdouble/Dklu_extract.java | 49 ++++++---- .../java/edu/ufl/cise/klu/tdouble/Dklu_factor.java | 29 +++--- .../java/edu/ufl/cise/klu/tdouble/Dklu_memory.java | 5 +- .../edu/ufl/cise/klu/tdouble/Dklu_refactor.java | 104 +++++++++++++-------- .../java/edu/ufl/cise/klu/tdouble/Dklu_solve.java | 64 ++++++------- .../java/edu/ufl/cise/klu/tdouble/Dklu_tsolve.java | 66 ++++++------- .../edu/ufl/cise/klu/tdouble/Dklu_version.java | 14 +-- 9 files changed, 245 insertions(+), 200 deletions(-) diff --git a/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_diagnostics.java b/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_diagnostics.java index 31b38fe..da6d522 100644 --- a/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_diagnostics.java +++ b/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_diagnostics.java @@ -63,11 +63,10 @@ public class Dklu_diagnostics extends Dklu_internal double[] Aentry, Ux, Ukk ; double[] Rs ; int i, newrow, oldrow, k1, k2, nk, j, oldcol, k, pend ; - int LU_offset = 0, Uip_offset = 0, Ulen_offset = 0, Ukk_offset = 0 ; + int Uip_offset = 0, Ulen_offset = 0, Ukk_offset = 0 ; int[] len = new int[1] ; int[] Ui_offset = new int[1] ; int[] Ux_offset = new int[1] ; - int[] Ui ; /* ---------------------------------------------------------------------- */ /* check inputs */ @@ -112,8 +111,7 @@ public class Dklu_diagnostics extends Dklu_internal { continue ; /* skip singleton blocks */ } - LU = Numeric.LUbx ; - LU_offset = i ; + LU = Numeric.LUbx[i] ; Uip = Numeric.Uip ; Uip_offset += k1 ; Ulen = Numeric.Ulen ; @@ -138,8 +136,8 @@ public class Dklu_diagnostics extends Dklu_internal ASSERT (newrow < k2) ; if (Rs != null) { - aik = Aentry [k] / Rs [newrow] ; //SCALE_DIV_ASSIGN (aik, Aentry [k], Rs [newrow]) ; + aik = Aentry [k] / Rs [newrow] ; } else { @@ -153,15 +151,11 @@ public class Dklu_diagnostics extends Dklu_internal } } - GET_POINTER (LU, LU_offset, - Uip, Uip_offset, - Ulen, Ulen_offset, - Ui, Ui_offset, - Ux, Ux_offset, - j, len) ; + Ux = GET_POINTER (LU, Uip, Uip_offset, Ulen, Ulen_offset, + Ui_offset, Ux_offset, j, len) ; for (k = 0 ; k < len[0] ; k++) { - temp = ABS (Ux [k]) ; + temp = ABS (Ux [Ux_offset[0] + k]) ; //ABS (temp, Ux [k]) ; if (temp > max_ui) { @@ -169,7 +163,7 @@ public class Dklu_diagnostics extends Dklu_internal } } /* consider the diagonal element */ - temp = ABS (Ukk [j]) ; + temp = ABS (Ukk [Ukk_offset + j]) ; //ABS (temp, Ukk [j]) ; if (temp > max_ui) { @@ -215,8 +209,7 @@ public class Dklu_diagnostics extends Dklu_internal { double xj, Xmax, csum, anorm, ainv_norm, est_old, est_new, abs_value ; double[] Udiag, Aentry, X, S ; - int[] R ; - int nblocks, i, j, jmax, jnew, pend, n ; + int i, j, jmax, jnew, pend, n ; int unchanged ; /* ---------------------------------------------------------------------- */ @@ -247,8 +240,6 @@ public class Dklu_diagnostics extends Dklu_internal /* ---------------------------------------------------------------------- */ n = Symbolic.n ; - nblocks = Symbolic.nblocks ; - R = Symbolic.R ; Udiag = Numeric.Udiag ; /* ---------------------------------------------------------------------- */ @@ -257,8 +248,8 @@ public class Dklu_diagnostics extends Dklu_internal for (i = 0 ; i < n ; i++) { - abs_value = ABS (Udiag [i]) ; //ABS (abs_value, Udiag [i]) ; + abs_value = ABS (Udiag [i]) ; if (SCALAR_IS_ZERO (abs_value)) { Common.condest = 1 / abs_value ; @@ -279,8 +270,8 @@ public class Dklu_diagnostics extends Dklu_internal csum = 0.0 ; for (j = Ap [i] ; j < pend ; j++) { - abs_value = ABS (Aentry [j]) ; //ABS (abs_value, Aentry [j]) ; + abs_value = ABS (Aentry [j]) ; csum += abs_value ; } if (csum > anorm) @@ -296,16 +287,17 @@ public class Dklu_diagnostics extends Dklu_internal /* get workspace (size 2*n double's) */ X = Numeric.Xwork ; /* size n space used in KLU_solve, tsolve */ //X += n ; /* X is size n */ - X = new double [n] ; + int X_offset = n ; //S = X + n ; /* S is size n */ - S = new double [n] ; + S = X ; + int S_offset = 2*n ; for (i = 0 ; i < n ; i++) { - CLEAR (S, i) ; - CLEAR (X, i) ; - X [i] = 1.0 / ((double) n); + CLEAR (S, S_offset + i) ; + CLEAR (X, X_offset + i) ; //REAL (X [i]) = 1.0 / ((double) n) ; + X [X_offset + i] = 1.0 / ((double) n); } jmax = 0 ; @@ -317,21 +309,21 @@ public class Dklu_diagnostics extends Dklu_internal /* X [jmax] is the largest entry in X */ for (j = 0 ; j < n ; j++) { - CLEAR (X, j) ; + CLEAR (X, X_offset + j) ; } - X [jmax] = 1 ; //REAL (X [jmax]) = 1 ; + X [X_offset + jmax] = 1 ; } - Dklu_solve.klu_solve (Symbolic, Numeric, n, 1, (double[]) X, Common) ; + klu_solve (Symbolic, Numeric, n, 1, (double[]) X, X_offset, Common) ; est_old = ainv_norm ; ainv_norm = 0.0 ; for (j = 0 ; j < n ; j++) { /* ainv_norm += ABS (X [j]) ;*/ - abs_value = ABS (X [j]) ; //ABS (abs_value, X [j]) ; + abs_value = ABS (X [X_offset + j]) ; ainv_norm += abs_value ; } @@ -339,10 +331,10 @@ public class Dklu_diagnostics extends Dklu_internal for (j = 0 ; j < n ; j++) { - double s = (X [j] >= 0) ? 1 : -1 ; - if (s != S [j]) // s != REAL (S [j]) + double s = (X [X_offset + j] >= 0) ? 1 : -1 ; + if (s != S [S_offset + j]) // s != REAL (S [j]) { - S [j] = s ; + S [S_offset + j] = s ; unchanged = FALSE ; } } @@ -354,19 +346,19 @@ public class Dklu_diagnostics extends Dklu_internal for (j = 0 ; j < n ; j++) { - X [j] = S [j] ; + X [j] = S [S_offset + j] ; } /* do a transpose solve */ - klu_tsolve (Symbolic, Numeric, n, 1, X, Common) ; + klu_tsolve (Symbolic, Numeric, n, 1, X, X_offset, Common) ; /* jnew = the position of the largest entry in X */ jnew = 0 ; Xmax = 0 ; for (j = 0 ; j < n ; j++) { - xj = ABS (X [j]) ; //ABS (xj, X [j]) ; + xj = ABS (X [X_offset + j]) ; if (xj > Xmax) { Xmax = xj ; @@ -388,27 +380,27 @@ public class Dklu_diagnostics extends Dklu_internal for (j = 0 ; j < n ; j++) { - CLEAR (X, j) ; + CLEAR (X, X_offset + j) ; if (j % 2 != 0) { - X [j] = 1 + ((double) j) / ((double) (n-1)) ; //REAL (X [j]) = 1 + ((double) j) / ((double) (n-1)) ; + X [X_offset + j] = 1 + ((double) j) / ((double) (n-1)) ; } else { - X [j] = -1 - ((double) j) / ((double) (n-1)) ; //REAL (X [j]) = -1 - ((double) j) / ((double) (n-1)) ; + X [X_offset + j] = -1 - ((double) j) / ((double) (n-1)) ; } } - klu_solve (Symbolic, Numeric, n, 1, (double[]) X, Common) ; + klu_solve (Symbolic, Numeric, n, 1, (double[]) X, X_offset, Common) ; est_new = 0.0 ; for (j = 0 ; j < n ; j++) { /* est_new += ABS (X [j]) ;*/ - abs_value = ABS (X [j]) ; //ABS (abs_value, X [j]) ; + abs_value = ABS (X [X_offset + j]) ; est_new += abs_value ; } est_new = 2 * est_new / (3 * n) ; @@ -434,10 +426,12 @@ public class Dklu_diagnostics extends Dklu_internal KLU_common Common) { double flops = 0 ; - int[] R, Ui, Uip, Llen, Ulen ; + int[] R, Uip, Llen, Ulen ; + /*int[]*/double[] Ui ; double[][] LUbx ; double[] LU ; - int k, ulen, p, n, nk, block, nblocks, k1 ; + int k, ulen, p, nk, block, nblocks, k1 ; + int Llen_offset = 0, Uip_offset = 0, Ulen_offset = 0 ; /* ---------------------------------------------------------------------- */ /* check inputs */ @@ -459,7 +453,6 @@ public class Dklu_diagnostics extends Dklu_internal /* get the contents of the Symbolic object */ /* ---------------------------------------------------------------------- */ - n = Symbolic.n ; R = Symbolic.R ; nblocks = Symbolic.nblocks ; @@ -479,21 +472,25 @@ public class Dklu_diagnostics extends Dklu_internal nk = R [block+1] - k1 ; if (nk > 1) { - Llen = Numeric.Llen + k1 ; - Uip = Numeric.Uip + k1 ; - Ulen = Numeric.Ulen + k1 ; + Llen = Numeric.Llen ; + Llen_offset += k1 ; + Uip = Numeric.Uip ; + Uip_offset += k1 ; + Ulen = Numeric.Ulen ; + Ulen_offset += k1 ; LU = LUbx [block] ; + int[] Ui_offset = new int[1] ; for (k = 0 ; k < nk ; k++) { /* compute kth column of U, and update kth column of A */ - GET_I_POINTER (LU, Uip, Ui, k) ; - ulen = Ulen [k] ; + Ui = GET_I_POINTER (LU, Uip, Uip_offset, Ui_offset, k) ; + ulen = Ulen [Ulen_offset + k] ; for (p = 0 ; p < ulen ; p++) { - flops += 2 * Llen [Ui [p]] ; + flops += 2 * Llen [Llen_offset + (int) Ui [Ui_offset[0] + p]] ; } /* gather and divide by pivot to get kth column of L */ - flops += Llen [k] ; + flops += Llen [Llen_offset + k] ; } } } @@ -548,8 +545,8 @@ public class Dklu_diagnostics extends Dklu_internal for (j = 0 ; j < n ; j++) { /* get the magnitude of the pivot */ - ukk = ABS (Udiag [j]) ; //ABS (ukk, Udiag [j]) ; + ukk = ABS (Udiag [j]) ; if (SCALAR_IS_NAN (ukk) || SCALAR_IS_ZERO (ukk)) { /* if NaN, or zero, the rcond is zero */ diff --git a/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_dump.java b/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_dump.java index 55f985a..fe931e9 100644 --- a/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_dump.java +++ b/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_dump.java @@ -31,7 +31,6 @@ package edu.ufl.cise.klu.tdouble; public class Dklu_dump extends Dklu_internal { - /** * Check if a column-form matrix is valid or not. The matrix A is * n-by-n. The row indices of entries in column j are in @@ -104,9 +103,12 @@ public class Dklu_dump extends Dklu_internal int[] Xip, int Xip_offset, int[] Xlen, int Xlen_offset, double[] LU) { - int[] Xi ; + /*int[]*/double[] Xi ; double[] Xx ; - int j, p1, p2, i, p, len ; + int j, p1, p2, i, p ; + int[] len = new int[1] ; + int[] Xi_offset = new int[1] ; + int[] Xx_offset = new int[1] ; PRINTF ("\ncolumn oriented matrix, n = %d\n", n) ; if (n <= 0) @@ -132,10 +134,11 @@ public class Dklu_dump extends Dklu_internal PRINTF ("column %d pointer bad\n", j) ; return (FALSE) ; } - GET_POINTER (LU, Xip, Xip_offset, Xlen, Xlen_offset, Xi, Xx, j, len) ; - for (p = 0 ; p < len ; p++) + Xi = Xx = GET_POINTER (LU, Xip, Xip_offset, Xlen, Xlen_offset, + Xi_offset, Xx_offset, j, len) ; + for (p = 0 ; p < len[0] ; p++) { - i = Xi [p] ; + i = (int) Xi [Xi_offset[0] + p] ; PRINTF ("row: %d", i) ; if (i < 0 || i >= n) { @@ -145,7 +148,7 @@ public class Dklu_dump extends Dklu_internal } if (Xx != null) { - PRINT_ENTRY (Xx [p]) ; + PRINT_ENTRY (Xx [Xx_offset[0] + p]) ; } PRINTF ("\n") ; } diff --git a/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_extract.java b/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_extract.java index 651523f..3b9262d 100644 --- a/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_extract.java +++ b/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_extract.java @@ -63,10 +63,18 @@ public class Dklu_extract extends Dklu_internal int[] Fp, int[] Fi, double[] Fx, int[] P, int[] Q, double[] Rs, int[] R, KLU_common Common) { - int[] Lip, Llen, Uip, Ulen, Li2, Ui2 ; + int[] Lip, Llen, Uip, Ulen ; + /*int[]*/double[] Li2, Ui2 ; double[] LU ; double[] Lx2, Ux2, Ukk ; - int i, k, block, nblocks, n, nz, k1, k2, nk, len, kk, p ; + int i, k, block, nblocks, n, nz, k1, k2, nk, kk, p ; + int Lip_offset = 0, Llen_offset = 0, Ukk_offset = 0 ; + int Uip_offset = 0, Ulen_offset = 0 ; + int[] len = new int[1] ; + int[] Li2_offset = new int[1] ; + int[] Lx2_offset = new int[1] ; + int[] Ui2_offset = new int[1] ; + int[] Ux2_offset = new int[1] ; if (Common == null) { @@ -166,8 +174,10 @@ public class Dklu_extract extends Dklu_internal { /* non-singleton block */ LU = Numeric.LUbx [block] ; - Lip = Numeric.Lip + k1 ; - Llen = Numeric.Llen + k1 ; + Lip = Numeric.Lip ; + Lip_offset += k1 ; + Llen = Numeric.Llen ; + Llen_offset += k1 ; for (kk = 0 ; kk < nk ; kk++) { Lp [k1+kk] = nz ; @@ -175,11 +185,13 @@ public class Dklu_extract extends Dklu_internal Li [nz] = k1 + kk ; Lx [nz] = 1 ; nz++ ; - GET_POINTER (LU, Lip, Llen, Li2, Lx2, kk, len) ; - for (p = 0 ; p < len ; p++) + Li2 = Lx2 = GET_POINTER (LU, Lip, Lip_offset, + Llen, Llen_offset, + Li2_offset, Lx2_offset, kk, len) ; + for (p = 0 ; p < len[0] ; p++) { - Li [nz] = k1 + Li2 [p] ; - Lx [nz] = Lx2 [p] ; //REAL (Lx2 [p]) ; + Li [nz] = k1 + (int) Li2 [Li2_offset[0] + p] ; + Lx [nz] = Lx2 [Lx2_offset[0] + p] ; //REAL (Lx2 [p]) ; nz++ ; } } @@ -201,34 +213,37 @@ public class Dklu_extract extends Dklu_internal k1 = Symbolic.R [block] ; k2 = Symbolic.R [block+1] ; nk = k2 - k1 ; - Ukk = ((double[]) Numeric.Udiag) + k1 ; + Ukk = Numeric.Udiag ; + Ukk_offset += k1 ; if (nk == 1) { /* singleton block */ Up [k1] = nz ; Ui [nz] = k1 ; - Ux [nz] = Ukk [0] ; //REAL (Ukk [0]) ; + Ux [nz] = Ukk [Ukk_offset + 0] ; //REAL (Ukk [0]) ; nz++ ; } else { /* non-singleton block */ LU = Numeric.LUbx [block] ; - Uip = Numeric.Uip + k1 ; - Ulen = Numeric.Ulen + k1 ; + Uip = Numeric.Uip ; + Uip_offset += k1 ; + Ulen = Numeric.Ulen ; + Ulen_offset += k1 ; for (kk = 0 ; kk < nk ; kk++) { Up [k1+kk] = nz ; - GET_POINTER (LU, Uip, Ulen, Ui2, Ux2, kk, len) ; - for (p = 0 ; p < len ; p++) + Ui2 = Ux2 = GET_POINTER (LU, Uip, Uip_offset, Ulen, Ulen_offset, Ui2_offset, Ux2_offset, kk, len) ; + for (p = 0 ; p < len[0] ; p++) { - Ui [nz] = k1 + Ui2 [p] ; - Ux [nz] = Ux2 [p] ; //REAL (Ux2 [p]) ; + Ui [nz] = k1 + (int) Ui2 [Ui2_offset[0] + p] ; + Ux [nz] = Ux2 [Ux2_offset[0] + p] ; //REAL (Ux2 [p]) ; nz++ ; } /* add the diagonal entry */ Ui [nz] = k1 + kk ; - Ux [nz] = Ukk [kk] ; //REAL (Ukk [kk]) ; + Ux [nz] = Ukk [Ukk_offset+ kk] ; //REAL (Ukk [kk]) ; nz++ ; } } diff --git a/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_factor.java b/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_factor.java index a00c657..5a310ff 100644 --- a/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_factor.java +++ b/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_factor.java @@ -405,17 +405,20 @@ public class Dklu_factor extends Dklu_internal } else { -// int[] Lip, Uip, Llen, Ulen ; double[] LU ; -// Lip = Numeric.Lip + k1 ; -// Llen = Numeric.Llen + k1 ; -// LU = (double[]) Numeric.LUbx [block] ; -// PRINTF ("\n---- L block %d\n", block); -// ASSERT (klu_valid_LU (nk, TRUE, Lip, Llen, LU)) ; -// Uip = Numeric.Uip + k1 ; -// Ulen = Numeric.Ulen + k1 ; -// PRINTF ("\n---- U block %d\n", block) ; -// ASSERT (klu_valid_LU (nk, FALSE, Uip, Ulen, LU)) ; + Lip = Numeric.Lip ; + int Lip_offset = k1 ; + Llen = Numeric.Llen ; + int Llen_offset = k1 ; + LU = (double[]) Numeric.LUbx [block] ; + PRINTF ("\n---- L block %d\n", block); + ASSERT (klu_valid_LU (nk, TRUE, Lip, Lip_offset, Llen, Llen_offset, LU)) ; + Uip = Numeric.Uip ; + int Uip_offset = k1 ; + Ulen = Numeric.Ulen ; + int Ulen_offset = k1 ; + PRINTF ("\n---- U block %d\n", block) ; + ASSERT (klu_valid_LU (nk, FALSE, Uip, Uip_offset, Ulen, Ulen_offset, LU)) ; } } } @@ -435,7 +438,6 @@ public class Dklu_factor extends Dklu_internal { int n, nzoff, nblocks, maxblock, k ; int[] ok = new int [] {TRUE} ; - int[] R ; KLU_numeric Numeric ; int n1, nzoff1, s, b6, n3 ; @@ -462,7 +464,6 @@ public class Dklu_factor extends Dklu_internal nzoff = Symbolic.nzoff ; nblocks = Symbolic.nblocks ; maxblock = Symbolic.maxblock ; - R = Symbolic.R ; PRINTF ("KLU_factor: n %d nzoff %d nblocks %d maxblock %d\n", n, nzoff, nblocks, maxblock) ; @@ -579,7 +580,7 @@ public class Dklu_factor extends Dklu_internal if (Common.status < KLU_OK) { /* out of memory or inputs invalid */ -// klu_free_numeric (Numeric, Common) ; + //klu_free_numeric (Numeric, Common) ; Numeric = null ; } else if (Common.status == KLU_SINGULAR) @@ -589,8 +590,8 @@ public class Dklu_factor extends Dklu_internal /* Matrix is singular, and the Numeric object is only partially * defined because we halted early. This is the default case for * a singular matrix. */ + //klu_free_numeric (Numeric, Common) ; Numeric = null ; -// klu_free_numeric (Numeric, Common) ; } } diff --git a/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_memory.java b/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_memory.java index 8d5f4b3..f6a8c48 100644 --- a/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_memory.java +++ b/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_memory.java @@ -148,7 +148,8 @@ public class Dklu_memory extends Dklu_internal { double[] p, KLU_common Common) { double[] pnew ; - int snew, sold ; + int snew ; +// int sold ; if (Common == null) { @@ -169,7 +170,7 @@ public class Dklu_memory extends Dklu_internal { /* The object exists, and is changing to some other nonzero size. */ /* call realloc, or its equivalent */ snew = MAX (1, nnew) ; - sold = MAX (1, nold) ; +// sold = MAX (1, nold) ; try { pnew = new double[snew] ; diff --git a/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_refactor.java b/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_refactor.java index 7f3308e..3b9089c 100644 --- a/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_refactor.java +++ b/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_refactor.java @@ -57,11 +57,18 @@ public class Dklu_refactor extends Dklu_internal { double ukk, ujk, s ; double[] Offx, Lx, Ux, X, Az, Udiag ; double[] Rs ; - int[] P, Q, R, Pnum, Offp, Offi, Ui, Li, Pinv, Lip, Uip, Llen, Ulen ; + int[] Q, R, Pnum, Offp, Offi, Pinv, Lip, Uip, Llen, Ulen ; + /*int[]*/double[] Ui, Li ; double[][] LUbx ; double[] LU ; int k1, k2, nk, k, block, oldcol, pend, oldrow, n, p, newrow, scale, - nblocks, poff, i, j, up, ulen, llen, maxblock, nzoff ; + nblocks, poff, i, j, up, maxblock, nzoff ; + int[] ulen = new int[1] ; + int[] Ui_offset = new int[1] ; + int[] Ux_offset = new int[1] ; + int[] llen = new int[1] ; + int[] Li_offset = new int[1] ; + int[] Lx_offset = new int[1] ; /* ---------------------------------------------------------------------- */ /* check inputs */ @@ -90,7 +97,6 @@ public class Dklu_refactor extends Dklu_internal { /* ---------------------------------------------------------------------- */ n = Symbolic.n ; - P = Symbolic.P ; Q = Symbolic.Q ; R = Symbolic.R ; nblocks = Symbolic.nblocks ; @@ -219,10 +225,14 @@ public class Dklu_refactor extends Dklu_internal { /* construct and factor the kth block */ /* ---------------------------------------------------------- */ - Lip = Numeric.Lip + k1 ; - Llen = Numeric.Llen + k1 ; - Uip = Numeric.Uip + k1 ; - Ulen = Numeric.Ulen + k1 ; + Lip = Numeric.Lip ; + int Lip_offset = k1 ; + Llen = Numeric.Llen ; + int Llen_offset = k1 ; + Uip = Numeric.Uip ; + int Uip_offset = k1 ; + Ulen = Numeric.Ulen ; + int Ulen_offset = k1 ; LU = LUbx [block] ; for (k = 0 ; k < nk ; k++) @@ -254,19 +264,21 @@ public class Dklu_refactor extends Dklu_internal { /* compute kth column of U, and update kth column of A */ /* ------------------------------------------------------ */ - GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, ulen) ; - for (up = 0 ; up < ulen ; up++) + Ui = Ux = GET_POINTER (LU, Uip, Uip_offset, Ulen, Ulen_offset, + Ui_offset, Ux_offset, k, ulen) ; + for (up = 0 ; up < ulen[0] ; up++) { - j = Ui [up] ; + j = (int) Ui [Ui_offset[0] + up] ; ujk = X [j] ; /* X [j] = 0 ; */ CLEAR (X, j) ; - Ux [up] = ujk ; - GET_POINTER (LU, Lip, Llen, Li, Lx, j, llen) ; - for (p = 0 ; p < llen ; p++) + Ux [Ux_offset[0] + up] = ujk ; + Li = Lx = GET_POINTER (LU, Lip, Lip_offset, Llen, Llen_offset, + Li_offset, Lx_offset, j, llen) ; + for (p = 0 ; p < llen[0] ; p++) { - X [Li [p]] -= Lx [p] * ujk ; //MULT_SUB (X [Li [p]], Lx [p], ujk) ; + X [(int) Li [Li_offset[0] + p]] -= Lx [Lx_offset[0] + p] * ujk ; } } /* get the diagonal entry of U */ @@ -291,12 +303,13 @@ public class Dklu_refactor extends Dklu_internal { } Udiag [k+k1] = ukk ; /* gather and divide by pivot to get kth column of L */ - GET_POINTER (LU, Lip, Llen, Li, Lx, k, llen) ; - for (p = 0 ; p < llen ; p++) + Li = Lx = GET_POINTER (LU, Lip, Lip_offset, Llen, Llen_offset, + Li_offset, Lx_offset, k, llen) ; + for (p = 0 ; p < llen[0] ; p++) { - i = Li [p] ; - Lx [p] = X [i] / ukk ; + i = (int) Li [Li_offset[0] + p] ; //DIV (Lx [p], X [i], ukk) ; + Lx [Lx_offset[0] + p] = X [i] / ukk ; CLEAR (X, i) ; } @@ -361,10 +374,14 @@ public class Dklu_refactor extends Dklu_internal { /* construct and factor the kth block */ /* ---------------------------------------------------------- */ - Lip = Numeric.Lip + k1 ; - Llen = Numeric.Llen + k1 ; - Uip = Numeric.Uip + k1 ; - Ulen = Numeric.Ulen + k1 ; + Lip = Numeric.Lip ; + int Lip_offset = k1 ; + Llen = Numeric.Llen ; + int Llen_offset = k1 ; + Uip = Numeric.Uip ; + int Uip_offset = k1 ; + Ulen = Numeric.Ulen ; + int Ulen_offset = k1 ; LU = LUbx [block] ; for (k = 0 ; k < nk ; k++) @@ -399,19 +416,22 @@ public class Dklu_refactor extends Dklu_internal { /* compute kth column of U, and update kth column of A */ /* ------------------------------------------------------ */ - GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, ulen) ; - for (up = 0 ; up < ulen ; up++) + Ui = Ux = GET_POINTER (LU, Uip, Uip_offset, Ulen, Ulen_offset, + Ui_offset, Ux_offset, k, ulen) ; +// GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, ulen) ; + for (up = 0 ; up < ulen[0] ; up++) { - j = Ui [up] ; + j = (int) Ui [Ui_offset[0] + up] ; ujk = X [j] ; /* X [j] = 0 ; */ CLEAR (X, j) ; - Ux [up] = ujk ; - GET_POINTER (LU, Lip, Llen, Li, Lx, j, llen) ; - for (p = 0 ; p < llen ; p++) + Ux [Ux_offset[0] + up] = ujk ; + Li = Lx = GET_POINTER (LU, Lip, Lip_offset, Llen, Llen_offset, + Li_offset, Lx_offset, j, llen) ; + for (p = 0 ; p < llen[0] ; p++) { - X [Li [p]] -= Lx [p] * ujk ; //MULT_SUB (X [Li [p]], Lx [p], ujk) ; + X [(int) Li [Li_offset[0] + p]] -= Lx [Lx_offset[0] + p] * ujk ; } } /* get the diagonal entry of U */ @@ -436,12 +456,14 @@ public class Dklu_refactor extends Dklu_internal { } Udiag [k+k1] = ukk ; /* gather and divide by pivot to get kth column of L */ - GET_POINTER (LU, Lip, Llen, Li, Lx, k, llen) ; - for (p = 0 ; p < llen ; p++) + Li = Lx = GET_POINTER (LU, Lip, Lip_offset, Llen, Llen_offset, + Li_offset, Lx_offset, k, llen) ; +// GET_POINTER (LU, Lip, Llen, Li, Lx, k, llen) ; + for (p = 0 ; p < llen[0] ; p++) { - i = Li [p] ; - Lx [p] = X [i] / ukk ; + i = (int) Li [Li_offset[0] + p] ; //DIV (Lx [p], X [i], ukk) ; + Lx [Lx_offset[0] + p] = X [i] / ukk ; CLEAR (X, i) ; } } @@ -492,15 +514,19 @@ public class Dklu_refactor extends Dklu_internal { } else { - Lip = Numeric.Lip + k1 ; - Llen = Numeric.Llen + k1 ; + Lip = Numeric.Lip ; + int Lip_offset = k1 ; + Llen = Numeric.Llen ; + int Llen_offset = k1 ; LU = (double[]) Numeric.LUbx [block] ; PRINTF ("\n---- L block %d\n", block) ; - ASSERT (klu_valid_LU (nk, TRUE, Lip, Llen, LU)) ; - Uip = Numeric.Uip + k1 ; - Ulen = Numeric.Ulen + k1 ; + ASSERT (klu_valid_LU (nk, TRUE, Lip, Lip_offset, Llen, Llen_offset, LU)) ; + Uip = Numeric.Uip ; + int Uip_offset = k1 ; + Ulen = Numeric.Ulen ; + int Ulen_offset = k1 ; PRINTF ("\n---- U block %d\n", block) ; - ASSERT (klu_valid_LU (nk, FALSE, Uip, Ulen, LU)) ; + ASSERT (klu_valid_LU (nk, FALSE, Uip, Uip_offset, Ulen, Ulen_offset, LU)) ; } } } diff --git a/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_solve.java b/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_solve.java index 9efc315..97ffa26 100644 --- a/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_solve.java +++ b/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_solve.java @@ -51,7 +51,7 @@ public class Dklu_solve extends Dklu_internal { * @return */ public static int klu_solve(KLU_symbolic Symbolic, KLU_numeric Numeric, - int d, int nrhs, double[] B, KLU_common Common) + int d, int nrhs, double[] B, int B_offset, KLU_common Common) { double offik, s ; double[] x = new double[4] ; @@ -137,7 +137,7 @@ public class Dklu_solve extends Dklu_internal { for (k = 0 ; k < n ; k++) { - X [k] = Bz [Pnum [k]] ; + X [k] = Bz [B_offset + Pnum [k]] ; } break ; @@ -146,8 +146,8 @@ public class Dklu_solve extends Dklu_internal { for (k = 0 ; k < n ; k++) { i = Pnum [k] ; - X [2*k ] = Bz [i ] ; - X [2*k + 1] = Bz [i + d ] ; + X [2*k ] = Bz [B_offset + i] ; + X [2*k + 1] = Bz [B_offset + i + d] ; } break ; @@ -156,9 +156,9 @@ public class Dklu_solve extends Dklu_internal { for (k = 0 ; k < n ; k++) { i = Pnum [k] ; - X [3*k ] = Bz [i ] ; - X [3*k + 1] = Bz [i + d ] ; - X [3*k + 2] = Bz [i + d*2] ; + X [3*k ] = Bz [B_offset + i ] ; + X [3*k + 1] = Bz [B_offset + i + d ] ; + X [3*k + 2] = Bz [B_offset + i + d*2] ; } break ; @@ -167,10 +167,10 @@ public class Dklu_solve extends Dklu_internal { for (k = 0 ; k < n ; k++) { i = Pnum [k] ; - X [4*k ] = Bz [i ] ; - X [4*k + 1] = Bz [i + d ] ; - X [4*k + 2] = Bz [i + d*2] ; - X [4*k + 3] = Bz [i + d*3] ; + X [4*k ] = Bz [B_offset + i ] ; + X [4*k + 1] = Bz [B_offset + i + d ] ; + X [4*k + 2] = Bz [B_offset + i + d*2] ; + X [4*k + 3] = Bz [B_offset + i + d*3] ; } break ; } @@ -186,7 +186,7 @@ public class Dklu_solve extends Dklu_internal { for (k = 0 ; k < n ; k++) { - X [k] = Bz [Pnum [k]] / Rs [k] ; + X [k] = Bz [B_offset + Pnum [k]] / Rs [k] ; //SCALE_DIV_ASSIGN (X [k], Bz [Pnum [k]], Rs [k]) ; } break ; @@ -197,9 +197,9 @@ public class Dklu_solve extends Dklu_internal { { i = Pnum [k] ; rs = Rs [k] ; - X [2*k] = Bz [i] / rs ; + X [2*k] = Bz [B_offset + i] / rs ; //SCALE_DIV_ASSIGN (X [2*k], Bz [i], rs) ; - X [2*k + 1] = Bz [i + d] / rs ; + X [2*k + 1] = Bz [B_offset + i + d] / rs ; //SCALE_DIV_ASSIGN (X [2*k + 1], Bz [i + d], rs) ; } break ; @@ -210,11 +210,11 @@ public class Dklu_solve extends Dklu_internal { { i = Pnum [k] ; rs = Rs [k] ; - X [3*k] = Bz [i] / rs ; + X [3*k] = Bz [B_offset + i] / rs ; //SCALE_DIV_ASSIGN (X [3*k], Bz [i], rs) ; - X [3*k + 1] = Bz [i + d] / rs ; + X [3*k + 1] = Bz [B_offset + i + d] / rs ; //SCALE_DIV_ASSIGN (X [3*k + 1], Bz [i + d], rs) ; - X [3*k + 2] = Bz [i + d*2] / rs ; + X [3*k + 2] = Bz [B_offset + i + d*2] / rs ; //SCALE_DIV_ASSIGN (X [3*k + 2], Bz [i + d*2], rs) ; } break ; @@ -225,13 +225,13 @@ public class Dklu_solve extends Dklu_internal { { i = Pnum [k] ; rs = Rs [k] ; - X [4*k] = Bz [i] / rs ; + X [4*k] = Bz [B_offset + i] / rs ; //SCALE_DIV_ASSIGN (X [4*k], Bz [i], rs) ; - X [4*k + 1] = Bz [i + d] / rs ; + X [4*k + 1] = Bz [B_offset + i + d] / rs ; //SCALE_DIV_ASSIGN (X [4*k + 1], Bz [i + d], rs) ; - X [4*k + 2] = Bz [i + d*2] / rs ; + X [4*k + 2] = Bz [B_offset + i + d*2] / rs ; //SCALE_DIV_ASSIGN (X [4*k + 2], Bz [i + d*2], rs) ; - X [4*k + 3] = Bz [i + d*3] / rs ; + X [4*k + 3] = Bz [B_offset + i + d*3] / rs ; //SCALE_DIV_ASSIGN (X [4*k + 3], Bz [i + d*3], rs) ; } break ; @@ -406,7 +406,7 @@ public class Dklu_solve extends Dklu_internal { for (k = 0 ; k < n ; k++) { - Bz [Q [k]] = X [k] ; + Bz [B_offset + Q [k]] = X [k] ; } break ; @@ -415,8 +415,8 @@ public class Dklu_solve extends Dklu_internal { for (k = 0 ; k < n ; k++) { i = Q [k] ; - Bz [i ] = X [2*k ] ; - Bz [i + d ] = X [2*k + 1] ; + Bz [B_offset + i ] = X [2*k ] ; + Bz [B_offset + i + d ] = X [2*k + 1] ; } break ; @@ -425,9 +425,9 @@ public class Dklu_solve extends Dklu_internal { for (k = 0 ; k < n ; k++) { i = Q [k] ; - Bz [i ] = X [3*k ] ; - Bz [i + d ] = X [3*k + 1] ; - Bz [i + d*2] = X [3*k + 2] ; + Bz [B_offset + i ] = X [3*k ] ; + Bz [B_offset + i + d ] = X [3*k + 1] ; + Bz [B_offset + i + d*2] = X [3*k + 2] ; } break ; @@ -436,10 +436,10 @@ public class Dklu_solve extends Dklu_internal { for (k = 0 ; k < n ; k++) { i = Q [k] ; - Bz [i ] = X [4*k ] ; - Bz [i + d ] = X [4*k + 1] ; - Bz [i + d*2] = X [4*k + 2] ; - Bz [i + d*3] = X [4*k + 3] ; + Bz [B_offset + i ] = X [4*k ] ; + Bz [B_offset + i + d ] = X [4*k + 1] ; + Bz [B_offset + i + d*2] = X [4*k + 2] ; + Bz [B_offset + i + d*3] = X [4*k + 3] ; } break ; } @@ -448,7 +448,7 @@ public class Dklu_solve extends Dklu_internal { /* go to the next chunk of B */ /* ------------------------------------------------------------------ */ - Bz += d*4 ; + B_offset += d*4 ; } return (TRUE) ; } diff --git a/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_tsolve.java b/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_tsolve.java index 37edd36..63ef36b 100644 --- a/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_tsolve.java +++ b/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_tsolve.java @@ -50,8 +50,8 @@ public class Dklu_tsolve extends Dklu_internal { * @return */ public static int klu_tsolve(KLU_symbolic Symbolic, - KLU_numeric Numeric, int d, int nrhs, double[] B, - KLU_common Common) + KLU_numeric Numeric, int d, int nrhs, + double[] B, int B_offset, KLU_common Common) { double[] x = new double[4] ; double offik, s ; @@ -133,7 +133,7 @@ public class Dklu_tsolve extends Dklu_internal { for (k = 0 ; k < n ; k++) { - X [k] = Bz [Q [k]] ; + X [k] = Bz [B_offset + Q [k]] ; } break ; @@ -142,8 +142,8 @@ public class Dklu_tsolve extends Dklu_internal { for (k = 0 ; k < n ; k++) { i = Q [k] ; - X [2*k ] = Bz [i ] ; - X [2*k + 1] = Bz [i + d ] ; + X [2*k ] = Bz [B_offset + i ] ; + X [2*k + 1] = Bz [B_offset + i + d ] ; } break ; @@ -152,9 +152,9 @@ public class Dklu_tsolve extends Dklu_internal { for (k = 0 ; k < n ; k++) { i = Q [k] ; - X [3*k ] = Bz [i ] ; - X [3*k + 1] = Bz [i + d ] ; - X [3*k + 2] = Bz [i + d*2] ; + X [3*k ] = Bz [B_offset + i ] ; + X [3*k + 1] = Bz [B_offset + i + d ] ; + X [3*k + 2] = Bz [B_offset + i + d*2] ; } break ; @@ -163,10 +163,10 @@ public class Dklu_tsolve extends Dklu_internal { for (k = 0 ; k < n ; k++) { i = Q [k] ; - X [4*k ] = Bz [i ] ; - X [4*k + 1] = Bz [i + d ] ; - X [4*k + 2] = Bz [i + d*2] ; - X [4*k + 3] = Bz [i + d*3] ; + X [4*k ] = Bz [B_offset + i ] ; + X [4*k + 1] = Bz [B_offset + i + d ] ; + X [4*k + 2] = Bz [B_offset + i + d*2] ; + X [4*k + 3] = Bz [B_offset + i + d*3] ; } break ; @@ -365,7 +365,7 @@ public class Dklu_tsolve extends Dklu_internal { for (k = 0 ; k < n ; k++) { - Bz [Pnum [k]] = X [k] ; + Bz [B_offset + Pnum [k]] = X [k] ; } break ; @@ -374,8 +374,8 @@ public class Dklu_tsolve extends Dklu_internal { for (k = 0 ; k < n ; k++) { i = Pnum [k] ; - Bz [i ] = X [2*k ] ; - Bz [i + d ] = X [2*k + 1] ; + Bz [B_offset + i ] = X [2*k ] ; + Bz [B_offset + i + d ] = X [2*k + 1] ; } break ; @@ -384,9 +384,9 @@ public class Dklu_tsolve extends Dklu_internal { for (k = 0 ; k < n ; k++) { i = Pnum [k] ; - Bz [i ] = X [3*k ] ; - Bz [i + d ] = X [3*k + 1] ; - Bz [i + d*2] = X [3*k + 2] ; + Bz [B_offset + i ] = X [3*k ] ; + Bz [B_offset + i + d ] = X [3*k + 1] ; + Bz [B_offset + i + d*2] = X [3*k + 2] ; } break ; @@ -395,10 +395,10 @@ public class Dklu_tsolve extends Dklu_internal { for (k = 0 ; k < n ; k++) { i = Pnum [k] ; - Bz [i ] = X [4*k ] ; - Bz [i + d ] = X [4*k + 1] ; - Bz [i + d*2] = X [4*k + 2] ; - Bz [i + d*3] = X [4*k + 3] ; + Bz [B_offset + i ] = X [4*k ] ; + Bz [B_offset + i + d ] = X [4*k + 1] ; + Bz [B_offset + i + d*2] = X [4*k + 2] ; + Bz [B_offset + i + d*3] = X [4*k + 3] ; } break ; } @@ -414,7 +414,7 @@ public class Dklu_tsolve extends Dklu_internal { for (k = 0 ; k < n ; k++) { - Bz [Pnum [k]] = X [k] / Rs [k] ; + Bz [B_offset + Pnum [k]] = X [k] / Rs [k] ; //SCALE_DIV_ASSIGN (Bz [Pnum [k]], X [k], Rs [k]) ; } break ; @@ -425,9 +425,9 @@ public class Dklu_tsolve extends Dklu_internal { { i = Pnum [k] ; rs = Rs [k] ; - Bz [i] = X [2*k] / rs ; + Bz [B_offset + i] = X [2*k] / rs ; //SCALE_DIV_ASSIGN (Bz [i], X [2*k], rs) ; - Bz [i + d] = X [2*k + 1] / rs ; + Bz [B_offset + i + d] = X [2*k + 1] / rs ; //SCALE_DIV_ASSIGN (Bz [i + d], X [2*k + 1], rs) ; } break ; @@ -438,11 +438,11 @@ public class Dklu_tsolve extends Dklu_internal { { i = Pnum [k] ; rs = Rs [k] ; - Bz [i] = X [3*k] / rs ; + Bz [B_offset + i] = X [3*k] / rs ; //SCALE_DIV_ASSIGN (Bz [i], X [3*k], rs) ; - Bz [i + d] = X [3*k + 1] / rs ; + Bz [B_offset + i + d] = X [3*k + 1] / rs ; //SCALE_DIV_ASSIGN (Bz [i + d], X [3*k + 1], rs) ; - Bz [i + d*2] = X [3*k + 2] / rs ; + Bz [B_offset + i + d*2] = X [3*k + 2] / rs ; //SCALE_DIV_ASSIGN (Bz [i + d*2], X [3*k + 2], rs) ; } break ; @@ -453,13 +453,13 @@ public class Dklu_tsolve extends Dklu_internal { { i = Pnum [k] ; rs = Rs [k] ; - Bz [i] = X [4*k] / rs ; + Bz [B_offset + i] = X [4*k] / rs ; //SCALE_DIV_ASSIGN (Bz [i], X [4*k], rs) ; - Bz [i + d] = X [4*k + 1] / rs ; + Bz [B_offset + i + d] = X [4*k + 1] / rs ; //SCALE_DIV_ASSIGN (Bz [i + d], X [4*k + 1], rs) ; - Bz [i + d*2] = X [4*k + 2] / rs ; + Bz [B_offset + i + d*2] = X [4*k + 2] / rs ; //SCALE_DIV_ASSIGN (Bz [i + d*2], X [4*k + 2], rs) ; - Bz [i + d*3] = X [4*k + 3] / rs ; + Bz [B_offset + i + d*3] = X [4*k + 3] / rs ; //SCALE_DIV_ASSIGN (Bz [i + d*3], X [4*k + 3], rs) ; } break ; @@ -470,7 +470,7 @@ public class Dklu_tsolve extends Dklu_internal { /* go to the next chunk of B */ /* ------------------------------------------------------------------ */ - Bz += d*4 ; + B_offset += d*4 ; } return (TRUE) ; } diff --git a/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_version.java b/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_version.java index 7234f30..d6f41a4 100644 --- a/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_version.java +++ b/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_version.java @@ -60,12 +60,14 @@ public abstract class Dklu_version { // { // return Math.ceil(BYTES (type, (double) n) / sizeof (double)) ; // } -// -// protected static void GET_I_POINTER(double[] LU, int[] Xip, int[] Xi, int k) -// { -// Xi = (Int[]) (LU + Xip [k]) ; -// } -// + + protected static double[] GET_I_POINTER(double[] LU, int[] Xip, + int Xip_offset, int[] Xi_offset, int k) + { + Xi_offset[0] = Xip [Xip_offset + k] ; + return LU ; + } + // protected static void GET_X_POINTER(double[] LU, int[] Xip, int Xlen, // double[] Xx, int k) // { -- 2.11.4.GIT