alpar@9: /* glpini01.c */ alpar@9: alpar@9: /*********************************************************************** alpar@9: * This code is part of GLPK (GNU Linear Programming Kit). alpar@9: * alpar@9: * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, alpar@9: * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics, alpar@9: * Moscow Aviation Institute, Moscow, Russia. All rights reserved. alpar@9: * E-mail: . alpar@9: * alpar@9: * GLPK is free software: you can redistribute it and/or modify it alpar@9: * under the terms of the GNU General Public License as published by alpar@9: * the Free Software Foundation, either version 3 of the License, or alpar@9: * (at your option) any later version. alpar@9: * alpar@9: * GLPK is distributed in the hope that it will be useful, but WITHOUT alpar@9: * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY alpar@9: * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public alpar@9: * License for more details. alpar@9: * alpar@9: * You should have received a copy of the GNU General Public License alpar@9: * along with GLPK. If not, see . alpar@9: ***********************************************************************/ alpar@9: alpar@9: #include "glpapi.h" alpar@9: alpar@9: /*---------------------------------------------------------------------- alpar@9: -- triang - find maximal triangular part of a rectangular matrix. alpar@9: -- alpar@9: -- *Synopsis* alpar@9: -- alpar@9: -- int triang(int m, int n, alpar@9: -- void *info, int (*mat)(void *info, int k, int ndx[]), alpar@9: -- int rn[], int cn[]); alpar@9: -- alpar@9: -- *Description* alpar@9: -- alpar@9: -- For a given rectangular (sparse) matrix A with m rows and n columns alpar@9: -- the routine triang tries to find such permutation matrices P and Q alpar@9: -- that the first rows and columns of the matrix B = P*A*Q form a lower alpar@9: -- triangular submatrix of as greatest size as possible: alpar@9: -- alpar@9: -- 1 n alpar@9: -- 1 * . . . . . . x x x x x x alpar@9: -- * * . . . . . x x x x x x alpar@9: -- * * * . . . . x x x x x x alpar@9: -- * * * * . . . x x x x x x alpar@9: -- B = P*A*Q = * * * * * . . x x x x x x alpar@9: -- * * * * * * . x x x x x x alpar@9: -- * * * * * * * x x x x x x alpar@9: -- x x x x x x x x x x x x x alpar@9: -- x x x x x x x x x x x x x alpar@9: -- m x x x x x x x x x x x x x alpar@9: -- alpar@9: -- where: '*' - elements of the lower triangular part, '.' - structural alpar@9: -- zeros, 'x' - other (either non-zero or zero) elements. alpar@9: -- alpar@9: -- The parameter info is a transit pointer passed to the formal routine alpar@9: -- mat (see below). alpar@9: -- alpar@9: -- The formal routine mat specifies the given matrix A in both row- and alpar@9: -- column-wise formats. In order to obtain an i-th row of the matrix A alpar@9: -- the routine triang calls the routine mat with the parameter k = +i, alpar@9: -- 1 <= i <= m. In response the routine mat should store column indices alpar@9: -- of (non-zero) elements of the i-th row to the locations ndx[1], ..., alpar@9: -- ndx[len], where len is number of non-zeros in the i-th row returned alpar@9: -- on exit. Analogously, in order to obtain a j-th column of the matrix alpar@9: -- A, the routine mat is called with the parameter k = -j, 1 <= j <= n, alpar@9: -- and should return pattern of the j-th column in the same way as for alpar@9: -- row patterns. Note that the routine mat may be called more than once alpar@9: -- for the same rows and columns. alpar@9: -- alpar@9: -- On exit the routine computes two resultant arrays rn and cn, which alpar@9: -- define the permutation matrices P and Q, respectively. The array rn alpar@9: -- should have at least 1+m locations, where rn[i] = i' (1 <= i <= m) alpar@9: -- means that i-th row of the original matrix A corresponds to i'-th row alpar@9: -- of the matrix B = P*A*Q. Similarly, the array cn should have at least alpar@9: -- 1+n locations, where cn[j] = j' (1 <= j <= n) means that j-th column alpar@9: -- of the matrix A corresponds to j'-th column of the matrix B. alpar@9: -- alpar@9: -- *Returns* alpar@9: -- alpar@9: -- The routine triang returns the size of the lower tringular part of alpar@9: -- the matrix B = P*A*Q (see the figure above). alpar@9: -- alpar@9: -- *Complexity* alpar@9: -- alpar@9: -- The time complexity of the routine triang is O(nnz), where nnz is alpar@9: -- number of non-zeros in the given matrix A. alpar@9: -- alpar@9: -- *Algorithm* alpar@9: -- alpar@9: -- The routine triang starts from the matrix B = P*Q*A, where P and Q alpar@9: -- are unity matrices, so initially B = A. alpar@9: -- alpar@9: -- Before the next iteration B = (B1 | B2 | B3), where B1 is partially alpar@9: -- built a lower triangular submatrix, B2 is the active submatrix, and alpar@9: -- B3 is a submatrix that contains rejected columns. Thus, the current alpar@9: -- matrix B looks like follows (initially k1 = 1 and k2 = n): alpar@9: -- alpar@9: -- 1 k1 k2 n alpar@9: -- 1 x . . . . . . . . . . . . . # # # alpar@9: -- x x . . . . . . . . . . . . # # # alpar@9: -- x x x . . . . . . . . . . # # # # alpar@9: -- x x x x . . . . . . . . . # # # # alpar@9: -- x x x x x . . . . . . . # # # # # alpar@9: -- k1 x x x x x * * * * * * * # # # # # alpar@9: -- x x x x x * * * * * * * # # # # # alpar@9: -- x x x x x * * * * * * * # # # # # alpar@9: -- x x x x x * * * * * * * # # # # # alpar@9: -- m x x x x x * * * * * * * # # # # # alpar@9: -- <--B1---> <----B2-----> <---B3--> alpar@9: -- alpar@9: -- On each iteartion the routine looks for a singleton row, i.e. some alpar@9: -- row that has the only non-zero in the active submatrix B2. If such alpar@9: -- row exists and the corresponding non-zero is b[i,j], where (by the alpar@9: -- definition) k1 <= i <= m and k1 <= j <= k2, the routine permutes alpar@9: -- k1-th and i-th rows and k1-th and j-th columns of the matrix B (in alpar@9: -- order to place the element in the position b[k1,k1]), removes the alpar@9: -- k1-th column from the active submatrix B2, and adds this column to alpar@9: -- the submatrix B1. If no row singletons exist, but B2 is not empty alpar@9: -- yet, the routine chooses a j-th column, which has maximal number of alpar@9: -- non-zeros among other columns of B2, removes this column from B2 and alpar@9: -- adds it to the submatrix B3 in the hope that new row singletons will alpar@9: -- appear in the active submatrix. */ alpar@9: alpar@9: static int triang(int m, int n, alpar@9: void *info, int (*mat)(void *info, int k, int ndx[]), alpar@9: int rn[], int cn[]) alpar@9: { int *ndx; /* int ndx[1+max(m,n)]; */ alpar@9: /* this array is used for querying row and column patterns of the alpar@9: given matrix A (the third parameter to the routine mat) */ alpar@9: int *rs_len; /* int rs_len[1+m]; */ alpar@9: /* rs_len[0] is not used; alpar@9: rs_len[i], 1 <= i <= m, is number of non-zeros in the i-th row alpar@9: of the matrix A, which (non-zeros) belong to the current active alpar@9: submatrix */ alpar@9: int *rs_head; /* int rs_head[1+n]; */ alpar@9: /* rs_head[len], 0 <= len <= n, is the number i of the first row alpar@9: of the matrix A, for which rs_len[i] = len */ alpar@9: int *rs_prev; /* int rs_prev[1+m]; */ alpar@9: /* rs_prev[0] is not used; alpar@9: rs_prev[i], 1 <= i <= m, is a number i' of the previous row of alpar@9: the matrix A, for which rs_len[i] = rs_len[i'] (zero marks the alpar@9: end of this linked list) */ alpar@9: int *rs_next; /* int rs_next[1+m]; */ alpar@9: /* rs_next[0] is not used; alpar@9: rs_next[i], 1 <= i <= m, is a number i' of the next row of the alpar@9: matrix A, for which rs_len[i] = rs_len[i'] (zero marks the end alpar@9: this linked list) */ alpar@9: int cs_head; alpar@9: /* is a number j of the first column of the matrix A, which has alpar@9: maximal number of non-zeros among other columns */ alpar@9: int *cs_prev; /* cs_prev[1+n]; */ alpar@9: /* cs_prev[0] is not used; alpar@9: cs_prev[j], 1 <= j <= n, is a number of the previous column of alpar@9: the matrix A with the same or greater number of non-zeros than alpar@9: in the j-th column (zero marks the end of this linked list) */ alpar@9: int *cs_next; /* cs_next[1+n]; */ alpar@9: /* cs_next[0] is not used; alpar@9: cs_next[j], 1 <= j <= n, is a number of the next column of alpar@9: the matrix A with the same or lesser number of non-zeros than alpar@9: in the j-th column (zero marks the end of this linked list) */ alpar@9: int i, j, ii, jj, k1, k2, len, t, size = 0; alpar@9: int *head, *rn_inv, *cn_inv; alpar@9: if (!(m > 0 && n > 0)) alpar@9: xerror("triang: m = %d; n = %d; invalid dimension\n", m, n); alpar@9: /* allocate working arrays */ alpar@9: ndx = xcalloc(1+(m >= n ? m : n), sizeof(int)); alpar@9: rs_len = xcalloc(1+m, sizeof(int)); alpar@9: rs_head = xcalloc(1+n, sizeof(int)); alpar@9: rs_prev = xcalloc(1+m, sizeof(int)); alpar@9: rs_next = xcalloc(1+m, sizeof(int)); alpar@9: cs_prev = xcalloc(1+n, sizeof(int)); alpar@9: cs_next = xcalloc(1+n, sizeof(int)); alpar@9: /* build linked lists of columns of the matrix A with the same alpar@9: number of non-zeros */ alpar@9: head = rs_len; /* currently rs_len is used as working array */ alpar@9: for (len = 0; len <= m; len ++) head[len] = 0; alpar@9: for (j = 1; j <= n; j++) alpar@9: { /* obtain length of the j-th column */ alpar@9: len = mat(info, -j, ndx); alpar@9: xassert(0 <= len && len <= m); alpar@9: /* include the j-th column in the corresponding linked list */ alpar@9: cs_prev[j] = head[len]; alpar@9: head[len] = j; alpar@9: } alpar@9: /* merge all linked lists of columns in one linked list, where alpar@9: columns are ordered by descending of their lengths */ alpar@9: cs_head = 0; alpar@9: for (len = 0; len <= m; len++) alpar@9: { for (j = head[len]; j != 0; j = cs_prev[j]) alpar@9: { cs_next[j] = cs_head; alpar@9: cs_head = j; alpar@9: } alpar@9: } alpar@9: jj = 0; alpar@9: for (j = cs_head; j != 0; j = cs_next[j]) alpar@9: { cs_prev[j] = jj; alpar@9: jj = j; alpar@9: } alpar@9: /* build initial doubly linked lists of rows of the matrix A with alpar@9: the same number of non-zeros */ alpar@9: for (len = 0; len <= n; len++) rs_head[len] = 0; alpar@9: for (i = 1; i <= m; i++) alpar@9: { /* obtain length of the i-th row */ alpar@9: rs_len[i] = len = mat(info, +i, ndx); alpar@9: xassert(0 <= len && len <= n); alpar@9: /* include the i-th row in the correspondng linked list */ alpar@9: rs_prev[i] = 0; alpar@9: rs_next[i] = rs_head[len]; alpar@9: if (rs_next[i] != 0) rs_prev[rs_next[i]] = i; alpar@9: rs_head[len] = i; alpar@9: } alpar@9: /* initially all rows and columns of the matrix A are active */ alpar@9: for (i = 1; i <= m; i++) rn[i] = 0; alpar@9: for (j = 1; j <= n; j++) cn[j] = 0; alpar@9: /* set initial bounds of the active submatrix */ alpar@9: k1 = 1, k2 = n; alpar@9: /* main loop starts here */ alpar@9: while (k1 <= k2) alpar@9: { i = rs_head[1]; alpar@9: if (i != 0) alpar@9: { /* the i-th row of the matrix A is a row singleton, since alpar@9: it has the only non-zero in the active submatrix */ alpar@9: xassert(rs_len[i] == 1); alpar@9: /* determine the number j of an active column of the matrix alpar@9: A, in which this non-zero is placed */ alpar@9: j = 0; alpar@9: t = mat(info, +i, ndx); alpar@9: xassert(0 <= t && t <= n); alpar@9: for (t = t; t >= 1; t--) alpar@9: { jj = ndx[t]; alpar@9: xassert(1 <= jj && jj <= n); alpar@9: if (cn[jj] == 0) alpar@9: { xassert(j == 0); alpar@9: j = jj; alpar@9: } alpar@9: } alpar@9: xassert(j != 0); alpar@9: /* the singleton is a[i,j]; move a[i,j] to the position alpar@9: b[k1,k1] of the matrix B */ alpar@9: rn[i] = cn[j] = k1; alpar@9: /* shift the left bound of the active submatrix */ alpar@9: k1++; alpar@9: /* increase the size of the lower triangular part */ alpar@9: size++; alpar@9: } alpar@9: else alpar@9: { /* the current active submatrix has no row singletons */ alpar@9: /* remove an active column with maximal number of non-zeros alpar@9: from the active submatrix */ alpar@9: j = cs_head; alpar@9: xassert(j != 0); alpar@9: cn[j] = k2; alpar@9: /* shift the right bound of the active submatrix */ alpar@9: k2--; alpar@9: } alpar@9: /* the j-th column of the matrix A has been removed from the alpar@9: active submatrix */ alpar@9: /* remove the j-th column from the linked list */ alpar@9: if (cs_prev[j] == 0) alpar@9: cs_head = cs_next[j]; alpar@9: else alpar@9: cs_next[cs_prev[j]] = cs_next[j]; alpar@9: if (cs_next[j] == 0) alpar@9: /* nop */; alpar@9: else alpar@9: cs_prev[cs_next[j]] = cs_prev[j]; alpar@9: /* go through non-zeros of the j-th columns and update active alpar@9: lengths of the corresponding rows */ alpar@9: t = mat(info, -j, ndx); alpar@9: xassert(0 <= t && t <= m); alpar@9: for (t = t; t >= 1; t--) alpar@9: { i = ndx[t]; alpar@9: xassert(1 <= i && i <= m); alpar@9: /* the non-zero a[i,j] has left the active submatrix */ alpar@9: len = rs_len[i]; alpar@9: xassert(len >= 1); alpar@9: /* remove the i-th row from the linked list of rows with alpar@9: active length len */ alpar@9: if (rs_prev[i] == 0) alpar@9: rs_head[len] = rs_next[i]; alpar@9: else alpar@9: rs_next[rs_prev[i]] = rs_next[i]; alpar@9: if (rs_next[i] == 0) alpar@9: /* nop */; alpar@9: else alpar@9: rs_prev[rs_next[i]] = rs_prev[i]; alpar@9: /* decrease the active length of the i-th row */ alpar@9: rs_len[i] = --len; alpar@9: /* return the i-th row to the corresponding linked list */ alpar@9: rs_prev[i] = 0; alpar@9: rs_next[i] = rs_head[len]; alpar@9: if (rs_next[i] != 0) rs_prev[rs_next[i]] = i; alpar@9: rs_head[len] = i; alpar@9: } alpar@9: } alpar@9: /* other rows of the matrix A, which are still active, correspond alpar@9: to rows k1, ..., m of the matrix B (in arbitrary order) */ alpar@9: for (i = 1; i <= m; i++) if (rn[i] == 0) rn[i] = k1++; alpar@9: /* but for columns this is not needed, because now the submatrix alpar@9: B2 has no columns */ alpar@9: for (j = 1; j <= n; j++) xassert(cn[j] != 0); alpar@9: /* perform some optional checks */ alpar@9: /* make sure that rn is a permutation of {1, ..., m} and cn is a alpar@9: permutation of {1, ..., n} */ alpar@9: rn_inv = rs_len; /* used as working array */ alpar@9: for (ii = 1; ii <= m; ii++) rn_inv[ii] = 0; alpar@9: for (i = 1; i <= m; i++) alpar@9: { ii = rn[i]; alpar@9: xassert(1 <= ii && ii <= m); alpar@9: xassert(rn_inv[ii] == 0); alpar@9: rn_inv[ii] = i; alpar@9: } alpar@9: cn_inv = rs_head; /* used as working array */ alpar@9: for (jj = 1; jj <= n; jj++) cn_inv[jj] = 0; alpar@9: for (j = 1; j <= n; j++) alpar@9: { jj = cn[j]; alpar@9: xassert(1 <= jj && jj <= n); alpar@9: xassert(cn_inv[jj] == 0); alpar@9: cn_inv[jj] = j; alpar@9: } alpar@9: /* make sure that the matrix B = P*A*Q really has the form, which alpar@9: was declared */ alpar@9: for (ii = 1; ii <= size; ii++) alpar@9: { int diag = 0; alpar@9: i = rn_inv[ii]; alpar@9: t = mat(info, +i, ndx); alpar@9: xassert(0 <= t && t <= n); alpar@9: for (t = t; t >= 1; t--) alpar@9: { j = ndx[t]; alpar@9: xassert(1 <= j && j <= n); alpar@9: jj = cn[j]; alpar@9: if (jj <= size) xassert(jj <= ii); alpar@9: if (jj == ii) alpar@9: { xassert(!diag); alpar@9: diag = 1; alpar@9: } alpar@9: } alpar@9: xassert(diag); alpar@9: } alpar@9: /* free working arrays */ alpar@9: xfree(ndx); alpar@9: xfree(rs_len); alpar@9: xfree(rs_head); alpar@9: xfree(rs_prev); alpar@9: xfree(rs_next); alpar@9: xfree(cs_prev); alpar@9: xfree(cs_next); alpar@9: /* return to the calling program */ alpar@9: return size; alpar@9: } alpar@9: alpar@9: /*---------------------------------------------------------------------- alpar@9: -- adv_basis - construct advanced initial LP basis. alpar@9: -- alpar@9: -- *Synopsis* alpar@9: -- alpar@9: -- #include "glpini.h" alpar@9: -- void adv_basis(glp_prob *lp); alpar@9: -- alpar@9: -- *Description* alpar@9: -- alpar@9: -- The routine adv_basis constructs an advanced initial basis for an LP alpar@9: -- problem object, which the parameter lp points to. alpar@9: -- alpar@9: -- In order to build the initial basis the routine does the following: alpar@9: -- alpar@9: -- 1) includes in the basis all non-fixed auxiliary variables; alpar@9: -- alpar@9: -- 2) includes in the basis as many as possible non-fixed structural alpar@9: -- variables preserving triangular form of the basis matrix; alpar@9: -- alpar@9: -- 3) includes in the basis appropriate (fixed) auxiliary variables alpar@9: -- in order to complete the basis. alpar@9: -- alpar@9: -- As a result the initial basis has minimum of fixed variables and the alpar@9: -- corresponding basis matrix is triangular. */ alpar@9: alpar@9: static int mat(void *info, int k, int ndx[]) alpar@9: { /* this auxiliary routine returns the pattern of a given row or alpar@9: a given column of the augmented constraint matrix A~ = (I|-A), alpar@9: in which columns of fixed variables are implicitly cleared */ alpar@9: LPX *lp = info; alpar@9: int m = lpx_get_num_rows(lp); alpar@9: int n = lpx_get_num_cols(lp); alpar@9: int typx, i, j, lll, len = 0; alpar@9: if (k > 0) alpar@9: { /* the pattern of the i-th row is required */ alpar@9: i = +k; alpar@9: xassert(1 <= i && i <= m); alpar@9: #if 0 /* 22/XII-2003 */ alpar@9: /* if the auxiliary variable x[i] is non-fixed, include its alpar@9: element (placed in the i-th column) in the pattern */ alpar@9: lpx_get_row_bnds(lp, i, &typx, NULL, NULL); alpar@9: if (typx != LPX_FX) ndx[++len] = i; alpar@9: /* include in the pattern elements placed in columns, which alpar@9: correspond to non-fixed structural varables */ alpar@9: i_beg = aa_ptr[i]; alpar@9: i_end = i_beg + aa_len[i] - 1; alpar@9: for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++) alpar@9: { j = m + sv_ndx[i_ptr]; alpar@9: lpx_get_col_bnds(lp, j-m, &typx, NULL, NULL); alpar@9: if (typx != LPX_FX) ndx[++len] = j; alpar@9: } alpar@9: #else alpar@9: lll = lpx_get_mat_row(lp, i, ndx, NULL); alpar@9: for (k = 1; k <= lll; k++) alpar@9: { lpx_get_col_bnds(lp, ndx[k], &typx, NULL, NULL); alpar@9: if (typx != LPX_FX) ndx[++len] = m + ndx[k]; alpar@9: } alpar@9: lpx_get_row_bnds(lp, i, &typx, NULL, NULL); alpar@9: if (typx != LPX_FX) ndx[++len] = i; alpar@9: #endif alpar@9: } alpar@9: else alpar@9: { /* the pattern of the j-th column is required */ alpar@9: j = -k; alpar@9: xassert(1 <= j && j <= m+n); alpar@9: /* if the (auxiliary or structural) variable x[j] is fixed, alpar@9: the pattern of its column is empty */ alpar@9: if (j <= m) alpar@9: lpx_get_row_bnds(lp, j, &typx, NULL, NULL); alpar@9: else alpar@9: lpx_get_col_bnds(lp, j-m, &typx, NULL, NULL); alpar@9: if (typx != LPX_FX) alpar@9: { if (j <= m) alpar@9: { /* x[j] is non-fixed auxiliary variable */ alpar@9: ndx[++len] = j; alpar@9: } alpar@9: else alpar@9: { /* x[j] is non-fixed structural variables */ alpar@9: #if 0 /* 22/XII-2003 */ alpar@9: j_beg = aa_ptr[j]; alpar@9: j_end = j_beg + aa_len[j] - 1; alpar@9: for (j_ptr = j_beg; j_ptr <= j_end; j_ptr++) alpar@9: ndx[++len] = sv_ndx[j_ptr]; alpar@9: #else alpar@9: len = lpx_get_mat_col(lp, j-m, ndx, NULL); alpar@9: #endif alpar@9: } alpar@9: } alpar@9: } alpar@9: /* return the length of the row/column pattern */ alpar@9: return len; alpar@9: } alpar@9: alpar@9: static void adv_basis(glp_prob *lp) alpar@9: { int m = lpx_get_num_rows(lp); alpar@9: int n = lpx_get_num_cols(lp); alpar@9: int i, j, jj, k, size; alpar@9: int *rn, *cn, *rn_inv, *cn_inv; alpar@9: int typx, *tagx = xcalloc(1+m+n, sizeof(int)); alpar@9: double lb, ub; alpar@9: xprintf("Constructing initial basis...\n"); alpar@9: #if 0 /* 13/V-2009 */ alpar@9: if (m == 0) alpar@9: xerror("glp_adv_basis: problem has no rows\n"); alpar@9: if (n == 0) alpar@9: xerror("glp_adv_basis: problem has no columns\n"); alpar@9: #else alpar@9: if (m == 0 || n == 0) alpar@9: { glp_std_basis(lp); alpar@9: return; alpar@9: } alpar@9: #endif alpar@9: /* use the routine triang (see above) to find maximal triangular alpar@9: part of the augmented constraint matrix A~ = (I|-A); in order alpar@9: to prevent columns of fixed variables to be included in the alpar@9: triangular part, such columns are implictly removed from the alpar@9: matrix A~ by the routine adv_mat */ alpar@9: rn = xcalloc(1+m, sizeof(int)); alpar@9: cn = xcalloc(1+m+n, sizeof(int)); alpar@9: size = triang(m, m+n, lp, mat, rn, cn); alpar@9: if (lpx_get_int_parm(lp, LPX_K_MSGLEV) >= 3) alpar@9: xprintf("Size of triangular part = %d\n", size); alpar@9: /* the first size rows and columns of the matrix P*A~*Q (where alpar@9: P and Q are permutation matrices defined by the arrays rn and alpar@9: cn) form a lower triangular matrix; build the arrays (rn_inv alpar@9: and cn_inv), which define the matrices inv(P) and inv(Q) */ alpar@9: rn_inv = xcalloc(1+m, sizeof(int)); alpar@9: cn_inv = xcalloc(1+m+n, sizeof(int)); alpar@9: for (i = 1; i <= m; i++) rn_inv[rn[i]] = i; alpar@9: for (j = 1; j <= m+n; j++) cn_inv[cn[j]] = j; alpar@9: /* include the columns of the matrix A~, which correspond to the alpar@9: first size columns of the matrix P*A~*Q, in the basis */ alpar@9: for (k = 1; k <= m+n; k++) tagx[k] = -1; alpar@9: for (jj = 1; jj <= size; jj++) alpar@9: { j = cn_inv[jj]; alpar@9: /* the j-th column of A~ is the jj-th column of P*A~*Q */ alpar@9: tagx[j] = LPX_BS; alpar@9: } alpar@9: /* if size < m, we need to add appropriate columns of auxiliary alpar@9: variables to the basis */ alpar@9: for (jj = size + 1; jj <= m; jj++) alpar@9: { /* the jj-th column of P*A~*Q should be replaced by the column alpar@9: of the auxiliary variable, for which the only unity element alpar@9: is placed in the position [jj,jj] */ alpar@9: i = rn_inv[jj]; alpar@9: /* the jj-th row of P*A~*Q is the i-th row of A~, but in the alpar@9: i-th row of A~ the unity element belongs to the i-th column alpar@9: of A~; therefore the disired column corresponds to the i-th alpar@9: auxiliary variable (note that this column doesn't belong to alpar@9: the triangular part found by the routine triang) */ alpar@9: xassert(1 <= i && i <= m); alpar@9: xassert(cn[i] > size); alpar@9: tagx[i] = LPX_BS; alpar@9: } alpar@9: /* free working arrays */ alpar@9: xfree(rn); alpar@9: xfree(cn); alpar@9: xfree(rn_inv); alpar@9: xfree(cn_inv); alpar@9: /* build tags of non-basic variables */ alpar@9: for (k = 1; k <= m+n; k++) alpar@9: { if (tagx[k] != LPX_BS) alpar@9: { if (k <= m) alpar@9: lpx_get_row_bnds(lp, k, &typx, &lb, &ub); alpar@9: else alpar@9: lpx_get_col_bnds(lp, k-m, &typx, &lb, &ub); alpar@9: switch (typx) alpar@9: { case LPX_FR: alpar@9: tagx[k] = LPX_NF; break; alpar@9: case LPX_LO: alpar@9: tagx[k] = LPX_NL; break; alpar@9: case LPX_UP: alpar@9: tagx[k] = LPX_NU; break; alpar@9: case LPX_DB: alpar@9: tagx[k] = alpar@9: (fabs(lb) <= fabs(ub) ? LPX_NL : LPX_NU); alpar@9: break; alpar@9: case LPX_FX: alpar@9: tagx[k] = LPX_NS; break; alpar@9: default: alpar@9: xassert(typx != typx); alpar@9: } alpar@9: } alpar@9: } alpar@9: for (k = 1; k <= m+n; k++) alpar@9: { if (k <= m) alpar@9: lpx_set_row_stat(lp, k, tagx[k]); alpar@9: else alpar@9: lpx_set_col_stat(lp, k-m, tagx[k]); alpar@9: } alpar@9: xfree(tagx); alpar@9: return; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * glp_adv_basis - construct advanced initial LP basis alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * void glp_adv_basis(glp_prob *lp, int flags); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine glp_adv_basis constructs an advanced initial basis for alpar@9: * the specified problem object. alpar@9: * alpar@9: * The parameter flags is reserved for use in the future and must be alpar@9: * specified as zero. */ alpar@9: alpar@9: void glp_adv_basis(glp_prob *lp, int flags) alpar@9: { if (flags != 0) alpar@9: xerror("glp_adv_basis: flags = %d; invalid flags\n", flags); alpar@9: if (lp->m == 0 || lp->n == 0) alpar@9: glp_std_basis(lp); alpar@9: else alpar@9: adv_basis(lp); alpar@9: return; alpar@9: } alpar@9: alpar@9: /* eof */