[2563] | 1 | /*************************************************************************
|
---|
| 2 | Copyright (c) 1992-2007 The University of Tennessee. All rights reserved.
|
---|
| 3 |
|
---|
| 4 | Contributors:
|
---|
| 5 | * Sergey Bochkanov (ALGLIB project). Translation from FORTRAN to
|
---|
| 6 | pseudocode.
|
---|
| 7 |
|
---|
| 8 | See subroutines comments for additional copyrights.
|
---|
| 9 |
|
---|
| 10 | >>> SOURCE LICENSE >>>
|
---|
| 11 | This program is free software; you can redistribute it and/or modify
|
---|
| 12 | it under the terms of the GNU General Public License as published by
|
---|
| 13 | the Free Software Foundation (www.fsf.org); either version 2 of the
|
---|
| 14 | License, or (at your option) any later version.
|
---|
| 15 |
|
---|
| 16 | This program is distributed in the hope that it will be useful,
|
---|
| 17 | but WITHOUT ANY WARRANTY; without even the implied warranty of
|
---|
| 18 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
---|
| 19 | GNU General Public License for more details.
|
---|
| 20 |
|
---|
| 21 | A copy of the GNU General Public License is available at
|
---|
| 22 | http://www.fsf.org/licensing/licenses
|
---|
| 23 |
|
---|
| 24 | >>> END OF LICENSE >>>
|
---|
| 25 | *************************************************************************/
|
---|
| 26 |
|
---|
| 27 | using System;
|
---|
| 28 |
|
---|
| 29 | namespace alglib
|
---|
| 30 | {
|
---|
| 31 | public class rcond
|
---|
| 32 | {
|
---|
| 33 | /*************************************************************************
|
---|
| 34 | Estimate of a matrix condition number (1-norm)
|
---|
| 35 |
|
---|
| 36 | The algorithm calculates a lower bound of the condition number. In this case,
|
---|
| 37 | the algorithm does not return a lower bound of the condition number, but an
|
---|
| 38 | inverse number (to avoid an overflow in case of a singular matrix).
|
---|
| 39 |
|
---|
| 40 | Input parameters:
|
---|
| 41 | A - matrix. Array whose indexes range within [0..N-1, 0..N-1].
|
---|
| 42 | N - size of matrix A.
|
---|
| 43 |
|
---|
| 44 | Result: 1/LowerBound(cond(A))
|
---|
| 45 | *************************************************************************/
|
---|
| 46 | public static double rmatrixrcond1(ref double[,] a,
|
---|
| 47 | int n)
|
---|
| 48 | {
|
---|
| 49 | double result = 0;
|
---|
| 50 | int i = 0;
|
---|
| 51 | double[,] a1 = new double[0,0];
|
---|
| 52 | int i_ = 0;
|
---|
| 53 | int i1_ = 0;
|
---|
| 54 |
|
---|
| 55 | System.Diagnostics.Debug.Assert(n>=1, "RMatrixRCond1: N<1!");
|
---|
| 56 | a1 = new double[n+1, n+1];
|
---|
| 57 | for(i=1; i<=n; i++)
|
---|
| 58 | {
|
---|
| 59 | i1_ = (0) - (1);
|
---|
| 60 | for(i_=1; i_<=n;i_++)
|
---|
| 61 | {
|
---|
| 62 | a1[i,i_] = a[i-1,i_+i1_];
|
---|
| 63 | }
|
---|
| 64 | }
|
---|
| 65 | result = rcond1(a1, n);
|
---|
| 66 | return result;
|
---|
| 67 | }
|
---|
| 68 |
|
---|
| 69 |
|
---|
| 70 | /*************************************************************************
|
---|
| 71 | Estimate of the condition number of a matrix given by its LU decomposition (1-norm)
|
---|
| 72 |
|
---|
| 73 | The algorithm calculates a lower bound of the condition number. In this case,
|
---|
| 74 | the algorithm does not return a lower bound of the condition number, but an
|
---|
| 75 | inverse number (to avoid an overflow in case of a singular matrix).
|
---|
| 76 |
|
---|
| 77 | Input parameters:
|
---|
| 78 | LUDcmp - LU decomposition of a matrix in compact form. Output of
|
---|
| 79 | the RMatrixLU subroutine.
|
---|
| 80 | N - size of matrix A.
|
---|
| 81 |
|
---|
| 82 | Result: 1/LowerBound(cond(A))
|
---|
| 83 | *************************************************************************/
|
---|
| 84 | public static double rmatrixlurcond1(ref double[,] ludcmp,
|
---|
| 85 | int n)
|
---|
| 86 | {
|
---|
| 87 | double result = 0;
|
---|
| 88 | int i = 0;
|
---|
| 89 | double[,] a1 = new double[0,0];
|
---|
| 90 | int i_ = 0;
|
---|
| 91 | int i1_ = 0;
|
---|
| 92 |
|
---|
| 93 | System.Diagnostics.Debug.Assert(n>=1, "RMatrixLURCond1: N<1!");
|
---|
| 94 | a1 = new double[n+1, n+1];
|
---|
| 95 | for(i=1; i<=n; i++)
|
---|
| 96 | {
|
---|
| 97 | i1_ = (0) - (1);
|
---|
| 98 | for(i_=1; i_<=n;i_++)
|
---|
| 99 | {
|
---|
| 100 | a1[i,i_] = ludcmp[i-1,i_+i1_];
|
---|
| 101 | }
|
---|
| 102 | }
|
---|
| 103 | result = rcond1lu(ref a1, n);
|
---|
| 104 | return result;
|
---|
| 105 | }
|
---|
| 106 |
|
---|
| 107 |
|
---|
| 108 | /*************************************************************************
|
---|
| 109 | Estimate of a matrix condition number (infinity-norm).
|
---|
| 110 |
|
---|
| 111 | The algorithm calculates a lower bound of the condition number. In this case,
|
---|
| 112 | the algorithm does not return a lower bound of the condition number, but an
|
---|
| 113 | inverse number (to avoid an overflow in case of a singular matrix).
|
---|
| 114 |
|
---|
| 115 | Input parameters:
|
---|
| 116 | A - matrix. Array whose indexes range within [0..N-1, 0..N-1].
|
---|
| 117 | N - size of matrix A.
|
---|
| 118 |
|
---|
| 119 | Result: 1/LowerBound(cond(A))
|
---|
| 120 | *************************************************************************/
|
---|
| 121 | public static double rmatrixrcondinf(ref double[,] a,
|
---|
| 122 | int n)
|
---|
| 123 | {
|
---|
| 124 | double result = 0;
|
---|
| 125 | int i = 0;
|
---|
| 126 | double[,] a1 = new double[0,0];
|
---|
| 127 | int i_ = 0;
|
---|
| 128 | int i1_ = 0;
|
---|
| 129 |
|
---|
| 130 | System.Diagnostics.Debug.Assert(n>=1, "RMatrixRCondInf: N<1!");
|
---|
| 131 | a1 = new double[n+1, n+1];
|
---|
| 132 | for(i=1; i<=n; i++)
|
---|
| 133 | {
|
---|
| 134 | i1_ = (0) - (1);
|
---|
| 135 | for(i_=1; i_<=n;i_++)
|
---|
| 136 | {
|
---|
| 137 | a1[i,i_] = a[i-1,i_+i1_];
|
---|
| 138 | }
|
---|
| 139 | }
|
---|
| 140 | result = rcondinf(a1, n);
|
---|
| 141 | return result;
|
---|
| 142 | }
|
---|
| 143 |
|
---|
| 144 |
|
---|
| 145 | /*************************************************************************
|
---|
| 146 | Estimate of the condition number of a matrix given by its LU decomposition
|
---|
| 147 | (infinity norm).
|
---|
| 148 |
|
---|
| 149 | The algorithm calculates a lower bound of the condition number. In this case,
|
---|
| 150 | the algorithm does not return a lower bound of the condition number, but an
|
---|
| 151 | inverse number (to avoid an overflow in case of a singular matrix).
|
---|
| 152 |
|
---|
| 153 | Input parameters:
|
---|
| 154 | LUDcmp - LU decomposition of a matrix in compact form. Output of
|
---|
| 155 | the RMatrixLU subroutine.
|
---|
| 156 | N - size of matrix A.
|
---|
| 157 |
|
---|
| 158 | Result: 1/LowerBound(cond(A))
|
---|
| 159 | *************************************************************************/
|
---|
| 160 | public static double rmatrixlurcondinf(ref double[,] ludcmp,
|
---|
| 161 | int n)
|
---|
| 162 | {
|
---|
| 163 | double result = 0;
|
---|
| 164 | int i = 0;
|
---|
| 165 | double[,] a1 = new double[0,0];
|
---|
| 166 | int i_ = 0;
|
---|
| 167 | int i1_ = 0;
|
---|
| 168 |
|
---|
| 169 | System.Diagnostics.Debug.Assert(n>=1, "RMatrixLURCondInf: N<1!");
|
---|
| 170 | a1 = new double[n+1, n+1];
|
---|
| 171 | for(i=1; i<=n; i++)
|
---|
| 172 | {
|
---|
| 173 | i1_ = (0) - (1);
|
---|
| 174 | for(i_=1; i_<=n;i_++)
|
---|
| 175 | {
|
---|
| 176 | a1[i,i_] = ludcmp[i-1,i_+i1_];
|
---|
| 177 | }
|
---|
| 178 | }
|
---|
| 179 | result = rcondinflu(ref a1, n);
|
---|
| 180 | return result;
|
---|
| 181 | }
|
---|
| 182 |
|
---|
| 183 |
|
---|
| 184 | public static double rcond1(double[,] a,
|
---|
| 185 | int n)
|
---|
| 186 | {
|
---|
| 187 | double result = 0;
|
---|
| 188 | int i = 0;
|
---|
| 189 | int j = 0;
|
---|
| 190 | double v = 0;
|
---|
| 191 | double nrm = 0;
|
---|
| 192 | int[] pivots = new int[0];
|
---|
| 193 |
|
---|
| 194 | a = (double[,])a.Clone();
|
---|
| 195 |
|
---|
| 196 | nrm = 0;
|
---|
| 197 | for(j=1; j<=n; j++)
|
---|
| 198 | {
|
---|
| 199 | v = 0;
|
---|
| 200 | for(i=1; i<=n; i++)
|
---|
| 201 | {
|
---|
| 202 | v = v+Math.Abs(a[i,j]);
|
---|
| 203 | }
|
---|
| 204 | nrm = Math.Max(nrm, v);
|
---|
| 205 | }
|
---|
| 206 | lu.ludecomposition(ref a, n, n, ref pivots);
|
---|
| 207 | internalestimatercondlu(ref a, n, true, true, nrm, ref v);
|
---|
| 208 | result = v;
|
---|
| 209 | return result;
|
---|
| 210 | }
|
---|
| 211 |
|
---|
| 212 |
|
---|
| 213 | public static double rcond1lu(ref double[,] ludcmp,
|
---|
| 214 | int n)
|
---|
| 215 | {
|
---|
| 216 | double result = 0;
|
---|
| 217 | double v = 0;
|
---|
| 218 |
|
---|
| 219 | internalestimatercondlu(ref ludcmp, n, true, false, 0, ref v);
|
---|
| 220 | result = v;
|
---|
| 221 | return result;
|
---|
| 222 | }
|
---|
| 223 |
|
---|
| 224 |
|
---|
| 225 | public static double rcondinf(double[,] a,
|
---|
| 226 | int n)
|
---|
| 227 | {
|
---|
| 228 | double result = 0;
|
---|
| 229 | int i = 0;
|
---|
| 230 | int j = 0;
|
---|
| 231 | double v = 0;
|
---|
| 232 | double nrm = 0;
|
---|
| 233 | int[] pivots = new int[0];
|
---|
| 234 |
|
---|
| 235 | a = (double[,])a.Clone();
|
---|
| 236 |
|
---|
| 237 | nrm = 0;
|
---|
| 238 | for(i=1; i<=n; i++)
|
---|
| 239 | {
|
---|
| 240 | v = 0;
|
---|
| 241 | for(j=1; j<=n; j++)
|
---|
| 242 | {
|
---|
| 243 | v = v+Math.Abs(a[i,j]);
|
---|
| 244 | }
|
---|
| 245 | nrm = Math.Max(nrm, v);
|
---|
| 246 | }
|
---|
| 247 | lu.ludecomposition(ref a, n, n, ref pivots);
|
---|
| 248 | internalestimatercondlu(ref a, n, false, true, nrm, ref v);
|
---|
| 249 | result = v;
|
---|
| 250 | return result;
|
---|
| 251 | }
|
---|
| 252 |
|
---|
| 253 |
|
---|
| 254 | public static double rcondinflu(ref double[,] ludcmp,
|
---|
| 255 | int n)
|
---|
| 256 | {
|
---|
| 257 | double result = 0;
|
---|
| 258 | double v = 0;
|
---|
| 259 |
|
---|
| 260 | internalestimatercondlu(ref ludcmp, n, false, false, 0, ref v);
|
---|
| 261 | result = v;
|
---|
| 262 | return result;
|
---|
| 263 | }
|
---|
| 264 |
|
---|
| 265 |
|
---|
| 266 | private static void internalestimatercondlu(ref double[,] ludcmp,
|
---|
| 267 | int n,
|
---|
| 268 | bool onenorm,
|
---|
| 269 | bool isanormprovided,
|
---|
| 270 | double anorm,
|
---|
| 271 | ref double rc)
|
---|
| 272 | {
|
---|
| 273 | double[] work0 = new double[0];
|
---|
| 274 | double[] work1 = new double[0];
|
---|
| 275 | double[] work2 = new double[0];
|
---|
| 276 | double[] work3 = new double[0];
|
---|
| 277 | int[] iwork = new int[0];
|
---|
| 278 | double v = 0;
|
---|
| 279 | bool normin = new bool();
|
---|
| 280 | int i = 0;
|
---|
| 281 | int im1 = 0;
|
---|
| 282 | int ip1 = 0;
|
---|
| 283 | int ix = 0;
|
---|
| 284 | int kase = 0;
|
---|
| 285 | int kase1 = 0;
|
---|
| 286 | double ainvnm = 0;
|
---|
| 287 | double ascale = 0;
|
---|
| 288 | double sl = 0;
|
---|
| 289 | double smlnum = 0;
|
---|
| 290 | double su = 0;
|
---|
| 291 | bool mupper = new bool();
|
---|
| 292 | bool mtrans = new bool();
|
---|
| 293 | bool munit = new bool();
|
---|
| 294 | int i_ = 0;
|
---|
| 295 |
|
---|
| 296 |
|
---|
| 297 | //
|
---|
| 298 | // Quick return if possible
|
---|
| 299 | //
|
---|
| 300 | if( n==0 )
|
---|
| 301 | {
|
---|
| 302 | rc = 1;
|
---|
| 303 | return;
|
---|
| 304 | }
|
---|
| 305 |
|
---|
| 306 | //
|
---|
| 307 | // init
|
---|
| 308 | //
|
---|
| 309 | if( onenorm )
|
---|
| 310 | {
|
---|
| 311 | kase1 = 1;
|
---|
| 312 | }
|
---|
| 313 | else
|
---|
| 314 | {
|
---|
| 315 | kase1 = 2;
|
---|
| 316 | }
|
---|
| 317 | mupper = true;
|
---|
| 318 | mtrans = true;
|
---|
| 319 | munit = true;
|
---|
| 320 | work0 = new double[n+1];
|
---|
| 321 | work1 = new double[n+1];
|
---|
| 322 | work2 = new double[n+1];
|
---|
| 323 | work3 = new double[n+1];
|
---|
| 324 | iwork = new int[n+1];
|
---|
| 325 |
|
---|
| 326 | //
|
---|
| 327 | // Estimate the norm of A.
|
---|
| 328 | //
|
---|
| 329 | if( !isanormprovided )
|
---|
| 330 | {
|
---|
| 331 | kase = 0;
|
---|
| 332 | anorm = 0;
|
---|
| 333 | while( true )
|
---|
| 334 | {
|
---|
| 335 | internalestimatenorm(n, ref work1, ref work0, ref iwork, ref anorm, ref kase);
|
---|
| 336 | if( kase==0 )
|
---|
| 337 | {
|
---|
| 338 | break;
|
---|
| 339 | }
|
---|
| 340 | if( kase==kase1 )
|
---|
| 341 | {
|
---|
| 342 |
|
---|
| 343 | //
|
---|
| 344 | // Multiply by U
|
---|
| 345 | //
|
---|
| 346 | for(i=1; i<=n; i++)
|
---|
| 347 | {
|
---|
| 348 | v = 0.0;
|
---|
| 349 | for(i_=i; i_<=n;i_++)
|
---|
| 350 | {
|
---|
| 351 | v += ludcmp[i,i_]*work0[i_];
|
---|
| 352 | }
|
---|
| 353 | work0[i] = v;
|
---|
| 354 | }
|
---|
| 355 |
|
---|
| 356 | //
|
---|
| 357 | // Multiply by L
|
---|
| 358 | //
|
---|
| 359 | for(i=n; i>=1; i--)
|
---|
| 360 | {
|
---|
| 361 | im1 = i-1;
|
---|
| 362 | if( i>1 )
|
---|
| 363 | {
|
---|
| 364 | v = 0.0;
|
---|
| 365 | for(i_=1; i_<=im1;i_++)
|
---|
| 366 | {
|
---|
| 367 | v += ludcmp[i,i_]*work0[i_];
|
---|
| 368 | }
|
---|
| 369 | }
|
---|
| 370 | else
|
---|
| 371 | {
|
---|
| 372 | v = 0;
|
---|
| 373 | }
|
---|
| 374 | work0[i] = work0[i]+v;
|
---|
| 375 | }
|
---|
| 376 | }
|
---|
| 377 | else
|
---|
| 378 | {
|
---|
| 379 |
|
---|
| 380 | //
|
---|
| 381 | // Multiply by L'
|
---|
| 382 | //
|
---|
| 383 | for(i=1; i<=n; i++)
|
---|
| 384 | {
|
---|
| 385 | ip1 = i+1;
|
---|
| 386 | v = 0.0;
|
---|
| 387 | for(i_=ip1; i_<=n;i_++)
|
---|
| 388 | {
|
---|
| 389 | v += ludcmp[i_,i]*work0[i_];
|
---|
| 390 | }
|
---|
| 391 | work0[i] = work0[i]+v;
|
---|
| 392 | }
|
---|
| 393 |
|
---|
| 394 | //
|
---|
| 395 | // Multiply by U'
|
---|
| 396 | //
|
---|
| 397 | for(i=n; i>=1; i--)
|
---|
| 398 | {
|
---|
| 399 | v = 0.0;
|
---|
| 400 | for(i_=1; i_<=i;i_++)
|
---|
| 401 | {
|
---|
| 402 | v += ludcmp[i_,i]*work0[i_];
|
---|
| 403 | }
|
---|
| 404 | work0[i] = v;
|
---|
| 405 | }
|
---|
| 406 | }
|
---|
| 407 | }
|
---|
| 408 | }
|
---|
| 409 |
|
---|
| 410 | //
|
---|
| 411 | // Quick return if possible
|
---|
| 412 | //
|
---|
| 413 | rc = 0;
|
---|
| 414 | if( (double)(anorm)==(double)(0) )
|
---|
| 415 | {
|
---|
| 416 | return;
|
---|
| 417 | }
|
---|
| 418 |
|
---|
| 419 | //
|
---|
| 420 | // Estimate the norm of inv(A).
|
---|
| 421 | //
|
---|
| 422 | smlnum = AP.Math.MinRealNumber;
|
---|
| 423 | ainvnm = 0;
|
---|
| 424 | normin = false;
|
---|
| 425 | kase = 0;
|
---|
| 426 | while( true )
|
---|
| 427 | {
|
---|
| 428 | internalestimatenorm(n, ref work1, ref work0, ref iwork, ref ainvnm, ref kase);
|
---|
| 429 | if( kase==0 )
|
---|
| 430 | {
|
---|
| 431 | break;
|
---|
| 432 | }
|
---|
| 433 | if( kase==kase1 )
|
---|
| 434 | {
|
---|
| 435 |
|
---|
| 436 | //
|
---|
| 437 | // Multiply by inv(L).
|
---|
| 438 | //
|
---|
| 439 | trlinsolve.safesolvetriangular(ref ludcmp, n, ref work0, ref sl, !mupper, !mtrans, munit, normin, ref work2);
|
---|
| 440 |
|
---|
| 441 | //
|
---|
| 442 | // Multiply by inv(U).
|
---|
| 443 | //
|
---|
| 444 | trlinsolve.safesolvetriangular(ref ludcmp, n, ref work0, ref su, mupper, !mtrans, !munit, normin, ref work3);
|
---|
| 445 | }
|
---|
| 446 | else
|
---|
| 447 | {
|
---|
| 448 |
|
---|
| 449 | //
|
---|
| 450 | // Multiply by inv(U').
|
---|
| 451 | //
|
---|
| 452 | trlinsolve.safesolvetriangular(ref ludcmp, n, ref work0, ref su, mupper, mtrans, !munit, normin, ref work3);
|
---|
| 453 |
|
---|
| 454 | //
|
---|
| 455 | // Multiply by inv(L').
|
---|
| 456 | //
|
---|
| 457 | trlinsolve.safesolvetriangular(ref ludcmp, n, ref work0, ref sl, !mupper, mtrans, munit, normin, ref work2);
|
---|
| 458 | }
|
---|
| 459 |
|
---|
| 460 | //
|
---|
| 461 | // Divide X by 1/(SL*SU) if doing so will not cause overflow.
|
---|
| 462 | //
|
---|
| 463 | ascale = sl*su;
|
---|
| 464 | normin = true;
|
---|
| 465 | if( (double)(ascale)!=(double)(1) )
|
---|
| 466 | {
|
---|
| 467 | ix = 1;
|
---|
| 468 | for(i=2; i<=n; i++)
|
---|
| 469 | {
|
---|
| 470 | if( (double)(Math.Abs(work0[i]))>(double)(Math.Abs(work0[ix])) )
|
---|
| 471 | {
|
---|
| 472 | ix = i;
|
---|
| 473 | }
|
---|
| 474 | }
|
---|
| 475 | if( (double)(ascale)<(double)(Math.Abs(work0[ix])*smlnum) | (double)(ascale)==(double)(0) )
|
---|
| 476 | {
|
---|
| 477 | return;
|
---|
| 478 | }
|
---|
| 479 | for(i=1; i<=n; i++)
|
---|
| 480 | {
|
---|
| 481 | work0[i] = work0[i]/ascale;
|
---|
| 482 | }
|
---|
| 483 | }
|
---|
| 484 | }
|
---|
| 485 |
|
---|
| 486 | //
|
---|
| 487 | // Compute the estimate of the reciprocal condition number.
|
---|
| 488 | //
|
---|
| 489 | if( (double)(ainvnm)!=(double)(0) )
|
---|
| 490 | {
|
---|
| 491 | rc = 1/ainvnm;
|
---|
| 492 | rc = rc/anorm;
|
---|
| 493 | }
|
---|
| 494 | }
|
---|
| 495 |
|
---|
| 496 |
|
---|
| 497 | private static void internalestimatenorm(int n,
|
---|
| 498 | ref double[] v,
|
---|
| 499 | ref double[] x,
|
---|
| 500 | ref int[] isgn,
|
---|
| 501 | ref double est,
|
---|
| 502 | ref int kase)
|
---|
| 503 | {
|
---|
| 504 | int itmax = 0;
|
---|
| 505 | int i = 0;
|
---|
| 506 | double t = 0;
|
---|
| 507 | bool flg = new bool();
|
---|
| 508 | int positer = 0;
|
---|
| 509 | int posj = 0;
|
---|
| 510 | int posjlast = 0;
|
---|
| 511 | int posjump = 0;
|
---|
| 512 | int posaltsgn = 0;
|
---|
| 513 | int posestold = 0;
|
---|
| 514 | int postemp = 0;
|
---|
| 515 | int i_ = 0;
|
---|
| 516 |
|
---|
| 517 | itmax = 5;
|
---|
| 518 | posaltsgn = n+1;
|
---|
| 519 | posestold = n+2;
|
---|
| 520 | postemp = n+3;
|
---|
| 521 | positer = n+1;
|
---|
| 522 | posj = n+2;
|
---|
| 523 | posjlast = n+3;
|
---|
| 524 | posjump = n+4;
|
---|
| 525 | if( kase==0 )
|
---|
| 526 | {
|
---|
| 527 | v = new double[n+3+1];
|
---|
| 528 | x = new double[n+1];
|
---|
| 529 | isgn = new int[n+4+1];
|
---|
| 530 | t = (double)(1)/(double)(n);
|
---|
| 531 | for(i=1; i<=n; i++)
|
---|
| 532 | {
|
---|
| 533 | x[i] = t;
|
---|
| 534 | }
|
---|
| 535 | kase = 1;
|
---|
| 536 | isgn[posjump] = 1;
|
---|
| 537 | return;
|
---|
| 538 | }
|
---|
| 539 |
|
---|
| 540 | //
|
---|
| 541 | // ................ ENTRY (JUMP = 1)
|
---|
| 542 | // FIRST ITERATION. X HAS BEEN OVERWRITTEN BY A*X.
|
---|
| 543 | //
|
---|
| 544 | if( isgn[posjump]==1 )
|
---|
| 545 | {
|
---|
| 546 | if( n==1 )
|
---|
| 547 | {
|
---|
| 548 | v[1] = x[1];
|
---|
| 549 | est = Math.Abs(v[1]);
|
---|
| 550 | kase = 0;
|
---|
| 551 | return;
|
---|
| 552 | }
|
---|
| 553 | est = 0;
|
---|
| 554 | for(i=1; i<=n; i++)
|
---|
| 555 | {
|
---|
| 556 | est = est+Math.Abs(x[i]);
|
---|
| 557 | }
|
---|
| 558 | for(i=1; i<=n; i++)
|
---|
| 559 | {
|
---|
| 560 | if( (double)(x[i])>=(double)(0) )
|
---|
| 561 | {
|
---|
| 562 | x[i] = 1;
|
---|
| 563 | }
|
---|
| 564 | else
|
---|
| 565 | {
|
---|
| 566 | x[i] = -1;
|
---|
| 567 | }
|
---|
| 568 | isgn[i] = Math.Sign(x[i]);
|
---|
| 569 | }
|
---|
| 570 | kase = 2;
|
---|
| 571 | isgn[posjump] = 2;
|
---|
| 572 | return;
|
---|
| 573 | }
|
---|
| 574 |
|
---|
| 575 | //
|
---|
| 576 | // ................ ENTRY (JUMP = 2)
|
---|
| 577 | // FIRST ITERATION. X HAS BEEN OVERWRITTEN BY TRANDPOSE(A)*X.
|
---|
| 578 | //
|
---|
| 579 | if( isgn[posjump]==2 )
|
---|
| 580 | {
|
---|
| 581 | isgn[posj] = 1;
|
---|
| 582 | for(i=2; i<=n; i++)
|
---|
| 583 | {
|
---|
| 584 | if( (double)(Math.Abs(x[i]))>(double)(Math.Abs(x[isgn[posj]])) )
|
---|
| 585 | {
|
---|
| 586 | isgn[posj] = i;
|
---|
| 587 | }
|
---|
| 588 | }
|
---|
| 589 | isgn[positer] = 2;
|
---|
| 590 |
|
---|
| 591 | //
|
---|
| 592 | // MAIN LOOP - ITERATIONS 2,3,...,ITMAX.
|
---|
| 593 | //
|
---|
| 594 | for(i=1; i<=n; i++)
|
---|
| 595 | {
|
---|
| 596 | x[i] = 0;
|
---|
| 597 | }
|
---|
| 598 | x[isgn[posj]] = 1;
|
---|
| 599 | kase = 1;
|
---|
| 600 | isgn[posjump] = 3;
|
---|
| 601 | return;
|
---|
| 602 | }
|
---|
| 603 |
|
---|
| 604 | //
|
---|
| 605 | // ................ ENTRY (JUMP = 3)
|
---|
| 606 | // X HAS BEEN OVERWRITTEN BY A*X.
|
---|
| 607 | //
|
---|
| 608 | if( isgn[posjump]==3 )
|
---|
| 609 | {
|
---|
| 610 | for(i_=1; i_<=n;i_++)
|
---|
| 611 | {
|
---|
| 612 | v[i_] = x[i_];
|
---|
| 613 | }
|
---|
| 614 | v[posestold] = est;
|
---|
| 615 | est = 0;
|
---|
| 616 | for(i=1; i<=n; i++)
|
---|
| 617 | {
|
---|
| 618 | est = est+Math.Abs(v[i]);
|
---|
| 619 | }
|
---|
| 620 | flg = false;
|
---|
| 621 | for(i=1; i<=n; i++)
|
---|
| 622 | {
|
---|
| 623 | if( (double)(x[i])>=(double)(0) & isgn[i]<0 | (double)(x[i])<(double)(0) & isgn[i]>=0 )
|
---|
| 624 | {
|
---|
| 625 | flg = true;
|
---|
| 626 | }
|
---|
| 627 | }
|
---|
| 628 |
|
---|
| 629 | //
|
---|
| 630 | // REPEATED SIGN VECTOR DETECTED, HENCE ALGORITHM HAS CONVERGED.
|
---|
| 631 | // OR MAY BE CYCLING.
|
---|
| 632 | //
|
---|
| 633 | if( !flg | (double)(est)<=(double)(v[posestold]) )
|
---|
| 634 | {
|
---|
| 635 | v[posaltsgn] = 1;
|
---|
| 636 | for(i=1; i<=n; i++)
|
---|
| 637 | {
|
---|
| 638 | x[i] = v[posaltsgn]*(1+((double)(i-1))/((double)(n-1)));
|
---|
| 639 | v[posaltsgn] = -v[posaltsgn];
|
---|
| 640 | }
|
---|
| 641 | kase = 1;
|
---|
| 642 | isgn[posjump] = 5;
|
---|
| 643 | return;
|
---|
| 644 | }
|
---|
| 645 | for(i=1; i<=n; i++)
|
---|
| 646 | {
|
---|
| 647 | if( (double)(x[i])>=(double)(0) )
|
---|
| 648 | {
|
---|
| 649 | x[i] = 1;
|
---|
| 650 | isgn[i] = 1;
|
---|
| 651 | }
|
---|
| 652 | else
|
---|
| 653 | {
|
---|
| 654 | x[i] = -1;
|
---|
| 655 | isgn[i] = -1;
|
---|
| 656 | }
|
---|
| 657 | }
|
---|
| 658 | kase = 2;
|
---|
| 659 | isgn[posjump] = 4;
|
---|
| 660 | return;
|
---|
| 661 | }
|
---|
| 662 |
|
---|
| 663 | //
|
---|
| 664 | // ................ ENTRY (JUMP = 4)
|
---|
| 665 | // X HAS BEEN OVERWRITTEN BY TRANDPOSE(A)*X.
|
---|
| 666 | //
|
---|
| 667 | if( isgn[posjump]==4 )
|
---|
| 668 | {
|
---|
| 669 | isgn[posjlast] = isgn[posj];
|
---|
| 670 | isgn[posj] = 1;
|
---|
| 671 | for(i=2; i<=n; i++)
|
---|
| 672 | {
|
---|
| 673 | if( (double)(Math.Abs(x[i]))>(double)(Math.Abs(x[isgn[posj]])) )
|
---|
| 674 | {
|
---|
| 675 | isgn[posj] = i;
|
---|
| 676 | }
|
---|
| 677 | }
|
---|
| 678 | if( (double)(x[isgn[posjlast]])!=(double)(Math.Abs(x[isgn[posj]])) & isgn[positer]<itmax )
|
---|
| 679 | {
|
---|
| 680 | isgn[positer] = isgn[positer]+1;
|
---|
| 681 | for(i=1; i<=n; i++)
|
---|
| 682 | {
|
---|
| 683 | x[i] = 0;
|
---|
| 684 | }
|
---|
| 685 | x[isgn[posj]] = 1;
|
---|
| 686 | kase = 1;
|
---|
| 687 | isgn[posjump] = 3;
|
---|
| 688 | return;
|
---|
| 689 | }
|
---|
| 690 |
|
---|
| 691 | //
|
---|
| 692 | // ITERATION COMPLETE. FINAL STAGE.
|
---|
| 693 | //
|
---|
| 694 | v[posaltsgn] = 1;
|
---|
| 695 | for(i=1; i<=n; i++)
|
---|
| 696 | {
|
---|
| 697 | x[i] = v[posaltsgn]*(1+((double)(i-1))/((double)(n-1)));
|
---|
| 698 | v[posaltsgn] = -v[posaltsgn];
|
---|
| 699 | }
|
---|
| 700 | kase = 1;
|
---|
| 701 | isgn[posjump] = 5;
|
---|
| 702 | return;
|
---|
| 703 | }
|
---|
| 704 |
|
---|
| 705 | //
|
---|
| 706 | // ................ ENTRY (JUMP = 5)
|
---|
| 707 | // X HAS BEEN OVERWRITTEN BY A*X.
|
---|
| 708 | //
|
---|
| 709 | if( isgn[posjump]==5 )
|
---|
| 710 | {
|
---|
| 711 | v[postemp] = 0;
|
---|
| 712 | for(i=1; i<=n; i++)
|
---|
| 713 | {
|
---|
| 714 | v[postemp] = v[postemp]+Math.Abs(x[i]);
|
---|
| 715 | }
|
---|
| 716 | v[postemp] = 2*v[postemp]/(3*n);
|
---|
| 717 | if( (double)(v[postemp])>(double)(est) )
|
---|
| 718 | {
|
---|
| 719 | for(i_=1; i_<=n;i_++)
|
---|
| 720 | {
|
---|
| 721 | v[i_] = x[i_];
|
---|
| 722 | }
|
---|
| 723 | est = v[postemp];
|
---|
| 724 | }
|
---|
| 725 | kase = 0;
|
---|
| 726 | return;
|
---|
| 727 | }
|
---|
| 728 | }
|
---|
| 729 | }
|
---|
| 730 | }
|
---|