| 1 | /*************************************************************************
|
|---|
| 2 | This file is a part of ALGLIB project.
|
|---|
| 3 |
|
|---|
| 4 | >>> SOURCE LICENSE >>>
|
|---|
| 5 | This program is free software; you can redistribute it and/or modify
|
|---|
| 6 | it under the terms of the GNU General Public License as published by
|
|---|
| 7 | the Free Software Foundation (www.fsf.org); either version 2 of the
|
|---|
| 8 | License, or (at your option) any later version.
|
|---|
| 9 |
|
|---|
| 10 | This program is distributed in the hope that it will be useful,
|
|---|
| 11 | but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|---|
| 12 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|---|
| 13 | GNU General Public License for more details.
|
|---|
| 14 |
|
|---|
| 15 | A copy of the GNU General Public License is available at
|
|---|
| 16 | http://www.fsf.org/licensing/licenses
|
|---|
| 17 |
|
|---|
| 18 | >>> END OF LICENSE >>>
|
|---|
| 19 | *************************************************************************/
|
|---|
| 20 |
|
|---|
| 21 | using System;
|
|---|
| 22 |
|
|---|
| 23 | namespace alglib
|
|---|
| 24 | {
|
|---|
| 25 | public class trlinsolve
|
|---|
| 26 | {
|
|---|
| 27 | /*************************************************************************
|
|---|
| 28 | Utility subroutine performing the "safe" solution of system of linear
|
|---|
| 29 | equations with triangular coefficient matrices.
|
|---|
| 30 |
|
|---|
| 31 | The subroutine uses scaling and solves the scaled system A*x=s*b (where s
|
|---|
| 32 | is a scalar value) instead of A*x=b, choosing s so that x can be
|
|---|
| 33 | represented by a floating-point number. The closer the system gets to a
|
|---|
| 34 | singular, the less s is. If the system is singular, s=0 and x contains the
|
|---|
| 35 | non-trivial solution of equation A*x=0.
|
|---|
| 36 |
|
|---|
| 37 | The feature of an algorithm is that it could not cause an overflow or a
|
|---|
| 38 | division by zero regardless of the matrix used as the input.
|
|---|
| 39 |
|
|---|
| 40 | The algorithm can solve systems of equations with upper/lower triangular
|
|---|
| 41 | matrices, with/without unit diagonal, and systems of type A*x=b or A'*x=b
|
|---|
| 42 | (where A' is a transposed matrix A).
|
|---|
| 43 |
|
|---|
| 44 | Input parameters:
|
|---|
| 45 | A - system matrix. Array whose indexes range within [0..N-1, 0..N-1].
|
|---|
| 46 | N - size of matrix A.
|
|---|
| 47 | X - right-hand member of a system.
|
|---|
| 48 | Array whose index ranges within [0..N-1].
|
|---|
| 49 | IsUpper - matrix type. If it is True, the system matrix is the upper
|
|---|
| 50 | triangular and is located in the corresponding part of
|
|---|
| 51 | matrix A.
|
|---|
| 52 | Trans - problem type. If it is True, the problem to be solved is
|
|---|
| 53 | A'*x=b, otherwise it is A*x=b.
|
|---|
| 54 | Isunit - matrix type. If it is True, the system matrix has a unit
|
|---|
| 55 | diagonal (the elements on the main diagonal are not used
|
|---|
| 56 | in the calculation process), otherwise the matrix is considered
|
|---|
| 57 | to be a general triangular matrix.
|
|---|
| 58 |
|
|---|
| 59 | Output parameters:
|
|---|
| 60 | X - solution. Array whose index ranges within [0..N-1].
|
|---|
| 61 | S - scaling factor.
|
|---|
| 62 |
|
|---|
| 63 | -- LAPACK auxiliary routine (version 3.0) --
|
|---|
| 64 | Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
|
|---|
| 65 | Courant Institute, Argonne National Lab, and Rice University
|
|---|
| 66 | June 30, 1992
|
|---|
| 67 | *************************************************************************/
|
|---|
| 68 | public static void rmatrixtrsafesolve(ref double[,] a,
|
|---|
| 69 | int n,
|
|---|
| 70 | ref double[] x,
|
|---|
| 71 | ref double s,
|
|---|
| 72 | bool isupper,
|
|---|
| 73 | bool istrans,
|
|---|
| 74 | bool isunit)
|
|---|
| 75 | {
|
|---|
| 76 | bool normin = new bool();
|
|---|
| 77 | double[] cnorm = new double[0];
|
|---|
| 78 | double[,] a1 = new double[0,0];
|
|---|
| 79 | double[] x1 = new double[0];
|
|---|
| 80 | int i = 0;
|
|---|
| 81 | int i_ = 0;
|
|---|
| 82 | int i1_ = 0;
|
|---|
| 83 |
|
|---|
| 84 |
|
|---|
| 85 | //
|
|---|
| 86 | // From 0-based to 1-based
|
|---|
| 87 | //
|
|---|
| 88 | normin = false;
|
|---|
| 89 | a1 = new double[n+1, n+1];
|
|---|
| 90 | x1 = new double[n+1];
|
|---|
| 91 | for(i=1; i<=n; i++)
|
|---|
| 92 | {
|
|---|
| 93 | i1_ = (0) - (1);
|
|---|
| 94 | for(i_=1; i_<=n;i_++)
|
|---|
| 95 | {
|
|---|
| 96 | a1[i,i_] = a[i-1,i_+i1_];
|
|---|
| 97 | }
|
|---|
| 98 | }
|
|---|
| 99 | i1_ = (0) - (1);
|
|---|
| 100 | for(i_=1; i_<=n;i_++)
|
|---|
| 101 | {
|
|---|
| 102 | x1[i_] = x[i_+i1_];
|
|---|
| 103 | }
|
|---|
| 104 |
|
|---|
| 105 | //
|
|---|
| 106 | // Solve 1-based
|
|---|
| 107 | //
|
|---|
| 108 | safesolvetriangular(ref a1, n, ref x1, ref s, isupper, istrans, isunit, normin, ref cnorm);
|
|---|
| 109 |
|
|---|
| 110 | //
|
|---|
| 111 | // From 1-based to 0-based
|
|---|
| 112 | //
|
|---|
| 113 | i1_ = (1) - (0);
|
|---|
| 114 | for(i_=0; i_<=n-1;i_++)
|
|---|
| 115 | {
|
|---|
| 116 | x[i_] = x1[i_+i1_];
|
|---|
| 117 | }
|
|---|
| 118 | }
|
|---|
| 119 |
|
|---|
| 120 |
|
|---|
| 121 | /*************************************************************************
|
|---|
| 122 | Obsolete 1-based subroutine.
|
|---|
| 123 | See RMatrixTRSafeSolve for 0-based replacement.
|
|---|
| 124 | *************************************************************************/
|
|---|
| 125 | public static void safesolvetriangular(ref double[,] a,
|
|---|
| 126 | int n,
|
|---|
| 127 | ref double[] x,
|
|---|
| 128 | ref double s,
|
|---|
| 129 | bool isupper,
|
|---|
| 130 | bool istrans,
|
|---|
| 131 | bool isunit,
|
|---|
| 132 | bool normin,
|
|---|
| 133 | ref double[] cnorm)
|
|---|
| 134 | {
|
|---|
| 135 | int i = 0;
|
|---|
| 136 | int imax = 0;
|
|---|
| 137 | int j = 0;
|
|---|
| 138 | int jfirst = 0;
|
|---|
| 139 | int jinc = 0;
|
|---|
| 140 | int jlast = 0;
|
|---|
| 141 | int jm1 = 0;
|
|---|
| 142 | int jp1 = 0;
|
|---|
| 143 | int ip1 = 0;
|
|---|
| 144 | int im1 = 0;
|
|---|
| 145 | int k = 0;
|
|---|
| 146 | int flg = 0;
|
|---|
| 147 | double v = 0;
|
|---|
| 148 | double vd = 0;
|
|---|
| 149 | double bignum = 0;
|
|---|
| 150 | double grow = 0;
|
|---|
| 151 | double rec = 0;
|
|---|
| 152 | double smlnum = 0;
|
|---|
| 153 | double sumj = 0;
|
|---|
| 154 | double tjj = 0;
|
|---|
| 155 | double tjjs = 0;
|
|---|
| 156 | double tmax = 0;
|
|---|
| 157 | double tscal = 0;
|
|---|
| 158 | double uscal = 0;
|
|---|
| 159 | double xbnd = 0;
|
|---|
| 160 | double xj = 0;
|
|---|
| 161 | double xmax = 0;
|
|---|
| 162 | bool notran = new bool();
|
|---|
| 163 | bool upper = new bool();
|
|---|
| 164 | bool nounit = new bool();
|
|---|
| 165 | int i_ = 0;
|
|---|
| 166 |
|
|---|
| 167 | upper = isupper;
|
|---|
| 168 | notran = !istrans;
|
|---|
| 169 | nounit = !isunit;
|
|---|
| 170 |
|
|---|
| 171 | //
|
|---|
| 172 | // Quick return if possible
|
|---|
| 173 | //
|
|---|
| 174 | if( n==0 )
|
|---|
| 175 | {
|
|---|
| 176 | return;
|
|---|
| 177 | }
|
|---|
| 178 |
|
|---|
| 179 | //
|
|---|
| 180 | // Determine machine dependent parameters to control overflow.
|
|---|
| 181 | //
|
|---|
| 182 | smlnum = AP.Math.MinRealNumber/(AP.Math.MachineEpsilon*2);
|
|---|
| 183 | bignum = 1/smlnum;
|
|---|
| 184 | s = 1;
|
|---|
| 185 | if( !normin )
|
|---|
| 186 | {
|
|---|
| 187 | cnorm = new double[n+1];
|
|---|
| 188 |
|
|---|
| 189 | //
|
|---|
| 190 | // Compute the 1-norm of each column, not including the diagonal.
|
|---|
| 191 | //
|
|---|
| 192 | if( upper )
|
|---|
| 193 | {
|
|---|
| 194 |
|
|---|
| 195 | //
|
|---|
| 196 | // A is upper triangular.
|
|---|
| 197 | //
|
|---|
| 198 | for(j=1; j<=n; j++)
|
|---|
| 199 | {
|
|---|
| 200 | v = 0;
|
|---|
| 201 | for(k=1; k<=j-1; k++)
|
|---|
| 202 | {
|
|---|
| 203 | v = v+Math.Abs(a[k,j]);
|
|---|
| 204 | }
|
|---|
| 205 | cnorm[j] = v;
|
|---|
| 206 | }
|
|---|
| 207 | }
|
|---|
| 208 | else
|
|---|
| 209 | {
|
|---|
| 210 |
|
|---|
| 211 | //
|
|---|
| 212 | // A is lower triangular.
|
|---|
| 213 | //
|
|---|
| 214 | for(j=1; j<=n-1; j++)
|
|---|
| 215 | {
|
|---|
| 216 | v = 0;
|
|---|
| 217 | for(k=j+1; k<=n; k++)
|
|---|
| 218 | {
|
|---|
| 219 | v = v+Math.Abs(a[k,j]);
|
|---|
| 220 | }
|
|---|
| 221 | cnorm[j] = v;
|
|---|
| 222 | }
|
|---|
| 223 | cnorm[n] = 0;
|
|---|
| 224 | }
|
|---|
| 225 | }
|
|---|
| 226 |
|
|---|
| 227 | //
|
|---|
| 228 | // Scale the column norms by TSCAL if the maximum element in CNORM is
|
|---|
| 229 | // greater than BIGNUM.
|
|---|
| 230 | //
|
|---|
| 231 | imax = 1;
|
|---|
| 232 | for(k=2; k<=n; k++)
|
|---|
| 233 | {
|
|---|
| 234 | if( (double)(cnorm[k])>(double)(cnorm[imax]) )
|
|---|
| 235 | {
|
|---|
| 236 | imax = k;
|
|---|
| 237 | }
|
|---|
| 238 | }
|
|---|
| 239 | tmax = cnorm[imax];
|
|---|
| 240 | if( (double)(tmax)<=(double)(bignum) )
|
|---|
| 241 | {
|
|---|
| 242 | tscal = 1;
|
|---|
| 243 | }
|
|---|
| 244 | else
|
|---|
| 245 | {
|
|---|
| 246 | tscal = 1/(smlnum*tmax);
|
|---|
| 247 | for(i_=1; i_<=n;i_++)
|
|---|
| 248 | {
|
|---|
| 249 | cnorm[i_] = tscal*cnorm[i_];
|
|---|
| 250 | }
|
|---|
| 251 | }
|
|---|
| 252 |
|
|---|
| 253 | //
|
|---|
| 254 | // Compute a bound on the computed solution vector to see if the
|
|---|
| 255 | // Level 2 BLAS routine DTRSV can be used.
|
|---|
| 256 | //
|
|---|
| 257 | j = 1;
|
|---|
| 258 | for(k=2; k<=n; k++)
|
|---|
| 259 | {
|
|---|
| 260 | if( (double)(Math.Abs(x[k]))>(double)(Math.Abs(x[j])) )
|
|---|
| 261 | {
|
|---|
| 262 | j = k;
|
|---|
| 263 | }
|
|---|
| 264 | }
|
|---|
| 265 | xmax = Math.Abs(x[j]);
|
|---|
| 266 | xbnd = xmax;
|
|---|
| 267 | if( notran )
|
|---|
| 268 | {
|
|---|
| 269 |
|
|---|
| 270 | //
|
|---|
| 271 | // Compute the growth in A * x = b.
|
|---|
| 272 | //
|
|---|
| 273 | if( upper )
|
|---|
| 274 | {
|
|---|
| 275 | jfirst = n;
|
|---|
| 276 | jlast = 1;
|
|---|
| 277 | jinc = -1;
|
|---|
| 278 | }
|
|---|
| 279 | else
|
|---|
| 280 | {
|
|---|
| 281 | jfirst = 1;
|
|---|
| 282 | jlast = n;
|
|---|
| 283 | jinc = 1;
|
|---|
| 284 | }
|
|---|
| 285 | if( (double)(tscal)!=(double)(1) )
|
|---|
| 286 | {
|
|---|
| 287 | grow = 0;
|
|---|
| 288 | }
|
|---|
| 289 | else
|
|---|
| 290 | {
|
|---|
| 291 | if( nounit )
|
|---|
| 292 | {
|
|---|
| 293 |
|
|---|
| 294 | //
|
|---|
| 295 | // A is non-unit triangular.
|
|---|
| 296 | //
|
|---|
| 297 | // Compute GROW = 1/G(j) and XBND = 1/M(j).
|
|---|
| 298 | // Initially, G(0) = max{x(i), i=1,...,n}.
|
|---|
| 299 | //
|
|---|
| 300 | grow = 1/Math.Max(xbnd, smlnum);
|
|---|
| 301 | xbnd = grow;
|
|---|
| 302 | j = jfirst;
|
|---|
| 303 | while( jinc>0 & j<=jlast | jinc<0 & j>=jlast )
|
|---|
| 304 | {
|
|---|
| 305 |
|
|---|
| 306 | //
|
|---|
| 307 | // Exit the loop if the growth factor is too small.
|
|---|
| 308 | //
|
|---|
| 309 | if( (double)(grow)<=(double)(smlnum) )
|
|---|
| 310 | {
|
|---|
| 311 | break;
|
|---|
| 312 | }
|
|---|
| 313 |
|
|---|
| 314 | //
|
|---|
| 315 | // M(j) = G(j-1) / abs(A(j,j))
|
|---|
| 316 | //
|
|---|
| 317 | tjj = Math.Abs(a[j,j]);
|
|---|
| 318 | xbnd = Math.Min(xbnd, Math.Min(1, tjj)*grow);
|
|---|
| 319 | if( (double)(tjj+cnorm[j])>=(double)(smlnum) )
|
|---|
| 320 | {
|
|---|
| 321 |
|
|---|
| 322 | //
|
|---|
| 323 | // G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) )
|
|---|
| 324 | //
|
|---|
| 325 | grow = grow*(tjj/(tjj+cnorm[j]));
|
|---|
| 326 | }
|
|---|
| 327 | else
|
|---|
| 328 | {
|
|---|
| 329 |
|
|---|
| 330 | //
|
|---|
| 331 | // G(j) could overflow, set GROW to 0.
|
|---|
| 332 | //
|
|---|
| 333 | grow = 0;
|
|---|
| 334 | }
|
|---|
| 335 | if( j==jlast )
|
|---|
| 336 | {
|
|---|
| 337 | grow = xbnd;
|
|---|
| 338 | }
|
|---|
| 339 | j = j+jinc;
|
|---|
| 340 | }
|
|---|
| 341 | }
|
|---|
| 342 | else
|
|---|
| 343 | {
|
|---|
| 344 |
|
|---|
| 345 | //
|
|---|
| 346 | // A is unit triangular.
|
|---|
| 347 | //
|
|---|
| 348 | // Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
|
|---|
| 349 | //
|
|---|
| 350 | grow = Math.Min(1, 1/Math.Max(xbnd, smlnum));
|
|---|
| 351 | j = jfirst;
|
|---|
| 352 | while( jinc>0 & j<=jlast | jinc<0 & j>=jlast )
|
|---|
| 353 | {
|
|---|
| 354 |
|
|---|
| 355 | //
|
|---|
| 356 | // Exit the loop if the growth factor is too small.
|
|---|
| 357 | //
|
|---|
| 358 | if( (double)(grow)<=(double)(smlnum) )
|
|---|
| 359 | {
|
|---|
| 360 | break;
|
|---|
| 361 | }
|
|---|
| 362 |
|
|---|
| 363 | //
|
|---|
| 364 | // G(j) = G(j-1)*( 1 + CNORM(j) )
|
|---|
| 365 | //
|
|---|
| 366 | grow = grow*(1/(1+cnorm[j]));
|
|---|
| 367 | j = j+jinc;
|
|---|
| 368 | }
|
|---|
| 369 | }
|
|---|
| 370 | }
|
|---|
| 371 | }
|
|---|
| 372 | else
|
|---|
| 373 | {
|
|---|
| 374 |
|
|---|
| 375 | //
|
|---|
| 376 | // Compute the growth in A' * x = b.
|
|---|
| 377 | //
|
|---|
| 378 | if( upper )
|
|---|
| 379 | {
|
|---|
| 380 | jfirst = 1;
|
|---|
| 381 | jlast = n;
|
|---|
| 382 | jinc = 1;
|
|---|
| 383 | }
|
|---|
| 384 | else
|
|---|
| 385 | {
|
|---|
| 386 | jfirst = n;
|
|---|
| 387 | jlast = 1;
|
|---|
| 388 | jinc = -1;
|
|---|
| 389 | }
|
|---|
| 390 | if( (double)(tscal)!=(double)(1) )
|
|---|
| 391 | {
|
|---|
| 392 | grow = 0;
|
|---|
| 393 | }
|
|---|
| 394 | else
|
|---|
| 395 | {
|
|---|
| 396 | if( nounit )
|
|---|
| 397 | {
|
|---|
| 398 |
|
|---|
| 399 | //
|
|---|
| 400 | // A is non-unit triangular.
|
|---|
| 401 | //
|
|---|
| 402 | // Compute GROW = 1/G(j) and XBND = 1/M(j).
|
|---|
| 403 | // Initially, M(0) = max{x(i), i=1,...,n}.
|
|---|
| 404 | //
|
|---|
| 405 | grow = 1/Math.Max(xbnd, smlnum);
|
|---|
| 406 | xbnd = grow;
|
|---|
| 407 | j = jfirst;
|
|---|
| 408 | while( jinc>0 & j<=jlast | jinc<0 & j>=jlast )
|
|---|
| 409 | {
|
|---|
| 410 |
|
|---|
| 411 | //
|
|---|
| 412 | // Exit the loop if the growth factor is too small.
|
|---|
| 413 | //
|
|---|
| 414 | if( (double)(grow)<=(double)(smlnum) )
|
|---|
| 415 | {
|
|---|
| 416 | break;
|
|---|
| 417 | }
|
|---|
| 418 |
|
|---|
| 419 | //
|
|---|
| 420 | // G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) )
|
|---|
| 421 | //
|
|---|
| 422 | xj = 1+cnorm[j];
|
|---|
| 423 | grow = Math.Min(grow, xbnd/xj);
|
|---|
| 424 |
|
|---|
| 425 | //
|
|---|
| 426 | // M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j))
|
|---|
| 427 | //
|
|---|
| 428 | tjj = Math.Abs(a[j,j]);
|
|---|
| 429 | if( (double)(xj)>(double)(tjj) )
|
|---|
| 430 | {
|
|---|
| 431 | xbnd = xbnd*(tjj/xj);
|
|---|
| 432 | }
|
|---|
| 433 | if( j==jlast )
|
|---|
| 434 | {
|
|---|
| 435 | grow = Math.Min(grow, xbnd);
|
|---|
| 436 | }
|
|---|
| 437 | j = j+jinc;
|
|---|
| 438 | }
|
|---|
| 439 | }
|
|---|
| 440 | else
|
|---|
| 441 | {
|
|---|
| 442 |
|
|---|
| 443 | //
|
|---|
| 444 | // A is unit triangular.
|
|---|
| 445 | //
|
|---|
| 446 | // Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
|
|---|
| 447 | //
|
|---|
| 448 | grow = Math.Min(1, 1/Math.Max(xbnd, smlnum));
|
|---|
| 449 | j = jfirst;
|
|---|
| 450 | while( jinc>0 & j<=jlast | jinc<0 & j>=jlast )
|
|---|
| 451 | {
|
|---|
| 452 |
|
|---|
| 453 | //
|
|---|
| 454 | // Exit the loop if the growth factor is too small.
|
|---|
| 455 | //
|
|---|
| 456 | if( (double)(grow)<=(double)(smlnum) )
|
|---|
| 457 | {
|
|---|
| 458 | break;
|
|---|
| 459 | }
|
|---|
| 460 |
|
|---|
| 461 | //
|
|---|
| 462 | // G(j) = ( 1 + CNORM(j) )*G(j-1)
|
|---|
| 463 | //
|
|---|
| 464 | xj = 1+cnorm[j];
|
|---|
| 465 | grow = grow/xj;
|
|---|
| 466 | j = j+jinc;
|
|---|
| 467 | }
|
|---|
| 468 | }
|
|---|
| 469 | }
|
|---|
| 470 | }
|
|---|
| 471 | if( (double)(grow*tscal)>(double)(smlnum) )
|
|---|
| 472 | {
|
|---|
| 473 |
|
|---|
| 474 | //
|
|---|
| 475 | // Use the Level 2 BLAS solve if the reciprocal of the bound on
|
|---|
| 476 | // elements of X is not too small.
|
|---|
| 477 | //
|
|---|
| 478 | if( upper & notran | !upper & !notran )
|
|---|
| 479 | {
|
|---|
| 480 | if( nounit )
|
|---|
| 481 | {
|
|---|
| 482 | vd = a[n,n];
|
|---|
| 483 | }
|
|---|
| 484 | else
|
|---|
| 485 | {
|
|---|
| 486 | vd = 1;
|
|---|
| 487 | }
|
|---|
| 488 | x[n] = x[n]/vd;
|
|---|
| 489 | for(i=n-1; i>=1; i--)
|
|---|
| 490 | {
|
|---|
| 491 | ip1 = i+1;
|
|---|
| 492 | if( upper )
|
|---|
| 493 | {
|
|---|
| 494 | v = 0.0;
|
|---|
| 495 | for(i_=ip1; i_<=n;i_++)
|
|---|
| 496 | {
|
|---|
| 497 | v += a[i,i_]*x[i_];
|
|---|
| 498 | }
|
|---|
| 499 | }
|
|---|
| 500 | else
|
|---|
| 501 | {
|
|---|
| 502 | v = 0.0;
|
|---|
| 503 | for(i_=ip1; i_<=n;i_++)
|
|---|
| 504 | {
|
|---|
| 505 | v += a[i_,i]*x[i_];
|
|---|
| 506 | }
|
|---|
| 507 | }
|
|---|
| 508 | if( nounit )
|
|---|
| 509 | {
|
|---|
| 510 | vd = a[i,i];
|
|---|
| 511 | }
|
|---|
| 512 | else
|
|---|
| 513 | {
|
|---|
| 514 | vd = 1;
|
|---|
| 515 | }
|
|---|
| 516 | x[i] = (x[i]-v)/vd;
|
|---|
| 517 | }
|
|---|
| 518 | }
|
|---|
| 519 | else
|
|---|
| 520 | {
|
|---|
| 521 | if( nounit )
|
|---|
| 522 | {
|
|---|
| 523 | vd = a[1,1];
|
|---|
| 524 | }
|
|---|
| 525 | else
|
|---|
| 526 | {
|
|---|
| 527 | vd = 1;
|
|---|
| 528 | }
|
|---|
| 529 | x[1] = x[1]/vd;
|
|---|
| 530 | for(i=2; i<=n; i++)
|
|---|
| 531 | {
|
|---|
| 532 | im1 = i-1;
|
|---|
| 533 | if( upper )
|
|---|
| 534 | {
|
|---|
| 535 | v = 0.0;
|
|---|
| 536 | for(i_=1; i_<=im1;i_++)
|
|---|
| 537 | {
|
|---|
| 538 | v += a[i_,i]*x[i_];
|
|---|
| 539 | }
|
|---|
| 540 | }
|
|---|
| 541 | else
|
|---|
| 542 | {
|
|---|
| 543 | v = 0.0;
|
|---|
| 544 | for(i_=1; i_<=im1;i_++)
|
|---|
| 545 | {
|
|---|
| 546 | v += a[i,i_]*x[i_];
|
|---|
| 547 | }
|
|---|
| 548 | }
|
|---|
| 549 | if( nounit )
|
|---|
| 550 | {
|
|---|
| 551 | vd = a[i,i];
|
|---|
| 552 | }
|
|---|
| 553 | else
|
|---|
| 554 | {
|
|---|
| 555 | vd = 1;
|
|---|
| 556 | }
|
|---|
| 557 | x[i] = (x[i]-v)/vd;
|
|---|
| 558 | }
|
|---|
| 559 | }
|
|---|
| 560 | }
|
|---|
| 561 | else
|
|---|
| 562 | {
|
|---|
| 563 |
|
|---|
| 564 | //
|
|---|
| 565 | // Use a Level 1 BLAS solve, scaling intermediate results.
|
|---|
| 566 | //
|
|---|
| 567 | if( (double)(xmax)>(double)(bignum) )
|
|---|
| 568 | {
|
|---|
| 569 |
|
|---|
| 570 | //
|
|---|
| 571 | // Scale X so that its components are less than or equal to
|
|---|
| 572 | // BIGNUM in absolute value.
|
|---|
| 573 | //
|
|---|
| 574 | s = bignum/xmax;
|
|---|
| 575 | for(i_=1; i_<=n;i_++)
|
|---|
| 576 | {
|
|---|
| 577 | x[i_] = s*x[i_];
|
|---|
| 578 | }
|
|---|
| 579 | xmax = bignum;
|
|---|
| 580 | }
|
|---|
| 581 | if( notran )
|
|---|
| 582 | {
|
|---|
| 583 |
|
|---|
| 584 | //
|
|---|
| 585 | // Solve A * x = b
|
|---|
| 586 | //
|
|---|
| 587 | j = jfirst;
|
|---|
| 588 | while( jinc>0 & j<=jlast | jinc<0 & j>=jlast )
|
|---|
| 589 | {
|
|---|
| 590 |
|
|---|
| 591 | //
|
|---|
| 592 | // Compute x(j) = b(j) / A(j,j), scaling x if necessary.
|
|---|
| 593 | //
|
|---|
| 594 | xj = Math.Abs(x[j]);
|
|---|
| 595 | flg = 0;
|
|---|
| 596 | if( nounit )
|
|---|
| 597 | {
|
|---|
| 598 | tjjs = a[j,j]*tscal;
|
|---|
| 599 | }
|
|---|
| 600 | else
|
|---|
| 601 | {
|
|---|
| 602 | tjjs = tscal;
|
|---|
| 603 | if( (double)(tscal)==(double)(1) )
|
|---|
| 604 | {
|
|---|
| 605 | flg = 100;
|
|---|
| 606 | }
|
|---|
| 607 | }
|
|---|
| 608 | if( flg!=100 )
|
|---|
| 609 | {
|
|---|
| 610 | tjj = Math.Abs(tjjs);
|
|---|
| 611 | if( (double)(tjj)>(double)(smlnum) )
|
|---|
| 612 | {
|
|---|
| 613 |
|
|---|
| 614 | //
|
|---|
| 615 | // abs(A(j,j)) > SMLNUM:
|
|---|
| 616 | //
|
|---|
| 617 | if( (double)(tjj)<(double)(1) )
|
|---|
| 618 | {
|
|---|
| 619 | if( (double)(xj)>(double)(tjj*bignum) )
|
|---|
| 620 | {
|
|---|
| 621 |
|
|---|
| 622 | //
|
|---|
| 623 | // Scale x by 1/b(j).
|
|---|
| 624 | //
|
|---|
| 625 | rec = 1/xj;
|
|---|
| 626 | for(i_=1; i_<=n;i_++)
|
|---|
| 627 | {
|
|---|
| 628 | x[i_] = rec*x[i_];
|
|---|
| 629 | }
|
|---|
| 630 | s = s*rec;
|
|---|
| 631 | xmax = xmax*rec;
|
|---|
| 632 | }
|
|---|
| 633 | }
|
|---|
| 634 | x[j] = x[j]/tjjs;
|
|---|
| 635 | xj = Math.Abs(x[j]);
|
|---|
| 636 | }
|
|---|
| 637 | else
|
|---|
| 638 | {
|
|---|
| 639 | if( (double)(tjj)>(double)(0) )
|
|---|
| 640 | {
|
|---|
| 641 |
|
|---|
| 642 | //
|
|---|
| 643 | // 0 < abs(A(j,j)) <= SMLNUM:
|
|---|
| 644 | //
|
|---|
| 645 | if( (double)(xj)>(double)(tjj*bignum) )
|
|---|
| 646 | {
|
|---|
| 647 |
|
|---|
| 648 | //
|
|---|
| 649 | // Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM
|
|---|
| 650 | // to avoid overflow when dividing by A(j,j).
|
|---|
| 651 | //
|
|---|
| 652 | rec = tjj*bignum/xj;
|
|---|
| 653 | if( (double)(cnorm[j])>(double)(1) )
|
|---|
| 654 | {
|
|---|
| 655 |
|
|---|
| 656 | //
|
|---|
| 657 | // Scale by 1/CNORM(j) to avoid overflow when
|
|---|
| 658 | // multiplying x(j) times column j.
|
|---|
| 659 | //
|
|---|
| 660 | rec = rec/cnorm[j];
|
|---|
| 661 | }
|
|---|
| 662 | for(i_=1; i_<=n;i_++)
|
|---|
| 663 | {
|
|---|
| 664 | x[i_] = rec*x[i_];
|
|---|
| 665 | }
|
|---|
| 666 | s = s*rec;
|
|---|
| 667 | xmax = xmax*rec;
|
|---|
| 668 | }
|
|---|
| 669 | x[j] = x[j]/tjjs;
|
|---|
| 670 | xj = Math.Abs(x[j]);
|
|---|
| 671 | }
|
|---|
| 672 | else
|
|---|
| 673 | {
|
|---|
| 674 |
|
|---|
| 675 | //
|
|---|
| 676 | // A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
|
|---|
| 677 | // scale = 0, and compute a solution to A*x = 0.
|
|---|
| 678 | //
|
|---|
| 679 | for(i=1; i<=n; i++)
|
|---|
| 680 | {
|
|---|
| 681 | x[i] = 0;
|
|---|
| 682 | }
|
|---|
| 683 | x[j] = 1;
|
|---|
| 684 | xj = 1;
|
|---|
| 685 | s = 0;
|
|---|
| 686 | xmax = 0;
|
|---|
| 687 | }
|
|---|
| 688 | }
|
|---|
| 689 | }
|
|---|
| 690 |
|
|---|
| 691 | //
|
|---|
| 692 | // Scale x if necessary to avoid overflow when adding a
|
|---|
| 693 | // multiple of column j of A.
|
|---|
| 694 | //
|
|---|
| 695 | if( (double)(xj)>(double)(1) )
|
|---|
| 696 | {
|
|---|
| 697 | rec = 1/xj;
|
|---|
| 698 | if( (double)(cnorm[j])>(double)((bignum-xmax)*rec) )
|
|---|
| 699 | {
|
|---|
| 700 |
|
|---|
| 701 | //
|
|---|
| 702 | // Scale x by 1/(2*abs(x(j))).
|
|---|
| 703 | //
|
|---|
| 704 | rec = rec*0.5;
|
|---|
| 705 | for(i_=1; i_<=n;i_++)
|
|---|
| 706 | {
|
|---|
| 707 | x[i_] = rec*x[i_];
|
|---|
| 708 | }
|
|---|
| 709 | s = s*rec;
|
|---|
| 710 | }
|
|---|
| 711 | }
|
|---|
| 712 | else
|
|---|
| 713 | {
|
|---|
| 714 | if( (double)(xj*cnorm[j])>(double)(bignum-xmax) )
|
|---|
| 715 | {
|
|---|
| 716 |
|
|---|
| 717 | //
|
|---|
| 718 | // Scale x by 1/2.
|
|---|
| 719 | //
|
|---|
| 720 | for(i_=1; i_<=n;i_++)
|
|---|
| 721 | {
|
|---|
| 722 | x[i_] = 0.5*x[i_];
|
|---|
| 723 | }
|
|---|
| 724 | s = s*0.5;
|
|---|
| 725 | }
|
|---|
| 726 | }
|
|---|
| 727 | if( upper )
|
|---|
| 728 | {
|
|---|
| 729 | if( j>1 )
|
|---|
| 730 | {
|
|---|
| 731 |
|
|---|
| 732 | //
|
|---|
| 733 | // Compute the update
|
|---|
| 734 | // x(1:j-1) := x(1:j-1) - x(j) * A(1:j-1,j)
|
|---|
| 735 | //
|
|---|
| 736 | v = x[j]*tscal;
|
|---|
| 737 | jm1 = j-1;
|
|---|
| 738 | for(i_=1; i_<=jm1;i_++)
|
|---|
| 739 | {
|
|---|
| 740 | x[i_] = x[i_] - v*a[i_,j];
|
|---|
| 741 | }
|
|---|
| 742 | i = 1;
|
|---|
| 743 | for(k=2; k<=j-1; k++)
|
|---|
| 744 | {
|
|---|
| 745 | if( (double)(Math.Abs(x[k]))>(double)(Math.Abs(x[i])) )
|
|---|
| 746 | {
|
|---|
| 747 | i = k;
|
|---|
| 748 | }
|
|---|
| 749 | }
|
|---|
| 750 | xmax = Math.Abs(x[i]);
|
|---|
| 751 | }
|
|---|
| 752 | }
|
|---|
| 753 | else
|
|---|
| 754 | {
|
|---|
| 755 | if( j<n )
|
|---|
| 756 | {
|
|---|
| 757 |
|
|---|
| 758 | //
|
|---|
| 759 | // Compute the update
|
|---|
| 760 | // x(j+1:n) := x(j+1:n) - x(j) * A(j+1:n,j)
|
|---|
| 761 | //
|
|---|
| 762 | jp1 = j+1;
|
|---|
| 763 | v = x[j]*tscal;
|
|---|
| 764 | for(i_=jp1; i_<=n;i_++)
|
|---|
| 765 | {
|
|---|
| 766 | x[i_] = x[i_] - v*a[i_,j];
|
|---|
| 767 | }
|
|---|
| 768 | i = j+1;
|
|---|
| 769 | for(k=j+2; k<=n; k++)
|
|---|
| 770 | {
|
|---|
| 771 | if( (double)(Math.Abs(x[k]))>(double)(Math.Abs(x[i])) )
|
|---|
| 772 | {
|
|---|
| 773 | i = k;
|
|---|
| 774 | }
|
|---|
| 775 | }
|
|---|
| 776 | xmax = Math.Abs(x[i]);
|
|---|
| 777 | }
|
|---|
| 778 | }
|
|---|
| 779 | j = j+jinc;
|
|---|
| 780 | }
|
|---|
| 781 | }
|
|---|
| 782 | else
|
|---|
| 783 | {
|
|---|
| 784 |
|
|---|
| 785 | //
|
|---|
| 786 | // Solve A' * x = b
|
|---|
| 787 | //
|
|---|
| 788 | j = jfirst;
|
|---|
| 789 | while( jinc>0 & j<=jlast | jinc<0 & j>=jlast )
|
|---|
| 790 | {
|
|---|
| 791 |
|
|---|
| 792 | //
|
|---|
| 793 | // Compute x(j) = b(j) - sum A(k,j)*x(k).
|
|---|
| 794 | // k<>j
|
|---|
| 795 | //
|
|---|
| 796 | xj = Math.Abs(x[j]);
|
|---|
| 797 | uscal = tscal;
|
|---|
| 798 | rec = 1/Math.Max(xmax, 1);
|
|---|
| 799 | if( (double)(cnorm[j])>(double)((bignum-xj)*rec) )
|
|---|
| 800 | {
|
|---|
| 801 |
|
|---|
| 802 | //
|
|---|
| 803 | // If x(j) could overflow, scale x by 1/(2*XMAX).
|
|---|
| 804 | //
|
|---|
| 805 | rec = rec*0.5;
|
|---|
| 806 | if( nounit )
|
|---|
| 807 | {
|
|---|
| 808 | tjjs = a[j,j]*tscal;
|
|---|
| 809 | }
|
|---|
| 810 | else
|
|---|
| 811 | {
|
|---|
| 812 | tjjs = tscal;
|
|---|
| 813 | }
|
|---|
| 814 | tjj = Math.Abs(tjjs);
|
|---|
| 815 | if( (double)(tjj)>(double)(1) )
|
|---|
| 816 | {
|
|---|
| 817 |
|
|---|
| 818 | //
|
|---|
| 819 | // Divide by A(j,j) when scaling x if A(j,j) > 1.
|
|---|
| 820 | //
|
|---|
| 821 | rec = Math.Min(1, rec*tjj);
|
|---|
| 822 | uscal = uscal/tjjs;
|
|---|
| 823 | }
|
|---|
| 824 | if( (double)(rec)<(double)(1) )
|
|---|
| 825 | {
|
|---|
| 826 | for(i_=1; i_<=n;i_++)
|
|---|
| 827 | {
|
|---|
| 828 | x[i_] = rec*x[i_];
|
|---|
| 829 | }
|
|---|
| 830 | s = s*rec;
|
|---|
| 831 | xmax = xmax*rec;
|
|---|
| 832 | }
|
|---|
| 833 | }
|
|---|
| 834 | sumj = 0;
|
|---|
| 835 | if( (double)(uscal)==(double)(1) )
|
|---|
| 836 | {
|
|---|
| 837 |
|
|---|
| 838 | //
|
|---|
| 839 | // If the scaling needed for A in the dot product is 1,
|
|---|
| 840 | // call DDOT to perform the dot product.
|
|---|
| 841 | //
|
|---|
| 842 | if( upper )
|
|---|
| 843 | {
|
|---|
| 844 | if( j>1 )
|
|---|
| 845 | {
|
|---|
| 846 | jm1 = j-1;
|
|---|
| 847 | sumj = 0.0;
|
|---|
| 848 | for(i_=1; i_<=jm1;i_++)
|
|---|
| 849 | {
|
|---|
| 850 | sumj += a[i_,j]*x[i_];
|
|---|
| 851 | }
|
|---|
| 852 | }
|
|---|
| 853 | else
|
|---|
| 854 | {
|
|---|
| 855 | sumj = 0;
|
|---|
| 856 | }
|
|---|
| 857 | }
|
|---|
| 858 | else
|
|---|
| 859 | {
|
|---|
| 860 | if( j<n )
|
|---|
| 861 | {
|
|---|
| 862 | jp1 = j+1;
|
|---|
| 863 | sumj = 0.0;
|
|---|
| 864 | for(i_=jp1; i_<=n;i_++)
|
|---|
| 865 | {
|
|---|
| 866 | sumj += a[i_,j]*x[i_];
|
|---|
| 867 | }
|
|---|
| 868 | }
|
|---|
| 869 | }
|
|---|
| 870 | }
|
|---|
| 871 | else
|
|---|
| 872 | {
|
|---|
| 873 |
|
|---|
| 874 | //
|
|---|
| 875 | // Otherwise, use in-line code for the dot product.
|
|---|
| 876 | //
|
|---|
| 877 | if( upper )
|
|---|
| 878 | {
|
|---|
| 879 | for(i=1; i<=j-1; i++)
|
|---|
| 880 | {
|
|---|
| 881 | v = a[i,j]*uscal;
|
|---|
| 882 | sumj = sumj+v*x[i];
|
|---|
| 883 | }
|
|---|
| 884 | }
|
|---|
| 885 | else
|
|---|
| 886 | {
|
|---|
| 887 | if( j<n )
|
|---|
| 888 | {
|
|---|
| 889 | for(i=j+1; i<=n; i++)
|
|---|
| 890 | {
|
|---|
| 891 | v = a[i,j]*uscal;
|
|---|
| 892 | sumj = sumj+v*x[i];
|
|---|
| 893 | }
|
|---|
| 894 | }
|
|---|
| 895 | }
|
|---|
| 896 | }
|
|---|
| 897 | if( (double)(uscal)==(double)(tscal) )
|
|---|
| 898 | {
|
|---|
| 899 |
|
|---|
| 900 | //
|
|---|
| 901 | // Compute x(j) := ( x(j) - sumj ) / A(j,j) if 1/A(j,j)
|
|---|
| 902 | // was not used to scale the dotproduct.
|
|---|
| 903 | //
|
|---|
| 904 | x[j] = x[j]-sumj;
|
|---|
| 905 | xj = Math.Abs(x[j]);
|
|---|
| 906 | flg = 0;
|
|---|
| 907 | if( nounit )
|
|---|
| 908 | {
|
|---|
| 909 | tjjs = a[j,j]*tscal;
|
|---|
| 910 | }
|
|---|
| 911 | else
|
|---|
| 912 | {
|
|---|
| 913 | tjjs = tscal;
|
|---|
| 914 | if( (double)(tscal)==(double)(1) )
|
|---|
| 915 | {
|
|---|
| 916 | flg = 150;
|
|---|
| 917 | }
|
|---|
| 918 | }
|
|---|
| 919 |
|
|---|
| 920 | //
|
|---|
| 921 | // Compute x(j) = x(j) / A(j,j), scaling if necessary.
|
|---|
| 922 | //
|
|---|
| 923 | if( flg!=150 )
|
|---|
| 924 | {
|
|---|
| 925 | tjj = Math.Abs(tjjs);
|
|---|
| 926 | if( (double)(tjj)>(double)(smlnum) )
|
|---|
| 927 | {
|
|---|
| 928 |
|
|---|
| 929 | //
|
|---|
| 930 | // abs(A(j,j)) > SMLNUM:
|
|---|
| 931 | //
|
|---|
| 932 | if( (double)(tjj)<(double)(1) )
|
|---|
| 933 | {
|
|---|
| 934 | if( (double)(xj)>(double)(tjj*bignum) )
|
|---|
| 935 | {
|
|---|
| 936 |
|
|---|
| 937 | //
|
|---|
| 938 | // Scale X by 1/abs(x(j)).
|
|---|
| 939 | //
|
|---|
| 940 | rec = 1/xj;
|
|---|
| 941 | for(i_=1; i_<=n;i_++)
|
|---|
| 942 | {
|
|---|
| 943 | x[i_] = rec*x[i_];
|
|---|
| 944 | }
|
|---|
| 945 | s = s*rec;
|
|---|
| 946 | xmax = xmax*rec;
|
|---|
| 947 | }
|
|---|
| 948 | }
|
|---|
| 949 | x[j] = x[j]/tjjs;
|
|---|
| 950 | }
|
|---|
| 951 | else
|
|---|
| 952 | {
|
|---|
| 953 | if( (double)(tjj)>(double)(0) )
|
|---|
| 954 | {
|
|---|
| 955 |
|
|---|
| 956 | //
|
|---|
| 957 | // 0 < abs(A(j,j)) <= SMLNUM:
|
|---|
| 958 | //
|
|---|
| 959 | if( (double)(xj)>(double)(tjj*bignum) )
|
|---|
| 960 | {
|
|---|
| 961 |
|
|---|
| 962 | //
|
|---|
| 963 | // Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
|
|---|
| 964 | //
|
|---|
| 965 | rec = tjj*bignum/xj;
|
|---|
| 966 | for(i_=1; i_<=n;i_++)
|
|---|
| 967 | {
|
|---|
| 968 | x[i_] = rec*x[i_];
|
|---|
| 969 | }
|
|---|
| 970 | s = s*rec;
|
|---|
| 971 | xmax = xmax*rec;
|
|---|
| 972 | }
|
|---|
| 973 | x[j] = x[j]/tjjs;
|
|---|
| 974 | }
|
|---|
| 975 | else
|
|---|
| 976 | {
|
|---|
| 977 |
|
|---|
| 978 | //
|
|---|
| 979 | // A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
|
|---|
| 980 | // scale = 0, and compute a solution to A'*x = 0.
|
|---|
| 981 | //
|
|---|
| 982 | for(i=1; i<=n; i++)
|
|---|
| 983 | {
|
|---|
| 984 | x[i] = 0;
|
|---|
| 985 | }
|
|---|
| 986 | x[j] = 1;
|
|---|
| 987 | s = 0;
|
|---|
| 988 | xmax = 0;
|
|---|
| 989 | }
|
|---|
| 990 | }
|
|---|
| 991 | }
|
|---|
| 992 | }
|
|---|
| 993 | else
|
|---|
| 994 | {
|
|---|
| 995 |
|
|---|
| 996 | //
|
|---|
| 997 | // Compute x(j) := x(j) / A(j,j) - sumj if the dot
|
|---|
| 998 | // product has already been divided by 1/A(j,j).
|
|---|
| 999 | //
|
|---|
| 1000 | x[j] = x[j]/tjjs-sumj;
|
|---|
| 1001 | }
|
|---|
| 1002 | xmax = Math.Max(xmax, Math.Abs(x[j]));
|
|---|
| 1003 | j = j+jinc;
|
|---|
| 1004 | }
|
|---|
| 1005 | }
|
|---|
| 1006 | s = s/tscal;
|
|---|
| 1007 | }
|
|---|
| 1008 |
|
|---|
| 1009 | //
|
|---|
| 1010 | // Scale the column norms by 1/TSCAL for return.
|
|---|
| 1011 | //
|
|---|
| 1012 | if( (double)(tscal)!=(double)(1) )
|
|---|
| 1013 | {
|
|---|
| 1014 | v = 1/tscal;
|
|---|
| 1015 | for(i_=1; i_<=n;i_++)
|
|---|
| 1016 | {
|
|---|
| 1017 | cnorm[i_] = v*cnorm[i_];
|
|---|
| 1018 | }
|
|---|
| 1019 | }
|
|---|
| 1020 | }
|
|---|
| 1021 | }
|
|---|
| 1022 | }
|
|---|
