Cbc  2.9.9
 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 1582 2013-04-06 14:33:07Z 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 { CoinMemcpyN(marr1,n,marr2);}
63 inline void c_ekkdcpy(int n, const double *marr1,double *marr2)
64 { CoinMemcpyN(marr1,n,marr2);}
65 int c_ekk_IsSet(const int * array,int bit);
66 void c_ekk_Set(int * array,int bit);
67 void c_ekk_Unset(int * array,int bit);
68 
69 void c_ekkzero(int length, int n, void * array);
70 inline void c_ekkdzero(int n, double *marray)
71 {CoinZeroN(marray,n);}
72 inline void c_ekkizero(int n, int *marray)
73 {CoinZeroN(marray,n);}
74 inline 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)
82 void clp_setup_pointers(EKKfactinfo * fact);
83 void clp_memory(int type);
84 double * clp_double(int number_entries);
85 int * clp_int(int number_entries);
86 void * clp_malloc(int number_entries);
87 void 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, 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  int ifdens = 0;
196 #endif
197  int irtcod = 0;
198  const int nrow = fact->nrow;
199 
200  /* Parameter adjustments */
201  --maction;
202 
203  /* Function Body */
204  lstart = nnetas - nnentl + 1;
205  for (i = lstart; i <= nnetas; ++i) {
206  hrowi[i] = SHIFT_INDEX(hcoli[i]);
207  }
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) {
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) {
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 
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:136
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:140
void c_ekk_Set(int *array, int bit)
void CoinMemcpyN(register const T *from, const int size, register T *to)
This helper function copies an array to another location.
T CoinMin(register const T x1, register const T x2)
Return the smaller (according to operator&lt;() of the arguments.
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)
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:118
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:74
void c_ekkdcpy(int n, const double *marr1, double *marr2)
Definition: CoinOslC.h:63
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 CoinZeroN(register T *to, const int size)
This helper function fills an array with zero.
void c_ekkizero(int n, int *marray)
Definition: CoinOslC.h:72
void clp_memory(int type)
T CoinMax(register const T x1, register const T x2)
Return the larger (according to operator&lt;() of the arguments.
#define UNROLL_LOOP_BODY4(code)
Definition: CoinOslC.h:142
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:105
void c_ekkdzero(int n, double *marray)
Definition: CoinOslC.h:70