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