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