CoinUtils
trunk
|
00001 /* $Id$ */ 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, ifdens, 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 #endif 00196 int irtcod = 0; 00197 const int nrow = fact->nrow; 00198 00199 /* Parameter adjustments */ 00200 --maction; 00201 00202 /* Function Body */ 00203 lstart = nnetas - nnentl + 1; 00204 for (i = lstart; i <= nnetas; ++i) { 00205 hrowi[i] = SHIFT_INDEX(hcoli[i]); 00206 } 00207 ifdens = 0; 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