Changeset 4068 for trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.5.0/ALGLIB-2.5.0/inverseupdate.cs
- Timestamp:
- 07/22/10 00:44:01 (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.5.0/ALGLIB-2.5.0/inverseupdate.cs
r3839 r4068 19 19 *************************************************************************/ 20 20 21 using System; 22 23 namespace alglib 24 { 25 public class inverseupdate 26 { 27 /************************************************************************* 28 Inverse matrix update by the Sherman-Morrison formula 29 30 The algorithm updates matrix A^-1 when adding a number to an element 31 of matrix A. 32 33 Input parameters: 34 InvA - inverse of matrix A. 21 22 namespace alglib { 23 public class inverseupdate { 24 /************************************************************************* 25 Inverse matrix update by the Sherman-Morrison formula 26 27 The algorithm updates matrix A^-1 when adding a number to an element 28 of matrix A. 29 30 Input parameters: 31 InvA - inverse of matrix A. 32 Array whose indexes range within [0..N-1, 0..N-1]. 33 N - size of matrix A. 34 UpdRow - row where the element to be updated is stored. 35 UpdColumn - column where the element to be updated is stored. 36 UpdVal - a number to be added to the element. 37 38 39 Output parameters: 40 InvA - inverse of modified matrix A. 41 42 -- ALGLIB -- 43 Copyright 2005 by Bochkanov Sergey 44 *************************************************************************/ 45 public static void rmatrixinvupdatesimple(ref double[,] inva, 46 int n, 47 int updrow, 48 int updcolumn, 49 double updval) { 50 double[] t1 = new double[0]; 51 double[] t2 = new double[0]; 52 int i = 0; 53 double lambda = 0; 54 double vt = 0; 55 int i_ = 0; 56 57 System.Diagnostics.Debug.Assert(updrow >= 0 & updrow < n, "RMatrixInvUpdateSimple: incorrect UpdRow!"); 58 System.Diagnostics.Debug.Assert(updcolumn >= 0 & updcolumn < n, "RMatrixInvUpdateSimple: incorrect UpdColumn!"); 59 t1 = new double[n - 1 + 1]; 60 t2 = new double[n - 1 + 1]; 61 62 // 63 // T1 = InvA * U 64 // 65 for (i_ = 0; i_ <= n - 1; i_++) { 66 t1[i_] = inva[i_, updrow]; 67 } 68 69 // 70 // T2 = v*InvA 71 // 72 for (i_ = 0; i_ <= n - 1; i_++) { 73 t2[i_] = inva[updcolumn, i_]; 74 } 75 76 // 77 // Lambda = v * InvA * U 78 // 79 lambda = updval * inva[updcolumn, updrow]; 80 81 // 82 // InvA = InvA - correction 83 // 84 for (i = 0; i <= n - 1; i++) { 85 vt = updval * t1[i]; 86 vt = vt / (1 + lambda); 87 for (i_ = 0; i_ <= n - 1; i_++) { 88 inva[i, i_] = inva[i, i_] - vt * t2[i_]; 89 } 90 } 91 } 92 93 94 /************************************************************************* 95 Inverse matrix update by the Sherman-Morrison formula 96 97 The algorithm updates matrix A^-1 when adding a vector to a row 98 of matrix A. 99 100 Input parameters: 101 InvA - inverse of matrix A. 102 Array whose indexes range within [0..N-1, 0..N-1]. 103 N - size of matrix A. 104 UpdRow - the row of A whose vector V was added. 105 0 <= Row <= N-1 106 V - the vector to be added to a row. 107 Array whose index ranges within [0..N-1]. 108 109 Output parameters: 110 InvA - inverse of modified matrix A. 111 112 -- ALGLIB -- 113 Copyright 2005 by Bochkanov Sergey 114 *************************************************************************/ 115 public static void rmatrixinvupdaterow(ref double[,] inva, 116 int n, 117 int updrow, 118 ref double[] v) { 119 double[] t1 = new double[0]; 120 double[] t2 = new double[0]; 121 int i = 0; 122 int j = 0; 123 double lambda = 0; 124 double vt = 0; 125 int i_ = 0; 126 127 t1 = new double[n - 1 + 1]; 128 t2 = new double[n - 1 + 1]; 129 130 // 131 // T1 = InvA * U 132 // 133 for (i_ = 0; i_ <= n - 1; i_++) { 134 t1[i_] = inva[i_, updrow]; 135 } 136 137 // 138 // T2 = v*InvA 139 // Lambda = v * InvA * U 140 // 141 for (j = 0; j <= n - 1; j++) { 142 vt = 0.0; 143 for (i_ = 0; i_ <= n - 1; i_++) { 144 vt += v[i_] * inva[i_, j]; 145 } 146 t2[j] = vt; 147 } 148 lambda = t2[updrow]; 149 150 // 151 // InvA = InvA - correction 152 // 153 for (i = 0; i <= n - 1; i++) { 154 vt = t1[i] / (1 + lambda); 155 for (i_ = 0; i_ <= n - 1; i_++) { 156 inva[i, i_] = inva[i, i_] - vt * t2[i_]; 157 } 158 } 159 } 160 161 162 /************************************************************************* 163 Inverse matrix update by the Sherman-Morrison formula 164 165 The algorithm updates matrix A^-1 when adding a vector to a column 166 of matrix A. 167 168 Input parameters: 169 InvA - inverse of matrix A. 35 170 Array whose indexes range within [0..N-1, 0..N-1]. 36 N - size of matrix A. 37 UpdRow - row where the element to be updated is stored. 38 UpdColumn - column where the element to be updated is stored. 39 UpdVal - a number to be added to the element. 40 41 42 Output parameters: 43 InvA - inverse of modified matrix A. 44 45 -- ALGLIB -- 46 Copyright 2005 by Bochkanov Sergey 47 *************************************************************************/ 48 public static void rmatrixinvupdatesimple(ref double[,] inva, 49 int n, 50 int updrow, 51 int updcolumn, 52 double updval) 53 { 54 double[] t1 = new double[0]; 55 double[] t2 = new double[0]; 56 int i = 0; 57 double lambda = 0; 58 double vt = 0; 59 int i_ = 0; 60 61 System.Diagnostics.Debug.Assert(updrow>=0 & updrow<n, "RMatrixInvUpdateSimple: incorrect UpdRow!"); 62 System.Diagnostics.Debug.Assert(updcolumn>=0 & updcolumn<n, "RMatrixInvUpdateSimple: incorrect UpdColumn!"); 63 t1 = new double[n-1+1]; 64 t2 = new double[n-1+1]; 65 66 // 67 // T1 = InvA * U 68 // 69 for(i_=0; i_<=n-1;i_++) 70 { 71 t1[i_] = inva[i_,updrow]; 72 } 73 74 // 75 // T2 = v*InvA 76 // 77 for(i_=0; i_<=n-1;i_++) 78 { 79 t2[i_] = inva[updcolumn,i_]; 80 } 81 82 // 83 // Lambda = v * InvA * U 84 // 85 lambda = updval*inva[updcolumn,updrow]; 86 87 // 88 // InvA = InvA - correction 89 // 90 for(i=0; i<=n-1; i++) 91 { 92 vt = updval*t1[i]; 93 vt = vt/(1+lambda); 94 for(i_=0; i_<=n-1;i_++) 95 { 96 inva[i,i_] = inva[i,i_] - vt*t2[i_]; 97 } 98 } 99 } 100 101 102 /************************************************************************* 103 Inverse matrix update by the Sherman-Morrison formula 104 105 The algorithm updates matrix A^-1 when adding a vector to a row 106 of matrix A. 107 108 Input parameters: 109 InvA - inverse of matrix A. 110 Array whose indexes range within [0..N-1, 0..N-1]. 111 N - size of matrix A. 112 UpdRow - the row of A whose vector V was added. 113 0 <= Row <= N-1 114 V - the vector to be added to a row. 171 N - size of matrix A. 172 UpdColumn - the column of A whose vector U was added. 173 0 <= UpdColumn <= N-1 174 U - the vector to be added to a column. 115 175 Array whose index ranges within [0..N-1]. 116 176 117 Output parameters: 118 InvA - inverse of modified matrix A. 119 120 -- ALGLIB -- 121 Copyright 2005 by Bochkanov Sergey 122 *************************************************************************/ 123 public static void rmatrixinvupdaterow(ref double[,] inva, 124 int n, 125 int updrow, 126 ref double[] v) 127 { 128 double[] t1 = new double[0]; 129 double[] t2 = new double[0]; 130 int i = 0; 131 int j = 0; 132 double lambda = 0; 133 double vt = 0; 134 int i_ = 0; 135 136 t1 = new double[n-1+1]; 137 t2 = new double[n-1+1]; 138 139 // 140 // T1 = InvA * U 141 // 142 for(i_=0; i_<=n-1;i_++) 143 { 144 t1[i_] = inva[i_,updrow]; 145 } 146 147 // 148 // T2 = v*InvA 149 // Lambda = v * InvA * U 150 // 151 for(j=0; j<=n-1; j++) 152 { 153 vt = 0.0; 154 for(i_=0; i_<=n-1;i_++) 155 { 156 vt += v[i_]*inva[i_,j]; 157 } 158 t2[j] = vt; 159 } 160 lambda = t2[updrow]; 161 162 // 163 // InvA = InvA - correction 164 // 165 for(i=0; i<=n-1; i++) 166 { 167 vt = t1[i]/(1+lambda); 168 for(i_=0; i_<=n-1;i_++) 169 { 170 inva[i,i_] = inva[i,i_] - vt*t2[i_]; 171 } 172 } 173 } 174 175 176 /************************************************************************* 177 Inverse matrix update by the Sherman-Morrison formula 178 179 The algorithm updates matrix A^-1 when adding a vector to a column 180 of matrix A. 181 182 Input parameters: 183 InvA - inverse of matrix A. 184 Array whose indexes range within [0..N-1, 0..N-1]. 185 N - size of matrix A. 186 UpdColumn - the column of A whose vector U was added. 187 0 <= UpdColumn <= N-1 188 U - the vector to be added to a column. 189 Array whose index ranges within [0..N-1]. 190 191 Output parameters: 192 InvA - inverse of modified matrix A. 193 194 -- ALGLIB -- 195 Copyright 2005 by Bochkanov Sergey 196 *************************************************************************/ 197 public static void rmatrixinvupdatecolumn(ref double[,] inva, 198 int n, 199 int updcolumn, 200 ref double[] u) 201 { 202 double[] t1 = new double[0]; 203 double[] t2 = new double[0]; 204 int i = 0; 205 double lambda = 0; 206 double vt = 0; 207 int i_ = 0; 208 209 t1 = new double[n-1+1]; 210 t2 = new double[n-1+1]; 211 212 // 213 // T1 = InvA * U 214 // Lambda = v * InvA * U 215 // 216 for(i=0; i<=n-1; i++) 217 { 218 vt = 0.0; 219 for(i_=0; i_<=n-1;i_++) 220 { 221 vt += inva[i,i_]*u[i_]; 222 } 223 t1[i] = vt; 224 } 225 lambda = t1[updcolumn]; 226 227 // 228 // T2 = v*InvA 229 // 230 for(i_=0; i_<=n-1;i_++) 231 { 232 t2[i_] = inva[updcolumn,i_]; 233 } 234 235 // 236 // InvA = InvA - correction 237 // 238 for(i=0; i<=n-1; i++) 239 { 240 vt = t1[i]/(1+lambda); 241 for(i_=0; i_<=n-1;i_++) 242 { 243 inva[i,i_] = inva[i,i_] - vt*t2[i_]; 244 } 245 } 246 } 247 248 249 /************************************************************************* 250 Inverse matrix update by the Sherman-Morrison formula 251 252 The algorithm computes the inverse of matrix A+u*v by using the given matrix 253 A^-1 and the vectors u and v. 254 255 Input parameters: 256 InvA - inverse of matrix A. 257 Array whose indexes range within [0..N-1, 0..N-1]. 258 N - size of matrix A. 259 U - the vector modifying the matrix. 260 Array whose index ranges within [0..N-1]. 261 V - the vector modifying the matrix. 262 Array whose index ranges within [0..N-1]. 263 264 Output parameters: 265 InvA - inverse of matrix A + u*v'. 266 267 -- ALGLIB -- 268 Copyright 2005 by Bochkanov Sergey 269 *************************************************************************/ 270 public static void rmatrixinvupdateuv(ref double[,] inva, 271 int n, 272 ref double[] u, 273 ref double[] v) 274 { 275 double[] t1 = new double[0]; 276 double[] t2 = new double[0]; 277 int i = 0; 278 int j = 0; 279 double lambda = 0; 280 double vt = 0; 281 int i_ = 0; 282 283 t1 = new double[n-1+1]; 284 t2 = new double[n-1+1]; 285 286 // 287 // T1 = InvA * U 288 // Lambda = v * T1 289 // 290 for(i=0; i<=n-1; i++) 291 { 292 vt = 0.0; 293 for(i_=0; i_<=n-1;i_++) 294 { 295 vt += inva[i,i_]*u[i_]; 296 } 297 t1[i] = vt; 298 } 299 lambda = 0.0; 300 for(i_=0; i_<=n-1;i_++) 301 { 302 lambda += v[i_]*t1[i_]; 303 } 304 305 // 306 // T2 = v*InvA 307 // 308 for(j=0; j<=n-1; j++) 309 { 310 vt = 0.0; 311 for(i_=0; i_<=n-1;i_++) 312 { 313 vt += v[i_]*inva[i_,j]; 314 } 315 t2[j] = vt; 316 } 317 318 // 319 // InvA = InvA - correction 320 // 321 for(i=0; i<=n-1; i++) 322 { 323 vt = t1[i]/(1+lambda); 324 for(i_=0; i_<=n-1;i_++) 325 { 326 inva[i,i_] = inva[i,i_] - vt*t2[i_]; 327 } 328 } 329 } 330 331 332 public static void shermanmorrisonsimpleupdate(ref double[,] inva, 333 int n, 334 int updrow, 335 int updcolumn, 336 double updval) 337 { 338 double[] t1 = new double[0]; 339 double[] t2 = new double[0]; 340 int i = 0; 341 double lambda = 0; 342 double vt = 0; 343 int i_ = 0; 344 345 t1 = new double[n+1]; 346 t2 = new double[n+1]; 347 348 // 349 // T1 = InvA * U 350 // 351 for(i_=1; i_<=n;i_++) 352 { 353 t1[i_] = inva[i_,updrow]; 354 } 355 356 // 357 // T2 = v*InvA 358 // 359 for(i_=1; i_<=n;i_++) 360 { 361 t2[i_] = inva[updcolumn,i_]; 362 } 363 364 // 365 // Lambda = v * InvA * U 366 // 367 lambda = updval*inva[updcolumn,updrow]; 368 369 // 370 // InvA = InvA - correction 371 // 372 for(i=1; i<=n; i++) 373 { 374 vt = updval*t1[i]; 375 vt = vt/(1+lambda); 376 for(i_=1; i_<=n;i_++) 377 { 378 inva[i,i_] = inva[i,i_] - vt*t2[i_]; 379 } 380 } 381 } 382 383 384 public static void shermanmorrisonupdaterow(ref double[,] inva, 385 int n, 386 int updrow, 387 ref double[] v) 388 { 389 double[] t1 = new double[0]; 390 double[] t2 = new double[0]; 391 int i = 0; 392 int j = 0; 393 double lambda = 0; 394 double vt = 0; 395 int i_ = 0; 396 397 t1 = new double[n+1]; 398 t2 = new double[n+1]; 399 400 // 401 // T1 = InvA * U 402 // 403 for(i_=1; i_<=n;i_++) 404 { 405 t1[i_] = inva[i_,updrow]; 406 } 407 408 // 409 // T2 = v*InvA 410 // Lambda = v * InvA * U 411 // 412 for(j=1; j<=n; j++) 413 { 414 vt = 0.0; 415 for(i_=1; i_<=n;i_++) 416 { 417 vt += v[i_]*inva[i_,j]; 418 } 419 t2[j] = vt; 420 } 421 lambda = t2[updrow]; 422 423 // 424 // InvA = InvA - correction 425 // 426 for(i=1; i<=n; i++) 427 { 428 vt = t1[i]/(1+lambda); 429 for(i_=1; i_<=n;i_++) 430 { 431 inva[i,i_] = inva[i,i_] - vt*t2[i_]; 432 } 433 } 434 } 435 436 437 public static void shermanmorrisonupdatecolumn(ref double[,] inva, 438 int n, 439 int updcolumn, 440 ref double[] u) 441 { 442 double[] t1 = new double[0]; 443 double[] t2 = new double[0]; 444 int i = 0; 445 double lambda = 0; 446 double vt = 0; 447 int i_ = 0; 448 449 t1 = new double[n+1]; 450 t2 = new double[n+1]; 451 452 // 453 // T1 = InvA * U 454 // Lambda = v * InvA * U 455 // 456 for(i=1; i<=n; i++) 457 { 458 vt = 0.0; 459 for(i_=1; i_<=n;i_++) 460 { 461 vt += inva[i,i_]*u[i_]; 462 } 463 t1[i] = vt; 464 } 465 lambda = t1[updcolumn]; 466 467 // 468 // T2 = v*InvA 469 // 470 for(i_=1; i_<=n;i_++) 471 { 472 t2[i_] = inva[updcolumn,i_]; 473 } 474 475 // 476 // InvA = InvA - correction 477 // 478 for(i=1; i<=n; i++) 479 { 480 vt = t1[i]/(1+lambda); 481 for(i_=1; i_<=n;i_++) 482 { 483 inva[i,i_] = inva[i,i_] - vt*t2[i_]; 484 } 485 } 486 } 487 488 489 public static void shermanmorrisonupdateuv(ref double[,] inva, 490 int n, 491 ref double[] u, 492 ref double[] v) 493 { 494 double[] t1 = new double[0]; 495 double[] t2 = new double[0]; 496 int i = 0; 497 int j = 0; 498 double lambda = 0; 499 double vt = 0; 500 int i_ = 0; 501 502 t1 = new double[n+1]; 503 t2 = new double[n+1]; 504 505 // 506 // T1 = InvA * U 507 // Lambda = v * T1 508 // 509 for(i=1; i<=n; i++) 510 { 511 vt = 0.0; 512 for(i_=1; i_<=n;i_++) 513 { 514 vt += inva[i,i_]*u[i_]; 515 } 516 t1[i] = vt; 517 } 518 lambda = 0.0; 519 for(i_=1; i_<=n;i_++) 520 { 521 lambda += v[i_]*t1[i_]; 522 } 523 524 // 525 // T2 = v*InvA 526 // 527 for(j=1; j<=n; j++) 528 { 529 vt = 0.0; 530 for(i_=1; i_<=n;i_++) 531 { 532 vt += v[i_]*inva[i_,j]; 533 } 534 t2[j] = vt; 535 } 536 537 // 538 // InvA = InvA - correction 539 // 540 for(i=1; i<=n; i++) 541 { 542 vt = t1[i]/(1+lambda); 543 for(i_=1; i_<=n;i_++) 544 { 545 inva[i,i_] = inva[i,i_] - vt*t2[i_]; 546 } 547 } 548 } 549 } 177 Output parameters: 178 InvA - inverse of modified matrix A. 179 180 -- ALGLIB -- 181 Copyright 2005 by Bochkanov Sergey 182 *************************************************************************/ 183 public static void rmatrixinvupdatecolumn(ref double[,] inva, 184 int n, 185 int updcolumn, 186 ref double[] u) { 187 double[] t1 = new double[0]; 188 double[] t2 = new double[0]; 189 int i = 0; 190 double lambda = 0; 191 double vt = 0; 192 int i_ = 0; 193 194 t1 = new double[n - 1 + 1]; 195 t2 = new double[n - 1 + 1]; 196 197 // 198 // T1 = InvA * U 199 // Lambda = v * InvA * U 200 // 201 for (i = 0; i <= n - 1; i++) { 202 vt = 0.0; 203 for (i_ = 0; i_ <= n - 1; i_++) { 204 vt += inva[i, i_] * u[i_]; 205 } 206 t1[i] = vt; 207 } 208 lambda = t1[updcolumn]; 209 210 // 211 // T2 = v*InvA 212 // 213 for (i_ = 0; i_ <= n - 1; i_++) { 214 t2[i_] = inva[updcolumn, i_]; 215 } 216 217 // 218 // InvA = InvA - correction 219 // 220 for (i = 0; i <= n - 1; i++) { 221 vt = t1[i] / (1 + lambda); 222 for (i_ = 0; i_ <= n - 1; i_++) { 223 inva[i, i_] = inva[i, i_] - vt * t2[i_]; 224 } 225 } 226 } 227 228 229 /************************************************************************* 230 Inverse matrix update by the Sherman-Morrison formula 231 232 The algorithm computes the inverse of matrix A+u*v by using the given matrix 233 A^-1 and the vectors u and v. 234 235 Input parameters: 236 InvA - inverse of matrix A. 237 Array whose indexes range within [0..N-1, 0..N-1]. 238 N - size of matrix A. 239 U - the vector modifying the matrix. 240 Array whose index ranges within [0..N-1]. 241 V - the vector modifying the matrix. 242 Array whose index ranges within [0..N-1]. 243 244 Output parameters: 245 InvA - inverse of matrix A + u*v'. 246 247 -- ALGLIB -- 248 Copyright 2005 by Bochkanov Sergey 249 *************************************************************************/ 250 public static void rmatrixinvupdateuv(ref double[,] inva, 251 int n, 252 ref double[] u, 253 ref double[] v) { 254 double[] t1 = new double[0]; 255 double[] t2 = new double[0]; 256 int i = 0; 257 int j = 0; 258 double lambda = 0; 259 double vt = 0; 260 int i_ = 0; 261 262 t1 = new double[n - 1 + 1]; 263 t2 = new double[n - 1 + 1]; 264 265 // 266 // T1 = InvA * U 267 // Lambda = v * T1 268 // 269 for (i = 0; i <= n - 1; i++) { 270 vt = 0.0; 271 for (i_ = 0; i_ <= n - 1; i_++) { 272 vt += inva[i, i_] * u[i_]; 273 } 274 t1[i] = vt; 275 } 276 lambda = 0.0; 277 for (i_ = 0; i_ <= n - 1; i_++) { 278 lambda += v[i_] * t1[i_]; 279 } 280 281 // 282 // T2 = v*InvA 283 // 284 for (j = 0; j <= n - 1; j++) { 285 vt = 0.0; 286 for (i_ = 0; i_ <= n - 1; i_++) { 287 vt += v[i_] * inva[i_, j]; 288 } 289 t2[j] = vt; 290 } 291 292 // 293 // InvA = InvA - correction 294 // 295 for (i = 0; i <= n - 1; i++) { 296 vt = t1[i] / (1 + lambda); 297 for (i_ = 0; i_ <= n - 1; i_++) { 298 inva[i, i_] = inva[i, i_] - vt * t2[i_]; 299 } 300 } 301 } 302 303 304 public static void shermanmorrisonsimpleupdate(ref double[,] inva, 305 int n, 306 int updrow, 307 int updcolumn, 308 double updval) { 309 double[] t1 = new double[0]; 310 double[] t2 = new double[0]; 311 int i = 0; 312 double lambda = 0; 313 double vt = 0; 314 int i_ = 0; 315 316 t1 = new double[n + 1]; 317 t2 = new double[n + 1]; 318 319 // 320 // T1 = InvA * U 321 // 322 for (i_ = 1; i_ <= n; i_++) { 323 t1[i_] = inva[i_, updrow]; 324 } 325 326 // 327 // T2 = v*InvA 328 // 329 for (i_ = 1; i_ <= n; i_++) { 330 t2[i_] = inva[updcolumn, i_]; 331 } 332 333 // 334 // Lambda = v * InvA * U 335 // 336 lambda = updval * inva[updcolumn, updrow]; 337 338 // 339 // InvA = InvA - correction 340 // 341 for (i = 1; i <= n; i++) { 342 vt = updval * t1[i]; 343 vt = vt / (1 + lambda); 344 for (i_ = 1; i_ <= n; i_++) { 345 inva[i, i_] = inva[i, i_] - vt * t2[i_]; 346 } 347 } 348 } 349 350 351 public static void shermanmorrisonupdaterow(ref double[,] inva, 352 int n, 353 int updrow, 354 ref double[] v) { 355 double[] t1 = new double[0]; 356 double[] t2 = new double[0]; 357 int i = 0; 358 int j = 0; 359 double lambda = 0; 360 double vt = 0; 361 int i_ = 0; 362 363 t1 = new double[n + 1]; 364 t2 = new double[n + 1]; 365 366 // 367 // T1 = InvA * U 368 // 369 for (i_ = 1; i_ <= n; i_++) { 370 t1[i_] = inva[i_, updrow]; 371 } 372 373 // 374 // T2 = v*InvA 375 // Lambda = v * InvA * U 376 // 377 for (j = 1; j <= n; j++) { 378 vt = 0.0; 379 for (i_ = 1; i_ <= n; i_++) { 380 vt += v[i_] * inva[i_, j]; 381 } 382 t2[j] = vt; 383 } 384 lambda = t2[updrow]; 385 386 // 387 // InvA = InvA - correction 388 // 389 for (i = 1; i <= n; i++) { 390 vt = t1[i] / (1 + lambda); 391 for (i_ = 1; i_ <= n; i_++) { 392 inva[i, i_] = inva[i, i_] - vt * t2[i_]; 393 } 394 } 395 } 396 397 398 public static void shermanmorrisonupdatecolumn(ref double[,] inva, 399 int n, 400 int updcolumn, 401 ref double[] u) { 402 double[] t1 = new double[0]; 403 double[] t2 = new double[0]; 404 int i = 0; 405 double lambda = 0; 406 double vt = 0; 407 int i_ = 0; 408 409 t1 = new double[n + 1]; 410 t2 = new double[n + 1]; 411 412 // 413 // T1 = InvA * U 414 // Lambda = v * InvA * U 415 // 416 for (i = 1; i <= n; i++) { 417 vt = 0.0; 418 for (i_ = 1; i_ <= n; i_++) { 419 vt += inva[i, i_] * u[i_]; 420 } 421 t1[i] = vt; 422 } 423 lambda = t1[updcolumn]; 424 425 // 426 // T2 = v*InvA 427 // 428 for (i_ = 1; i_ <= n; i_++) { 429 t2[i_] = inva[updcolumn, i_]; 430 } 431 432 // 433 // InvA = InvA - correction 434 // 435 for (i = 1; i <= n; i++) { 436 vt = t1[i] / (1 + lambda); 437 for (i_ = 1; i_ <= n; i_++) { 438 inva[i, i_] = inva[i, i_] - vt * t2[i_]; 439 } 440 } 441 } 442 443 444 public static void shermanmorrisonupdateuv(ref double[,] inva, 445 int n, 446 ref double[] u, 447 ref double[] v) { 448 double[] t1 = new double[0]; 449 double[] t2 = new double[0]; 450 int i = 0; 451 int j = 0; 452 double lambda = 0; 453 double vt = 0; 454 int i_ = 0; 455 456 t1 = new double[n + 1]; 457 t2 = new double[n + 1]; 458 459 // 460 // T1 = InvA * U 461 // Lambda = v * T1 462 // 463 for (i = 1; i <= n; i++) { 464 vt = 0.0; 465 for (i_ = 1; i_ <= n; i_++) { 466 vt += inva[i, i_] * u[i_]; 467 } 468 t1[i] = vt; 469 } 470 lambda = 0.0; 471 for (i_ = 1; i_ <= n; i_++) { 472 lambda += v[i_] * t1[i_]; 473 } 474 475 // 476 // T2 = v*InvA 477 // 478 for (j = 1; j <= n; j++) { 479 vt = 0.0; 480 for (i_ = 1; i_ <= n; i_++) { 481 vt += v[i_] * inva[i_, j]; 482 } 483 t2[j] = vt; 484 } 485 486 // 487 // InvA = InvA - correction 488 // 489 for (i = 1; i <= n; i++) { 490 vt = t1[i] / (1 + lambda); 491 for (i_ = 1; i_ <= n; i_++) { 492 inva[i, i_] = inva[i, i_] - vt * t2[i_]; 493 } 494 } 495 } 496 } 550 497 }
Note: See TracChangeset
for help on using the changeset viewer.
