CoinOslC.h

Go to the documentation of this file.
00001 /* $Id: CoinOslC.h 1582 2013-04-06 14:33:07Z 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 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

Generated on 28 Aug 2016 for CoinUtils by  doxygen 1.6.1