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