1 /* $Id: CoinOslC.h 2259 2020-01-16 13:43:43Z stefan $ */
2 #ifndef COIN_OSL_C_INCLUDE
3 /*
4 Copyright (C) 1987, 2009, International Business Machines Corporation
5 and others. All Rights Reserved.
6
7 This code is licensed under the terms of the Eclipse Public License (EPL).
8 */
9 #define COIN_OSL_C_INCLUDE
10
11 #ifndef CLP_OSL
12 #define CLP_OSL 0
13 #endif
14 #define C_EKK_GO_SPARSE 200
15
16 #ifdef HAVE_ENDIAN_H
17 #include <endian.h>
18 #if __BYTE_ORDER == __LITTLE_ENDIAN
19 #define INTEL
20 #endif
21 #endif
22
23 #include <math.h>
24 #include <string.h>
25 #include <stdio.h>
26 #include <stdlib.h>
27
28 #define SPARSE_UPDATE
29 #define NO_SHIFT
30 #include "CoinHelperFunctions.hpp"
31
32 #include <stddef.h>
33 #ifdef __cplusplus
34 extern "C" {
35 #endif
36
37 int c_ekkbtrn(register const EKKfactinfo *fact,
38 double *dwork1,
39 int *mpt, int first_nonzero);
40 int c_ekkbtrn_ipivrw(register const EKKfactinfo *fact,
41 double *dwork1,
42 int *mpt, int ipivrw, int *spare);
43
44 int c_ekketsj(register /*const*/ EKKfactinfo *fact,
45 double *dwork1,
46 int *mpt2, double dalpha, int orig_nincol,
47 int npivot, int *nuspikp,
48 const int ipivrw, int *spare);
49 int c_ekkftrn(register const EKKfactinfo *fact,
50 double *dwork1,
51 double *dpermu, int *mpt, int numberNonZero);
52
53 int c_ekkftrn_ft(register EKKfactinfo *fact,
54 double *dwork1, int *mpt, int *nincolp);
55 void c_ekkftrn2(register EKKfactinfo *fact, double *dwork1,
56 double *dpermu1, int *mpt1, int *nincolp,
57 double *dwork1_ft, int *mpt_ft, int *nincolp_ft);
58
59 int c_ekklfct(register EKKfactinfo *fact);
60 int c_ekkslcf(register const EKKfactinfo *fact);
c_ekkscpy(int n,const int * marr1,int * marr2)61 inline void c_ekkscpy(int n, const int *marr1, int *marr2)
62 {
63 CoinMemcpyN(marr1, n, marr2);
64 }
c_ekkdcpy(int n,const double * marr1,double * marr2)65 inline void c_ekkdcpy(int n, const double *marr1, double *marr2)
66 {
67 CoinMemcpyN(marr1, n, marr2);
68 }
69 int c_ekk_IsSet(const int *array, int bit);
70 void c_ekk_Set(int *array, int bit);
71 void c_ekk_Unset(int *array, int bit);
72
73 void c_ekkzero(int length, int n, void *array);
c_ekkdzero(int n,double * marray)74 inline void c_ekkdzero(int n, double *marray)
75 {
76 CoinZeroN(marray, n);
77 }
c_ekkizero(int n,int * marray)78 inline void c_ekkizero(int n, int *marray)
79 {
80 CoinZeroN(marray, n);
81 }
c_ekkczero(int n,char * marray)82 inline void c_ekkczero(int n, char *marray)
83 {
84 CoinZeroN(marray, n);
85 }
86 #ifdef __cplusplus
87 }
88 #endif
89
90 #define c_ekkscpy_0_1(s, ival, array) CoinFillN(array, s, ival)
91 #define c_ekks1cpy(n, marr1, marr2) CoinMemcpyN(marr1, n, marr2)
92 void clp_setup_pointers(EKKfactinfo *fact);
93 void clp_memory(int type);
94 double *clp_double(int number_entries);
95 int *clp_int(int number_entries);
96 void *clp_malloc(int number_entries);
97 void clp_free(void *oldArray);
98
99 #define SLACK_VALUE -1.0
100 #define C_EKK_REMOVE_LINK(hpiv, hin, link, ipivot) \
101 { \
102 int ipre = link[ipivot].pre; \
103 int isuc = link[ipivot].suc; \
104 if (ipre > 0) { \
105 link[ipre].suc = isuc; \
106 } \
107 if (ipre <= 0) { \
108 hpiv[hin[ipivot]] = isuc; \
109 } \
110 if (isuc > 0) { \
111 link[isuc].pre = ipre; \
112 } \
113 }
114
115 #define C_EKK_ADD_LINK(hpiv, nzi, link, npr) \
116 { \
117 int ifiri = hpiv[nzi]; \
118 hpiv[nzi] = npr; \
119 link[npr].suc = ifiri; \
120 link[npr].pre = 0; \
121 if (ifiri != 0) { \
122 link[ifiri].pre = npr; \
123 } \
124 }
125 #include <assert.h>
126 #ifdef NO_SHIFT
127
128 #define SHIFT_INDEX(limit) (limit)
129 #define UNSHIFT_INDEX(limit) (limit)
130 #define SHIFT_REF(arr, ind) (arr)[ind]
131
132 #else
133
134 #define SHIFT_INDEX(limit) ((limit) << 3)
135 #define UNSHIFT_INDEX(limit) ((unsigned int)(limit) >> 3)
136 #define SHIFT_REF(arr, ind) (*(double *)((char *)(arr) + (ind)))
137
138 #endif
139
140 #ifdef INTEL
141 #define NOT_ZERO(x) (((*((reinterpret_cast< unsigned char * >(&x)) + 7)) & 0x7F) != 0)
142 #else
143 #define NOT_ZERO(x) ((x) != 0.0)
144 #endif
145
146 #define SWAP(type, _x, _y) \
147 { \
148 type _tmp = (_x); \
149 (_x) = (_y); \
150 (_y) = _tmp; \
151 }
152
153 #define UNROLL_LOOP_BODY1(code) \
154 { \
155 { \
156 code \
157 } \
158 }
159 #define UNROLL_LOOP_BODY2(code) \
160 { \
161 { code } { \
162 code \
163 } \
164 }
165 #define UNROLL_LOOP_BODY4(code) \
166 { \
167 { code } { code } { code } { \
168 code \
169 } \
170 }
171 #endif
172 #ifdef COIN_OSL_CMFC
173 /* Return codes in IRTCOD/IRTCOD are */
174 /* 4: numerical problems */
175 /* 5: not enough space in row file */
176 /* 6: not enough space in column file */
177 /* 23: system error at label 320 */
178 {
179 #if 1
180 int *hcoli = fact->xecadr;
181 double *dluval = fact->xeeadr;
182 double *dvalpv = fact->kw3adr;
183 int *mrstrt = fact->xrsadr;
184 int *hrowi = fact->xeradr;
185 int *mcstrt = fact->xcsadr;
186 int *hinrow = fact->xrnadr;
187 int *hincol = fact->xcnadr;
188 int *hpivro = fact->krpadr;
189 int *hpivco = fact->kcpadr;
190 #endif
191 int nnentl = fact->nnentl;
192 int nnentu = fact->nnentu;
193 int kmxeta = fact->kmxeta;
194 int xnewro = *xnewrop;
195 int ncompactions = *ncompactionsp;
196
197 MACTION_T *maction = reinterpret_cast< MACTION_T * >(maction_void);
198
199 int i, j, k;
200 double d1;
201 int j1, j2;
202 int jj, kk, kr, nz, jj1, jj2, kce, kcs, kqq, npr;
203 int fill, naft;
204 int enpr;
205 int nres, npre;
206 int knpr, irow, iadd32, ibase;
207 double pivot;
208 int count, nznpr;
209 int nlast, epivr1;
210 int kipis;
211 double dpivx;
212 int kipie, kcpiv, knprs, knpre;
213 bool cancel;
214 double multip, elemnt;
215 int ipivot, jpivot, epivro, epivco, lstart, nfirst;
216 int nzpivj, kfill, kstart;
217 int nmove, ileft;
218 #ifndef C_EKKCMFY
219 int iput, nspare;
220 int noRoomForDense = 0;
221 int if_sparse_update = fact->if_sparse_update;
222 int ifdens = 0;
223 #endif
224 int irtcod = 0;
225 const int nrow = fact->nrow;
226
227 /* Parameter adjustments */
228 --maction;
229
230 /* Function Body */
231 lstart = nnetas - nnentl + 1;
232 for (i = lstart; i <= nnetas; ++i) {
233 hrowi[i] = SHIFT_INDEX(hcoli[i]);
234 }
235
236 for (i = 1; i <= nrow; ++i) {
237 maction[i] = 0;
238 mwork[i].pre = i - 1;
239 mwork[i].suc = i + 1;
240 }
241
242 iadd32 = 0;
243 nlast = nrow;
244 nfirst = 1;
245 mwork[1].pre = nrow;
246 mwork[nrow].suc = 1;
247
248 for (count = 1; count <= nrow; ++count) {
249
250 /* Pick column singletons */
251 if (!(hpivco[1] <= 0)) {
252 int small_pivot = c_ekkcsin(fact,
253 rlink, clink,
254 nsingp);
255
256 if (small_pivot) {
257 irtcod = 7; /* pivot too small */
258 if (fact->invok >= 0) {
259 goto L1050;
260 }
261 }
262 if (fact->npivots >= nrow) {
263 goto L1050;
264 }
265 }
266
267 /* Pick row singletons */
268 if (!(hpivro[1] <= 0)) {
269 irtcod = c_ekkrsin(fact,
270 rlink, clink,
271 mwork, nfirst,
272 nsingp,
273
274 &xnewco, &xnewro,
275 &nnentu,
276 &kmxeta, &ncompactions,
277 &nnentl);
278 if (irtcod != 0) {
279 if (irtcod < 0 || fact->invok >= 0) {
280 /* -5 */
281 goto L1050;
282 }
283 /* ASSERT: irtcod == 7 - pivot too small */
284 /* why don't we return with an error? */
285 }
286 if (fact->npivots >= nrow) {
287 goto L1050;
288 }
289 lstart = nnetas - nnentl + 1;
290 }
291
292 /* Find a pivot element */
293 irtcod = c_ekkfpvt(fact,
294 rlink, clink,
295 nsingp, xrejctp, &ipivot, &jpivot);
296 if (irtcod != 0) {
297 /* irtcod == 10 */
298 goto L1050;
299 }
300 /* Update list structures and prepare for numerical phase */
301 c_ekkprpv(fact, rlink, clink,
302 *xrejctp, ipivot, jpivot);
303
304 epivco = hincol[jpivot];
305 ++fact->xnetal;
306 mcstrt[fact->xnetal] = lstart - 1;
307 hpivco[fact->xnetal] = ipivot;
308 epivro = hinrow[ipivot];
309 epivr1 = epivro - 1;
310 kipis = mrstrt[ipivot];
311 pivot = dluval[kipis];
312 dpivx = 1. / pivot;
313 kipie = kipis + epivr1;
314 ++kipis;
315 #ifndef C_EKKCMFY
316 {
317 double size = nrow - fact->npivots;
318 if (size > GO_DENSE && (nnentu - fact->nuspike) * GO_DENSE_RATIO > size * size) {
319 /* say going to dense coding */
320 if (*nsingp == 0) {
321 ifdens = 1;
322 }
323 }
324 }
325 #endif
326 /* copy the pivot row entries into dvalpv */
327 /* the maction array tells us the index into dvalpv for a given row */
328 /* the alternative would be using a large array of doubles */
329 for (k = kipis; k <= kipie; ++k) {
330 irow = hcoli[k];
331 dvalpv[k - kipis + 1] = dluval[k];
332 maction[irow] = static_cast< MACTION_T >(k - kipis + 1);
333 }
334
335 /* Loop over nonzeros in pivot column */
336 kcpiv = mcstrt[jpivot] - 1;
337 for (nzpivj = 1; nzpivj <= epivco; ++nzpivj) {
338 ++kcpiv;
339 npr = hrowi[kcpiv];
340 hrowi[kcpiv] = 0; /* zero out for possible compaction later on */
341
342 --hincol[jpivot];
343
344 ++mcstrt[jpivot];
345 /* loop invariant: kcpiv == mcstrt[jpivot] - 1 */
346
347 --hinrow[npr];
348 enpr = hinrow[npr];
349 knprs = mrstrt[npr];
350 knpre = knprs + enpr;
351
352 /* Search for element to be eliminated */
353 knpr = knprs;
354 while (1) {
355 UNROLL_LOOP_BODY4({
356 if (jpivot == hcoli[knpr]) {
357 break;
358 }
359 knpr++;
360 });
361 }
362
363 multip = -dluval[knpr] * dpivx;
364
365 /* swap last entry with pivot */
366 dluval[knpr] = dluval[knpre];
367 hcoli[knpr] = hcoli[knpre];
368 --knpre;
369
370 #if 1
371 /* MONSTER_UNROLLED_CODE - see below */
372 kfill = epivr1 - (knpre - knprs + 1);
373 nres = ((knpre - knprs + 1) & 1) + knprs;
374 cancel = false;
375 d1 = 1e33;
376 j1 = hcoli[nres];
377
378 if (nres != knprs) {
379 j = hcoli[knprs];
380 if (maction[j] == 0) {
381 ++kfill;
382 } else {
383 jj = maction[j];
384 maction[j] = static_cast< MACTION_T >(-maction[j]);
385 dluval[knprs] += multip * dvalpv[jj];
386 d1 = fabs(dluval[knprs]);
387 }
388 }
389 j2 = hcoli[nres + 1];
390 jj1 = maction[j1];
391 for (kr = nres; kr < knpre; kr += 2) {
392 jj2 = maction[j2];
393 if (jj1 == 0) {
394 ++kfill;
395 } else {
396 maction[j1] = static_cast< MACTION_T >(-maction[j1]);
397 dluval[kr] += multip * dvalpv[jj1];
398 cancel = cancel || !(fact->zeroTolerance < d1);
399 d1 = fabs(dluval[kr]);
400 }
401 j1 = hcoli[kr + 2];
402 if (jj2 == 0) {
403 ++kfill;
404 } else {
405 maction[j2] = static_cast< MACTION_T >(-maction[j2]);
406 dluval[kr + 1] += multip * dvalpv[jj2];
407 cancel = cancel || !(fact->zeroTolerance < d1);
408 d1 = fabs(dluval[kr + 1]);
409 }
410 jj1 = maction[j1];
411 j2 = hcoli[kr + 3];
412 }
413 cancel = cancel || !(fact->zeroTolerance < d1);
414 #else
415 /*
416 * This is apparently what the above code does.
417 * In addition to being unrolled, the assignments to j[12] and jj[12]
418 * are shifted so that the result of dereferencing maction doesn't
419 * have to be used immediately afterwards for the branch test.
420 * This would would cause a pipeline delay. (The apparent dereference
421 * of hcoli will be removed by the compiler using strength reduction).
422 *
423 * loop through the entries in the row being processed,
424 * flipping the sign of the maction entries as we go along.
425 * Afterwards, we look for positive entries to see what pivot
426 * row entries will cause fill-in. We count the number of fill-ins, too.
427 * "cancel" says if the result of combining the pivot row with this one
428 * causes an entry to get too small; if so, we discard those entries.
429 */
430 kfill = epivr1 - (knpre - knprs + 1);
431 cancel = false;
432
433 for (kr = knprs; kr <= knpre; kr++) {
434 j1 = hcoli[kr];
435 jj1 = maction[j1];
436 if ((jj1 == 0)) {
437 /* no entry - this pivot column entry will have to be added */
438 ++kfill;
439 } else {
440 /* there is an entry for this column in the pivot row */
441 maction[j1] = -maction[j1];
442 dluval[kr] += multip * dvalpv[jj1];
443 d1 = fabs(dluval[kr]);
444 cancel = cancel || !(fact->zeroTolerance < d1);
445 }
446 }
447 #endif
448 kstart = knpre;
449 fill = kfill;
450
451 if (cancel) {
452 /* KSTART is used as a stack pointer for nonzeros in factored row */
453 kstart = knprs - 1;
454 for (kr = knprs; kr <= knpre; ++kr) {
455 j = hcoli[kr];
456 if (fabs(dluval[kr]) > fact->zeroTolerance) {
457 ++kstart;
458 dluval[kstart] = dluval[kr];
459 hcoli[kstart] = j;
460 } else {
461 /* Remove element from column file */
462 --nnentu;
463 --hincol[j];
464 --enpr;
465 kcs = mcstrt[j];
466 kce = kcs + hincol[j];
467 for (kk = kcs; kk <= kce; ++kk) {
468 if (hrowi[kk] == npr) {
469 hrowi[kk] = hrowi[kce];
470 hrowi[kce] = 0;
471 break;
472 }
473 }
474 /* ASSERT !(kk>kce) */
475 }
476 }
477 knpre = kstart;
478 }
479 /* Fill contains an upper bound on the amount of fill-in */
480 if (fill == 0) {
481 for (k = kipis; k <= kipie; ++k) {
482 maction[hcoli[k]] = static_cast< MACTION_T >(-maction[hcoli[k]]);
483 }
484 } else {
485 naft = mwork[npr].suc;
486 kqq = mrstrt[naft] - knpre - 1;
487
488 if (fill > kqq) {
489 /* Fill-in exceeds space left. Check if there is enough */
490 /* space in row file for the new row. */
491 nznpr = enpr + fill;
492 if (!(xnewro + nznpr + 1 < lstart)) {
493 if (!(nnentu + nznpr + 1 < lstart)) {
494 irtcod = -5;
495 goto L1050;
496 }
497 /* idea 1 is to compress every time xnewro increases by x thousand */
498 /* idea 2 is to copy nucleus rows with a reasonable gap */
499 /* then copy each row down when used */
500 /* compressions would just be 1 remainder which eventually will */
501 /* fit in cache */
502 {
503 int iput = c_ekkrwcs(fact, dluval, hcoli, mrstrt, hinrow, mwork, nfirst);
504 kmxeta += xnewro - iput;
505 xnewro = iput - 1;
506 ++ncompactions;
507 }
508
509 kipis = mrstrt[ipivot] + 1;
510 kipie = kipis + epivr1 - 1;
511 knprs = mrstrt[npr];
512 }
513
514 /* I think this assignment should be inside the previous if-stmt */
515 /* otherwise, it does nothing */
516 /*assert(knpre == knprs + enpr - 1);*/
517 knpre = knprs + enpr - 1;
518
519 /*
520 * copy this row to the end of the row file and adjust its links.
521 * The links keep track of the order of rows in memory.
522 * Rows are only moved from the middle all the way to the end.
523 */
524 if (npr != nlast) {
525 npre = mwork[npr].pre;
526 if (npr == nfirst) {
527 nfirst = naft;
528 }
529 /* take out of chain */
530 mwork[naft].pre = npre;
531 mwork[npre].suc = naft;
532 /* and put in at end */
533 mwork[nfirst].pre = npr;
534 mwork[nlast].suc = npr;
535 mwork[npr].pre = nlast;
536 mwork[npr].suc = nfirst;
537 nlast = npr;
538 kstart = xnewro;
539 mrstrt[npr] = kstart + 1;
540 nmove = knpre - knprs + 1;
541 ibase = kstart + 1 - knprs;
542 for (kr = knprs; kr <= knpre; ++kr) {
543 dluval[ibase + kr] = dluval[kr];
544 hcoli[ibase + kr] = hcoli[kr];
545 }
546 kstart += nmove;
547 } else {
548 kstart = knpre;
549 }
550
551 /* extra space ? */
552 /*
553 * The mystery of iadd32.
554 * This code assigns to xnewro, possibly using iadd32.
555 * However, in that case xnewro is assigned to just after
556 * the for-loop below, and there is no intervening reference.
557 * Therefore, I believe that this code can be entirely eliminated;
558 * it is the leftover of an interrupted or dropped experiment.
559 * Presumably, this was trying to implement the ideas about
560 * padding expressed above.
561 */
562 if (iadd32 != 0) {
563 xnewro += iadd32;
564 } else {
565 if (kstart + (nrow << 1) + 100 < lstart) {
566 ileft = ((nrow - fact->npivots + 32) & -32);
567 if (kstart + ileft * ileft + 32 < lstart) {
568 iadd32 = ileft;
569 xnewro = CoinMax(kstart, xnewro);
570 xnewro = (xnewro & -32) + ileft;
571 } else {
572 xnewro = ((kstart + 31) & -32);
573 }
574 } else {
575 xnewro = kstart;
576 }
577 }
578
579 hinrow[npr] = enpr;
580 } else if (!(nnentu + kqq + 2 < lstart)) {
581 irtcod = -5;
582 goto L1050;
583 }
584 /* Scan pivot row again to generate fill in. */
585 for (kr = kipis; kr <= kipie; ++kr) {
586 j = hcoli[kr];
587 jj = maction[j];
588 if (jj > 0) {
589 elemnt = multip * dvalpv[jj];
590 if (fabs(elemnt) > fact->zeroTolerance) {
591 ++kstart;
592 dluval[kstart] = elemnt;
593 //printf("pivot %d at %d col %d el %g\n",
594 // npr,kstart,j,elemnt);
595 hcoli[kstart] = j;
596 ++nnentu;
597 nz = hincol[j];
598 kcs = mcstrt[j];
599 kce = kcs + nz - 1;
600 if (kce == xnewco) {
601 if (xnewco + 1 >= lstart) {
602 if (xnewco + nz + 1 >= lstart) {
603 /* Compress column file */
604 if (nnentu + nz + 1 < lstart) {
605 xnewco = c_ekkclco(fact, hrowi, mcstrt, hincol, xnewco);
606 ++ncompactions;
607
608 kcpiv = mcstrt[jpivot] - 1;
609 kcs = mcstrt[j];
610 /* HINCOL MAY HAVE CHANGED? (JJHF) ??? */
611 nz = hincol[j];
612 kce = kcs + nz - 1;
613 } else {
614 irtcod = -5;
615 goto L1050;
616 }
617 }
618 /* Copy column */
619 mcstrt[j] = xnewco + 1;
620 ibase = mcstrt[j] - kcs;
621 for (kk = kcs; kk <= kce; ++kk) {
622 hrowi[ibase + kk] = hrowi[kk];
623 hrowi[kk] = 0;
624 }
625 kce = xnewco + kce - kcs + 1;
626 xnewco = kce + 1;
627 } else {
628 ++xnewco;
629 }
630 } else if (hrowi[kce + 1] != 0) {
631 /* here we use the fact that hrowi entries not "in use" are zeroed */
632 if (xnewco + nz + 1 >= lstart) {
633 /* Compress column file */
634 if (nnentu + nz + 1 < lstart) {
635 xnewco = c_ekkclco(fact, hrowi, mcstrt, hincol, xnewco);
636 ++ncompactions;
637
638 kcpiv = mcstrt[jpivot] - 1;
639 kcs = mcstrt[j];
640 /* HINCOL MAY HAVE CHANGED? (JJHF) ??? */
641 nz = hincol[j];
642 kce = kcs + nz - 1;
643 } else {
644 irtcod = -5;
645 goto L1050;
646 }
647 }
648 /* move the column to the end of the column file */
649 mcstrt[j] = xnewco + 1;
650 ibase = mcstrt[j] - kcs;
651 for (kk = kcs; kk <= kce; ++kk) {
652 hrowi[ibase + kk] = hrowi[kk];
653 hrowi[kk] = 0;
654 }
655 kce = xnewco + kce - kcs + 1;
656 xnewco = kce + 1;
657 }
658 /* store element */
659 hrowi[kce + 1] = npr;
660 hincol[j] = nz + 1;
661 }
662 } else {
663 maction[j] = static_cast< MACTION_T >(-maction[j]);
664 }
665 }
666 if (fill > kqq) {
667 xnewro = kstart;
668 }
669 }
670 hinrow[npr] = kstart - mrstrt[npr] + 1;
671 /* Check if row or column file needs compression */
672 if (!(xnewco + 1 < lstart)) {
673 xnewco = c_ekkclco(fact, hrowi, mcstrt, hincol, xnewco);
674 ++ncompactions;
675
676 kcpiv = mcstrt[jpivot] - 1;
677 }
678 if (!(xnewro + 1 < lstart)) {
679 int iput = c_ekkrwcs(fact, dluval, hcoli, mrstrt, hinrow, mwork, nfirst);
680 kmxeta += xnewro - iput;
681 xnewro = iput - 1;
682 ++ncompactions;
683
684 kipis = mrstrt[ipivot] + 1;
685 kipie = kipis + epivr1 - 1;
686 }
687 /* Store elementary row transformation */
688 ++nnentl;
689 --nnentu;
690 --lstart;
691 dluval[lstart] = multip;
692
693 hrowi[lstart] = SHIFT_INDEX(npr);
694 #define INLINE_AFPV 3
695 /* We could do this while computing values but
696 it makes it much more complex. At least we should get
697 reasonable cache behavior by doing it each row */
698 #if INLINE_AFPV
699 {
700 int j;
701 int nel, krs;
702 int koff;
703 int *index;
704 double *els;
705 nel = hinrow[npr];
706 krs = mrstrt[npr];
707 index = &hcoli[krs];
708 els = &dluval[krs];
709 #if INLINE_AFPV < 3
710 #if INLINE_AFPV == 1
711 double maxaij = 0.0;
712 koff = 0;
713 j = 0;
714 while (j < nel) {
715 double d = fabs(els[j]);
716 if (maxaij < d) {
717 maxaij = d;
718 koff = j;
719 }
720 j++;
721 }
722 #else
723 assert(nel);
724 koff = 0;
725 double maxaij = fabs(els[0]);
726 for (j = 1; j < nel; j++) {
727 double d = fabs(els[j]);
728 if (maxaij < d) {
729 maxaij = d;
730 koff = j;
731 }
732 }
733 #endif
734 #else
735 double maxaij = 0.0;
736 koff = 0;
737 j = 0;
738 if ((nel & 1) != 0) {
739 maxaij = fabs(els[0]);
740 j = 1;
741 }
742
743 while (j < nel) {
744 UNROLL_LOOP_BODY2({
745 double d = fabs(els[j]);
746 if (maxaij < d) {
747 maxaij = d;
748 koff = j;
749 }
750 j++;
751 });
752 }
753 #endif
754 SWAP(int, index[koff], index[0]);
755 SWAP(double, els[koff], els[0]);
756 }
757 #endif
758
759 {
760 int nzi = hinrow[npr];
761 if (nzi > 0) {
762 C_EKK_ADD_LINK(hpivro, nzi, rlink, npr);
763 }
764 }
765 }
766
767 /* after pivot move biggest to first in each row */
768 #if INLINE_AFPV == 0
769 int nn = mcstrt[fact->xnetal] - lstart + 1;
770 c_ekkafpv(hrowi + lstart, hcoli, dluval, mrstrt, hinrow, nn);
771 #endif
772
773 /* Restore work array */
774 for (k = kipis; k <= kipie; ++k) {
775 maction[hcoli[k]] = 0;
776 }
777
778 if (*xrejctp > 0) {
779 for (k = kipis; k <= kipie; ++k) {
780 int j = hcoli[k];
781 int nzj = hincol[j];
782 if (!(nzj <= 0) && !((clink[j].pre > nrow && nzj != 1))) {
783 C_EKK_ADD_LINK(hpivco, nzj, clink, j);
784 }
785 }
786 } else {
787 for (k = kipis; k <= kipie; ++k) {
788 int j = hcoli[k];
789 int nzj = hincol[j];
790 if (!(nzj <= 0)) {
791 C_EKK_ADD_LINK(hpivco, nzj, clink, j);
792 }
793 }
794 }
795 fact->nuspike += hinrow[ipivot];
796
797 /* Go to dense coding if appropriate */
798 #ifndef C_EKKCMFY
799 if (ifdens != 0) {
800 int ndense = nrow - fact->npivots;
801 if (!(xnewro + ndense * ndense >= lstart)) {
802
803 /* set up sort order in MACTION */
804 c_ekkizero(nrow, reinterpret_cast< int * >(maction + 1));
805 iput = 0;
806 for (i = 1; i <= nrow; ++i) {
807 if (clink[i].pre >= 0) {
808 ++iput;
809 maction[i] = static_cast< short int >(iput);
810 }
811 }
812 /* and get number spare needed */
813 nspare = 0;
814 for (i = 1; i <= nrow; ++i) {
815 if (rlink[i].pre >= 0) {
816 nspare = nspare + ndense - hinrow[i];
817 }
818 }
819 if (iput != nrow - fact->npivots) {
820 /* must be singular */
821 c_ekkizero(nrow, reinterpret_cast< int * >(maction + 1));
822 } else {
823 /* pack down then back up */
824 int iput = c_ekkrwcs(fact, dluval, hcoli, mrstrt, hinrow, mwork, nfirst);
825 kmxeta += xnewro - iput;
826 xnewro = iput - 1;
827 ++ncompactions;
828
829 --ncompactions;
830 if (xnewro + nspare + ndense * ndense >= lstart) {
831 c_ekkizero(nrow, reinterpret_cast< int * >(maction + 1));
832 } else {
833 xnewro += nspare;
834 c_ekkrwct(fact, dluval, hcoli, mrstrt, hinrow, mwork,
835 rlink, maction, dvalpv,
836 nlast, xnewro);
837 kmxeta += xnewro;
838 if (nnentu + nnentl > nrow * 5 && (ndense * ndense) > (nnentu + nnentl) >> 2 && !if_sparse_update) {
839 fact->ndenuc = ndense;
840 }
841 irtcod = c_ekkcmfd(fact,
842 (reinterpret_cast< int * >(dvalpv) + 1),
843 rlink, clink,
844 (reinterpret_cast< int * >(maction + 1)) + 1,
845 nnetas,
846 &nnentl, &nnentu,
847 nsingp);
848 /* irtcod == 0 || irtcod == 10 */
849 /* 10 == found 0.0 pivot */
850 goto L1050;
851 }
852 }
853 } else {
854 /* say not enough room */
855 /*printf("no room %d\n",ndense);*/
856 if (1) {
857 /* return and increase size of etas if possible */
858 if (!noRoomForDense) {
859 int etasize = CoinMax(4 * fact->nnentu + (nnetas - fact->nnentl) + 1000, fact->eta_size);
860 noRoomForDense = ndense;
861 fact->eta_size = CoinMin(static_cast< int >(1.2 * fact->eta_size), etasize);
862 if (fact->maxNNetas > 0 && fact->eta_size > fact->maxNNetas) {
863 fact->eta_size = fact->maxNNetas;
864 }
865 }
866 }
867 }
868 }
869 #endif /* C_EKKCMFY */
870 }
871
872 L1050 : {
873 int iput = c_ekkrwcs(fact, dluval, hcoli, mrstrt, hinrow, mwork, nfirst);
874 kmxeta += xnewro - iput;
875 xnewro = iput - 1;
876 ++ncompactions;
877 }
878
879 nnentu = xnewro;
880 /* save order of row copy for c_ekkshfv */
881 mwork[nrow + 1].pre = nfirst;
882 mwork[nrow + 1].suc = nlast;
883
884 fact->nnentl = nnentl;
885 fact->nnentu = nnentu;
886 fact->kmxeta = kmxeta;
887 *xnewrop = xnewro;
888 *ncompactionsp = ncompactions;
889
890 return (irtcod);
891 } /* c_ekkcmfc */
892 #endif
893
894 /* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
895 */
896