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