/* 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: . * * 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 . ***********************************************************************/ #include "glpenv.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; } /*********************************************************************** * 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 */