1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
|
/* spxat.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 "spxat.h"
/***********************************************************************
* spx_alloc_at - allocate constraint matrix in sparse row-wise format
*
* This routine allocates the memory for arrays needed to represent the
* constraint matrix in sparse row-wise format. */
void spx_alloc_at(SPXLP *lp, SPXAT *at)
{ int m = lp->m;
int n = lp->n;
int nnz = lp->nnz;
at->ptr = talloc(1+m+1, int);
at->ind = talloc(1+nnz, int);
at->val = talloc(1+nnz, double);
at->work = talloc(1+n, double);
return;
}
/***********************************************************************
* spx_build_at - build constraint matrix in sparse row-wise format
*
* This routine builds sparse row-wise representation of the constraint
* matrix A using its sparse column-wise representation stored in the
* lp object, and stores the result in the at object. */
void spx_build_at(SPXLP *lp, SPXAT *at)
{ 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 *AT_ptr = at->ptr;
int *AT_ind = at->ind;
double *AT_val = at->val;
int i, k, ptr, end, pos;
/* calculate AT_ptr[i] = number of non-zeros in i-th row */
memset(&AT_ptr[1], 0, m * sizeof(int));
for (k = 1; k <= n; k++)
{ ptr = A_ptr[k];
end = A_ptr[k+1];
for (; ptr < end; ptr++)
AT_ptr[A_ind[ptr]]++;
}
/* set AT_ptr[i] to position after last element in i-th row */
AT_ptr[1]++;
for (i = 2; i <= m; i++)
AT_ptr[i] += AT_ptr[i-1];
xassert(AT_ptr[m] == nnz+1);
AT_ptr[m+1] = nnz+1;
/* build row-wise representation and re-arrange AT_ptr[i] */
for (k = n; k >= 1; k--)
{ /* copy elements from k-th column to corresponding rows */
ptr = A_ptr[k];
end = A_ptr[k+1];
for (; ptr < end; ptr++)
{ pos = --AT_ptr[A_ind[ptr]];
AT_ind[pos] = k;
AT_val[pos] = A_val[ptr];
}
}
xassert(AT_ptr[1] == 1);
return;
}
/***********************************************************************
* spx_at_prod - compute product y := y + s * A'* x
*
* This routine computes the product:
*
* y := y + s * A'* x,
*
* where A' is a matrix transposed to the mxn-matrix A of constraint
* coefficients, x is a m-vector, s is a scalar, y is a n-vector.
*
* The routine uses the row-wise representation of the matrix A and
* computes the product as a linear combination:
*
* y := y + s * (A'[1] * x[1] + ... + A'[m] * x[m]),
*
* where A'[i] is i-th row of A, 1 <= i <= m. */
void spx_at_prod(SPXLP *lp, SPXAT *at, double y[/*1+n*/], double s,
const double x[/*1+m*/])
{ int m = lp->m;
int *AT_ptr = at->ptr;
int *AT_ind = at->ind;
double *AT_val = at->val;
int i, ptr, end;
double t;
for (i = 1; i <= m; i++)
{ if (x[i] != 0.0)
{ /* y := y + s * (i-th row of A) * x[i] */
t = s * x[i];
ptr = AT_ptr[i];
end = AT_ptr[i+1];
for (; ptr < end; ptr++)
y[AT_ind[ptr]] += AT_val[ptr] * t;
}
}
return;
}
/***********************************************************************
* spx_nt_prod1 - 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. */
void spx_nt_prod1(SPXLP *lp, SPXAT *at, double y[/*1+n-m*/], int ign,
double s, const double x[/*1+m*/])
{ int m = lp->m;
int n = lp->n;
int *head = lp->head;
double *work = at->work;
int j, k;
for (k = 1; k <= n; k++)
work[k] = 0.0;
if (!ign)
{ for (j = 1; j <= n-m; j++)
work[head[m+j]] = y[j];
}
spx_at_prod(lp, at, work, s, x);
for (j = 1; j <= n-m; j++)
y[j] = work[head[m+j]];
return;
}
/***********************************************************************
* spx_eval_trow1 - compute i-th row of simplex table
*
* This routine computes i-th row of the current simplex table
* T = (T[i,j]) = - inv(B) * N, 1 <= i <= m, using representation of
* the constraint matrix A in row-wise format.
*
* The vector rho = (rho[j]), which is i-th row of the basis inverse
* inv(B), should be previously computed with the routine spx_eval_rho.
* It is assumed that elements of this vector are stored in the array
* locations rho[1], ..., rho[m].
*
* There exist two ways to compute the simplex table row.
*
* 1. T[i,j], j = 1,...,n-m, is computed as inner product:
*
* m
* T[i,j] = - sum a[i,k] * rho[i],
* i=1
*
* where N[j] = A[k] is a column of the constraint matrix corresponding
* to non-basic variable xN[j]. The estimated number of operations in
* this case is:
*
* n1 = (n - m) * (nnz(A) / n),
*
* (n - m) is the number of columns of N, nnz(A) / n is the average
* number of non-zeros in one column of A and, therefore, of N.
*
* 2. The simplex table row is computed as part of a linear combination
* of rows of A with coefficients rho[i] != 0. The estimated number
* of operations in this case is:
*
* n2 = nnz(rho) * (nnz(A) / m),
*
* where nnz(rho) is the number of non-zeros in the vector rho,
* nnz(A) / m is the average number of non-zeros in one row of A.
*
* If n1 < n2, the routine computes the simples table row using the
* first way (like the routine spx_eval_trow). Otherwise, the routine
* uses the second way calling the routine spx_nt_prod1.
*
* On exit components of the simplex table row are stored in the array
* locations trow[1], ... trow[n-m]. */
void spx_eval_trow1(SPXLP *lp, SPXAT *at, const double rho[/*1+m*/],
double trow[/*1+n-m*/])
{ int m = lp->m;
int n = lp->n;
int nnz = lp->nnz;
int i, j, nnz_rho;
double cnt1, cnt2;
/* determine nnz(rho) */
nnz_rho = 0;
for (i = 1; i <= m; i++)
{ if (rho[i] != 0.0)
nnz_rho++;
}
/* estimate the number of operations for both ways */
cnt1 = (double)(n - m) * ((double)nnz / (double)n);
cnt2 = (double)nnz_rho * ((double)nnz / (double)m);
/* compute i-th row of simplex table */
if (cnt1 < cnt2)
{ /* as inner products */
int *A_ptr = lp->A_ptr;
int *A_ind = lp->A_ind;
double *A_val = lp->A_val;
int *head = lp->head;
int k, ptr, end;
double tij;
for (j = 1; j <= n-m; j++)
{ k = head[m+j]; /* x[k] = xN[j] */
/* compute t[i,j] = - N'[j] * pi */
tij = 0.0;
ptr = A_ptr[k];
end = A_ptr[k+1];
for (; ptr < end; ptr++)
tij -= A_val[ptr] * rho[A_ind[ptr]];
trow[j] = tij;
}
}
else
{ /* as linear combination */
spx_nt_prod1(lp, at, trow, 1, -1.0, rho);
}
return;
}
/***********************************************************************
* spx_free_at - deallocate constraint matrix in sparse row-wise format
*
* This routine deallocates the memory used for arrays of the program
* object at. */
void spx_free_at(SPXLP *lp, SPXAT *at)
{ xassert(lp == lp);
tfree(at->ptr);
tfree(at->ind);
tfree(at->val);
tfree(at->work);
return;
}
/* eof */
|
