lp_solve / bfp / bfp_LUSOL / LUSOL / lusol8a.c

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279``` ``` /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ File lusol8a lu8rpc Sparse LU update: Replace Column LUSOL's sparse implementation of the Bartels-Golub update. 01 May 2002: Derived from LUSOL's original lu8a.f file. 01 May 2002: Current version of lusol8a.f. ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ /* ================================================================== lu8rpc updates the LU factorization A = L*U when column jrep is replaced by some vector a(new). lu8rpc is an implementation of the Bartels-Golub update, designed for the case where A is rectangular and/or singular. L is a product of stabilized eliminations (m x m, nonsingular). P U Q is upper trapezoidal (m x n, rank nrank). If mode1 = 0, the old column is taken to be zero (so it does not have to be removed from U). If mode1 = 1, the old column need not have been zero. If mode2 = 0, the new column is taken to be zero. v(*) is not used or altered. If mode2 = 1, v(*) must contain the new column a(new). On exit, v(*) will satisfy L*v = a(new). If mode2 = 2, v(*) must satisfy L*v = a(new). The array w(*) is not used or altered. On entry, all elements of locc are assumed to be zero. On a successful exit (inform != 7), this will again be true. On exit: inform = -1 if the rank of U decreased by 1. inform = 0 if the rank of U stayed the same. inform = 1 if the rank of U increased by 1. inform = 2 if the update seemed to be unstable (diag much bigger than vnorm). inform = 7 if the update was not completed (lack of storage). inform = 8 if jrep is not between 1 and n. ------------------------------------------------------------------ -- Jan 1985: Original F66 version. -- Jul 1987: Modified to maintain U in trapezoidal form. 10 May 1988: First f77 version. 16 Oct 2000: Added test for instability (inform = 2). ================================================================== */ void LU8RPC(LUSOLrec *LUSOL, int MODE1, int MODE2, int JREP, REAL V[], REAL W[], int *INFORM, REAL *DIAG, REAL *VNORM) { MYBOOL SINGLR; int LPRINT, NRANK, LENL, LENU, LROW, NRANK0, KREP, KLAST, IW, L1, J1, JSING; REAL UTOL1, UTOL2; LPRINT = LUSOL->luparm[LUSOL_IP_PRINTLEVEL]; NRANK = LUSOL->luparm[LUSOL_IP_RANK_U]; LENL = LUSOL->luparm[LUSOL_IP_NONZEROS_L]; LENU = LUSOL->luparm[LUSOL_IP_NONZEROS_U]; LROW = LUSOL->luparm[LUSOL_IP_NONZEROS_ROW]; UTOL1 = LUSOL->parmlu[LUSOL_RP_SMALLDIAG_U]; UTOL2 = LUSOL->parmlu[LUSOL_RP_EPSDIAG_U]; NRANK0 = NRANK; *DIAG = ZERO; *VNORM = ZERO; if(JREP<1) goto x980; if(JREP>LUSOL->n) goto x980; /* ------------------------------------------------------------------ If mode1 = 0, there are no elements to be removed from U but we still have to set krep (using a backward loop). Otherwise, use lu7zap to remove column jrep from U and set krep at the same time. ------------------------------------------------------------------ */ if(MODE1==LUSOL_UPDATE_OLDEMPTY) { KREP = LUSOL->n+1; x10: KREP--; if(LUSOL->iq[KREP]!=JREP) goto x10; } else LU7ZAP(LUSOL, JREP,&KREP,&LENU,&LROW,NRANK); /* ------------------------------------------------------------------ Insert a new column of u and find klast. ------------------------------------------------------------------ */ if(MODE2==LUSOL_UPDATE_NEWEMPTY) { KLAST = 0; } else { if(MODE2==LUSOL_UPDATE_NEWNONEMPTY) { /* Transform v = a(new) to satisfy L*v = a(new). */ LU6SOL(LUSOL, LUSOL_SOLVE_Lv_v, V,W, NULL, INFORM); } else if(V==NULL) /* Otherwise, the V vector is taken to satisfy this already, or stored earlier. */ V=LUSOL->vLU6L; /* Insert into U any nonzeros in the top of v. row ip(klast) will contain the last nonzero in pivotal order. Note that klast will be in the range ( 0, nrank ). */ LU7ADD(LUSOL, JREP,V,LENL,&LENU,&LROW,NRANK,INFORM,&KLAST,VNORM); if(*INFORM==LUSOL_INFORM_ANEEDMEM) goto x970; } /* ------------------------------------------------------------------ In general, the new column causes U to look like this: krep n krep n ....a......... ..........a... . a . . a . . a . . a . .a . . a . P U Q = a . or . a . b. . . a . b . . . a . b . . . a . b ...... ..a... nrank c c c c c c m klast points to the last nonzero "a" or "b". klast = 0 means all "a" and "b" entries are zero. ------------------------------------------------------------------ */ if(MODE2==LUSOL_UPDATE_NEWEMPTY) { if(KREP>NRANK) goto x900; } else if(NRANKm) { /* Eliminate any "c"s (in either case). Row nrank + 1 may end up containing one nonzero. */ LU7ELM(LUSOL, JREP,V,&LENL,&LROW,NRANK,INFORM,DIAG); if(*INFORM==LUSOL_INFORM_ANEEDMEM) goto x970; if(*INFORM==LUSOL_INFORM_LUSINGULAR) { /* The nonzero is apparently significant. Increase nrank by 1 and make klast point to the bottom. */ NRANK++; KLAST = NRANK; } } if(NRANKn) { /* The column rank is low. In the first case, we want the new column to end up in position nrank, so the trapezoidal columns will have a chance later on (in lu7rnk) to pivot in that position. Otherwise the new column is not part of the triangle. We swap it into position nrank so we can judge it for singularity. lu7rnk might choose some other trapezoidal column later. */ if(KREPiq[KREP] = LUSOL->iq[NRANK]; LUSOL->iq[NRANK] = JREP; KREP = NRANK; } } /* ------------------------------------------------------------------ If krep .lt. klast, there are some "b"s to eliminate: krep ....a......... . a . . a . .a . P U Q = a . krep b. . b . . b . . b ...... nrank If krep .eq. klast, there are no "b"s, but the last "a" still has to be moved to the front of row krep (by lu7for). ------------------------------------------------------------------ */ if(KREP<=KLAST) { /* Perform a cyclic permutation on the current pivotal order, and eliminate the resulting row spike. krep becomes klast. The final diagonal (if any) will be correctly positioned at the front of the new krep-th row. nrank stays the same. */ LU7CYC(LUSOL, KREP,KLAST,LUSOL->ip); LU7CYC(LUSOL, KREP,KLAST,LUSOL->iq); LU7FOR(LUSOL, KREP,KLAST,&LENL,&LENU,&LROW,INFORM,DIAG); if(*INFORM==LUSOL_INFORM_ANEEDMEM) goto x970; KREP = KLAST; /* Test for instability (diag much bigger than vnorm). */ SINGLR = (MYBOOL) ((*VNORM)ip[KREP]; SINGLR = (MYBOOL) (LUSOL->lenr[IW]==0); if(!SINGLR) { L1 = LUSOL->locr[IW]; J1 = LUSOL->indr[L1]; SINGLR = (MYBOOL) (J1!=JREP); if(!SINGLR) { *DIAG = LUSOL->a[L1]; SINGLR = (MYBOOL) (fabs(*DIAG)<=UTOL1 || fabs(*DIAG)<=UTOL2*(*VNORM)); } } if(SINGLR && (KREPip); LU7CYC(LUSOL, KREP,LUSOL->n,LUSOL->iq); LU7FOR(LUSOL, KREP,NRANK,&LENL,&LENU,&LROW,INFORM,DIAG); if(*INFORM==LUSOL_INFORM_ANEEDMEM) goto x970; } /* Find the best column to be in position nrank. If singlr, it can't be the new column, jrep. If nothing satisfactory exists, nrank will be decreased. */ if(SINGLR || (NRANKn)) { JSING = 0; if(SINGLR) JSING = JREP; LU7RNK(LUSOL, JSING,&LENU,&LROW,&NRANK,INFORM,DIAG); } /* ------------------------------------------------------------------ Update indeces of optional row-based version of L0. ------------------------------------------------------------------ */ #if 0 if(LUSOL->L0 != NULL) LU1L0UPD(LUSOL, INFORM); #endif /* ------------------------------------------------------------------ Set inform for exit. ------------------------------------------------------------------ */ x900: if(NRANK==NRANK0) *INFORM = LUSOL_INFORM_LUSUCCESS; else if(NRANKn) { if(LPRINT>=LUSOL_MSG_SINGULARITY) LUSOL_report(LUSOL, 0, "lu8rpc warning...\nSingularity after replacing column. jrep=%8d diag=%g\n", JREP,DIAG); } } else *INFORM = LUSOL_INFORM_LUSINGULAR; goto x990; /* Instability. */ x920: *INFORM = LUSOL_INFORM_LUUNSTABLE; if(LPRINT>=LUSOL_MSG_SINGULARITY) LUSOL_report(LUSOL, 0, "lu8rpc warning...\nInstability after replacing column. jrep=%8d diag=%g\n", JREP,DIAG); goto x990; /* Not enough storage. */ x970: *INFORM = LUSOL_INFORM_ANEEDMEM; if(LPRINT>=LUSOL_MSG_SINGULARITY) LUSOL_report(LUSOL, 0, "lu8rpc error...\nInsufficient memory. lena=%8d\n", LUSOL->lena); goto x990; /* jrep is out of range. */ x980: *INFORM = LUSOL_INFORM_FATALERR; if(LPRINT>=LUSOL_MSG_SINGULARITY) LUSOL_report(LUSOL, 0, "lu8rpc error...\njrep is out of range. m=%8d n=%8d jrep=%8d\n", LUSOL->m,LUSOL->n,JREP); /* Exit. */ x990: LUSOL->luparm[LUSOL_IP_UPDATECOUNT]++; LUSOL->luparm[LUSOL_IP_RANK_U] = NRANK; LUSOL->luparm[LUSOL_IP_NONZEROS_L] = LENL; LUSOL->luparm[LUSOL_IP_NONZEROS_U] = LENU; LUSOL->luparm[LUSOL_IP_NONZEROS_ROW] = LROW; LUSOL->luparm[LUSOL_IP_INFORM] = *INFORM; } ```