/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 

Generated on Tue Jan 19 03:07:47 2010 for coin-Cgl by  doxygen 1.4.7