From 03579c6350c886b4e9a4f4748e437cff565c6316 Mon Sep 17 00:00:00 2001 From: Richard Lincoln Date: Sun, 31 Jul 2011 21:20:34 +0100 Subject: [PATCH] Re-uppercasing macros. --- .../java/edu/ufl/cise/klu/common/KLU_common.java | 58 ++--- .../java/edu/ufl/cise/klu/common/KLU_numeric.java | 48 ++-- .../java/edu/ufl/cise/klu/common/KLU_symbolic.java | 23 +- src/main/java/edu/ufl/cise/klu/tdouble/Dklu.java | 184 +++++++-------- .../edu/ufl/cise/klu/tdouble/Dklu_internal.java | 32 +-- .../java/edu/ufl/cise/klu/tdouble/Dklu_kernel.java | 256 +++++++++++---------- .../java/edu/ufl/cise/klu/tdouble/Dklu_memory.java | 44 ++-- .../edu/ufl/cise/klu/tdouble/Dklu_refactor.java | 116 +++++----- .../java/edu/ufl/cise/klu/tdouble/Dklu_scale.java | 40 ++-- .../java/edu/ufl/cise/klu/tdouble/Dklu_solve.java | 86 +++---- .../java/edu/ufl/cise/klu/tdouble/Dklu_sort.java | 26 ++- .../java/edu/ufl/cise/klu/tdouble/Dklu_tsolve.java | 90 ++++---- .../edu/ufl/cise/klu/tdouble/Dklu_version.java | 128 +++++------ 13 files changed, 573 insertions(+), 558 deletions(-) diff --git a/src/main/java/edu/ufl/cise/klu/common/KLU_common.java b/src/main/java/edu/ufl/cise/klu/common/KLU_common.java index d908600..6c98115 100644 --- a/src/main/java/edu/ufl/cise/klu/common/KLU_common.java +++ b/src/main/java/edu/ufl/cise/klu/common/KLU_common.java @@ -30,29 +30,29 @@ package edu.ufl.cise.klu.common; public class KLU_common { - static final int KLU_OK = 0; + public static final int KLU_OK = 0; /* status > 0 is a warning, not an error */ - static final int KLU_SINGULAR = 1; - static final int KLU_OUT_OF_MEMORY = -2; - static final int KLU_INVALID = -3; + public static final int KLU_SINGULAR = 1; + public static final int KLU_OUT_OF_MEMORY = -2; + public static final int KLU_INVALID = -3; /* integer overflow has occured */ - static final int KLU_TOO_LARGE = -4; + public static final int KLU_TOO_LARGE = -4; /* ---------------------------------------------------------------------- */ /* parameters */ /* ---------------------------------------------------------------------- */ - double tol; /* pivot tolerance for diagonal preference */ - double memgrow; /* realloc memory growth size for LU factors */ - double initmem_amd; /* init. memory size with AMD: c*nnz(L) + n */ - double initmem; /* init. memory size: c*nnz(A) + n */ - double maxwork; /* maxwork for BTF, <= 0 if no limit */ + public double tol; /* pivot tolerance for diagonal preference */ + public double memgrow; /* realloc memory growth size for LU factors */ + public double initmem_amd; /* init. memory size with AMD: c*nnz(L) + n */ + public double initmem; /* init. memory size: c*nnz(A) + n */ + public double maxwork; /* maxwork for BTF, <= 0 if no limit */ - int btf; /* use BTF pre-ordering, or not */ - int ordering; /* 0: AMD, 1: COLAMD, 2: user P and Q, - * 3: user function */ - int scale; /* row scaling: -1: none (and no error check), - * 0: none, 1: sum, 2: max */ + public int btf; /* use BTF pre-ordering, or not */ + public int ordering; /* 0: AMD, 1: COLAMD, 2: user P and Q, + * 3: user function */ + public int scale; /* row scaling: -1: none (and no error check), + * 0: none, 1: sum, 2: max */ /* memory management routines */ // void *(*malloc_memory) (size_t); /* pointer to malloc */ @@ -68,7 +68,7 @@ public class KLU_common * information). */ Object user_data; - int halt_if_singular; /* how to handle a singular matrix: + public int halt_if_singular; /* how to handle a singular matrix: * FALSE: keep going. Return a Numeric object with a zero U(k,k). A * divide-by-zero may occur when computing L(:,k). The Numeric object * can be passed to klu_solve (a divide-by-zero will occur). It can @@ -81,34 +81,34 @@ public class KLU_common /* statistics */ /* ---------------------------------------------------------------------- */ - int status; /* KLU_OK if OK, < 0 if error */ - int nrealloc; /* # of reallocations of L and U */ + public int status; /* KLU_OK if OK, < 0 if error */ + public int nrealloc; /* # of reallocations of L and U */ - int structural_rank; /* 0 to n-1 if the matrix is structurally rank + public int structural_rank; /* 0 to n-1 if the matrix is structurally rank * deficient (as determined by maxtrans). -1 if not computed. n if the * matrix has full structural rank. This is computed by klu_analyze * if a BTF preordering is requested. */ - int numerical_rank; /* First k for which a zero U(k,k) was found, + public int numerical_rank; /* First k for which a zero U(k,k) was found, * if the matrix was singular (in the range 0 to n-1). n if the matrix * has full rank. This is not a true rank-estimation. It just reports * where the first zero pivot was found. -1 if not computed. * Computed by klu_factor and klu_refactor. */ - int singular_col; /* n if the matrix is not singular. If in the + public int singular_col; /* n if the matrix is not singular. If in the * range 0 to n-1, this is the column index of the original matrix A that * corresponds to the column of U that contains a zero diagonal entry. * -1 if not computed. Computed by klu_factor and klu_refactor. */ - int noffdiag; /* # of off-diagonal pivots, -1 if not computed */ + public int noffdiag; /* # of off-diagonal pivots, -1 if not computed */ - double flops; /* actual factorization flop count, from klu_flops */ - double rcond; /* crude reciprocal condition est., from klu_rcond */ - double condest; /* accurate condition est., from klu_condest */ - double rgrowth; /* reciprocal pivot rgrowth, from klu_rgrowth */ - double work; /* actual work done in BTF, in klu_analyze */ + public double flops; /* actual factorization flop count, from klu_flops */ + public double rcond; /* crude reciprocal condition est., from klu_rcond */ + public double condest; /* accurate condition est., from klu_condest */ + public double rgrowth; /* reciprocal pivot rgrowth, from klu_rgrowth */ + public double work; /* actual work done in BTF, in klu_analyze */ - size_t memusage; /* current memory usage, in bytes */ - size_t mempeak; /* peak memory usage, in bytes */ + public size_t memusage; /* current memory usage, in bytes */ + public size_t mempeak; /* peak memory usage, in bytes */ } diff --git a/src/main/java/edu/ufl/cise/klu/common/KLU_numeric.java b/src/main/java/edu/ufl/cise/klu/common/KLU_numeric.java index 24076b6..5087dc4 100644 --- a/src/main/java/edu/ufl/cise/klu/common/KLU_numeric.java +++ b/src/main/java/edu/ufl/cise/klu/common/KLU_numeric.java @@ -33,37 +33,37 @@ public class KLU_numeric /* LU factors of each block, the pivot row permutation, and the * entries in the off-diagonal blocks */ - int n; /* A is n-by-n */ - int nblocks; /* number of diagonal blocks */ - int lnz; /* actual nz in L, including diagonal */ - int unz; /* actual nz in U, including diagonal */ - int max_lnz_block; /* max actual nz in L in any one block, incl. diag */ - int max_unz_block; /* max actual nz in U in any one block, incl. diag */ - int[] Pnum; /* size n. final pivot permutation */ - int[] Pinv; /* size n. inverse of final pivot permutation */ + public int n; /* A is n-by-n */ + public int nblocks; /* number of diagonal blocks */ + public int lnz; /* actual nz in L, including diagonal */ + public int unz; /* actual nz in U, including diagonal */ + public int max_lnz_block; /* max actual nz in L in any one block, incl. diag */ + public int max_unz_block; /* max actual nz in U in any one block, incl. diag */ + public int[] Pnum; /* size n. final pivot permutation */ + public int[] Pinv; /* size n. inverse of final pivot permutation */ /* LU factors of each block */ - int[] Lip; /* size n. pointers into LUbx[block] for L */ - int[] Uip; /* size n. pointers into LUbx[block] for U */ - int[] Llen; /* size n. Llen [k] = # of entries in kth column of L */ - int[] Ulen; /* size n. Ulen [k] = # of entries in kth column of U */ - void[][] LUbx; /* L and U indices and entries (excl. diagonal of U) */ - size_t[] LUsize; /* size of each LUbx [block], in sizeof (Unit) */ - void[] Udiag; /* diagonal of U */ + public int[] Lip; /* size n. pointers into LUbx[block] for L */ + public int[] Uip; /* size n. pointers into LUbx[block] for U */ + public int[] Llen; /* size n. Llen [k] = # of entries in kth column of L */ + public int[] Ulen; /* size n. Ulen [k] = # of entries in kth column of U */ + public void[][] LUbx; /* L and U indices and entries (excl. diagonal of U) */ + public size_t[] LUsize; /* size of each LUbx [block], in sizeof (Unit) */ + public void[] Udiag; /* diagonal of U */ /* scale factors; can be NULL if no scaling */ - double[] Rs; /* size n. Rs [i] is scale factor for row i */ + public double[] Rs; /* size n. Rs [i] is scale factor for row i */ /* permanent workspace for factorization and solve */ - size_t[] worksize; /* size (in bytes) of Work */ - void[] Work; /* workspace */ - void[] Xwork; /* alias into Numeric->Work */ - int[] Iwork; /* alias into Numeric->Work */ + public size_t[] worksize; /* size (in bytes) of Work */ + public void[] Work; /* workspace */ + public void[] Xwork; /* alias into Numeric->Work */ + public int[] Iwork; /* alias into Numeric->Work */ /* off-diagonal entries in a conventional compressed-column sparse matrix */ - int[] Offp; /* size n+1, column pointers */ - int[] Offi; /* size nzoff, row indices */ - void[] Offx; /* size nzoff, numerical values */ - int nzoff; + public int[] Offp; /* size n+1, column pointers */ + public int[] Offi; /* size nzoff, row indices */ + public void[] Offx; /* size nzoff, numerical values */ + public int nzoff; } diff --git a/src/main/java/edu/ufl/cise/klu/common/KLU_symbolic.java b/src/main/java/edu/ufl/cise/klu/common/KLU_symbolic.java index 6633ee8..1365aed 100644 --- a/src/main/java/edu/ufl/cise/klu/common/KLU_symbolic.java +++ b/src/main/java/edu/ufl/cise/klu/common/KLU_symbolic.java @@ -36,28 +36,23 @@ public class KLU_symbolic */ /* only computed if the AMD ordering is chosen: */ - double symmetry; /* symmetry of largest block */ - double est_flops; /* est. factorization flop count */ - double lnz, unz; /* estimated nz in L and U, including diagonals */ - double[] Lnz; /* size n, but only Lnz [0..nblocks-1] is used */ + public double symmetry; /* symmetry of largest block */ + public double est_flops; /* est. factorization flop count */ + public double lnz, unz; /* estimated nz in L and U, including diagonals */ + public double[] Lnz; /* size n, but only Lnz [0..nblocks-1] is used */ /* computed for all orderings: */ - int - n, /* input matrix A is n-by-n */ - nz, /* # entries in input matrix */ - nzoff, /* nz in off-diagonal blocks */ - nblocks, /* number of blocks */ - maxblock, /* size of largest block */ - ordering, /* ordering used (AMD, COLAMD, or GIVEN) */ - do_btf; /* whether or not BTF preordering was requested */ + public int + n; /* whether or not BTF preordering was requested */ + public int nz, nzoff, nblocks, maxblock, ordering, do_btf; - int[] + public int[] P, /* size n */ Q, /* size n */ R; /* size n+1, but only R [0..nblocks] is used */ /* only computed if BTF preordering requested */ - int structural_rank; /* 0 to n-1 if the matrix is structurally rank + public int structural_rank; /* 0 to n-1 if the matrix is structurally rank * deficient. -1 if not computed. n if the matrix has * full structural rank */ diff --git a/src/main/java/edu/ufl/cise/klu/tdouble/Dklu.java b/src/main/java/edu/ufl/cise/klu/tdouble/Dklu.java index 78e7674..8900270 100644 --- a/src/main/java/edu/ufl/cise/klu/tdouble/Dklu.java +++ b/src/main/java/edu/ufl/cise/klu/tdouble/Dklu.java @@ -128,13 +128,13 @@ public class Dklu extends Dklu_internal { /* get control parameters, or use defaults */ /* ------------------------------------------------------------------ */ - n = max(1, n); + n = MAX (1, n); anz = Ap[n+k1] - Ap[k1]; if (Lsize <= 0) { Lsize = -Lsize; - Lsize = max(Lsize, 1.0); + Lsize = MAX (Lsize, 1.0); lsize = Lsize * anz + n; } else @@ -144,15 +144,15 @@ public class Dklu extends Dklu_internal { usize = lsize; - lsize = max(n+1, lsize); - usize = max(n+1, usize); + lsize = MAX (n+1, lsize); + usize = MAX (n+1, usize); maxlnz = (((double) n) * ((double) n) + ((double) n)) / 2.; - maxlnz = min(maxlnz, ((double) INT_MAX)); - lsize = min(maxlnz, lsize); - usize = min(maxlnz, usize); + maxlnz = MIN (maxlnz, ((double) INT_MAX)); + lsize = MIN (maxlnz, lsize); + usize = MIN (maxlnz, usize); - printf("Welcome to klu: n %d anz %d k1 %d lsize %d usize %d maxlnz %g\n", + PRINTF ("Welcome to klu: n %d anz %d k1 %d lsize %d usize %d maxlnz %g\n", n, anz, k1, lsize, usize, maxlnz); /* ------------------------------------------------------------------ */ @@ -170,15 +170,15 @@ public class Dklu extends Dklu_internal { Lpend = (int[]) W; W += n; Ap_pos = (int[]) W; W += n; - dunits= dunits(Int, lsize) + dunits(Entry, lsize) + - dunits(Int, usize) + dunits(Entry, usize); - lusize = (size_t) dunits; - ok = !int_overflow(dunits); + dunits= DUNITS (Int, lsize) + DUNITS (Entry, lsize) + + DUNITS (Int, usize) + DUNITS (Entry, usize); + lusize = (size_t) DUNITS ; + ok = !INT_OVERFLOW (dunits); LU = ok != 0 ? Dklu_malloc.klu_malloc(lusize, sizeof(Unit), Common) : null; if (LU == null) { /* out of memory, or problem too large */ - Common.status = KLU_OUT_OF_MEMORY; + Common.status = KLU_common.KLU_OUT_OF_MEMORY; lusize = 0; return (lusize); } @@ -197,13 +197,13 @@ public class Dklu extends Dklu_internal { /* return LU factors, or return nothing if an error occurred */ /* ------------------------------------------------------------------ */ - if (Common.status < KLU_OK) + if (Common.status < KLU_common.KLU_OK) { - LU = Dklu_free.klu_free(LU, lusize, sizeof(Unit), Common); + LU = Dklu_free.klu_free(LU, lusize, sizeof (Unit), Common); lusize = 0; } p_LU = LU; - printf(" in klu noffdiag %d\n", Common.noffdiag); + PRINTF (" in klu noffdiag %d\n", Common.noffdiag); return lusize; } @@ -237,12 +237,12 @@ public class Dklu extends Dklu_internal { for (k = 0; k < n; k++) { x[0] = X[k]; - get_pointer(LU, Lip, Llen, Li, Lx, k, len); + GET_POINTER (LU, Lip, Llen, Li, Lx, k, len); /* unit diagonal of L is not stored*/ for (p = 0; p < len; p++) { /* X[Li[p]] -= Lx[p] * x[0]; */ - mult_sub(X[Li[p]], Lx[p], x[0]); + MULT_SUB (X[Li[p]], Lx[p], x[0]); } } break; @@ -253,13 +253,13 @@ public class Dklu extends Dklu_internal { { x[0] = X[2*k ]; x[1] = X[2*k + 1]; - get_pointer(LU, Lip, Llen, Li, Lx, k, len); + GET_POINTER (LU, Lip, Llen, Li, Lx, k, len); for (p = 0; p < len; p++) { i = Li[p]; lik = Lx[p]; - mult_sub(X[2*i], lik, x[0]); - mult_sub(X[2*i + 1], lik, x[1]); + MULT_SUB (X[2*i], lik, x[0]); + MULT_SUB (X[2*i + 1], lik, x[1]); } } break; @@ -271,14 +271,14 @@ public class Dklu extends Dklu_internal { x[0] = X[3*k ]; x[1] = X[3*k + 1]; x[2] = X[3*k + 2]; - get_pointer(LU, Lip, Llen, Li, Lx, k, len); + GET_POINTER (LU, Lip, Llen, Li, Lx, k, len); for (p = 0; p < len; p++) { i = Li[p]; lik = Lx[p]; - mult_sub(X[3*i], lik, x[0]); - mult_sub(X[3*i + 1], lik, x[1]); - mult_sub(X[3*i + 2], lik, x[2]); + MULT_SUB (X[3*i], lik, x[0]); + MULT_SUB (X[3*i + 1], lik, x[1]); + MULT_SUB (X[3*i + 2], lik, x[2]); } } break; @@ -291,15 +291,15 @@ public class Dklu extends Dklu_internal { x[1] = X[4*k + 1]; x[2] = X[4*k + 2]; x[3] = X[4*k + 3]; - get_pointer(LU, Lip, Llen, Li, Lx, k, len); + GET_POINTER (LU, Lip, Llen, Li, Lx, k, len); for (p = 0; p < len; p++) { i = Li[p]; lik = Lx[p]; - mult_sub(X[4*i], lik, x[0]); - mult_sub(X[4*i + 1], lik, x[1]); - mult_sub(X[4*i + 2], lik, x[2]); - mult_sub(X[4*i + 3], lik, x[3]); + MULT_SUB (X[4*i], lik, x[0]); + MULT_SUB (X[4*i + 1], lik, x[1]); + MULT_SUB (X[4*i + 2], lik, x[2]); + MULT_SUB (X[4*i + 3], lik, x[3]); } } break; @@ -338,14 +338,14 @@ public class Dklu extends Dklu_internal { for (k = n-1; k >= 0; k--) { - get_pointer(LU, Uip, Ulen, Ui, Ux, k, len); + GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, len); /* x[0] = X[k] / Udiag[k]; */ - div(x[0], X[k], Udiag[k]); + DIV (x[0], X[k], Udiag[k]); X[k] = x[0]; for (p = 0; p < len; p++) { /* X[Ui[p]] -= Ux[p] * x[0]; */ - mult_sub(X[Ui[p]], Ux[p], x[0]); + MULT_SUB (X[Ui[p]], Ux[p], x[0]); } } @@ -356,12 +356,12 @@ public class Dklu extends Dklu_internal { for (k = n-1; k >= 0; k--) { - get_pointer(LU, Uip, Ulen, Ui, Ux, k, len); + GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, len); ukk = Udiag[k]; /* x[0] = X[2*k ] / ukk; x[1] = X[2*k + 1] / ukk; */ - div(x[0], X[2*k], ukk); - div(x[1], X[2*k + 1], ukk); + DIV (x[0], X[2*k], ukk); + DIV (x[1], X[2*k + 1], ukk); X[2*k ] = x[0]; X[2*k + 1] = x[1]; @@ -371,8 +371,8 @@ public class Dklu extends Dklu_internal { uik = Ux[p]; /* X[2*i ] -= uik * x[0]; X[2*i + 1] -= uik * x[1]; */ - mult_sub(X[2*i], uik, x[0]); - mult_sub(X[2*i + 1], uik, x[1]); + MULT_SUB (X[2*i], uik, x[0]); + MULT_SUB (X[2*i + 1], uik, x[1]); } } @@ -382,12 +382,12 @@ public class Dklu extends Dklu_internal { for (k = n-1; k >= 0; k--) { - get_pointer(LU, Uip, Ulen, Ui, Ux, k, len); + GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, len); ukk = Udiag[k]; - div(x[0], X[3*k], ukk); - div(x[1], X[3*k + 1], ukk); - div(x[2], X[3*k + 2], ukk); + DIV (x[0], X[3*k], ukk); + DIV (x[1], X[3*k + 1], ukk); + DIV (x[2], X[3*k + 2], ukk); X[3*k ] = x[0]; X[3*k + 1] = x[1]; @@ -396,9 +396,9 @@ public class Dklu extends Dklu_internal { { i = Ui[p]; uik = Ux[p]; - mult_sub(X[3*i], uik, x[0]); - mult_sub(X[3*i + 1], uik, x[1]); - mult_sub(X[3*i + 2], uik, x[2]); + MULT_SUB (X[3*i], uik, x[0]); + MULT_SUB (X[3*i + 1], uik, x[1]); + MULT_SUB (X[3*i + 2], uik, x[2]); } } @@ -408,13 +408,13 @@ public class Dklu extends Dklu_internal { for (k = n-1; k >= 0; k--) { - get_pointer(LU, Uip, Ulen, Ui, Ux, k, len); + GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, len); ukk = Udiag[k]; - div(x[0], X[4*k], ukk); - div(x[1], X[4*k + 1], ukk); - div(x[2], X[4*k + 2], ukk); - div(x[3], X[4*k + 3], ukk); + DIV (x[0], X[4*k], ukk); + DIV (x[1], X[4*k + 1], ukk); + DIV (x[2], X[4*k + 2], ukk); + DIV (x[3], X[4*k + 3], ukk); X[4*k ] = x[0]; X[4*k + 1] = x[1]; @@ -425,10 +425,10 @@ public class Dklu extends Dklu_internal { i = Ui[p]; uik = Ux[p]; - mult_sub(X[4*i], uik, x[0]); - mult_sub(X[4*i + 1], uik, x[1]); - mult_sub(X[4*i + 2], uik, x[2]); - mult_sub(X[4*i + 3], uik, x[3]); + MULT_SUB (X[4*i], uik, x[0]); + MULT_SUB (X[4*i + 1], uik, x[1]); + MULT_SUB (X[4*i + 2], uik, x[2]); + MULT_SUB (X[4*i + 3], uik, x[3]); } } @@ -467,12 +467,12 @@ public class Dklu extends Dklu_internal { for (k = n-1; k >= 0; k--) { - get_pointer(LU, Lip, Llen, Li, Lx, k, len); + GET_POINTER (LU, Lip, Llen, Li, Lx, k, len); x[0] = X[k]; for (p = 0; p < len; p++) { /*x[0] -= Lx[p] * X[Li[p]];*/ - mult_sub(x[0], Lx[p], X[Li[p]]); + MULT_SUB (x[0], Lx[p], X[Li[p]]); } X[k] = x[0]; } @@ -484,13 +484,13 @@ public class Dklu extends Dklu_internal { { x[0] = X[2*k ]; x[1] = X[2*k + 1]; - get_pointer(LU, Lip, Llen, Li, Lx, k, len); + GET_POINTER (LU, Lip, Llen, Li, Lx, k, len); for (p = 0; p < len; p++) { i = Li[p]; lik = Lx[p]; - mult_sub(x[0], lik, X[2*i]); - mult_sub(x[1], lik, X[2*i + 1]); + MULT_SUB (x[0], lik, X[2*i]); + MULT_SUB (x[1], lik, X[2*i + 1]); } X[2*k ] = x[0]; X[2*k + 1] = x[1]; @@ -504,14 +504,14 @@ public class Dklu extends Dklu_internal { x[0] = X[3*k ]; x[1] = X[3*k + 1]; x[2] = X[3*k + 2]; - get_pointer(LU, Lip, Llen, Li, Lx, k, len); + GET_POINTER (LU, Lip, Llen, Li, Lx, k, len); for (p = 0; p < len; p++) { i = Li[p]; lik = Lx[p]; - mult_sub(x[0], lik, X[3*i]); - mult_sub(x[1], lik, X[3*i + 1]); - mult_sub(x[2], lik, X[3*i + 2]); + MULT_SUB (x[0], lik, X[3*i]); + MULT_SUB (x[1], lik, X[3*i + 1]); + MULT_SUB (x[2], lik, X[3*i + 2]); } X[3*k ] = x[0]; X[3*k + 1] = x[1]; @@ -527,15 +527,15 @@ public class Dklu extends Dklu_internal { x[1] = X[4*k + 1]; x[2] = X[4*k + 2]; x[3] = X[4*k + 3]; - get_pointer(LU, Lip, Llen, Li, Lx, k, len); + GET_POINTER (LU, Lip, Llen, Li, Lx, k, len); for (p = 0; p < len; p++) { i = Li[p]; lik = Lx[p]; - mult_sub(x[0], lik, X[4*i]); - mult_sub(x[1], lik, X[4*i + 1]); - mult_sub(x[2], lik, X[4*i + 2]); - mult_sub(x[3], lik, X[4*i + 3]); + MULT_SUB (x[0], lik, X[4*i]); + MULT_SUB (x[1], lik, X[4*i + 1]); + MULT_SUB (x[2], lik, X[4*i + 2]); + MULT_SUB (x[3], lik, X[4*i + 3]); } X[4*k ] = x[0]; X[4*k + 1] = x[1]; @@ -576,15 +576,15 @@ public class Dklu extends Dklu_internal { for (k = 0; k < n; k++) { - get_pointer(LU, Uip, Ulen, Ui, Ux, k, len); + GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, len); x[0] = X[k]; for (p = 0; p < len; p++) { /* x[0] -= Ux[p] * X[Ui[p]]; */ - mult_sub(x[0], Ux[p], X[Ui[p]]); + MULT_SUB (x[0], Ux[p], X[Ui[p]]); } ukk = Udiag[k]; - div(X[k], x[0], ukk); + DIV (X[k], x[0], ukk); } break; @@ -592,19 +592,19 @@ public class Dklu extends Dklu_internal { for (k = 0; k < n; k++) { - get_pointer(LU, Uip, Ulen, Ui, Ux, k, len); + GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, len); x[0] = X[2*k ]; x[1] = X[2*k + 1]; for (p = 0; p < len; p++) { i = Ui[p]; uik = Ux[p]; - mult_sub(x[0], uik, X[2*i]); - mult_sub(x[1], uik, X[2*i + 1]); + MULT_SUB (x[0], uik, X[2*i]); + MULT_SUB (x[1], uik, X[2*i + 1]); } ukk = Udiag[k]; - div(X[2*k], x[0], ukk); - div(X[2*k + 1], x[1], ukk); + DIV (X[2*k], x[0], ukk); + DIV (X[2*k + 1], x[1], ukk); } break; @@ -612,7 +612,7 @@ public class Dklu extends Dklu_internal { for (k = 0; k < n; k++) { - get_pointer(LU, Uip, Ulen, Ui, Ux, k, len); + GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, len); x[0] = X[3*k ]; x[1] = X[3*k + 1]; x[2] = X[3*k + 2]; @@ -620,14 +620,14 @@ public class Dklu extends Dklu_internal { { i = Ui[p]; uik = Ux[p]; - mult_sub(x[0], uik, X[3*i]); - mult_sub(x[1], uik, X[3*i + 1]); - mult_sub(x[2], uik, X[3*i + 2]); + MULT_SUB (x[0], uik, X[3*i]); + MULT_SUB (x[1], uik, X[3*i + 1]); + MULT_SUB (x[2], uik, X[3*i + 2]); } ukk = Udiag[k]; - div(X[3*k], x[0], ukk); - div(X[3*k + 1], x[1], ukk); - div(X[3*k + 2], x[2], ukk); + DIV (X[3*k], x[0], ukk); + DIV (X[3*k + 1], x[1], ukk); + DIV (X[3*k + 2], x[2], ukk); } break; @@ -635,7 +635,7 @@ public class Dklu extends Dklu_internal { for (k = 0; k < n; k++) { - get_pointer(LU, Uip, Ulen, Ui, Ux, k, len); + GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, len); x[0] = X[4*k ]; x[1] = X[4*k + 1]; x[2] = X[4*k + 2]; @@ -644,16 +644,16 @@ public class Dklu extends Dklu_internal { { i = Ui[p]; uik = Ux[p]; - mult_sub(x[0], uik, X[4*i]); - mult_sub(x[1], uik, X[4*i + 1]); - mult_sub(x[2], uik, X[4*i + 2]); - mult_sub(x[3], uik, X[4*i + 3]); + MULT_SUB (x[0], uik, X[4*i]); + MULT_SUB (x[1], uik, X[4*i + 1]); + MULT_SUB (x[2], uik, X[4*i + 2]); + MULT_SUB (x[3], uik, X[4*i + 3]); } ukk = Udiag[k]; - div(X[4*k], x[0], ukk); - div(X[4*k + 1], x[1], ukk); - div(X[4*k + 2], x[2], ukk); - div(X[4*k + 3], x[3], ukk); + DIV (X[4*k], x[0], ukk); + DIV (X[4*k + 1], x[1], ukk); + DIV (X[4*k + 2], x[2], ukk); + DIV (X[4*k + 3], x[3], ukk); } break; } diff --git a/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_internal.java b/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_internal.java index 7b3b907..603523b 100644 --- a/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_internal.java +++ b/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_internal.java @@ -29,13 +29,13 @@ public abstract class Dklu_internal extends Dklu_version { /** * enable debugging and assertions */ - static boolean NDEBUG = true; + protected static boolean NDEBUG = true ; - static void ASSERT (boolean a) + protected static void ASSERT (boolean a) { if (!NDEBUG) { - assert a; + assert a ; } } @@ -43,39 +43,39 @@ public abstract class Dklu_internal extends Dklu_version { * @return true if an integer (stored in double x) would overflow (or if * x is NaN) */ - static boolean INT_OVERFLOW (double x) + protected static boolean INT_OVERFLOW (double x) { return ((!(x * (1.0+1e-8) <= (double) INT_MAX)) - || SCALAR_IS_NAN (x)); + || SCALAR_IS_NAN (x)) ; } - static final int TRUE = 1; - static final int FALSE = 0; + protected static final int TRUE = 1 ; + protected static final int FALSE = 0 ; - static double MAX (double a, double b) + protected static double MAX (double a, double b) { - return a > b ? a : b; + return a > b ? a : b ; } - static double MIN (double a, double b) + protected static double MIN (double a, double b) { - return a < b ? a : b; + return a < b ? a : b ; } /* FLIP is a "negation about -1", and is used to mark an integer i that is * normally non-negative. FLIP (EMPTY) is EMPTY. FLIP of a number > EMPTY * is negative, and FLIP of a number < EMTPY is positive. FLIP (FLIP (i)) = i * for all integers i. UNFLIP (i) is >= EMPTY. */ - static final int EMPTY = -1; + protected static final int EMPTY = -1 ; - static double FLIP (double i) + protected static double FLIP (double i) { - return -i - 2; + return -i - 2 ; } - static double UNFLIP (double i) + protected static double UNFLIP (double i) { - return (i < EMPTY) ? FLIP (i) : i; + return (i < EMPTY) ? FLIP (i) : i ; } } diff --git a/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_kernel.java b/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_kernel.java index c7002d5..d851dd5 100644 --- a/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_kernel.java +++ b/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_kernel.java @@ -24,12 +24,14 @@ package edu.ufl.cise.klu.tdouble; +import edu.ufl.cise.klu.common.KLU_common; + /** * Sparse left-looking LU factorization, with partial pivoting. Based on * Gilbert & Peierl's method, with a non-recursive DFS and with Eisenstat & * Liu's symmetric pruning. No user-callable routines are in this file. */ -public class Dklu_kernel { +public class Dklu_kernel extends Dklu_internal { /** * Does a depth-first-search, starting at node j. @@ -61,19 +63,19 @@ public class Dklu_kernel { head = 0; Stack[0] = j; - assert (Flag[j] != k); + ASSERT (Flag[j] != k); while (head >= 0) { j = Stack[head]; jnew = Pinv[j]; - assert (jnew >= 0 && jnew < k); /* j is pivotal */ + ASSERT (jnew >= 0 && jnew < k); /* j is pivotal */ if (Flag[j] != k) /* a node is not yet visited */ { /* first time that j has been visited */ Flag[j] = k; - printf("[ start dfs at %d : new %d\n", j, jnew); + PRINTF ("[ start dfs at %d : new %d\n", j, jnew); /* set Ap_pos[head] to one past the last entry in col j to scan */ Ap_pos[head] = (Lpend[jnew] == EMPTY) ? Llen[jnew] : Lpend[jnew]; @@ -118,11 +120,11 @@ public class Dklu_kernel { * recursive stack and push j onto output stack */ head--; Stack[--top] = j; - printf(" end dfs at %d ] head : %d\n", j, head); + PRINTF (" end dfs at %d ] head : %d\n", j, head); } } - *plength = l_length; + plength = l_length; return top; } @@ -175,7 +177,7 @@ public class Dklu_kernel { if (i < 0) continue; /* skip entry outside the block */ /* (i,k) is an entry in the block. start a DFS at node i */ - printf("\n ===== DFS at node %d in b, inew: %d\n", i, Pinv[i]); + PRINTF ("\n ===== DFS at node %d in b, inew: %d\n", i, Pinv[i]); if (Flag[i] != k) { if (Pinv[i] >= 0) @@ -263,7 +265,7 @@ public class Dklu_kernel { oldrow = Ai[p]; i = PSinv[oldrow] - k1; aik = Ax[p]; - scale_div(aik, Rs[oldrow]); + SCALE_DIV (aik, Rs[oldrow]); if (i < 0) { /* this is an entry in the off-diagonal part */ @@ -313,14 +315,14 @@ public class Dklu_kernel { /* forward solve with column j of L */ j = Stack[s]; jnew = Pinv[j]; - assert (jnew >= 0); + ASSERT (jnew >= 0); xj = X[j]; - get_pointer(LU, Lip, Llen, Li, Lx, jnew, len); - assert (Lip[jnew] <= Lip[jnew+1]); + GET_POINTER (LU, Lip, Llen, Li, Lx, jnew, len); + ASSERT (Lip[jnew] <= Lip[jnew+1]); for (p = 0; p < len; p++) { /*X[Li[p]] -= Lx[p] * xj; */ - mult_sub(X[Li[p]], Lx[p], xj); + MULT_SUB (X[Li[p]], Lx[p], xj); } } } @@ -345,7 +347,7 @@ public class Dklu_kernel { * @param Common * @return */ - public static boolean lpivot(int diagrow, int[] p_pivrow, Entry[] p_pivot, + public static int lpivot(int diagrow, int[] p_pivrow, Entry[] p_pivot, double[] p_abs_pivot, double tol, Entry[] X, Unit[] LU, int[] Lip, int[] Llen, int k, int n, int[] Pinv , int[] p_firstrow, KLU_common Common) @@ -360,40 +362,40 @@ public class Dklu_kernel { if (Llen[k] == 0) { /* matrix is structurally singular */ - if (Common.halt_if_singular) + if (Common.halt_if_singular == 1) { - return false; + return FALSE; } for (firstrow = p_firstrow; firstrow < n; firstrow++) { - printf("check %d\n", firstrow); + PRINTF ("check %d\n", firstrow); if (Pinv[firstrow] < 0) { /* found the lowest-numbered non-pivotal row. Pick it. */ pivrow = firstrow; - printf("Got pivotal row: %d\n", pivrow); + PRINTF ("Got pivotal row: %d\n", pivrow); break; } } - assert (pivrow >= 0 && pivrow < n); - clear(pivot); + ASSERT (pivrow >= 0 && pivrow < n); + CLEAR (pivot); p_pivrow = pivrow; p_pivot = pivot; p_abs_pivot = 0; p_firstrow = firstrow; - return false; + return FALSE; } pdiag = EMPTY; ppivrow = EMPTY; abs_pivot = EMPTY; i = Llen[k] - 1; - get_pointer(LU, Lip, Llen, Li, Lx, k, len); + GET_POINTER (LU, Lip, Llen, Li, Lx, k, len); last_row_index = Li[i]; /* decrement the length by 1 */ Llen[k] = i; - get_pointer(LU, Lip, Llen, Li, Lx, k, len); + GET_POINTER (LU, Lip, Llen, Li, Lx, k, len); /* look in Li[0 ..Llen[k] - 1 ] for a pivot row */ for (p = 0; p < len; p++) @@ -401,11 +403,11 @@ public class Dklu_kernel { /* gather the entry from X and store in L */ i = Li[p]; x = X[i]; - clear(X[i]); + CLEAR (X[i]); Lx[p] = x; - /* xabs = abs(x); */ - abs(xabs, x); + /* xabs = ABS (x); */ + ABS (xabs, x); /* find the diagonal */ if (i == diagrow) @@ -421,8 +423,8 @@ public class Dklu_kernel { } } - /* xabs = abs(X[last_row_index]);*/ - abs(xabs, X[last_row_index]); + /* xabs = ABS (X[last_row_index]);*/ + ABS (xabs, X[last_row_index]); if (xabs > abs_pivot) { abs_pivot = xabs; @@ -440,8 +442,8 @@ public class Dklu_kernel { } else if (pdiag != EMPTY) { - /* xabs = abs(Lx[pdiag]);*/ - abs(xabs, Lx[pdiag]); + /* xabs = ABS (Lx[pdiag]);*/ + ABS (xabs, Lx[pdiag]); if (xabs >= tol * abs_pivot) { /* the diagonal is large enough */ @@ -463,27 +465,27 @@ public class Dklu_kernel { pivrow = last_row_index; pivot = X[last_row_index]; } - clear(X[last_row_index]); + CLEAR (X[last_row_index]); p_pivrow = pivrow; p_pivot = pivot; p_abs_pivot = abs_pivot; - assert (pivrow >= 0 && pivrow < n); + ASSERT (pivrow >= 0 && pivrow < n); - if (is_zero(pivot) && Common.halt_if_singular) + if (IS_ZERO (pivot) && Common.halt_if_singular) { /* numerically singular case */ - return false; + return FALSE; } /* divide L by the pivot value */ for (p = 0; p < Llen[k]; p++) { /* Lx[p] /= pivot; */ - div(Lx[p], Lx[p], pivot); + DIV (Lx[p], Lx[p], pivot); } - return true; + return TRUE; } /** @@ -509,29 +511,29 @@ public class Dklu_kernel { int p, i, j, p2, phead, ptail, llen, ulen; /* check to see if any column of L can be pruned */ - get_pointer(LU, Uip, Ulen, Ui, Ux, k, ulen); + GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, ulen); for (p = 0; p < ulen; p++) { j = Ui[p]; - assert (j < k); - printf(("%d is pruned: %d. Lpend[j] %d Lip[j+1] %d\n", - j, Lpend[j] != EMPTY, Lpend[j], Lip[j+1])); + ASSERT (j < k); + PRINTF ("%d is pruned: %d. Lpend[j] %d Lip[j+1] %d\n", + j, Lpend[j] != EMPTY, Lpend[j], Lip[j+1]); if (Lpend[j] == EMPTY) { /* scan column j of L for the pivot row */ - get_pointer(LU, Lip, Llen, Li, Lx, j, llen); + GET_POINTER (LU, Lip, Llen, Li, Lx, j, llen); for (p2 = 0; p2 < llen; p2++) { if (pivrow == Li[p2]) { /* found it! This column can be pruned */ - if (debug) + if (!NDEBUG) { - printf(("==== PRUNE: col j %d of L\n", j)); + PRINTF ("==== PRUNE: col j %d of L\n", j); int p3; for (p3 = 0; p3 < Llen[j]; p3++) { - printf("before: %i pivotal: %d\n", Li[p3], + PRINTF ("before: %i pivotal: %d\n", Li[p3], Pinv[Li[p3]] >= 0); } } @@ -568,13 +570,13 @@ public class Dklu_kernel { * column j as pruned. */ Lpend[j] = ptail; - if (debug) + if (!NDEBUG) { int p3; for (p3 = 0; p3 < Llen[j]; p3++) { - if (p3 == Lpend[j]) printf("----\n"); - printf("after: %i pivotal: %d\n", Li[p3], + if (p3 == Lpend[j]) PRINTF ("----\n"); + PRINTF ("after: %i pivotal: %d\n", Li[p3], Pinv[Li[p3]] >= 0); } } @@ -635,25 +637,25 @@ public class Dklu_kernel { int k, p, i, j, pivrow = 0, kbar, diagrow, firstrow, lup, top, scale, len; size_t newlusize; - if (debug) + if (!NDEBUG) { - Entry *Lx; + Entry Lx; } - assert (Common != null); + ASSERT (Common != null); scale = Common.scale; tol = Common.tol; memgrow = Common.memgrow; lnz = 0; unz = 0; - clear(pivot); + CLEAR (pivot); /* ---------------------------------------------------------------------- */ /* get initial Li, Lx, Ui, and Ux */ /* ---------------------------------------------------------------------- */ - printf("input: lusize %d \n", lusize); - assert (lusize > 0); + PRINTF ("input: lusize %d \n", lusize); + ASSERT (lusize > 0); LU = p_LU; /* ---------------------------------------------------------------------- */ @@ -666,7 +668,7 @@ public class Dklu_kernel { for (k = 0; k < n; k++) { /* X[k] = 0; */ - clear(X[k]); + CLEAR (X[k]); Flag[k] = EMPTY; Lpend[k] = EMPTY; /* flag k as not pruned */ } @@ -679,21 +681,21 @@ public class Dklu_kernel { for (k = 0; k < n; k++) { P[k] = k; - Pinv[k] = flip(k); /* mark all rows as non-pivotal */ + Pinv[k] = FLIP (k); /* mark all rows as non-pivotal */ } /* initialize the construction of the off-diagonal matrix */ Offp[0] = 0; - /* P[k] = row means that unflip(Pinv[row]) = k, and visa versa. + /* P[k] = row means that UNFLIP (Pinv[row]) = k, and visa versa. * If row is pivotal, then Pinv[row] >= 0. A row is initially "flipped" * (Pinv[k] < EMPTY), and then marked "unflipped" when it becomes * pivotal. */ - if (debug) + if (!NDEBUG) { for (k = 0; k < n; k++) { - printf("Initial P[%d] = %d\n", k, P[k]); + PRINTF ("Initial P[%d] = %d\n", k, P[k]); } } @@ -704,28 +706,28 @@ public class Dklu_kernel { for (k = 0; k < n; k++) { - printf("\n\n==================================== k: %d\n", k); + PRINTF ("\n\n==================================== k: %d\n", k); /* ------------------------------------------------------------------ */ /* determine if LU factors have grown too big */ /* ------------------------------------------------------------------ */ /* (n - k) entries for L and k entries for U */ - nunits = Dunits(Int, n - k) + Dunits(Int, k) + - Dunits(Entry, n - k) + Dunits(Entry, k); + nunits = DUNITS (Int, n - k) + DUNITS (Int, k) + + DUNITS (Entry, n - k) + DUNITS (Entry, k); /* LU can grow by at most 'nunits' entries if the column is dense */ - printf("lup %d lusize %g lup+nunits: %g\n", lup, (double) lusize, + PRINTF ("lup %d lusize %g lup+nunits: %g\n", lup, (double) lusize, lup+nunits); xsize = ((double) lup) + nunits; if (xsize > (double) lusize) { /* check here how much to grow */ xsize = (memgrow * ((double) lusize) + 4*n + 1); - if (int_overflow(xsize)) + if (INT_OVERFLOW (xsize)) { - printf("Matrix is too large (int overflow)\n"); - Common.status = KLU_TOO_LARGE; + PRINTF ("Matrix is too large (int overflow)\n"); + Common.status = KLU_common.KLU_TOO_LARGE; return (lusize); } newlusize = memgrow * lusize + 2*n + 1; @@ -733,13 +735,13 @@ public class Dklu_kernel { LU = Dklu_realloc.klu_realloc(newlusize, lusize, sizeof (Unit), LU, Common); Common.nrealloc++; p_LU = LU; - if (Common.status == KLU_OUT_OF_MEMORY) + if (Common.status == KLU_common.KLU_OUT_OF_MEMORY) { - printf(("Matrix is too large (LU)\n")); + PRINTF (("Matrix is too large (LU)\n")); return (lusize); } lusize = newlusize; - printf("inc LU to %d done\n", lusize); + PRINTF ("inc LU to %d done\n", lusize); } /* ------------------------------------------------------------------ */ @@ -752,40 +754,40 @@ public class Dklu_kernel { /* compute the nonzero pattern of the kth column of L and U */ /* ------------------------------------------------------------------ */ - if (debug) + if (!NDEBUG) { for (i = 0; i < n; i++) { - assert (Flag[i] < k); - /* assert (X[i] == 0); */ - assert (is_zero(X[i])); + ASSERT (Flag[i] < k); + /* ASSERT (X[i] == 0); */ + ASSERT (IS_ZERO (X[i])); } } - top = lsolve_symbolic(n, k, Ap, Ai, Q, Pinv, Stack, Flag, + top = lsolve_symbolic (n, k, Ap, Ai, Q, Pinv, Stack, Flag, Lpend, Ap_pos, LU, lup, Llen, Lip, k1, PSinv); - if (debug) + if (!NDEBUG) { - printf(("--- in U:\n")); + PRINTF (("--- in U:\n")); for (p = top; p < n; p++) { - printf(("pattern of X for U: %d : %d pivot row: %d\n", - p, Stack[p], Pinv[Stack[p]])); - assert (Flag[Stack[p]] == k); + PRINTF ("pattern of X for U: %d : %d pivot row: %d\n", + p, Stack[p], Pinv[Stack[p]]); + ASSERT (Flag[Stack[p]] == k); } - printf(("--- in L:\n")); - Li = (int *) (LU + Lip[k]); + PRINTF ("--- in L:\n"); + Li = (int[]) (LU + Lip[k]); for (p = 0; p < Llen[k]; p++) { - printf(("pattern of X in L: %d : %d pivot row: %d\n", - p, Li[p], Pinv[Li[p]])); - assert (Flag[Li[p]] == k); + PRINTF ("pattern of X in L: %d : %d pivot row: %d\n", + p, Li[p], Pinv[Li[p]]); + ASSERT (Flag[Li[p]] == k); } p = 0; for (i = 0; i < n; i++) { - assert (Flag[i] <= k); + ASSERT (Flag[i] <= k); if (Flag[i] == k) p++; } } @@ -794,27 +796,27 @@ public class Dklu_kernel { /* get the column of the matrix to factorize and scatter into X */ /* ------------------------------------------------------------------ */ - construct_column(k, Ap, Ai, Ax, Q, X, + construct_column (k, Ap, Ai, Ax, Q, X, k1, PSinv, Rs, scale, Offp, Offi, Offx); /* ------------------------------------------------------------------ */ /* compute the numerical values of the kth column (s = L \ A (:,k)) */ /* ------------------------------------------------------------------ */ - lsolve_numeric(Pinv, LU, Stack, Lip, top, n, Llen, X); + lsolve_numeric (Pinv, LU, Stack, Lip, top, n, Llen, X); - if (debug) + if (!NDEBUG) { for (p = top; p < n; p++) { - printf("X for U %d : ", Stack[p]); - print_entry(X[Stack[p]]); + PRINTF ("X for U %d : ", Stack[p]); + PRINT_ENTRY (X[Stack[p]]); } Li = (int[]) (LU + Lip[k]); for (p = 0; p < Llen[k]; p++) { - printf("X for L %d : ", Li[p]); - print_entry(X[Li[p]]); + PRINTF ("X for L %d : ", Li[p]); + PRINT_ENTRY (X[Li[p]]); } } @@ -824,21 +826,21 @@ public class Dklu_kernel { /* determine what the "diagonal" is */ diagrow = P[k]; /* might already be pivotal */ - printf("k %d, diagrow = %d, unflip(diagrow) = %d\n", - k, diagrow, unflip(diagrow)); + PRINTF ("k %d, diagrow = %d, UNFLIP (diagrow) = %d\n", + k, diagrow, UNFLIP (diagrow)); /* find a pivot and scale the pivot column */ - if (!lpivot(diagrow, pivrow, pivot, abs_pivot, tol, X, LU, Lip, + if (!lpivot (diagrow, pivrow, pivot, abs_pivot, tol, X, LU, Lip, Llen, k, n, Pinv, firstrow, Common)) { /* matrix is structurally or numerically singular */ - Common.status = KLU_SINGULAR; + Common.status = KLU_common.KLU_SINGULAR; if (Common.numerical_rank == EMPTY) { Common.numerical_rank = k+k1; Common.singular_col = Q[k+k1]; } - if (Common.halt_if_singular) + if (Common.halt_if_singular == 1) { /* do not continue the factorization */ return (lusize); @@ -847,32 +849,32 @@ public class Dklu_kernel { /* we now have a valid pivot row, even if the column has NaN's or * has no entries on or below the diagonal at all. */ - printf(("\nk %d : Pivot row %d : ", k, pivrow)); - print_entry(pivot); - assert (pivrow >= 0 && pivrow < n); - assert (Pinv[pivrow] < 0); + PRINTF ("\nk %d : Pivot row %d : ", k, pivrow); + PRINT_ENTRY (pivot); + ASSERT (pivrow >= 0 && pivrow < n); + ASSERT (Pinv[pivrow] < 0); /* set the Uip pointer */ - Uip[k] = Lip[k] + units(int, Llen[k]) + units(Entry, Llen[k]); + Uip[k] = Lip[k] + UNITS (int, Llen[k]) + UNITS (Entry, Llen[k]); /* move the lup pointer to the position where indices of U * should be stored */ - lup += units(int, Llen[k]) + units(Entry, Llen[k]); + lup += UNITS (int, Llen[k]) + UNITS (Entry, Llen[k]); Ulen[k] = n - top; /* extract Stack[top..n-1] to Ui and the values to Ux and clear X */ - get_pointer(LU, Uip, Ulen, Ui, Ux, k, len); + GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, len); for (p = top, i = 0; p < n; p++, i++) { j = Stack[p]; Ui[i] = Pinv[j]; Ux[i] = X[j]; - clear(X[j]); + CLEAR (X[j]); } /* position the lu index at the starting point for next column */ - lup += units(int, Ulen[k]) + units(Entry, Ulen[k]); + lup += UNITS (int, Ulen[k]) + UNITS (Entry, Ulen[k]); /* U(k,k) = pivot */ Udiag[k] = pivot; @@ -881,43 +883,43 @@ public class Dklu_kernel { /* log the pivot permutation */ /* ------------------------------------------------------------------ */ - assert (unflip(Pinv[diagrow]) < n); - assert (P[unflip(Pinv[diagrow])] == diagrow); + ASSERT (UNFLIP (Pinv[diagrow]) < n); + ASSERT (P[UNFLIP (Pinv[diagrow])] == diagrow); if (pivrow != diagrow) { /* an off-diagonal pivot has been chosen */ Common.noffdiag++; - printf((">>>>>>>>>>>>>>>>> pivrow %d k %d off-diagonal\n", - pivrow, k)); + PRINTF (">>>>>>>>>>>>>>>>> pivrow %d k %d off-diagonal\n", + pivrow, k); if (Pinv[diagrow] < 0) { /* the former diagonal row index, diagrow, has not yet been * chosen as a pivot row. Log this diagrow as the "diagonal" * entry in the column kbar for which the chosen pivot row, * pivrow, was originally logged as the "diagonal" */ - kbar = flip(Pinv[pivrow]); + kbar = FLIP (Pinv[pivrow]); P[kbar] = diagrow; - Pinv[diagrow] = flip(kbar); + Pinv[diagrow] = FLIP (kbar); } } P[k] = pivrow; Pinv[pivrow] = k; - if (debug) + if (!NDEBUG) { - for (i = 0; i < n; i++) { assert (is_zero(X[i]));} - get_pointer(LU, Uip, Ulen, Ui, Ux, k, len); + for (i = 0; i < n; i++) { ASSERT (IS_ZERO (X[i]));} + GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, len); for (p = 0; p < len; p++) { - printf(("Column %d of U: %d : ", k, Ui[p])); - print_entry(Ux[p]); + PRINTF ("Column %d of U: %d : ", k, Ui[p]); + PRINT_ENTRY (Ux[p]); } - get_pointer(LU, Lip, Llen, Li, Lx, k, len); + GET_POINTER (LU, Lip, Llen, Li, Lx, k, len); for (p = 0; p < len; p++) { - printf(("Column %d of L: %d : ", k, Li[p])); - print_entry(Lx[p]); + PRINTF ("Column %d of L: %d : ", k, Li[p]); + PRINT_ENTRY (Lx[p]); } } @@ -925,7 +927,7 @@ public class Dklu_kernel { /* symmetric pruning */ /* ------------------------------------------------------------------ */ - prune(Lpend, Pinv, k, pivrow, LU, Uip, Lip, Ulen, Llen); + prune (Lpend, Pinv, k, pivrow, LU, Uip, Lip, Ulen, Llen); lnz += Llen[k] + 1; /* 1 added to lnz for diagonal */ unz += Ulen[k] + 1; /* 1 added to unz for diagonal */ @@ -937,25 +939,25 @@ public class Dklu_kernel { for (p = 0; p < n; p++) { - Li = (int *) (LU + Lip[p]); + Li = (int[]) (LU + Lip[p]); for (i = 0; i < Llen[p]; i++) { Li[i] = Pinv[Li[i]]; } } - if (debug) + if (!NDEBUG) { for (i = 0; i < n; i++) { - printf(("P[%d] = %d Pinv[%d] = %d\n", i, P[i], i, Pinv[i])); + PRINTF ("P[%d] = %d Pinv[%d] = %d\n", i, P[i], i, Pinv[i]); } for (i = 0; i < n; i++) { - assert (Pinv[i] >= 0 && Pinv[i] < n); - assert (P[i] >= 0 && P[i] < n); - assert (P[Pinv[i]] == i); - assert (is_zero(X[i])); + ASSERT (Pinv[i] >= 0 && Pinv[i] < n); + ASSERT (P[i] >= 0 && P[i] < n); + ASSERT (P[Pinv[i]] == i); + ASSERT (IS_ZERO (X[i])); } } @@ -964,7 +966,7 @@ public class Dklu_kernel { /* ---------------------------------------------------------------------- */ newlusize = lup; - assert ((size_t) newlusize <= lusize); + ASSERT ((size_t) newlusize <= lusize); /* this cannot fail, since the block is descreasing in size */ LU = Dklu_realloc.klu_realloc(newlusize, lusize, sizeof (Unit), LU, 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 db40cf9..571ba60 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 @@ -24,17 +24,19 @@ package edu.ufl.cise.klu.tdouble; +import edu.ufl.cise.klu.common.KLU_common; + /** * KLU memory management routines. */ -public class Dklu_memory { +public class Dklu_memory extends Dklu_internal { /** * Safely compute a+b, and check for size_t overflow. */ public static size_t klu_add_size_t(size_t a, size_t b, int ok) { - ok = ok && ((a + b) >= max(a, b)); + ok = ok && ((a + b) >= MAX (a, b)); return ok ? (a + b) : (size_t) -1; } @@ -71,7 +73,7 @@ public class Dklu_memory { { Object p ; size_t s ; - Int ok = true ; + int ok = TRUE ; if (Common == null) { @@ -80,30 +82,30 @@ public class Dklu_memory { else if (size == 0) { /* size must be > 0 */ - Common.status = KLU_INVALID ; + Common.status = KLU_common.KLU_INVALID ; p = null ; } else if (n >= INT_MAX) { /* object is too big to allocate; p[i] where i is an Int will not * be enough. */ - Common.status = KLU_TOO_LARGE ; + Common.status = KLU_common.KLU_TOO_LARGE ; p = null ; } else { /* call malloc, or its equivalent */ - s = klu_mult_size_t(max(1, n), size, ok) ; - p = ok ? ((Common.malloc_memory) (s)) : null ; + s = klu_mult_size_t (MAX (1, n), size, ok) ; + p = ok == 1 ? ((Common.malloc_memory) (s)) : null ; if (p == null) { /* failure: out of memory */ - Common.status = KLU_OUT_OF_MEMORY ; + Common.status = KLU_common.KLU_OUT_OF_MEMORY ; } else { Common.memusage += s ; - Common.mempeak = max(Common.mempeak, Common.memusage) ; + Common.mempeak = MAX (Common.mempeak, Common.memusage) ; } } return (p) ; @@ -120,18 +122,18 @@ public class Dklu_memory { public static void klu_free(Object p, size_t n, size_t size, KLU_common Common) { size_t s ; - Int ok = true ; + int ok = TRUE ; if (p != null && Common != null) { /* only free the object if the pointer is not null */ /* call free, or its equivalent */ (Common.free_memory) (p) ; - s = klu_mult_size_t(max(1, n), size, ok) ; + s = klu_mult_size_t (MAX (1, n), size, ok) ; Common.memusage -= s ; } /* return null, and the caller should assign this to p. This avoids * freeing the same pointer twice. */ - return (null) ; + //return (null) ; } /** @@ -163,7 +165,7 @@ public class Dklu_memory { { Object pnew ; size_t snew, sold ; - Int ok = true ; + int ok = TRUE ; if (Common == null) { @@ -172,36 +174,36 @@ public class Dklu_memory { else if (size == 0) { /* size must be > 0 */ - Common.status = KLU_INVALID ; + Common.status = KLU_common.KLU_INVALID ; p = null ; } else if (p == null) { /* A fresh object is being allocated. */ - p = KLU_malloc (nnew, size, Common) ; + p = klu_malloc (nnew, size, Common) ; } else if (nnew >= INT_MAX) { /* failure: nnew is too big. Do not change p */ - Common.status = KLU_TOO_LARGE ; + Common.status = KLU_common.KLU_TOO_LARGE ; } else { /* The object exists, and is changing to some other nonzero size. */ /* call realloc, or its equivalent */ - snew = klu_mult_size_t(max(1,nnew), size, ok) ; - sold = klu_mult_size_t(max(1,nold), size, ok) ; - pnew = ok ? ((Common.realloc_memory) (p, snew)) : null ; + snew = klu_mult_size_t (MAX (1,nnew), size, ok) ; + sold = klu_mult_size_t (MAX (1,nold), size, ok) ; + pnew = ok == 1 ? ((Common.realloc_memory) (p, snew)) : null ; if (pnew == null) { /* Do not change p, since it still points to allocated memory */ - Common.status = KLU_OUT_OF_MEMORY ; + Common.status = KLU_common.KLU_OUT_OF_MEMORY ; } else { /* success: return the new p and change the size of the block */ Common.memusage += (snew - sold) ; - Common.mempeak = max(Common.mempeak, Common.memusage) ; + Common.mempeak = MAX (Common.mempeak, Common.memusage) ; p = pnew ; } } 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 70bb02b..a9a9a51 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 @@ -24,13 +24,17 @@ package edu.ufl.cise.klu.tdouble; +import edu.ufl.cise.klu.common.KLU_common; +import edu.ufl.cise.klu.common.KLU_numeric; +import edu.ufl.cise.klu.common.KLU_symbolic; + /** * Factor the matrix, after ordering and analyzing it with KLU_analyze, and * factoring it once with KLU_factor. This routine cannot do any numerical * pivoting. The pattern of the input matrix (Ap, Ai) must be identical to * the pattern given to KLU_factor. */ -public class Dklu_refactor { +public class Dklu_refactor extends Dklu_internal { /** * @@ -42,7 +46,7 @@ public class Dklu_refactor { * @param Common * @return true if successful, false otherwise */ - public static boolean klu_refactor(int[] Ap, int[] Ai, double[] Ax, + public static int klu_refactor(int[] Ap, int[] Ai, double[] Ax, KLU_symbolic Symbolic, KLU_numeric Numeric, KLU_common Common) { Entry ukk, ujk, s; @@ -60,15 +64,15 @@ public class Dklu_refactor { if (Common == null) { - return false; + return FALSE; } - Common.status = KLU_OK; + Common.status = KLU_common.KLU_OK; if (Numeric == null) { /* invalid Numeric object */ - Common.status = KLU_INVALID; - return false; + Common.status = KLU_common.KLU_INVALID; + return FALSE; } Common.numerical_rank = EMPTY; @@ -105,10 +109,10 @@ public class Dklu_refactor { if (Numeric.Rs == null) { Numeric.Rs = KLU_malloc (n, sizeof(Double), Common); - if (Common.status < KLU_OK) + if (Common.status < KLU_common.KLU_OK) { - Common.status = KLU_OUT_OF_MEMORY; - return false; + Common.status = KLU_common.KLU_OUT_OF_MEMORY; + return FALSE; } } } @@ -116,7 +120,7 @@ public class Dklu_refactor { { /* no scaling for refactorization; ensure Numeric.Rs is freed. This * does nothing if Numeric.Rs is already null. */ - Numeric.Rs = KLU_free (Numeric.Rs, n, sizeof(Double), Common); + Numeric.Rs = Dklu_free.klu_free (Numeric.Rs, n, sizeof(Double), Common); } Rs = Numeric.Rs; @@ -134,9 +138,9 @@ public class Dklu_refactor { if (scale >= 0) { /* check for out-of-range indices, but do not check for duplicates */ - if (!Dklu_scale.klu_scale(scale, n, Ap, Ai, Ax, Rs, null, Common)) + if (Dklu_scale.klu_scale (scale, n, Ap, Ai, Ax, Rs, null, Common) == 0) { - return false; + return FALSE; } } @@ -147,7 +151,7 @@ public class Dklu_refactor { for (k = 0; k < maxblock; k++) { /* X[k] = 0 */ - clear(X[k]); + CLEAR (X[k]); } poff = 0; @@ -183,7 +187,7 @@ public class Dklu_refactor { oldcol = Q[k1]; pend = Ap[oldcol+1]; - clear(s); + CLEAR (s); for (p = Ap[oldcol]; p < pend; p++) { newrow = Pinv[Ai[p]] - k1; @@ -244,27 +248,27 @@ public class Dklu_refactor { /* compute kth column of U, and update kth column of A */ /* -------------------------------------------------- */ - get_pointer(LU, Uip, Ulen, Ui, Ux, k, ulen); + GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, ulen); for (up = 0; up < ulen; up++) { j = Ui[up]; ujk = X[j]; /* X[j] = 0 */ - clear(X[j]); + CLEAR (X[j]); Ux[up] = ujk; - get_pointer(LU, Lip, Llen, Li, Lx, j, llen); + GET_POINTER (LU, Lip, Llen, Li, Lx, j, llen); for (p = 0; p < llen; p++) { /* X[Li[p]] -= Lx[p] * ujk */ - mult_sub(X[Li[p]], Lx[p], ujk); + MULT_SUB (X[Li[p]], Lx[p], ujk); } } /* get the diagonal entry of U */ ukk = X[k]; /* X[k] = 0 */ - clear(X[k]); + CLEAR (X[k]); /* singular case */ - if (is_zero(ukk)) + if (IS_ZERO (ukk)) { /* matrix is numerically singular */ Common.status = KLU_SINGULAR; @@ -276,17 +280,17 @@ public class Dklu_refactor { if (Common.halt_if_singular) { /* do not continue the factorization */ - return false; + return FALSE; } } Udiag[k+k1] = ukk; /* gather and divide by pivot to get kth column of L */ - get_pointer(LU, Lip, Llen, Li, Lx, k, llen); + GET_POINTER (LU, Lip, Llen, Li, Lx, k, llen); for (p = 0; p < llen; p++) { i = Li[p]; - div(Lx[p], X[i], ukk); - clear(X[i]); + DIV (Lx[p], X[i], ukk); + CLEAR (X[i]); } } @@ -321,7 +325,7 @@ public class Dklu_refactor { oldcol = Q[k1]; pend = Ap[oldcol+1]; - clear(s); + CLEAR (s); for (p = Ap[oldcol]; p < pend; p++) { oldrow = Ai[p]; @@ -330,14 +334,14 @@ public class Dklu_refactor { { /* entry in off-diagonal block */ /* Offx[poff] = Az[p] / Rs[oldrow] */ - scale_div_assign(Offx[poff], Az[p], Rs[oldrow]); + SCALE_DIV_ASSIGN (Offx[poff], Az[p], Rs[oldrow]); poff++; } else { /* singleton */ /* s = Az[p] / Rs[oldrow] */ - scale_div_assign(s, Az[p], Rs[oldrow]); + SCALE_DIV_ASSIGN (s, Az[p], Rs[oldrow]); } } Udiag[k1] = s; @@ -373,14 +377,14 @@ public class Dklu_refactor { { /* entry in off-diagonal part */ /* Offx[poff] = Az[p] / Rs[oldrow] */ - scale_div_assign(Offx[poff], Az[p], Rs[oldrow]); + SCALE_DIV_ASSIGN (Offx[poff], Az[p], Rs[oldrow]); poff++; } else { /* (newrow,k) is an entry in the block */ /* X[newrow] = Az[p] / Rs[oldrow] */ - scale_div_assign(X[newrow], Az[p], Rs[oldrow]); + SCALE_DIV_ASSIGN (X[newrow], Az[p], Rs[oldrow]); } } @@ -388,49 +392,49 @@ public class Dklu_refactor { /* compute kth column of U, and update kth column of A */ /* -------------------------------------------------- */ - get_pointer(LU, Uip, Ulen, Ui, Ux, k, ulen); + GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, ulen); for (up = 0; up < ulen; up++) { j = Ui[up]; ujk = X[j]; /* X[j] = 0 */ - clear(X[j]); + CLEAR (X[j]); Ux[up] = ujk; - get_pointer(LU, Lip, Llen, Li, Lx, j, llen); + GET_POINTER (LU, Lip, Llen, Li, Lx, j, llen); for (p = 0; p < llen; p++) { /* X[Li[p]] -= Lx[p] * ujk */ - mult_sub(X[Li[p]], Lx[p], ujk); + MULT_SUB (X[Li[p]], Lx[p], ujk); } } /* get the diagonal entry of U */ ukk = X[k]; /* X[k] = 0 */ - clear(X[k]); + CLEAR (X[k]); /* singular case */ if (is_zero(ukk)) { /* matrix is numerically singular */ - Common.status = KLU_SINGULAR; + Common.status = KLU_common.KLU_SINGULAR; if (Common.numerical_rank == EMPTY) { Common.numerical_rank = k+k1; Common.singular_col = Q[k+k1]; } - if (Common.halt_if_singular) + if (Common.halt_if_singular == 1) { /* do not continue the factorization */ - return false; + return FALSE; } } Udiag[k+k1] = ukk; /* gather and divide by pivot to get kth column of L */ - get_pointer(LU, Lip, Llen, Li, Lx, k, llen); + GET_POINTER (LU, Lip, Llen, Li, Lx, k, llen); for (p = 0; p < llen; p++) { i = Li[p]; - div(Lx[p], X[i], ukk); - clear(X[i]); + DIV (Lx[p], X[i], ukk); + CLEAR (X[i]); } } } @@ -445,53 +449,53 @@ public class Dklu_refactor { { for (k = 0; k < n; k++) { - real(X[k]) = Rs[Pnum[k]]; + REAL (X[k]) = Rs[Pnum[k]]; } for (k = 0; k < n; k++) { - Rs[k] = real(X[k]); + Rs[k] = REAL (X[k]); } } if (false) { - assert (Offp[n] == poff); - assert (Symbolic.nzoff == poff); - printf("\n------------------- Off diagonal entries, new:\n"); - assert Dklu_valid.klu_valid(n, Offp, Offi, Offx); - if (Common.status == KLU_OK) + ASSERT (Offp[n] == poff); + ASSERT (Symbolic.nzoff == poff); + PRINTF ("\n------------------- Off diagonal entries, new:\n"); + ASSERT (Dklu_valid.klu_valid(n, Offp, Offi, Offx)); + if (Common.status == KLU_common.KLU_OK) { - printf("\n ########### KLU_BTF_REFACTOR done, nblocks %d\n", + PRINTF ("\n ########### KLU_BTF_REFACTOR done, nblocks %d\n", nblocks); for (block = 0; block < nblocks; block++) { k1 = R[block]; k2 = R[block+1]; nk = k2 - k1; - printf( + PRINTF ( "\n================KLU_refactor output: k1 %d k2 %d nk %d\n", k1, k2, nk); if (nk == 1) { - printf("singleton "); - print_entry(Udiag[k1]); + PRINTF ("singleton "); + PRINT_ENTRY (Udiag[k1]); } else { Lip = Numeric.Lip + k1; Llen = Numeric.Llen + k1; LU = (Unit[]) Numeric.LUbx[block]; - printf("\n---- L block %d\n", block); - assert (Dklu_valid_LU.klu_valid_LU(nk, true, Lip, Llen, LU)); + PRINTF ("\n---- L block %d\n", block); + ASSERT (Dklu_valid_LU.klu_valid_LU(nk, true, Lip, Llen, LU)); Uip = Numeric.Uip + k1; Ulen = Numeric.Ulen + k1; - printf("\n---- U block %d\n", block); - assert (Dklu_valid_LU.klu_valid_LU(nk, false, Uip, Ulen, LU)); + PRINTF ("\n---- U block %d\n", block); + ASSERT (Dklu_valid_LU.klu_valid_LU(nk, false, Uip, Ulen, LU)); } } } } - return true; + return TRUE; } } diff --git a/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_scale.java b/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_scale.java index f999c96..d8ea791 100644 --- a/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_scale.java +++ b/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_scale.java @@ -24,6 +24,8 @@ package edu.ufl.cise.klu.tdouble; +import edu.ufl.cise.klu.common.KLU_common; + /** * Scale a matrix and check to see if it is valid. Can be called by the user. * This is called by KLU_factor and KLU_refactor. Returns true if the input @@ -36,7 +38,7 @@ package edu.ufl.cise.klu.tdouble; * 1: the scale factor for row i is sum (abs (A (i,:))) * 2 or more: the scale factor for row i is max(abs (A (i,:))) */ -public class Dklu_scale { +public class Dklu_scale extends Dklu_internal { /** * @@ -50,7 +52,7 @@ public class Dklu_scale { * @param Common * @return true if successful, false otherwise */ - public static boolean klu_scale(int scale, int n, int[] Ap, int[] Ai, + public static int klu_scale(int scale, int n, int[] Ap, int[] Ai, double[] Ax, double[] Rs, int[] W, KLU_common Common) { double a; @@ -64,15 +66,15 @@ public class Dklu_scale { if (Common == null) { - return false; + return FALSE; } - Common.status = KLU_OK; + Common.status = KLU_common.KLU_OK; if (scale < 0) { /* return without checking anything and without computing the * scale factors */ - return true; + return TRUE; } Az = (Entry[]) Ax; @@ -81,22 +83,22 @@ public class Dklu_scale { (scale > 0 && Rs == null)) { /* Ap, Ai, Ax and Rs must be present, and n must be > 0 */ - Common.status = KLU_INVALID; - return false; + Common.status = KLU_common.KLU_INVALID; + return FALSE; } if (Ap[0] != 0 || Ap[n] < 0) { /* nz = Ap[n] must be >= 0 and Ap[0] must equal zero */ - Common.status = KLU_INVALID; - return false; + Common.status = KLU_common.KLU_INVALID; + return FALSE; } for (col = 0; col < n; col++) { if (Ap[col] > Ap[col+1]) { /* column pointers must be non-decreasing */ - Common.status = KLU_INVALID; - return false; + Common.status = KLU_common.KLU_INVALID; + return FALSE; } } @@ -132,16 +134,16 @@ public class Dklu_scale { if (row < 0 || row >= n) { /* row index out of range, or duplicate entry */ - Common.status = KLU_INVALID; - return false; + Common.status = KLU_common.KLU_INVALID; + return FALSE; } if (check_duplicates) { if (W[row] == col) { /* duplicate entry */ - Common.status = KLU_INVALID; - return false; + Common.status = KLU_common.KLU_INVALID; + return FALSE; } /* flag row i as appearing in column col */ W[row] = col; @@ -156,7 +158,7 @@ public class Dklu_scale { else if (scale > 1) { /* find the maxabs. value in the row */ - Rs[row] = max(Rs[row], a); + Rs[row] = MAX (Rs[row], a); } } } @@ -167,17 +169,17 @@ public class Dklu_scale { for (row = 0; row < n; row++) { /* matrix is singular */ - printf("Rs[%d] = %g\n", row, Rs[row]); + PRINTF ("Rs[%d] = %g\n", row, Rs[row]); if (Rs[row] == 0.0) { - printf("Row %d of A is all zero\n", row); + PRINTF ("Row %d of A is all zero\n", row); Rs[row] = 1.0; } } } - return true; + return TRUE; } } 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 04a021e..a9d16ef 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 @@ -24,6 +24,10 @@ package edu.ufl.cise.klu.tdouble; +import edu.ufl.cise.klu.common.KLU_common; +import edu.ufl.cise.klu.common.KLU_numeric; +import edu.ufl.cise.klu.common.KLU_symbolic; + /** * Solve Ax=b using the symbolic and numeric objects from KLU_analyze * (or KLU_analyze_given) and KLU_factor. Note that no iterative refinement is @@ -31,7 +35,7 @@ package edu.ufl.cise.klu.tdouble; * of size 4n Entry's (note that columns 2 to 4 of Xwork overlap with * Numeric.Iwork). */ -public class Dklu_solve { +public class Dklu_solve extends Dklu_internal { /** * @@ -44,7 +48,7 @@ public class Dklu_solve { * @param Common * @return */ - public static boolean klu_solve(KLU_symbolic Symbolic, KLU_numeric Numeric, + public static int klu_solve(KLU_symbolic Symbolic, KLU_numeric Numeric, int d, int nrhs, double[] B, KLU_common Common) { Entry offik, s; @@ -61,15 +65,15 @@ public class Dklu_solve { if (Common == null) { - return false; + return FALSE; } if (Numeric == null || Symbolic == null || d < Symbolic.n || nrhs < 0 || B == null) { - Common.status = KLU_INVALID; - return false; + Common.status = KLU_common.KLU_INVALID; + return FALSE; } - Common.status = KLU_OK; + Common.status = KLU_common.KLU_OK; /* ------------------------------------------------------------------ */ /* get the contents of the Symbolic object */ @@ -85,7 +89,7 @@ public class Dklu_solve { /* get the contents of the Numeric object */ /* ------------------------------------------------------------------ */ - assert (nblocks == Numeric.nblocks); + ASSERT (nblocks == Numeric.nblocks); Pnum = Numeric.Pnum; Offp = Numeric.Offp; Offi = Numeric.Offi; @@ -101,7 +105,7 @@ public class Dklu_solve { Rs = Numeric.Rs; X = (Entry) Numeric.Xwork; - assert Dklu_valid.klu_valid(n, Offp, Offi, Offx); + ASSERT (Dklu_valid.klu_valid(n, Offp, Offi, Offx)); /* ------------------------------------------------------------------ */ /* solve in chunks of 4 columns at a time */ @@ -114,7 +118,7 @@ public class Dklu_solve { /* get the size of the current chunk */ /* -------------------------------------------------------------- */ - nr = min(nrhs - chunk, 4); + nr = MIN (nrhs - chunk, 4); /* -------------------------------------------------------------- */ /* scale and permute the right hand side, X = P*(R\B) */ @@ -180,7 +184,7 @@ public class Dklu_solve { for (k = 0; k < n; k++) { - scale_div_assign(X[k], Bz[Pnum[k]], Rs[k]); + SCALE_DIV_ASSIGN (X[k], Bz[Pnum[k]], Rs[k]); } break; @@ -190,8 +194,8 @@ public class Dklu_solve { { i = Pnum[k]; rs = Rs[k]; - scale_div_assign(X[2*k], Bz[i], rs); - scale_div_assign(X[2*k + 1], Bz[i + d], rs); + SCALE_DIV_ASSIGN (X[2*k], Bz[i], rs); + SCALE_DIV_ASSIGN (X[2*k + 1], Bz[i + d], rs); } break; @@ -201,9 +205,9 @@ public class Dklu_solve { { i = Pnum[k]; rs = Rs[k]; - scale_div_assign(X[3*k], Bz[i], rs); - scale_div_assign(X[3*k + 1], Bz[i + d], rs); - scale_div_assign(X[3*k + 2], Bz[i + d*2], rs); + SCALE_DIV_ASSIGN (X[3*k], Bz[i], rs); + SCALE_DIV_ASSIGN (X[3*k + 1], Bz[i + d], rs); + SCALE_DIV_ASSIGN (X[3*k + 2], Bz[i + d*2], rs); } break; @@ -213,10 +217,10 @@ public class Dklu_solve { { i = Pnum[k]; rs = Rs[k]; - scale_div_assign(X[4*k], Bz[i], rs); - scale_div_assign(X[4*k + 1], Bz[i + d], rs); - scale_div_assign(X[4*k + 2], Bz[i + d*2], rs); - scale_div_assign(X[4*k + 3], Bz[i + d*3], rs); + SCALE_DIV_ASSIGN (X[4*k], Bz[i], rs); + SCALE_DIV_ASSIGN (X[4*k + 1], Bz[i + d], rs); + SCALE_DIV_ASSIGN (X[4*k + 2], Bz[i + d*2], rs); + SCALE_DIV_ASSIGN (X[4*k + 3], Bz[i + d*3], rs); } break; } @@ -237,7 +241,7 @@ public class Dklu_solve { k1 = R[block]; k2 = R[block+1]; nk = k2 - k1; - printf("solve %d, k1 %d k2-1 %d nk %d\n", block, k1, k2-1, nk); + PRINTF ("solve %d, k1 %d k2-1 %d nk %d\n", block, k1, k2-1, nk); /* solve the block system */ if (nk == 1) @@ -247,25 +251,25 @@ public class Dklu_solve { { case 1: - div(X[k1], X[k1], s); + DIV (X[k1], X[k1], s); break; case 2: - div(X[2*k1], X[2*k1], s); - div(X[2*k1 + 1], X[2*k1 + 1], s); + DIV (X[2*k1], X[2*k1], s); + DIV (X[2*k1 + 1], X[2*k1 + 1], s); break; case 3: - div(X[3*k1], X[3*k1], s); - div(X[3*k1 + 1], X[3*k1 + 1], s); - div(X[3*k1 + 2], X[3*k1 + 2], s); + DIV (X[3*k1], X[3*k1], s); + DIV (X[3*k1 + 1], X[3*k1 + 1], s); + DIV (X[3*k1 + 2], X[3*k1 + 2], s); break; case 4: - div(X[4*k1], X[4*k1], s); - div(X[4*k1 + 1], X[4*k1 + 1], s); - div(X[4*k1 + 2], X[4*k1 + 2], s); - div(X[4*k1 + 3], X[4*k1 + 3], s); + DIV (X[4*k1], X[4*k1], s); + DIV (X[4*k1 + 1], X[4*k1 + 1], s); + DIV (X[4*k1 + 2], X[4*k1 + 2], s); + DIV (X[4*k1 + 3], X[4*k1 + 3], s); break; } @@ -295,7 +299,7 @@ public class Dklu_solve { x[0] = X[k]; for (p = Offp[k]; p < pend; p++) { - mult_sub(X[Offi[p]], Offx[p], x[0]); + MULT_SUB (X[Offi[p]], Offx[p], x[0]); } } break; @@ -311,8 +315,8 @@ public class Dklu_solve { { i = Offi[p]; offik = Offx[p]; - mult_sub(X[2*i], offik, x[0]); - mult_sub(X[2*i + 1], offik, x[1]); + MULT_SUB (X[2*i], offik, x[0]); + MULT_SUB (X[2*i + 1], offik, x[1]); } } break; @@ -329,9 +333,9 @@ public class Dklu_solve { { i = Offi[p]; offik = Offx[p]; - mult_sub(X[3*i], offik, x[0]); - mult_sub(X[3*i + 1], offik, x[1]); - mult_sub(X[3*i + 2], offik, x[2]); + MULT_SUB (X[3*i], offik, x[0]); + MULT_SUB (X[3*i + 1], offik, x[1]); + MULT_SUB (X[3*i + 2], offik, x[2]); } } break; @@ -349,10 +353,10 @@ public class Dklu_solve { { i = Offi[p]; offik = Offx[p]; - mult_sub(X[4*i], offik, x[0]); - mult_sub(X[4*i + 1], offik, x[1]); - mult_sub(X[4*i + 2], offik, x[2]); - mult_sub(X[4*i + 3], offik, x[3]); + MULT_SUB (X[4*i], offik, x[0]); + MULT_SUB (X[4*i + 1], offik, x[1]); + MULT_SUB (X[4*i + 2], offik, x[2]); + MULT_SUB (X[4*i + 3], offik, x[3]); } } break; @@ -415,7 +419,7 @@ public class Dklu_solve { Bz += d*4; } - return true; + return TRUE; } } diff --git a/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_sort.java b/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_sort.java index db1f92c..6ef24a6 100644 --- a/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_sort.java +++ b/src/main/java/edu/ufl/cise/klu/tdouble/Dklu_sort.java @@ -24,12 +24,14 @@ package edu.ufl.cise.klu.tdouble; +import edu.ufl.cise.klu.common.KLU_common; + /** * Sorts the columns of L and U so that the row indices appear in strictly * increasing order. * */ -public class Dklu_sort { +public class Dklu_sort extends Dklu_internal { /** * Sort L or U using a double-transpose. @@ -41,7 +43,7 @@ public class Dklu_sort { Entry Xx; int p, i, j, len, nz, tp, xlen, pend; - assert Dklu_valid_LU.klu_valid_LU(n, false, Xip, Xlen, LU); + ASSERT (Dklu_valid_LU.klu_valid_LU(n, false, Xip, Xlen, LU)); /* count the number of entries in each row of L or U */ for (i = 0; i < n; i++) @@ -50,7 +52,7 @@ public class Dklu_sort { } for (j = 0; j < n; j++) { - get_pointer(LU, Xip, Xlen, Xi, Xx, j, len); + GET_POINTER (LU, Xip, Xlen, Xi, Xx, j, len); for (p = 0; p < len; p++) { W[Xi[p]]++; @@ -73,7 +75,7 @@ public class Dklu_sort { /* transpose the matrix into Tp, Ti, Tx */ for (j = 0; j < n; j++) { - get_pointer(LU, Xip, Xlen, Xi, Xx, j, len); + GET_POINTER (LU, Xip, Xlen, Xi, Xx, j, len); for (p = 0; p < len; p++) { tp = W[Xi[p]]++; @@ -93,14 +95,14 @@ public class Dklu_sort { for (p = Tp[i]; p < pend; p++) { j = Tj[p]; - get_pointer(LU, Xip, Xlen, Xi, Xx, j, len); + GET_POINTER (LU, Xip, Xlen, Xi, Xx, j, len); xlen = W[j]++; Xi[xlen] = i; Xx[xlen] = Tx[p]; } } - assert Dklu_valid_LU.klu_valid_LU(n, false, Xip, Xlen, LU); + ASSERT (Dklu_valid_LU.klu_valid_LU(n, false, Xip, Xlen, LU)); } @@ -117,7 +119,7 @@ public class Dklu_sort { { return false; } - Common.status = KLU_OK; + Common.status = KLU_common.KLU_OK; n = Symbolic.n; R = Symbolic.R; @@ -139,9 +141,9 @@ public class Dklu_sort { Ti = Dklu_malloc.klu_malloc(nz, (Int) sizeof, Common); Tx = Dklu_malloc.klu_malloc(nz, (Entry) sizeof, Common); - printf("\n======================= Start sort:\n"); + PRINTF ("\n======================= Start sort:\n"); - if (Common.status == KLU_OK) + if (Common.status == KLU_common.KLU_OK) { /* sort each block of L and U */ for (block = 0; block < nblocks; block++) @@ -150,14 +152,14 @@ public class Dklu_sort { nk = R[block+1] - k1; if (nk > 1) { - printf("\n-------------------block: %d nk %d\n", block, nk); + PRINTF ("\n-------------------block: %d nk %d\n", block, nk); sort(nk, Lip + k1, Llen + k1, LUbx[block], Tp, Ti, Tx, W); sort(nk, Uip + k1, Ulen + k1, LUbx[block], Tp, Ti, Tx, W); } } } - printf("\n======================= sort done.\n"); + PRINTF ("\n======================= sort done.\n"); /* free workspace */ Dklu_free.klu_free(W, maxblock, (Int) sizeof, Common); @@ -165,7 +167,7 @@ public class Dklu_sort { Dklu_free.klu_free(Ti, nz, (Int) sizeof, Common); Dklu_free.klu_free(Tx, nz, (Entry) sizeof, Common); - return Common.status == KLU_OK; + return Common.status == KLU_common.KLU_OK; } } 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 7b840b2..8839b5e 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 @@ -24,6 +24,8 @@ package edu.ufl.cise.klu.tdouble; +import edu.ufl.cise.klu.common.KLU_common; + /** * Solve A'x=b using the symbolic and numeric objects from KLU_analyze * (or KLU_analyze_given) and KLU_factor. Note that no iterative refinement is @@ -31,7 +33,7 @@ package edu.ufl.cise.klu.tdouble; * of size 4n Entry's (note that columns 2 to 4 of Xwork overlap with * Numeric.Iwork). */ -public class Dklu_tsolve { +public class Dklu_tsolve extends Dklu_internal { /** * @@ -43,7 +45,7 @@ public class Dklu_tsolve { * output. Size n*nrhs, in column-oriented form, with leading dimension d. * @return */ - public static boolean klu_tsolve(KLU_symbolic Symbolic, + public static int klu_tsolve(KLU_symbolic Symbolic, KLU_numeric Numeric, int d, int nrhs, double[] B, KLU_common Common) { @@ -62,15 +64,15 @@ public class Dklu_tsolve { if (Common == null) { - return (false); + return (FALSE); } if (Numeric == null || Symbolic == null || d < Symbolic.n || nrhs < 0 || B == null) { - Common.status = KLU_INVALID; - return false; + Common.status = KLU_common.KLU_INVALID; + return (FALSE); } - Common.status = KLU_OK; + Common.status = KLU_common.KLU_OK; /* ------------------------------------------------------------------ */ /* get the contents of the Symbolic object */ @@ -86,7 +88,7 @@ public class Dklu_tsolve { /* get the contents of the Numeric object */ /* ------------------------------------------------------------------ */ - assert (nblocks == Numeric.nblocks); + ASSERT (nblocks == Numeric.nblocks); Pnum = Numeric.Pnum; Offp = Numeric.Offp; Offi = Numeric.Offi; @@ -101,7 +103,7 @@ public class Dklu_tsolve { Rs = Numeric.Rs; X = (Entry) Numeric.Xwork; - assert Dklu_valid.klu_valid(n, Offp, Offi, Offx); + ASSERT (Dklu_valid.klu_valid(n, Offp, Offi, Offx)); /* ------------------------------------------------------------------ */ /* solve in chunks of 4 columns at a time */ @@ -114,7 +116,7 @@ public class Dklu_tsolve { /* get the size of the current chunk */ /* -------------------------------------------------------------- */ - nr = min(nrhs - chunk, 4); + nr = MIN (nrhs - chunk, 4); /* -------------------------------------------------------------- */ /* permute the right hand side, X = Q'*B */ @@ -180,7 +182,7 @@ public class Dklu_tsolve { k1 = R[block]; k2 = R[block+1]; nk = k2 - k1; - printf("tsolve %d, k1 %d k2-1 %d nk %d\n", block, k1, k2-1, nk); + PRINTF ("tsolve %d, k1 %d k2-1 %d nk %d\n", block, k1, k2-1, nk); /* ---------------------------------------------------------- */ /* block back-substitution for the off-diagonal-block entries */ @@ -198,7 +200,7 @@ public class Dklu_tsolve { pend = Offp[k+1]; for (p = Offp[k]; p < pend; p++) { - mult_sub(X[k], Offx[p], X[Offi[p]]); + MULT_SUB (X[k], Offx[p], X[Offi[p]]); } } break; @@ -214,8 +216,8 @@ public class Dklu_tsolve { { i = Offi[p]; offik = Offx[p]; - mult_sub(x[0], offik, X[2*i]); - mult_sub(x[1], offik, X[2*i + 1]); + MULT_SUB (x[0], offik, X[2*i]); + MULT_SUB (x[1], offik, X[2*i + 1]); } X[2*k ] = x[0]; X[2*k + 1] = x[1]; @@ -236,9 +238,9 @@ public class Dklu_tsolve { { offik = Offx[p]; } - mult_sub(x[0], offik, X[3*i]); - mult_sub(x[1], offik, X[3*i + 1]); - mult_sub(x[2], offik, X[3*i + 2]); + MULT_SUB (x[0], offik, X[3*i]); + MULT_SUB (x[1], offik, X[3*i + 1]); + MULT_SUB (x[2], offik, X[3*i + 2]); } X[3*k ] = x[0]; X[3*k + 1] = x[1]; @@ -259,10 +261,10 @@ public class Dklu_tsolve { { i = Offi[p]; offik = Offx[p]; - mult_sub(x[0], offik, X[4*i]); - mult_sub(x[1], offik, X[4*i + 1]); - mult_sub(x[2], offik, X[4*i + 2]); - mult_sub(x[3], offik, X[4*i + 3]); + MULT_SUB (x[0], offik, X[4*i]); + MULT_SUB (x[1], offik, X[4*i + 1]); + MULT_SUB (x[2], offik, X[4*i + 2]); + MULT_SUB (x[3], offik, X[4*i + 3]); } X[4*k ] = x[0]; X[4*k + 1] = x[1]; @@ -284,35 +286,37 @@ public class Dklu_tsolve { { case 1: - div(X[k1], X[k1], s); + DIV (X[k1], X[k1], s); break; case 2: - div(X[2*k1], X[2*k1], s); - div(X[2*k1 + 1], X[2*k1 + 1], s); + DIV (X[2*k1], X[2*k1], s); + DIV (X[2*k1 + 1], X[2*k1 + 1], s); break; case 3: - div(X[3*k1], X[3*k1], s); - div(X[3*k1 + 1], X[3*k1 + 1], s); - div(X[3*k1 + 2], X[3*k1 + 2], s); + DIV (X[3*k1], X[3*k1], s); + DIV (X[3*k1 + 1], X[3*k1 + 1], s); + DIV (X[3*k1 + 2], X[3*k1 + 2], s); break; case 4: - div(X[4*k1], X[4*k1], s); - div(X[4*k1 + 1], X[4*k1 + 1], s); - div(X[4*k1 + 2], X[4*k1 + 2], s); - div(X[4*k1 + 3], X[4*k1 + 3], s); + DIV (X[4*k1], X[4*k1], s); + DIV (X[4*k1 + 1], X[4*k1 + 1], s); + DIV (X[4*k1 + 2], X[4*k1 + 2], s); + DIV (X[4*k1 + 3], X[4*k1 + 3], s); break; } } else { - KLU_utsolve (nk, Uip + k1, Ulen + k1, LUbx[block], + Dklu_utsolve.klu_utsolve (nk, Uip + k1, Ulen + k1, + LUbx[block], Udiag + k1, nr, X + nr*k1); - KLU_ltsolve (nk, Lip + k1, Llen + k1, LUbx[block], nr, + Dklu_ltsolve.klu_ltsolve (nk, Lip + k1, Llen + k1, + LUbx[block], nr, X + nr*k1); } } @@ -381,7 +385,7 @@ public class Dklu_tsolve { for (k = 0; k < n; k++) { - scale_div_assign(Bz[Pnum[k]], X[k], Rs[k]); + SCALE_DIV_ASSIGN (Bz[Pnum[k]], X[k], Rs[k]); } break; @@ -391,8 +395,8 @@ public class Dklu_tsolve { { i = Pnum[k]; rs = Rs[k]; - scale_div_assign(Bz[i], X[2*k], rs); - scale_div_assign(Bz[i + d], X[2*k + 1], rs); + SCALE_DIV_ASSIGN (Bz[i], X[2*k], rs); + SCALE_DIV_ASSIGN (Bz[i + d], X[2*k + 1], rs); } break; @@ -402,9 +406,9 @@ public class Dklu_tsolve { { i = Pnum[k]; rs = Rs[k]; - scale_div_assign(Bz[i], X[3*k], rs); - scale_div_assign(Bz[i + d], X[3*k + 1], rs); - scale_div_assign(Bz[i + d*2], X[3*k + 2], rs); + SCALE_DIV_ASSIGN (Bz[i], X[3*k], rs); + SCALE_DIV_ASSIGN (Bz[i + d], X[3*k + 1], rs); + SCALE_DIV_ASSIGN (Bz[i + d*2], X[3*k + 2], rs); } break; @@ -414,10 +418,10 @@ public class Dklu_tsolve { { i = Pnum[k]; rs = Rs[k]; - scale_div_assign(Bz[i], X[4*k], rs); - scale_div_assign(Bz[i + d], X[4*k + 1], rs); - scale_div_assign(Bz[i + d*2], X[4*k + 2], rs); - scale_div_assign(Bz[i + d*3], X[4*k + 3], rs); + SCALE_DIV_ASSIGN (Bz[i], X[4*k], rs); + SCALE_DIV_ASSIGN (Bz[i + d], X[4*k + 1], rs); + SCALE_DIV_ASSIGN (Bz[i + d*2], X[4*k + 2], rs); + SCALE_DIV_ASSIGN (Bz[i + d*3], X[4*k + 3], rs); } break; } @@ -429,7 +433,7 @@ public class Dklu_tsolve { Bz += d*4; } - return true; + 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 d649650..765b1c6 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 @@ -27,43 +27,43 @@ package edu.ufl.cise.klu.tdouble; public class Dklu_version { /** enable diagnostic printing */ - static boolean NPRINT = true; + protected static boolean NPRINT = true ; - static final int INT_MAX = 0x7fffffff; + protected static final int INT_MAX = 0x7fffffff ; - static final String INT_ID = "%d"; + protected static final String INT_ID = "%d" ; - static int BYTES (Object type, double n) + protected static int BYTES (Object type, double n) { - return sizeof (type * n); + return sizeof (type * n) ; } - static double CEILING (double b, double u) + protected static double CEILING (double b, double u) { - return (b+u-1) / u; + return (b+u-1) / u ; } - static double UNITS (Object type, double n) + protected static double UNITS (Object type, double n) { - return CEILING (BYTES (type, n), sizeof (Unit)); + return CEILING (BYTES (type, n), sizeof (Unit)) ; } - static double DUNITS (Object type, int n) + protected static double DUNITS (Object type, int n) { - return Math.ceil(BYTES (type, (double) n) / sizeof (Unit)); + return Math.ceil(BYTES (type, (double) n) / sizeof (Unit)) ; } - static void GET_I_POINTER(LU, Xip, Xi, k) + protected static void GET_I_POINTER(LU, Xip, Xi, k) { Xi = (Int[]) (LU + Xip [k]) ; } - static void GET_X_POINTER(LU, Xip, Xlen, Xx, k) + protected static void GET_X_POINTER(LU, Xip, Xlen, Xx, k) { Xx = (Entry[]) (LU + Xip [k] + UNITS (Int, Xlen [k])) ; } - static void GET_POINTER(LU, Xip, Xlen, Xi, Xx, k, xlen) + protected static void GET_POINTER(LU, Xip, Xlen, Xi, Xx, k, xlen) { Unit xp = LU + Xip [k] ; xlen = Xlen [k] ; @@ -71,38 +71,38 @@ public class Dklu_version { Xx = (Entry[]) (xp + UNITS (Int, xlen)) ; } - static boolean SCALAR_IS_NAN (double x) + protected static boolean SCALAR_IS_NAN (double x) { - return x != x; + return x != x ; } - static boolean SCALAR_IS_ZERO (double x) + protected static boolean SCALAR_IS_ZERO (double x) { - return x == 0.0; + return x == 0.0 ; } - static boolean SCALAR_IS_NONZERO (double x) + protected static boolean SCALAR_IS_NONZERO (double x) { - return x != 0.0; + return x != 0.0 ; } - static boolean SCALAR_IS_LTZERO (double x) + protected static boolean SCALAR_IS_LTZERO (double x) { - return x < 0.0; + return x < 0.0 ; } /* scalar absolute value macro. If x is NaN, the result is NaN */ - static double SCALAR_ABS (double x) + protected static double SCALAR_ABS (double x) { - return SCALAR_IS_LTZERO (x) ? -x : x; + return SCALAR_IS_LTZERO (x) ? -x : x ; } - static void PRINTF (String format, Object... args) + protected static void PRINTF (String format, Object... args) { - System.out.printf(format, args); + System.out.printf(format, args) ; } - static double PRINT_SCALAR (double a) + protected static double PRINT_SCALAR (double a) { if (!NPRINT) { @@ -122,15 +122,15 @@ public class Dklu_version { /* Real floating-point arithmetic */ /* ---------------------------------------------------------------------- */ - static final Class Unit = Double.class; + protected static final Class Unit = Double.class ; - static final Class Entry = Double.class; + protected static final Class Entry = Double.class ; /** * @return TRUE if a complex number is in split form, FALSE if in packed * form. */ - static int SPLIT (double s) + protected static int SPLIT (double s) { return 1 ; } @@ -138,7 +138,7 @@ public class Dklu_version { /** * @return real part of c */ - static double REAL (double c) + protected static double REAL (double c) { return c ; } @@ -146,7 +146,7 @@ public class Dklu_version { /** * @return imag part of c */ - static double IMAG (double c) + protected static double IMAG (double c) { return 0.0 ; } @@ -154,27 +154,27 @@ public class Dklu_version { /** * c = (s1) + (s2)*i */ - static void ASSIGN (Double c, double[] s1, double[] s2, int p, + protected static void ASSIGN (Double c, double[] s1, double[] s2, int p, boolean split) { c = s1[p] ; } - static void CLEAR (Double c) + protected static void CLEAR (Double c) { c = 0.0 ; } - static void CLEAR_AND_INCREMENT (Double p) + protected static void CLEAR_AND_INCREMENT (Double p) { p = 0.0 ; - p++; + p++ ; } /** * @return True if a is NaN */ - static boolean IS_NAN (double a) + protected static boolean IS_NAN (double a) { return SCALAR_IS_NAN (a) ; } @@ -182,7 +182,7 @@ public class Dklu_version { /** * @return True if a == 0 */ - static boolean IS_ZERO (double a) + protected static boolean IS_ZERO (double a) { return SCALAR_IS_ZERO (a) ; } @@ -190,7 +190,7 @@ public class Dklu_version { /** * @return True if a != 0 */ - static boolean IS_NONZERO (double a) + protected static boolean IS_NONZERO (double a) { return SCALAR_IS_NONZERO (a) ; } @@ -198,7 +198,7 @@ public class Dklu_version { /** * c /= s */ - static void SCALE_DIV (Double c, double s) + protected static void SCALE_DIV (Double c, double s) { c /= s ; } @@ -206,7 +206,7 @@ public class Dklu_version { /** * a = c/s */ - static void SCALE_DIV_ASSIGN (Double a, double c, double s) + protected static void SCALE_DIV_ASSIGN (Double a, double c, double s) { a = c / s ; } @@ -214,7 +214,7 @@ public class Dklu_version { /** * c *= s */ - static void SCALE (Double c, double s) + protected static void SCALE (Double c, double s) { c *= s ; } @@ -222,7 +222,7 @@ public class Dklu_version { /** * c += a */ - static void ASSEMBLE (Double c, double a) + protected static void ASSEMBLE (Double c, double a) { c += a ; } @@ -230,7 +230,7 @@ public class Dklu_version { /** * c += *p++ */ - static void ASSEMBLE_AND_INCREMENT (Double c, double p) + protected static void ASSEMBLE_AND_INCREMENT (Double c, double p) { c += p++ ; } @@ -238,7 +238,7 @@ public class Dklu_version { /** * c -= a */ - static void DECREMENT (Double c, double a) + protected static void DECREMENT (Double c, double a) { c -= a ; } @@ -246,7 +246,7 @@ public class Dklu_version { /** * c = a*b */ - static void MULT (Double c, double a, double b) + protected static void MULT (Double c, double a, double b) { c = a * b ; } @@ -254,7 +254,7 @@ public class Dklu_version { /** * c = a*conjugate(b) */ - static void MULT_CONJ (Double c, double a, double b) + protected static void MULT_CONJ (Double c, double a, double b) { c = a * b ; } @@ -262,7 +262,7 @@ public class Dklu_version { /** * c -= a*b */ - static void MULT_SUB (Double c, double a, double b) + protected static void MULT_SUB (Double c, double a, double b) { c -= a * b ; } @@ -270,7 +270,7 @@ public class Dklu_version { /** * c -= a*conjugate(b) */ - static void MULT_SUB_CONJ (Double c, double a, double b) + protected static void MULT_SUB_CONJ (Double c, double a, double b) { c -= a * b ; } @@ -278,7 +278,7 @@ public class Dklu_version { /** * c = a/b */ - static void DIV (Double c, double a, double b) + protected static void DIV (Double c, double a, double b) { c = a / b ; } @@ -286,7 +286,7 @@ public class Dklu_version { /** * c = 1/c */ - static void RECIPROCAL (Double c) + protected static void RECIPROCAL (Double c) { c = 1.0 / c ; } @@ -294,7 +294,7 @@ public class Dklu_version { /** * c = a/conjugate(b) */ - static void DIV_CONJ (Double c, double a, double b) + protected static void DIV_CONJ (Double c, double a, double b) { c = a / b ; } @@ -302,7 +302,7 @@ public class Dklu_version { /** * approximate absolute value, s = |r|+|i| */ - static void APPROX_ABS (Double s, double a) + protected static void APPROX_ABS (Double s, double a) { s = SCALAR_ABS (a) ; } @@ -310,28 +310,28 @@ public class Dklu_version { /** * exact absolute value, s = sqrt (a.real^2 + amag^2) */ - static void ABS (Double s, double a) + protected static void ABS (Double s, double a) { s = SCALAR_ABS (a) ; } - static void PRINT_ENTRY (double a) + protected static void PRINT_ENTRY (double a) { PRINT_SCALAR (a) ; } - static void CONJ (Double a, double x) + protected static void CONJ (Double a, double x) { - a = x; + a = x ; } /* for flop counts */ - static final double MULTSUB_FLOPS = 2.0 ; /* c -= a*b */ - static final double DIV_FLOPS = 1.0 ; /* c = a/b */ - static final double ABS_FLOPS = 0.0 ; /* c = abs (a) */ - static final double ASSEMBLE_FLOPS = 1.0 ; /* c += a */ - static final double DECREMENT_FLOPS = 1.0 ; /* c -= a */ - static final double MULT_FLOPS = 1.0 ; /* c = a*b */ - static final double SCALE_FLOPS = 1.0 ; /* c = a/s */ + protected static final double MULTSUB_FLOPS = 2.0 ; /* c -= a*b */ + protected static final double DIV_FLOPS = 1.0 ; /* c = a/b */ + protected static final double ABS_FLOPS = 0.0 ; /* c = abs(a) */ + protected static final double ASSEMBLE_FLOPS = 1.0 ; /* c += a */ + protected static final double DECREMENT_FLOPS = 1.0 ; /* c -= a */ + protected static final double MULT_FLOPS = 1.0 ; /* c = a*b */ + protected static final double SCALE_FLOPS = 1.0 ; /* c = a/s */ } -- 2.11.4.GIT