1 /* $Id: CoinPresolveHelperFunctions.cpp 2083 2019-01-06 19:38:09Z unxusr $ */
2 // Copyright (C) 2002, International Business Machines
3 // Corporation and others. All Rights Reserved.
4 // This code is licensed under the terms of the Eclipse Public License (EPL).
5
6 /*! \file
7
8 This file contains helper functions for CoinPresolve. The declarations needed
9 for use are included in CoinPresolveMatrix.hpp.
10 */
11
12 #include <stdio.h>
13
14 #include <cassert>
15 #include <iostream>
16
17 #include "CoinHelperFunctions.hpp"
18 #include "CoinPresolveMatrix.hpp"
19
20 /*! \defgroup PMMDVX Packed Matrix Major Dimension Vector Expansion
21 \brief Functions to help with major-dimension vector expansion in a
22 packed matrix structure.
23
24 This next block of functions handles the problems associated with expanding
25 a column in a column-major representation or a row in a row-major
26 representation.
27
28 We need to be able to answer the questions:
29 * Is there room to expand a major vector in place?
30 * Is there sufficient free space at the end of the element and minor
31 index storage areas (bulk storage) to hold the major vector?
32
33 When the answer to the first question is `no', we need to be able to move
34 the major vector to the free space at the end of bulk storage.
35
36 When the answer to the second question is `no', we need to be able to
37 compact the major vectors in the bulk storage area in order to regain a
38 block of useable space at the end.
39
40 presolve_make_memlists initialises a linked list that tracks the position of
41 major vectors in the bulk storage area. It's used to locate physically
42 adjacent vectors.
43
44 presolve_expand deals with adding a coefficient to a major vector, either
45 in-place or by moving it to free space at the end of the storage areas.
46 There are two inline wrappers, presolve_expand_col and presolve_expand_row,
47 defined in CoinPresolveMatrix.hpp.
48
49 compact_rep compacts the major vectors in the storage areas to
50 leave a single block of free space at the end.
51 */
52 //@{
53
54 /*
55 This first function doesn't need to be known outside of this file.
56 */
57 namespace {
58
59 /*
60 compact_rep
61
62 This routine compacts the major vectors in the bulk storage area,
63 leaving a single block of free space at the end. The vectors are not
64 reordered, just shifted down to remove gaps.
65 */
66
compact_rep(double * elems,int * indices,CoinBigIndex * starts,const int * lengths,int n,const presolvehlink * link)67 void compact_rep(double *elems, int *indices,
68 CoinBigIndex *starts, const int *lengths, int n,
69 const presolvehlink *link)
70 {
71 #if PRESOLVE_SUMMARY
72 printf("****COMPACTING****\n");
73 #endif
74
75 // for now, just look for the first element of the list
76 int i = n;
77 while (link[i].pre != NO_LINK)
78 i = link[i].pre;
79
80 CoinBigIndex j = 0;
81 for (; i != n; i = link[i].suc) {
82 CoinBigIndex s = starts[i];
83 CoinBigIndex e = starts[i] + lengths[i];
84
85 // because of the way link is organized, j <= s
86 starts[i] = j;
87 for (CoinBigIndex k = s; k < e; k++) {
88 elems[j] = elems[k];
89 indices[j] = indices[k];
90 j++;
91 }
92 }
93 }
94
95 } /* end unnamed namespace */
96
97 /*
98 \brief Initialise linked list for major vector order in bulk storage
99
100 Initialise the linked list that will track the order of major vectors in
101 the element and row index bulk storage arrays. When finished, link[j].pre
102 contains the index of the previous non-empty vector in the storage arrays
103 and link[j].suc contains the index of the next non-empty vector.
104
105 For an empty vector j, link[j].pre = link[j].suc = NO_LINK.
106
107 If n is the number of major-dimension vectors, link[n] is valid;
108 link[n].pre is the index of the last non-empty vector, and
109 link[n].suc = NO_LINK.
110
111 This routine makes the implicit assumption that the order of vectors in the
112 storage areas matches the order in starts. (I.e., there's no check that
113 starts[j] > starts[i] for j < i.)
114 */
115
presolve_make_memlists(int * lengths,presolvehlink * link,int n)116 void presolve_make_memlists(/*CoinBigIndex *starts,*/ int *lengths,
117 presolvehlink *link, int n)
118 {
119 int i;
120 int pre = NO_LINK;
121
122 for (i = 0; i < n; i++) {
123 if (lengths[i]) {
124 link[i].pre = pre;
125 if (pre != NO_LINK)
126 link[pre].suc = i;
127 pre = i;
128 } else {
129 link[i].pre = NO_LINK;
130 link[i].suc = NO_LINK;
131 }
132 }
133 if (pre != NO_LINK)
134 link[pre].suc = n;
135
136 // (1) Arbitrarily place the last non-empty entry in link[n].pre
137 link[n].pre = pre;
138
139 link[n].suc = NO_LINK;
140 }
141
142 /*
143 presolve_expand_major
144
145 The routine looks at the space currently occupied by major-dimension vector
146 k and makes sure that there's room to add one coefficient.
147
148 This may require moving the vector to the vacant area at the end of the
149 bulk storage array. If there's no room left at the end of the array, an
150 attempt is made to compact the existing vectors to make space.
151
152 Returns true for failure, false for success.
153 */
154
presolve_expand_major(CoinBigIndex * majstrts,double * els,int * minndxs,int * majlens,presolvehlink * majlinks,int nmaj,int k)155 bool presolve_expand_major(CoinBigIndex *majstrts, double *els,
156 int *minndxs, int *majlens,
157 presolvehlink *majlinks, int nmaj, int k)
158
159 {
160 const CoinBigIndex bulkCap = majstrts[nmaj];
161
162 /*
163 Get the start and end of column k, and the index of the column which
164 follows in the bulk storage.
165 */
166 CoinBigIndex kcsx = majstrts[k];
167 CoinBigIndex kcex = kcsx + majlens[k];
168 int nextcol = majlinks[k].suc;
169 /*
170 Do we have room to add one coefficient in place?
171 */
172 if (kcex + 1 < majstrts[nextcol]) { /* no action required */
173 }
174 /*
175 Is k the last non-empty column? In that case, attempt to compact the
176 bulk storage. This will move k, so update the column start and end.
177 If we still have no space, it's a fatal error.
178 */
179 else if (nextcol == nmaj) {
180 compact_rep(els, minndxs, majstrts, majlens, nmaj, majlinks);
181 kcsx = majstrts[k];
182 kcex = kcsx + majlens[k];
183 if (kcex + 1 >= bulkCap) {
184 return (true);
185 }
186 }
187 /*
188 The most complicated case --- we need to move k from its current location
189 to empty space at the end of the bulk storage. And we may need to make that!
190 Compaction is identical to the above case.
191 */
192 else {
193 int lastcol = majlinks[nmaj].pre;
194 CoinBigIndex newkcsx = majstrts[lastcol] + majlens[lastcol];
195 CoinBigIndex newkcex = newkcsx + majlens[k];
196
197 if (newkcex + 1 >= bulkCap) {
198 compact_rep(els, minndxs, majstrts, majlens, nmaj, majlinks);
199 kcsx = majstrts[k];
200 kcex = kcsx + majlens[k];
201 newkcsx = majstrts[lastcol] + majlens[lastcol];
202 newkcex = newkcsx + majlens[k];
203 }
204 /*
205 It is possible (probably on very small problems) that
206 there is not room for two copies of modified column -
207 so allow temporary overflow
208 */
209 /*
210 Moving the vector requires three actions. First we move the data, then
211 update the packed matrix vector start, then relink the storage order list,
212 */
213 memcpy(reinterpret_cast< void * >(&minndxs[newkcsx]),
214 reinterpret_cast< void * >(&minndxs[kcsx]), majlens[k] * sizeof(int));
215 memcpy(reinterpret_cast< void * >(&els[newkcsx]),
216 reinterpret_cast< void * >(&els[kcsx]), majlens[k] * sizeof(double));
217 majstrts[k] = newkcsx;
218 PRESOLVE_REMOVE_LINK(majlinks, k);
219 PRESOLVE_INSERT_LINK(majlinks, k, lastcol);
220 if (newkcex + 1 >= bulkCap) {
221 // compact - faking extra one
222 //majlens[k]++;
223 //majstrts[nmaj]=newkcex;
224 compact_rep(els, minndxs, majstrts, majlens, nmaj, majlinks);
225 //majlens[k]--;
226 //majstrts[nmaj]=bulkCap;
227 kcsx = majstrts[k];
228 kcex = kcsx + majlens[k];
229 if (kcex > bulkCap) {
230 // still no room
231 return (true);
232 }
233 }
234 }
235 /*
236 Success --- the vector has room for one more coefficient.
237 */
238 return (false);
239 }
240
241 //@}
242
243 /*
244 Helper function to duplicate a major-dimension vector.
245 */
246
247 /*
248 A major-dimension vector is composed of paired element and minor index
249 arrays. We want to duplicate length entries from both arrays, starting at
250 offset.
251
252 If tgt > 0, we'll run a more complicated copy loop which will elide the
253 entry with minor index == tgt. In this case, we want to reduce the size of
254 the allocated array by 1.
255
256 Pigs will fly before sizeof(int) > sizeof(double), but if it ever
257 happens this code will fail.
258 */
259
presolve_dupmajor(const double * elems,const int * indices,int length,CoinBigIndex offset,int tgt)260 double *presolve_dupmajor(const double *elems, const int *indices,
261 int length, CoinBigIndex offset, int tgt)
262
263 {
264 int n;
265
266 if (tgt >= 0)
267 length--;
268
269 if (2 * sizeof(int) <= sizeof(double))
270 n = (3 * length + 1) >> 1;
271 else
272 n = 2 * length;
273
274 double *dArray = new double[n];
275 int *iArray = reinterpret_cast< int * >(dArray + length);
276
277 if (tgt < 0) {
278 memcpy(dArray, elems + offset, length * sizeof(double));
279 memcpy(iArray, indices + offset, length * sizeof(int));
280 } else {
281 int korig;
282 int kcopy = 0;
283 indices += offset;
284 elems += offset;
285 for (korig = 0; korig <= length; korig++) {
286 int i = indices[korig];
287 if (i != tgt) {
288 dArray[kcopy] = elems[korig];
289 iArray[kcopy++] = indices[korig];
290 }
291 }
292 }
293
294 return (dArray);
295 }
296
297 /*
298 Routines to find the position of the entry for a given minor index in a
299 major vector. Inline wrappers with column-major and row-major parameter
300 names are defined in CoinPresolveMatrix.hpp. The threaded matrix used in
301 postsolve exists only as a column-major form, so only one wrapper is
302 defined.
303 */
304
305 /*
306 presolve_find_minor
307
308 Find the position (k) of the entry for a given minor index (tgt) within
309 the range of entries for a major vector (ks, ke).
310
311 Print a tag and abort (DIE) if there's no entry for tgt.
312 */
313 #if 0
314 CoinBigIndex presolve_find_minor (int tgt, CoinBigIndex ks,
315 CoinBigIndex ke, const int *minndxs)
316
317 { CoinBigIndex k ;
318 for (k = ks ; k < ke ; k++)
319 { if (minndxs[k] == tgt)
320 return (k) ; }
321 DIE("FIND_MINOR") ;
322
323 abort () ; return -1; }
324 #endif
325 /*
326 As presolve_find_minor, but return a position one past the end of
327 the major vector when the entry is not already present.
328 */
presolve_find_minor1(int tgt,CoinBigIndex ks,CoinBigIndex ke,const int * minndxs)329 CoinBigIndex presolve_find_minor1(int tgt, CoinBigIndex ks,
330 CoinBigIndex ke, const int *minndxs)
331 {
332 CoinBigIndex k;
333 for (k = ks; k < ke; k++) {
334 if (minndxs[k] == tgt)
335 return (k);
336 }
337
338 return (k);
339 }
340
341 /*
342 In a threaded matrix, the major vector does not occupy a contiguous block
343 in the bulk storage area. For example, in a threaded column-major matrix,
344 if a<i,p> is in pos'n kp of hrow, the next coefficient a<i,q> will be
345 in pos'n kq = link[kp]. Abort if we don't find it.
346 */
presolve_find_minor2(int tgt,CoinBigIndex ks,int majlen,const int * minndxs,const CoinBigIndex * majlinks)347 CoinBigIndex presolve_find_minor2(int tgt, CoinBigIndex ks,
348 int majlen, const int *minndxs,
349 const CoinBigIndex *majlinks)
350
351 {
352 for (int i = 0; i < majlen; ++i) {
353 if (minndxs[ks] == tgt)
354 return (ks);
355 ks = majlinks[ks];
356 }
357 DIE("FIND_MINOR2");
358
359 abort();
360 return -1;
361 }
362
363 /*
364 As presolve_find_minor2, but return -1 if the entry is missing
365 */
presolve_find_minor3(int tgt,CoinBigIndex ks,int majlen,const int * minndxs,const CoinBigIndex * majlinks)366 CoinBigIndex presolve_find_minor3(int tgt, CoinBigIndex ks,
367 int majlen, const int *minndxs,
368 const CoinBigIndex *majlinks)
369
370 {
371 for (int i = 0; i < majlen; ++i) {
372 if (minndxs[ks] == tgt)
373 return (ks);
374 ks = majlinks[ks];
375 }
376
377 return (-1);
378 }
379
380 /*
381 Delete the entry for a minor index from a major vector. The last entry in
382 the major vector is moved into the hole left by the deleted entry. This
383 leaves some space between the end of this major vector and the start of the
384 next in the bulk storage areas (this is termed loosely packed). Inline
385 wrappers with column-major and row-major parameter names are defined in
386 CoinPresolveMatrix.hpp. The threaded matrix used in postsolve exists only
387 as a column-major form, so only one wrapper is defined.
388 */
389
390 #if 0
391 void presolve_delete_from_major (int majndx, int minndx,
392 const CoinBigIndex *majstrts,
393 int *majlens, int *minndxs, double *els)
394
395 { CoinBigIndex ks = majstrts[majndx] ;
396 CoinBigIndex ke = ks + majlens[majndx] ;
397
398 CoinBigIndex kmi = presolve_find_minor(minndx,ks,ke,minndxs) ;
399
400 minndxs[kmi] = minndxs[ke-1] ;
401 els[kmi] = els[ke-1] ;
402 majlens[majndx]-- ;
403
404 return ; }
405 // Delete all marked and zero marked
406 void presolve_delete_many_from_major (int majndx, char * marked,
407 const CoinBigIndex *majstrts,
408 int *majlens, int *minndxs, double *els)
409
410 {
411 CoinBigIndex ks = majstrts[majndx] ;
412 CoinBigIndex ke = ks + majlens[majndx] ;
413 CoinBigIndex put=ks;
414 for (CoinBigIndex k=ks;k<ke;k++) {
415 int iMinor = minndxs[k];
416 if (!marked[iMinor]) {
417 minndxs[put]=iMinor;
418 els[put++]=els[k];
419 } else {
420 marked[iMinor]=0;
421 }
422 }
423 majlens[majndx] = put-ks ;
424 return ;
425 }
426 #endif
427
428 /*
429 Delete the entry for a minor index from a major vector in a threaded matrix.
430
431 This involves properly relinking the free list.
432 */
presolve_delete_from_major2(int majndx,int minndx,CoinBigIndex * majstrts,int * majlens,int * minndxs,CoinBigIndex * majlinks,CoinBigIndex * free_listp)433 void presolve_delete_from_major2(int majndx, int minndx,
434 CoinBigIndex *majstrts, int *majlens,
435 int *minndxs, CoinBigIndex *majlinks,
436 CoinBigIndex *free_listp)
437
438 {
439 CoinBigIndex k = majstrts[majndx];
440
441 /*
442 Desired entry is the first in its major vector. We need to touch up majstrts
443 to point to the next entry and link the deleted entry to the front of the
444 free list.
445 */
446 if (minndxs[k] == minndx) {
447 majstrts[majndx] = majlinks[k];
448 majlinks[k] = *free_listp;
449 *free_listp = k;
450 majlens[majndx]--;
451 } else {
452 /*
453 Desired entry is somewhere in the major vector. We need to relink around
454 it and then link it on the front of the free list.
455
456 The loop runs over elements 1 .. len-1; we've already ruled out element 0.
457 */
458 int len = majlens[majndx];
459 CoinBigIndex kpre = k;
460 k = majlinks[k];
461 for (int i = 1; i < len; ++i) {
462 if (minndxs[k] == minndx) {
463 majlinks[kpre] = majlinks[k];
464 majlinks[k] = *free_listp;
465 *free_listp = k;
466 majlens[majndx]--;
467 return;
468 }
469 kpre = k;
470 k = majlinks[k];
471 }
472 DIE("DELETE_FROM_MAJOR2");
473 }
474
475 assert(*free_listp >= 0);
476
477 return;
478 }
479
480 /* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
481 */
482