- Main Page
- Modules
- Namespaces
- Classes
- Files
- Related Pages

- File List
- File Members

# /home/coin/SVN-release/Cgl-0.55.0/CoinUtils/src/CoinOslC.h

Go to the documentation of this file.

    00001 /* $Id: CoinOslC.h 1202 2009-09-25 15:27:15Z forrest $ */
    00002 #ifndef COIN_OSL_C_INCLUDE
    00003 /* Copyright (C) 1987, 2009, International Business Machines
    00004    Corporation and others.  All Rights Reserved. */
    00005 #define COIN_OSL_C_INCLUDE
    00006 typedef struct {int suc, pre;} EKKHlink;
    00007 typedef struct _EKKfactinfo {
    00008   double drtpiv;
    00009   double demark;
    00010   double zpivlu;
    00011   double zeroTolerance;
    00012   double areaFactor;
    00013   int *xrsadr;
    00014   int *xcsadr;
    00015   int *xrnadr;
    00016   int *xcnadr;
    00017   int *krpadr;
    00018   int *kcpadr;
    00019   int *mpermu;
    00020   int *bitArray;
    00021   int * back;
    00022   char * nonzero;
    00023   double * trueStart;
    00024   mutable double *kadrpm;
    00025   int *R_etas_index;
    00026   int *R_etas_start;
    00027   double *R_etas_element;
    00028
    00029   int *xecadr;
    00030   int *xeradr;
    00031   double *xeeadr;
    00032   double *xe2adr;
    00033   EKKHlink * kp1adr;
    00034   EKKHlink * kp2adr;
    00035   double * kw1adr;
    00036   double * kw2adr;
    00037   double * kw3adr;
    00038   int * hpivcoR;
    00039   int nrow;
    00040   int nrowmx;
    00041   int firstDoRow;
    00042   int firstLRow;
    00043   int maxinv;
    00044   int nnetas;
    00045   int iterin;
    00046   int iter0;
    00047   int invok;
    00048   int nbfinv;
    00049   int num_resets;
    00050   int nnentl;
    00051   int nnentu;
    00052   int ndenuc;
    00053   int npivots; /* use as xpivsq in factorization */
    00054   int kmxeta;
    00055   int xnetal;
    00056   int first_dense;
    00057   int last_dense;
    00058   int iterno;
    00059   int numberSlacks;
    00060   int lastSlack;
    00061   int firstNonSlack;
    00062   int xnetalval;
    00063   int lstart;
    00064   int if_sparse_update;
    00065   mutable int packedMode;
    00066   int switch_off_sparse_update;
    00067   int nuspike;
    00068   bool rows_ok; /* replaces test using mrstrt[1] */
    00069
    00070   int nR_etas;
    00071   int sortedEta; /* if vector for F-T is sorted */
    00072   int lastEtaCount;
    00073   int ifvsol;
    00074   int eta_size;
    00075   int last_eta_size;
    00076   int maxNNetas;
    00077 } EKKfactinfo;
    00078 #ifndef CLP_OSL
    00079 #define CLP_OSL 0
    00080 #endif
    00081 #define C_EKK_GO_SPARSE 200
    00082
    00083 #ifdef HAVE_ENDIAN_H
    00084 #include <endian.h>
    00085 #if __BYTE_ORDER == __LITTLE_ENDIAN
    00086 #define INTEL
    00087 #endif
    00088 #endif
    00089
    00090 #include <math.h>
    00091 #include <string.h>
    00092 #include <stdio.h>
    00093 #include <stdlib.h>
    00094
    00095 #define SPARSE_UPDATE
    00096 #define NO_SHIFT
    00097 #include "CoinHelperFunctions.hpp"
    00098
    00099 #include <stddef.h>
    00100 #ifdef __cplusplus
    00101 extern "C"{
    00102 #endif
    00103
    00104 int c_ekkbtrn( register const EKKfactinfo *fact,
    00105             double *dwork1,
    00106             int * mpt,int first_nonzero);
    00107 int c_ekkbtrn_ipivrw( register const EKKfactinfo *fact,
    00108                    double *dwork1,
    00109                    int * mpt, int ipivrw,int * spare);
    00110
    00111 int c_ekketsj( register /*const*/ EKKfactinfo *fact,
    00112             double *dwork1,
    00113             int *mpt2, double dalpha, int orig_nincol,
    00114             int npivot, int *nuspikp,
    00115             const int ipivrw, int * spare);
    00116 int c_ekkftrn( register const EKKfactinfo *fact,
    00117             double *dwork1,
    00118             double * dpermu,int * mpt, int numberNonZero);
    00119
    00120 int c_ekkftrn_ft( register EKKfactinfo *fact,
    00121                double *dwork1, int *mpt, int *nincolp);
    00122 void c_ekkftrn2( register EKKfactinfo *fact, double *dwork1,
    00123               double * dpermu1,int * mpt1, int *nincolp,
    00124              double *dwork1_ft, int *mpt_ft, int *nincolp_ft);
    00125
    00126 int c_ekklfct( register EKKfactinfo *fact);
    00127 int c_ekkslcf( register const EKKfactinfo *fact);
    00128 inline void c_ekkscpy(int n, const int *marr1,int *marr2)
    00129 { CoinMemcpyN(marr1,n,marr2);}
    00130 inline void c_ekkdcpy(int n, const double *marr1,double *marr2)
    00131 { CoinMemcpyN(marr1,n,marr2);}
    00132 int c_ekk_IsSet(const int * array,int bit);
    00133 void c_ekk_Set(int * array,int bit);
    00134 void c_ekk_Unset(int * array,int bit);
    00135
    00136 void c_ekkzero(int length, int n, void * array);
    00137 inline void c_ekkdzero(int n, double *marray)
    00138 {CoinZeroN(marray,n);}
    00139 inline void c_ekkizero(int n, int *marray)
    00140 {CoinZeroN(marray,n);}
    00141 inline void c_ekkczero(int n, char *marray)
    00142 {CoinZeroN(marray,n);}
    00143 #ifdef __cplusplus
    00144           }
    00145 #endif
    00146
    00147 #define c_ekkscpy_0_1(s,ival,array) CoinFillN(array,s,ival)
    00148 #define c_ekks1cpy( n,marr1,marr2)  CoinMemcpyN(marr1,n, marr2)
    00149 void clp_setup_pointers(EKKfactinfo * fact);
    00150 void clp_memory(int type);
    00151 double * clp_double(int number_entries);
    00152 int * clp_int(int number_entries);
    00153 void * clp_malloc(int number_entries);
    00154 void clp_free(void * oldArray);
    00155
    00156 #define SLACK_VALUE -1.0
    00157 #define C_EKK_REMOVE_LINK(hpiv,hin,link,ipivot) \
    00158   {                                             \
    00159     int ipre = link[ipivot].pre;                \
    00160     int isuc = link[ipivot].suc;                \
    00161     if (ipre > 0) {                             \
    00162       link[ipre].suc = isuc;                    \
    00163     }                                           \
    00164     if (ipre <= 0) {                            \
    00165       hpiv[hin[ipivot]] = isuc;                 \
    00166     }                                           \
    00167     if (isuc > 0) {                             \
    00168       link[isuc].pre = ipre;                    \
    00169     }                                           \
    00170   }
    00171
    00172 #define C_EKK_ADD_LINK(hpiv,nzi,link, npr)      \
    00173   {                                             \
    00174     int ifiri = hpiv[nzi];                      \
    00175     hpiv[nzi] = npr;                            \
    00176     link[npr].suc = ifiri;                      \
    00177     link[npr].pre = 0;                          \
    00178     if (ifiri != 0) {                           \
    00179       link[ifiri].pre = npr;                    \
    00180     }                                           \
    00181   }
    00182 #include <assert.h>
    00183 #ifdef  NO_SHIFT
    00184
    00185 #define SHIFT_INDEX(limit)      (limit)
    00186 #define UNSHIFT_INDEX(limit)    (limit)
    00187 #define SHIFT_REF(arr,ind)      (arr)[ind]
    00188
    00189 #else
    00190
    00191 #define SHIFT_INDEX(limit)      ((limit)<<3)
    00192 #define UNSHIFT_INDEX(limit)    ((unsigned int)(limit)>>3)
    00193 #define SHIFT_REF(arr,ind)      (*(double*)((char*)(arr) + (ind)))
    00194
    00195 #endif
    00196
    00197 #ifdef INTEL
    00198 #define NOT_ZERO(x)     (((*((reinterpret_cast<unsigned char *>(&x))+7)) & 0x7F) != 0)
    00199 #else
    00200 #define NOT_ZERO(x)     ((x) != 0.0)
    00201 #endif
    00202
    00203 #define SWAP(type,_x,_y)        { type _tmp = (_x); (_x) = (_y); (_y) = _tmp;}
    00204
    00205 #define UNROLL_LOOP_BODY1(code)                 \
    00206   {{code}}
    00207 #define UNROLL_LOOP_BODY2(code)                 \
    00208   {{code} {code}}
    00209 #define UNROLL_LOOP_BODY4(code)                 \
    00210   {{code} {code} {code} {code}}
    00211 #endif
    00212 #ifdef COIN_OSL_CMFC
    00213 /*     Return codes in IRTCOD/IRTCOD are */
    00214 /*     4: numerical problems */
    00215 /*     5: not enough space in row file */
    00216 /*     6: not enough space in column file */
    00217 /*    23: system error at label 320 */
    00218 {
    00219 #if 1
    00220   int *hcoli    = fact->xecadr;
    00221   double *dluval        = fact->xeeadr;
    00222   double *dvalpv = fact->kw3adr;
    00223   int *mrstrt   = fact->xrsadr;
    00224   int *hrowi    = fact->xeradr;
    00225   int *mcstrt   = fact->xcsadr;
    00226   int *hinrow   = fact->xrnadr;
    00227   int *hincol   = fact->xcnadr;
    00228   int *hpivro   = fact->krpadr;
    00229   int *hpivco   = fact->kcpadr;
    00230 #endif
    00231   int nnentl    = fact->nnentl;
    00232   int nnentu    = fact->nnentu;
    00233   int kmxeta    = fact->kmxeta;
    00234   int xnewro    = *xnewrop;
    00235   int ncompactions      = *ncompactionsp;
    00236
    00237   MACTION_T *maction = reinterpret_cast<MACTION_T*>(maction_void);
    00238
    00239   int i, j, k;
    00240   double d1;
    00241   int j1, j2;
    00242   int jj, kk, kr, nz, jj1, jj2, kce, kcs, kqq, npr;
    00243   int fill, naft;
    00244   int enpr;
    00245   int nres, npre;
    00246   int knpr, irow, iadd32, ibase;
    00247   double pivot;
    00248   int count, nznpr;
    00249   int nlast, epivr1;
    00250   int kipis;
    00251   double dpivx;
    00252   int kipie, kcpiv, knprs, knpre;
    00253   bool cancel;
    00254   double multip, elemnt;
    00255   int ipivot, jpivot, epivro, epivco, lstart, ifdens, nfirst;
    00256   int nzpivj, kfill, kstart;
    00257   int nmove, ileft;
    00258 #ifndef C_EKKCMFY
    00259   int iput, nspare;
    00260   int noRoomForDense=0;
    00261   int if_sparse_update=fact->if_sparse_update;
    00262 #endif
    00263   int irtcod    = 0;
    00264   const int nrow        = fact->nrow;
    00265
    00266   /* Parameter adjustments */
    00267   --maction;
    00268
    00269   /* Function Body */
    00270   lstart = nnetas - nnentl + 1;
    00271   for (i = lstart; i <= nnetas; ++i) {
    00272       hrowi[i] = SHIFT_INDEX(hcoli[i]);
    00273   }
    00274   ifdens = 0;
    00275
    00276   for (i = 1; i <= nrow; ++i) {
    00277     maction[i] = 0;
    00278     mwork[i].pre = i - 1;
    00279     mwork[i].suc = i + 1;
    00280   }
    00281
    00282   iadd32 = 0;
    00283   nlast = nrow;
    00284   nfirst = 1;
    00285   mwork[1].pre = nrow;
    00286   mwork[nrow].suc = 1;
    00287
    00288   for (count = 1; count <= nrow; ++count) {
    00289
    00290     /* Pick column singletons */
    00291     if (! (hpivco[1] <= 0)) {
    00292       int small_pivot = c_ekkcsin(fact,
    00293                                  rlink, clink,
    00294                                     nsingp);
    00295
    00296       if (small_pivot) {
    00297         irtcod = 7; /* pivot too small */
    00298         if (fact->invok >= 0) {
    00299           goto L1050;
    00300         }
    00301       }
    00302       if (fact->npivots >= nrow) {
    00303         goto L1050;
    00304       }
    00305     }
    00306
    00307     /* Pick row singletons */
    00308     if (! (hpivro[1] <= 0)) {
    00309       irtcod = c_ekkrsin(fact,
    00310                          rlink, clink,
    00311                          mwork,nfirst,
    00312                          nsingp,
    00313
    00314                      &xnewco, &xnewro,
    00315                      &nnentu,
    00316                      &kmxeta, &ncompactions,
    00317                          &nnentl);
    00318         if (irtcod != 0) {
    00319           if (irtcod < 0 || fact->invok >= 0) {
    00320             /* -5 */
    00321             goto L1050;
    00322           }
    00323           /* ASSERT:  irtcod == 7 - pivot too small */
    00324           /* why don't we return with an error? */
    00325         }
    00326         if (fact->npivots >= nrow) {
    00327             goto L1050;
    00328         }
    00329         lstart = nnetas - nnentl + 1;
    00330     }
    00331
    00332     /* Find a pivot element */
    00333     irtcod = c_ekkfpvt(fact,
    00334                       rlink, clink,
    00335                      nsingp, xrejctp, &ipivot, &jpivot);
    00336     if (irtcod != 0) {
    00337       /* irtcod == 10 */
    00338         goto L1050;
    00339     }
    00340     /*        Update list structures and prepare for numerical phase */
    00341     c_ekkprpv(fact, rlink, clink,
    00342                      *xrejctp, ipivot, jpivot);
    00343
    00344     epivco = hincol[jpivot];
    00345     ++fact->xnetal;
    00346     mcstrt[fact->xnetal] = lstart - 1;
    00347     hpivco[fact->xnetal] = ipivot;
    00348     epivro = hinrow[ipivot];
    00349     epivr1 = epivro - 1;
    00350     kipis = mrstrt[ipivot];
    00351     pivot = dluval[kipis];
    00352     dpivx = 1. / pivot;
    00353     kipie = kipis + epivr1;
    00354     ++kipis;
    00355 #ifndef C_EKKCMFY
    00356     {
    00357       double size = nrow - fact->npivots;
    00358       if (size > GO_DENSE && (nnentu - fact->nuspike) * GO_DENSE_RATIO > size * size) {
    00359         /* say going to dense coding */
    00360         if (*nsingp == 0) {
    00361           ifdens = 1;
    00362         }
    00363       }
    00364     }
    00365 #endif
    00366     /* copy the pivot row entries into dvalpv */
    00367     /* the maction array tells us the index into dvalpv for a given row */
    00368     /* the alternative would be using a large array of doubles */
    00369     for (k = kipis; k <= kipie; ++k) {
    00370       irow = hcoli[k];
    00371       dvalpv[k - kipis + 1] = dluval[k];
    00372       maction[irow] = static_cast<MACTION_T>(k - kipis + 1);
    00373     }
    00374
    00375     /* Loop over nonzeros in pivot column */
    00376     kcpiv = mcstrt[jpivot] - 1;
    00377     for (nzpivj = 1; nzpivj <= epivco; ++nzpivj) {
    00378       ++kcpiv;
    00379       npr = hrowi[kcpiv];
    00380       hrowi[kcpiv] = 0; /* zero out for possible compaction later on */
    00381
    00382       --hincol[jpivot];
    00383
    00384       ++mcstrt[jpivot];
    00385       /* loop invariant:  kcpiv == mcstrt[jpivot] - 1 */
    00386
    00387       --hinrow[npr];
    00388       enpr = hinrow[npr];
    00389       knprs = mrstrt[npr];
    00390       knpre = knprs + enpr;
    00391
    00392       /* Search for element to be eliminated */
    00393       knpr = knprs;
    00394       while (1) {
    00395           UNROLL_LOOP_BODY4({
    00396             if (jpivot == hcoli[knpr]) {
    00397               break;
    00398             }
    00399             knpr++;
    00400           });
    00401       }
    00402
    00403       multip = -dluval[knpr] * dpivx;
    00404
    00405       /* swap last entry with pivot */
    00406       dluval[knpr] = dluval[knpre];
    00407       hcoli[knpr] = hcoli[knpre];
    00408       --knpre;
    00409
    00410 #if     1
    00411       /* MONSTER_UNROLLED_CODE - see below */
    00412       kfill = epivr1 - (knpre - knprs + 1);
    00413       nres = ((knpre - knprs + 1) & 1) + knprs;
    00414       cancel = false;
    00415       d1 = 1e33;
    00416       j1 = hcoli[nres];
    00417
    00418       if (nres != knprs) {
    00419         j = hcoli[knprs];
    00420         if (maction[j] == 0) {
    00421           ++kfill;
    00422         } else {
    00423           jj = maction[j];
    00424           maction[j] = static_cast<MACTION_T>(-maction[j]);
    00425           dluval[knprs] += multip * dvalpv[jj];
    00426           d1 = fabs(dluval[knprs]);
    00427         }
    00428       }
    00429       j2 = hcoli[nres + 1];
    00430       jj1 = maction[j1];
    00431       for (kr = nres; kr < knpre; kr += 2) {
    00432         jj2 = maction[j2];
    00433         if ( (jj1 == 0)) {
    00434           ++kfill;
    00435         } else {
    00436           maction[j1] = static_cast<MACTION_T>(-maction[j1]);
    00437           dluval[kr] += multip * dvalpv[jj1];
    00438           cancel = cancel || ! (fact->zeroTolerance < d1);
    00439           d1 = fabs(dluval[kr]);
    00440         }
    00441         j1 = hcoli[kr + 2];
    00442         if ( (jj2 == 0)) {
    00443           ++kfill;
    00444         } else {
    00445           maction[j2] = static_cast<MACTION_T>(-maction[j2]);
    00446           dluval[kr + 1] += multip * dvalpv[jj2];
    00447           cancel = cancel || ! (fact->zeroTolerance < d1);
    00448           d1 = fabs(dluval[kr + 1]);
    00449         }
    00450         jj1 = maction[j1];
    00451         j2 = hcoli[kr + 3];
    00452       }
    00453       cancel = cancel || ! (fact->zeroTolerance < d1);
    00454 #else
    00455       /*
    00456        * This is apparently what the above code does.
    00457        * In addition to being unrolled, the assignments to j[12] and jj[12]
    00458        * are shifted so that the result of dereferencing maction doesn't
    00459        * have to be used immediately afterwards for the branch test.
    00460        * This would would cause a pipeline delay.  (The apparent dereference
    00461        * of hcoli will be removed by the compiler using strength reduction).
    00462        *
    00463        * loop through the entries in the row being processed,
    00464        * flipping the sign of the maction entries as we go along.
    00465        * Afterwards, we look for positive entries to see what pivot
    00466        * row entries will cause fill-in.  We count the number of fill-ins, too.
    00467        * "cancel" says if the result of combining the pivot row with this one
    00468        * causes an entry to get too small; if so, we discard those entries.
    00469        */
    00470       kfill = epivr1 - (knpre - knprs + 1);
    00471       cancel = false;
    00472
    00473       for (kr = knprs; kr <= knpre; kr++) {
    00474         j1 = hcoli[kr];
    00475         jj1 = maction[j1];
    00476         if ( (jj1 == 0)) {
    00477           /* no entry - this pivot column entry will have to be added */
    00478           ++kfill;
    00479         } else {
    00480           /* there is an entry for this column in the pivot row */
    00481           maction[j1] = -maction[j1];
    00482           dluval[kr] += multip * dvalpv[jj1];
    00483           d1 = fabs(dluval[kr]);
    00484           cancel = cancel || ! (fact->zeroTolerance < d1);
    00485         }
    00486       }
    00487 #endif
    00488       kstart = knpre;
    00489       fill = kfill;
    00490
    00491       if (cancel) {
    00492         /* KSTART is used as a stack pointer for nonzeros in factored row */
    00493         kstart = knprs - 1;
    00494         for (kr = knprs; kr <= knpre; ++kr) {
    00495           j = hcoli[kr];
    00496           if (fabs(dluval[kr]) > fact->zeroTolerance) {
    00497             ++kstart;
    00498             dluval[kstart] = dluval[kr];
    00499             hcoli[kstart] = j;
    00500           } else {
    00501             /* Remove element from column file */
    00502             --nnentu;
    00503             --hincol[j];
    00504             --enpr;
    00505             kcs = mcstrt[j];
    00506             kce = kcs + hincol[j];
    00507             for (kk = kcs; kk <= kce; ++kk) {
    00508               if (hrowi[kk] == npr) {
    00509                 hrowi[kk] = hrowi[kce];
    00510                 hrowi[kce] = 0;
    00511                 break;
    00512               }
    00513             }
    00514             /* ASSERT !(kk>kce) */
    00515           }
    00516         }
    00517         knpre = kstart;
    00518       }
    00519       /* Fill contains an upper bound on the amount of fill-in */
    00520       if (fill == 0) {
    00521         for (k = kipis; k <= kipie; ++k) {
    00522           maction[hcoli[k]] = static_cast<MACTION_T>(-maction[hcoli[k]]);
    00523         }
    00524       }
    00525       else {
    00526         naft = mwork[npr].suc;
    00527         kqq = mrstrt[naft] - knpre - 1;
    00528
    00529         if (fill > kqq) {
    00530           /* Fill-in exceeds space left. Check if there is enough */
    00531           /* space in row file for the new row. */
    00532           nznpr = enpr + fill;
    00533           if (! (xnewro + nznpr + 1 < lstart)) {
    00534             if (! (nnentu + nznpr + 1 < lstart)) {
    00535               irtcod = -5;
    00536               goto L1050;
    00537             }
    00538             /* idea 1 is to compress every time xnewro increases by x thousand */
    00539             /* idea 2 is to copy nucleus rows with a reasonable gap */
    00540             /* then copy each row down when used */
    00541             /* compressions would just be 1 remainder which eventually will */
    00542             /* fit in cache */
    00543             {
    00544               int iput = c_ekkrwcs(fact,dluval, hcoli, mrstrt, hinrow, mwork, nfirst);
    00545               kmxeta += xnewro - iput ;
    00546               xnewro = iput - 1;
    00547               ++ncompactions;
    00548             }
    00549
    00550             kipis = mrstrt[ipivot] + 1;
    00551             kipie = kipis + epivr1 - 1;
    00552             knprs = mrstrt[npr];
    00553           }
    00554
    00555           /* I think this assignment should be inside the previous if-stmt */
    00556           /* otherwise, it does nothing */
    00557           /*assert(knpre == knprs + enpr - 1);*/
    00558           knpre = knprs + enpr - 1;
    00559
    00560           /*
    00561            * copy this row to the end of the row file and adjust its links.
    00562            * The links keep track of the order of rows in memory.
    00563            * Rows are only moved from the middle all the way to the end.
    00564            */
    00565           if (npr != nlast) {
    00566             npre = mwork[npr].pre;
    00567             if (npr == nfirst) {
    00568               nfirst = naft;
    00569             }
    00570             /*             take out of chain */
    00571             mwork[naft].pre = npre;
    00572             mwork[npre].suc = naft;
    00573             /*             and put in at end */
    00574             mwork[nfirst].pre = npr;
    00575             mwork[nlast].suc = npr;
    00576             mwork[npr].pre = nlast;
    00577             mwork[npr].suc = nfirst;
    00578             nlast = npr;
    00579             kstart = xnewro;
    00580             mrstrt[npr] = kstart + 1;
    00581             nmove = knpre - knprs + 1;
    00582             ibase = kstart + 1 - knprs;
    00583             for (kr = knprs; kr <= knpre; ++kr) {
    00584               dluval[ibase + kr] = dluval[kr];
    00585               hcoli[ibase + kr] = hcoli[kr];
    00586             }
    00587             kstart += nmove;
    00588           } else {
    00589             kstart = knpre;
    00590           }
    00591
    00592           /* extra space ? */
    00593           /*
    00594            * The mystery of iadd32.
    00595            * This code assigns to xnewro, possibly using iadd32.
    00596            * However, in that case xnewro is assigned to just after
    00597            * the for-loop below, and there is no intervening reference.
    00598            * Therefore, I believe that this code can be entirely eliminated;
    00599            * it is the leftover of an interrupted or dropped experiment.
    00600            * Presumably, this was trying to implement the ideas about
    00601            * padding expressed above.
    00602            */
    00603           if (iadd32 != 0) {
    00604             xnewro += iadd32;
    00605           } else {
    00606             if (kstart + (nrow << 1) + 100 < lstart) {
    00607               ileft = ((nrow - fact->npivots + 32) & -32);
    00608               if (kstart + ileft * ileft + 32 < lstart) {
    00609                 iadd32 = ileft;
    00610                 xnewro = CoinMax(kstart,xnewro);
    00611                 xnewro = (xnewro & -32) + ileft;
    00612               } else {
    00613                 xnewro = ((kstart + 31) & -32);
    00614               }
    00615             } else {
    00616               xnewro = kstart;
    00617             }
    00618           }
    00619
    00620           hinrow[npr] = enpr;
    00621         } else if (! (nnentu + kqq + 2 < lstart)) {
    00622           irtcod = -5;
    00623           goto L1050;
    00624         }
    00625         /* Scan pivot row again to generate fill in. */
    00626         for (kr = kipis; kr <= kipie; ++kr) {
    00627           j = hcoli[kr];
    00628           jj = maction[j];
    00629           if (jj >0) {
    00630             elemnt = multip * dvalpv[jj];
    00631             if (fabs(elemnt) > fact->zeroTolerance) {
    00632               ++kstart;
    00633               dluval[kstart] = elemnt;
    00634               //printf("pivot %d at %d col %d el %g\n",
    00635               // npr,kstart,j,elemnt);
    00636               hcoli[kstart] = j;
    00637               ++nnentu;
    00638               nz = hincol[j];
    00639               kcs = mcstrt[j];
    00640               kce = kcs + nz - 1;
    00641               if (kce == xnewco) {
    00642                 if (xnewco + 1 >= lstart) {
    00643                   if (xnewco + nz + 1 >= lstart) {
    00644                     /*                  Compress column file */
    00645                     if (nnentu + nz + 1 < lstart) {
    00646                       xnewco = c_ekkclco(fact,hrowi, mcstrt, hincol, xnewco);
    00647                       ++ncompactions;
    00648
    00649                       kcpiv = mcstrt[jpivot] - 1;
    00650                       kcs = mcstrt[j];
    00651                       /*                  HINCOL MAY HAVE CHANGED? (JJHF) ??? */
    00652                       nz = hincol[j];
    00653                       kce = kcs + nz - 1;
    00654                     } else {
    00655                       irtcod = -5;
    00656                       goto L1050;
    00657                     }
    00658                   }
    00659                   /*              Copy column */
    00660                   mcstrt[j] = xnewco + 1;
    00661                   ibase = mcstrt[j] - kcs;
    00662                   for (kk = kcs; kk <= kce; ++kk) {
    00663                     hrowi[ibase + kk] = hrowi[kk];
    00664                     hrowi[kk] = 0;
    00665                   }
    00666                   kce = xnewco + kce - kcs + 1;
    00667                   xnewco = kce + 1;
    00668                 } else {
    00669                   ++xnewco;
    00670                 }
    00671               } else if (hrowi[kce + 1] != 0) {
    00672                 /* here we use the fact that hrowi entries not "in use" are zeroed */
    00673                 if (xnewco + nz + 1 >= lstart) {
    00674                   /* Compress column file */
    00675                   if (nnentu + nz + 1 < lstart) {
    00676                     xnewco = c_ekkclco(fact,hrowi, mcstrt, hincol, xnewco);
    00677                     ++ncompactions;
    00678
    00679                     kcpiv = mcstrt[jpivot] - 1;
    00680                     kcs = mcstrt[j];
    00681                     /*                  HINCOL MAY HAVE CHANGED? (JJHF) ??? */
    00682                     nz = hincol[j];
    00683                     kce = kcs + nz - 1;
    00684                   } else {
    00685                     irtcod = -5;
    00686                     goto L1050;
    00687                   }
    00688                 }
    00689                 /* move the column to the end of the column file */
    00690                 mcstrt[j] = xnewco + 1;
    00691                 ibase = mcstrt[j] - kcs;
    00692                 for (kk = kcs; kk <= kce; ++kk) {
    00693                   hrowi[ibase + kk] = hrowi[kk];
    00694                   hrowi[kk] = 0;
    00695                 }
    00696                 kce = xnewco + kce - kcs + 1;
    00697                 xnewco = kce + 1;
    00698               }
    00699               /* store element */
    00700               hrowi[kce + 1] = npr;
    00701               hincol[j] = nz + 1;
    00702             }
    00703           } else {
    00704             maction[j] = static_cast<MACTION_T>(-maction[j]);
    00705           }
    00706         }
    00707         if (fill > kqq) {
    00708           xnewro = kstart;
    00709         }
    00710       }
    00711       hinrow[npr] = kstart - mrstrt[npr] + 1;
    00712       /* Check if row or column file needs compression */
    00713       if (! (xnewco + 1 < lstart)) {
    00714         xnewco = c_ekkclco(fact,hrowi, mcstrt, hincol, xnewco);
    00715         ++ncompactions;
    00716
    00717         kcpiv = mcstrt[jpivot] - 1;
    00718       }
    00719       if (! (xnewro + 1 < lstart)) {
    00720         int iput = c_ekkrwcs(fact,dluval, hcoli, mrstrt, hinrow, mwork, nfirst);
    00721         kmxeta += xnewro - iput ;
    00722         xnewro = iput - 1;
    00723         ++ncompactions;
    00724
    00725         kipis = mrstrt[ipivot] + 1;
    00726         kipie = kipis + epivr1 - 1;
    00727       }
    00728       /* Store elementary row transformation */
    00729       ++nnentl;
    00730       --nnentu;
    00731       --lstart;
    00732       dluval[lstart] = multip;
    00733
    00734       hrowi[lstart] = SHIFT_INDEX(npr);
    00735 #define INLINE_AFPV 3
    00736       /* We could do this while computing values but
    00737          it makes it much more complex.  At least we should get
    00738          reasonable cache behavior by doing it each row */
    00739 #if INLINE_AFPV
    00740       {
    00741         int j;
    00742         int nel, krs;
    00743         int koff;
    00744         int * index;
    00745         double * els;
    00746         nel = hinrow[npr];
    00747         krs = mrstrt[npr];
    00748         index=&hcoli[krs];
    00749         els=&dluval[krs];
    00750 #if INLINE_AFPV<3
    00751 #if INLINE_AFPV==1
    00752         double maxaij = 0.0;
    00753         koff = 0;
    00754         j=0;
    00755         while (j<nel) {
    00756           double d = fabs(els[j]);
    00757           if (maxaij < d) {
    00758             maxaij = d;
    00759             koff=j;
    00760           }
    00761           j++;
    00762         }
    00763 #else
    00764         assert (nel);
    00765         koff=0;
    00766         double maxaij=fabs(els[0]);
    00767         for (j=1;j<nel;j++) {
    00768           double d = fabs(els[j]);
    00769           if (maxaij < d) {
    00770             maxaij = d;
    00771             koff=j;
    00772           }
    00773         }
    00774 #endif
    00775 #else
    00776         double maxaij = 0.0;
    00777         koff = 0;
    00778         j=0;
    00779         if ((nel&1)!=0) {
    00780           maxaij=fabs(els[0]);
    00781           j=1;
    00782         }
    00783
    00784         while (j<nel) {
    00785           UNROLL_LOOP_BODY2({
    00786               double d = fabs(els[j]);
    00787               if (maxaij < d) {
    00788                 maxaij = d;
    00789                 koff=j;
    00790               }
    00791               j++;
    00792             });
    00793         }
    00794 #endif
    00795         SWAP(int, index[koff], index[0]);
    00796         SWAP(double, els[koff], els[0]);
    00797       }
    00798 #endif
    00799
    00800       {
    00801         int nzi = hinrow[npr];
    00802         if (nzi > 0) {
    00803           C_EKK_ADD_LINK(hpivro, nzi, rlink, npr);
    00804         }
    00805       }
    00806     }
    00807
    00808     /* after pivot move biggest to first in each row */
    00809 #if INLINE_AFPV==0
    00810     int nn = mcstrt[fact->xnetal] - lstart + 1;
    00811     c_ekkafpv(hrowi+lstart, hcoli, dluval, mrstrt, hinrow, nn);
    00812 #endif
    00813
    00814     /* Restore work array */
    00815     for (k = kipis; k <= kipie; ++k) {
    00816       maction[hcoli[k]] = 0;
    00817     }
    00818
    00819     if (*xrejctp > 0) {
    00820       for (k = kipis; k <= kipie; ++k) {
    00821         int j = hcoli[k];
    00822         int nzj = hincol[j];
    00823         if (! (nzj <= 0) &&
    00824             ! ((clink[j].pre > nrow && nzj != 1))) {
    00825           C_EKK_ADD_LINK(hpivco, nzj, clink, j);
    00826         }
    00827       }
    00828     } else {
    00829       for (k = kipis; k <= kipie; ++k) {
    00830         int j = hcoli[k];
    00831         int nzj = hincol[j];
    00832         if (! (nzj <= 0)) {
    00833           C_EKK_ADD_LINK(hpivco, nzj, clink, j);
    00834         }
    00835       }
    00836     }
    00837     fact->nuspike += hinrow[ipivot];
    00838
    00839     /* Go to dense coding if appropriate */
    00840 #ifndef C_EKKCMFY
    00841     if (ifdens != 0) {
    00842       int ndense = nrow - fact->npivots;
    00843       if (! (xnewro + ndense * ndense >= lstart)) {
    00844
    00845         /* set up sort order in MACTION */
    00846         c_ekkizero( nrow, reinterpret_cast<int *> (maction+1));
    00847         iput = 0;
    00848         for (i = 1; i <= nrow; ++i) {
    00849           if (clink[i].pre >= 0) {
    00850             ++iput;
    00851             maction[i] = static_cast<short int>(iput);
    00852           }
    00853         }
    00854         /* and get number spare needed */
    00855         nspare = 0;
    00856         for (i = 1; i <= nrow; ++i) {
    00857           if (rlink[i].pre >= 0) {
    00858             nspare = nspare + ndense - hinrow[i];
    00859           }
    00860         }
    00861         if (iput != nrow - fact->npivots) {
    00862           /* must be singular */
    00863           c_ekkizero( nrow, reinterpret_cast<int *> (maction+1));
    00864         } else {
    00865           /* pack down then back up */
    00866           int iput = c_ekkrwcs(fact,dluval, hcoli, mrstrt, hinrow, mwork, nfirst);
    00867           kmxeta += xnewro - iput ;
    00868           xnewro = iput - 1;
    00869           ++ncompactions;
    00870
    00871           --ncompactions;
    00872           if (xnewro + nspare + ndense * ndense >= lstart) {
    00873             c_ekkizero( nrow, reinterpret_cast<int *> (maction+1));
    00874           }
    00875           else {
    00876             xnewro += nspare;
    00877             c_ekkrwct(fact,dluval, hcoli, mrstrt, hinrow, mwork,
    00878                     rlink, maction, dvalpv,
    00879                     nlast,  xnewro);
    00880             kmxeta += xnewro ;
    00881             if (nnentu + nnentl > nrow * 5 &&
    00882                 (ndense*ndense)>(nnentu+nnentl)>>2 &&
    00883                 !if_sparse_update) {
    00884               fact->ndenuc = ndense;
    00885             }
    00886             irtcod = c_ekkcmfd(fact,
    00887                              (reinterpret_cast<int*>(dvalpv)+1),
    00888                              rlink, clink,
    00889                              (reinterpret_cast<int*>(maction+1))+1,
    00890                              nnetas,
    00891                              &nnentl, &nnentu,
    00892                              nsingp);
    00893             /* irtcod == 0 || irtcod == 10 */
    00894             /* 10 == found 0.0 pivot */
    00895             goto L1050;
    00896           }
    00897         }
    00898       } else {
    00899         /* say not enough room */
    00900         /*printf("no room %d\n",ndense);*/
    00901         if (1) {
    00902           /* return and increase size of etas if possible */
    00903           if (!noRoomForDense) {
    00904             int etasize =CoinMax(4*fact->nnentu+(nnetas-fact->nnentl)+1000,fact->eta_size);
    00905             noRoomForDense=ndense;
    00906             fact->eta_size=CoinMin(static_cast<int>(1.2*fact->eta_size),etasize);
    00907             if (fact->maxNNetas>0&&fact->eta_size>
    00908                 fact->maxNNetas) {
    00909               fact->eta_size=fact->maxNNetas;
    00910             }
    00911           }
    00912         }
    00913       }
    00914     }
    00915 #endif  /* C_EKKCMFY */
    00916   }
    00917
    00918  L1050:
    00919   {
    00920     int iput = c_ekkrwcs(fact,dluval, hcoli, mrstrt, hinrow, mwork, nfirst);
    00921     kmxeta += xnewro - iput;
    00922     xnewro = iput - 1;
    00923     ++ncompactions;
    00924   }
    00925
    00926   nnentu = xnewro;
    00927   /* save order of row copy for c_ekkshfv */
    00928   mwork[nrow+1].pre = nfirst;
    00929   mwork[nrow+1].suc = nlast;
    00930
    00931   fact->nnentl = nnentl;
    00932   fact->nnentu = nnentu;
    00933   fact->kmxeta = kmxeta;
    00934   *xnewrop = xnewro;
    00935   *ncompactionsp = ncompactions;
    00936
    00937   return (irtcod);
    00938 } /* c_ekkcmfc */
    00939 #endif
    00940