• Main Page
  • Related Pages
  • Modules
  • Namespaces
  • Classes
  • Files
  • File List
  • File Members

/build/buildd/coinutils-2.6.4/CoinUtils/src/CoinOslC.h

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

Generated on Fri Oct 15 2010 18:21:02 by  doxygen 1.7.1