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

- File List
- File Members

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