aboutsummaryrefslogtreecommitdiffstats
path: root/test/monniaux/glpk-4.65/src/simplex/spxnt.c
diff options
context:
space:
mode:
Diffstat (limited to 'test/monniaux/glpk-4.65/src/simplex/spxnt.c')
-rw-r--r--test/monniaux/glpk-4.65/src/simplex/spxnt.c303
1 files changed, 303 insertions, 0 deletions
diff --git a/test/monniaux/glpk-4.65/src/simplex/spxnt.c b/test/monniaux/glpk-4.65/src/simplex/spxnt.c
new file mode 100644
index 00000000..7eaac852
--- /dev/null
+++ b/test/monniaux/glpk-4.65/src/simplex/spxnt.c
@@ -0,0 +1,303 @@
+/* spxnt.c */
+
+/***********************************************************************
+* This code is part of GLPK (GNU Linear Programming Kit).
+*
+* Copyright (C) 2015 Andrew Makhorin, Department for Applied
+* Informatics, Moscow Aviation Institute, Moscow, Russia. All rights
+* reserved. E-mail: <mao@gnu.org>.
+*
+* GLPK is free software: you can redistribute it and/or modify it
+* under the terms of the GNU General Public License as published by
+* the Free Software Foundation, either version 3 of the License, or
+* (at your option) any later version.
+*
+* GLPK is distributed in the hope that it will be useful, but WITHOUT
+* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
+* License for more details.
+*
+* You should have received a copy of the GNU General Public License
+* along with GLPK. If not, see <http://www.gnu.org/licenses/>.
+***********************************************************************/
+
+#include "env.h"
+#include "spxnt.h"
+
+/***********************************************************************
+* spx_alloc_nt - allocate matrix N in sparse row-wise format
+*
+* This routine allocates the memory for arrays needed to represent the
+* matrix N composed of non-basic columns of the constraint matrix A. */
+
+void spx_alloc_nt(SPXLP *lp, SPXNT *nt)
+{ int m = lp->m;
+ int nnz = lp->nnz;
+ nt->ptr = talloc(1+m, int);
+ nt->len = talloc(1+m, int);
+ nt->ind = talloc(1+nnz, int);
+ nt->val = talloc(1+nnz, double);
+ return;
+}
+
+/***********************************************************************
+* spx_init_nt - initialize row pointers for matrix N
+*
+* This routine initializes (sets up) row pointers for the matrix N
+* using column-wise representation of the constraint matrix A.
+*
+* This routine needs to be called only once. */
+
+void spx_init_nt(SPXLP *lp, SPXNT *nt)
+{ int m = lp->m;
+ int n = lp->n;
+ int nnz = lp->nnz;
+ int *A_ptr = lp->A_ptr;
+ int *A_ind = lp->A_ind;
+ int *NT_ptr = nt->ptr;
+ int *NT_len = nt->len;
+ int i, k, ptr, end;
+ /* calculate NT_len[i] = maximal number of non-zeros in i-th row
+ * of N = number of non-zeros in i-th row of A */
+ memset(&NT_len[1], 0, m * sizeof(int));
+ for (k = 1; k <= n; k++)
+ { ptr = A_ptr[k];
+ end = A_ptr[k+1];
+ for (; ptr < end; ptr++)
+ NT_len[A_ind[ptr]]++;
+ }
+ /* initialize row pointers NT_ptr[i], i = 1,...,n-m */
+ NT_ptr[1] = 1;
+ for (i = 2; i <= m; i++)
+ NT_ptr[i] = NT_ptr[i-1] + NT_len[i-1];
+ xassert(NT_ptr[m] + NT_len[m] == nnz+1);
+ return;
+}
+
+/***********************************************************************
+* spx_nt_add_col - add column N[j] = A[k] to matrix N
+*
+* This routine adds elements of column N[j] = A[k], 1 <= j <= n-m,
+* 1 <= k <= n, to the row-wise represntation of the matrix N. It is
+* assumed (with no check) that elements of the specified column are
+* missing in the row-wise represntation of N. */
+
+void spx_nt_add_col(SPXLP *lp, SPXNT *nt, int j, int k)
+{ int m = lp->m;
+ int n = lp->n;
+ int nnz = lp->nnz;
+ int *A_ptr = lp->A_ptr;
+ int *A_ind = lp->A_ind;
+ double *A_val = lp->A_val;
+ int *NT_ptr = nt->ptr;
+ int *NT_len = nt->len;
+ int *NT_ind = nt->ind;
+ double *NT_val = nt->val;
+ int i, ptr, end, pos;
+ xassert(1 <= j && j <= n-m);
+ xassert(1 <= k && k <= n);
+ ptr = A_ptr[k];
+ end = A_ptr[k+1];
+ for (; ptr < end; ptr++)
+ { i = A_ind[ptr];
+ /* add element N[i,j] = A[i,k] to i-th row of matrix N */
+ pos = NT_ptr[i] + (NT_len[i]++);
+ if (i < m)
+ xassert(pos < NT_ptr[i+1]);
+ else
+ xassert(pos <= nnz);
+ NT_ind[pos] = j;
+ NT_val[pos] = A_val[ptr];
+ }
+ return;
+}
+
+/***********************************************************************
+* spx_build_nt - build matrix N for current basis
+*
+* This routine builds the row-wise represntation of the matrix N
+* for the current basis by adding columns of the constraint matrix A
+* corresponding to non-basic variables. */
+
+void spx_build_nt(SPXLP *lp, SPXNT *nt)
+{ int m = lp->m;
+ int n = lp->n;
+ int *head = lp->head;
+ int *NT_len = nt->len;
+ int j, k;
+ /* N := 0 */
+ memset(&NT_len[1], 0, m * sizeof(int));
+ /* add non-basic columns N[j] = A[k] */
+ for (j = 1; j <= n-m; j++)
+ { k = head[m+j]; /* x[k] = xN[j] */
+ spx_nt_add_col(lp, nt, j, k);
+ }
+ return;
+}
+
+/***********************************************************************
+* spx_nt_del_col - remove column N[j] = A[k] from matrix N
+*
+* This routine removes elements of column N[j] = A[k], 1 <= j <= n-m,
+* 1 <= k <= n, from the row-wise representation of the matrix N. It is
+* assumed (with no check) that elements of the specified column are
+* present in the row-wise representation of N. */
+
+void spx_nt_del_col(SPXLP *lp, SPXNT *nt, int j, int k)
+{ int m = lp->m;
+ int n = lp->n;
+ int *A_ptr = lp->A_ptr;
+ int *A_ind = lp->A_ind;
+ int *NT_ptr = nt->ptr;
+ int *NT_len = nt->len;
+ int *NT_ind = nt->ind;
+ double *NT_val = nt->val;
+ int i, ptr, end, ptr1, end1;
+ xassert(1 <= j && j <= n-m);
+ xassert(1 <= k && k <= n);
+ ptr = A_ptr[k];
+ end = A_ptr[k+1];
+ for (; ptr < end; ptr++)
+ { i = A_ind[ptr];
+ /* find element N[i,j] = A[i,k] in i-th row of matrix N */
+ ptr1 = NT_ptr[i];
+ end1 = ptr1 + NT_len[i];
+ for (; NT_ind[ptr1] != j; ptr1++)
+ /* nop */;
+ xassert(ptr1 < end1);
+ /* and remove it from i-th row element list */
+ NT_len[i]--;
+ NT_ind[ptr1] = NT_ind[end1-1];
+ NT_val[ptr1] = NT_val[end1-1];
+ }
+ return;
+}
+
+/***********************************************************************
+* spx_update_nt - update matrix N for adjacent basis
+*
+* This routine updates the row-wise represntation of matrix N for
+* the adjacent basis, where column N[q], 1 <= q <= n-m, is replaced by
+* column B[p], 1 <= p <= m, of the current basis matrix B. */
+
+void spx_update_nt(SPXLP *lp, SPXNT *nt, int p, int q)
+{ int m = lp->m;
+ int n = lp->n;
+ int *head = lp->head;
+ xassert(1 <= p && p <= m);
+ xassert(1 <= q && q <= n-m);
+ /* remove old column N[q] corresponding to variable xN[q] */
+ spx_nt_del_col(lp, nt, q, head[m+q]);
+ /* add new column N[q] corresponding to variable xB[p] */
+ spx_nt_add_col(lp, nt, q, head[p]);
+ return;
+}
+
+/***********************************************************************
+* spx_nt_prod - compute product y := y + s * N'* x
+*
+* This routine computes the product:
+*
+* y := y + s * N'* x,
+*
+* where N' is a matrix transposed to the mx(n-m)-matrix N composed
+* from non-basic columns of the constraint matrix A, x is a m-vector,
+* s is a scalar, y is (n-m)-vector.
+*
+* If the flag ign is non-zero, the routine ignores the input content
+* of the array y assuming that y = 0.
+*
+* The routine uses the row-wise representation of the matrix N and
+* computes the product as a linear combination:
+*
+* y := y + s * (N'[1] * x[1] + ... + N'[m] * x[m]),
+*
+* where N'[i] is i-th row of N, 1 <= i <= m. */
+
+void spx_nt_prod(SPXLP *lp, SPXNT *nt, double y[/*1+n-m*/], int ign,
+ double s, const double x[/*1+m*/])
+{ int m = lp->m;
+ int n = lp->n;
+ int *NT_ptr = nt->ptr;
+ int *NT_len = nt->len;
+ int *NT_ind = nt->ind;
+ double *NT_val = nt->val;
+ int i, j, ptr, end;
+ double t;
+ if (ign)
+ { /* y := 0 */
+ for (j = 1; j <= n-m; j++)
+ y[j] = 0.0;
+ }
+ for (i = 1; i <= m; i++)
+ { if (x[i] != 0.0)
+ { /* y := y + s * (i-th row of N) * x[i] */
+ t = s * x[i];
+ ptr = NT_ptr[i];
+ end = ptr + NT_len[i];
+ for (; ptr < end; ptr++)
+ y[NT_ind[ptr]] += NT_val[ptr] * t;
+ }
+ }
+ return;
+}
+
+#if 1 /* 31/III-2016 */
+void spx_nt_prod_s(SPXLP *lp, SPXNT *nt, FVS *y, int ign, double s,
+ const FVS *x, double eps)
+{ /* sparse version of spx_nt_prod */
+ int *NT_ptr = nt->ptr;
+ int *NT_len = nt->len;
+ int *NT_ind = nt->ind;
+ double *NT_val = nt->val;
+ int *x_ind = x->ind;
+ double *x_vec = x->vec;
+ int *y_ind = y->ind;
+ double *y_vec = y->vec;
+ int i, j, k, nnz, ptr, end;
+ double t;
+ xassert(x->n == lp->m);
+ xassert(y->n == lp->n-lp->m);
+ if (ign)
+ { /* y := 0 */
+ fvs_clear_vec(y);
+ }
+ nnz = y->nnz;
+ for (k = x->nnz; k >= 1; k--)
+ { i = x_ind[k];
+ /* y := y + s * (i-th row of N) * x[i] */
+ t = s * x_vec[i];
+ ptr = NT_ptr[i];
+ end = ptr + NT_len[i];
+ for (; ptr < end; ptr++)
+ { j = NT_ind[ptr];
+ if (y_vec[j] == 0.0)
+ y_ind[++nnz] = j;
+ y_vec[j] += NT_val[ptr] * t;
+ /* don't forget about numeric cancellation */
+ if (y_vec[j] == 0.0)
+ y_vec[j] = DBL_MIN;
+ }
+ }
+ y->nnz = nnz;
+ fvs_adjust_vec(y, eps);
+ return;
+}
+#endif
+
+/***********************************************************************
+* spx_free_nt - deallocate matrix N in sparse row-wise format
+*
+* This routine deallocates the memory used for arrays of the program
+* object nt. */
+
+void spx_free_nt(SPXLP *lp, SPXNT *nt)
+{ xassert(lp == lp);
+ tfree(nt->ptr);
+ tfree(nt->len);
+ tfree(nt->ind);
+ tfree(nt->val);
+ return;
+}
+
+/* eof */
generated by cgit (git 2.34.1) at 2024-06-02 19:26:05 +0000