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

- File List
- File Members

# /home/coin/svn-release/Alps-1.4.6/CoinUtils/src/CoinOslC.h

Go to the documentation of this file.

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