Cgl  0.60.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CoinOslC.h
Go to the documentation of this file.
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);
61 inline void c_ekkscpy(int n, const int *marr1, int *marr2)
62 {
63  CoinMemcpyN(marr1, n, marr2);
64 }
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);
74 inline void c_ekkdzero(int n, double *marray)
75 {
76  CoinZeroN(marray, n);
77 }
78 inline void c_ekkizero(int n, int *marray)
79 {
80  CoinZeroN(marray, n);
81 }
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)
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) {
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) {
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 */
void clp_free(void *oldArray)
void c_ekkscpy(int n, const int *marr1, int *marr2)
Definition: CoinOslC.h:61
#define SWAP(type, _x, _y)
Definition: CoinOslC.h:146
int * clp_int(int number_entries)
int c_ekketsj(registerEKKfactinfo *fact, double *dwork1, int *mpt2, double dalpha, int orig_nincol, int npivot, int *nuspikp, const int ipivrw, int *spare)
#define UNROLL_LOOP_BODY2(code)
Definition: CoinOslC.h:159
void c_ekk_Set(int *array, int bit)
void clp_setup_pointers(EKKfactinfo *fact)
int c_ekkbtrn(register const EKKfactinfo *fact, double *dwork1, int *mpt, int first_nonzero)
int c_ekkftrn(register const EKKfactinfo *fact, double *dwork1, double *dpermu, int *mpt, int numberNonZero)
void CoinMemcpyN(const T *from, const CoinBigIndex size, T *to)
This helper function copies an array to another location.
double * clp_double(int number_entries)
int c_ekkbtrn_ipivrw(register const EKKfactinfo *fact, double *dwork1, int *mpt, int ipivrw, int *spare)
#define SHIFT_INDEX(limit)
Definition: CoinOslC.h:128
void * clp_malloc(int number_entries)
void c_ekk_Unset(int *array, int bit)
int c_ekklfct(register EKKfactinfo *fact)
void c_ekkzero(int length, int n, void *array)
void c_ekkczero(int n, char *marray)
Definition: CoinOslC.h:82
void CoinZeroN(T *to, const CoinBigIndex size)
This helper function fills an array with zero.
T CoinMax(const T x1, const T x2)
Return the larger (according to operator&lt;() of the arguments.
void c_ekkdcpy(int n, const double *marr1, double *marr2)
Definition: CoinOslC.h:65
void c_ekkftrn2(register EKKfactinfo *fact, double *dwork1, double *dpermu1, int *mpt1, int *nincolp, double *dwork1_ft, int *mpt_ft, int *nincolp_ft)
int c_ekkftrn_ft(register EKKfactinfo *fact, double *dwork1, int *mpt, int *nincolp)
void c_ekkizero(int n, int *marray)
Definition: CoinOslC.h:78
void clp_memory(int type)
T CoinMin(const T x1, const T x2)
Return the smaller (according to operator&lt;() of the arguments.
#define UNROLL_LOOP_BODY4(code)
Definition: CoinOslC.h:165
int c_ekk_IsSet(const int *array, int bit)
int c_ekkslcf(register const EKKfactinfo *fact)
#define C_EKK_ADD_LINK(hpiv, nzi, link, npr)
Definition: CoinOslC.h:115
void c_ekkdzero(int n, double *marray)
Definition: CoinOslC.h:74