`   1:  #region Translated by Jose Antonio De Santiago-Castillo.`
`   2:   `
`   3:  //Translated by Jose Antonio De Santiago-Castillo. `
`   4:  //E-mail:JAntonioDeSantiago@gmail.com`
`   5:  //Web: www.DotNumerics.com`
`   6:  //`
`   7:  //Fortran to C# Translation.`
`   8:  //Translated by:`
`   9:  //F2CSharp Version 0.71 (November 10, 2009)`
`  10:  //Code Optimizations: None`
`  11:  //`
`  12:  #endregion`
`  13:   `
`  14:  using System;`
`  15:  using DotNumerics.FortranLibrary;`
`  16:   `
`  17:  namespace DotNumerics.CSLapack`
`  18:  {`
`  19:      /// <summary>`
`  20:      /// -- LAPACK auxiliary routine (version 3.1) --`
`  21:      /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..`
`  22:      /// November 2006`
`  23:      /// Purpose`
`  24:      /// =======`
`  25:      /// `
`  26:      /// DLALN2 solves a system of the form  (ca A - w D ) X = s B`
`  27:      /// or (ca A' - w D) X = s B   with possible scaling ("s") and`
`  28:      /// perturbation of A.  (A' means A-transpose.)`
`  29:      /// `
`  30:      /// A is an NA x NA real matrix, ca is a real scalar, D is an NA x NA`
`  31:      /// real diagonal matrix, w is a real or complex value, and X and B are`
`  32:      /// NA x 1 matrices -- real if w is real, complex if w is complex.  NA`
`  33:      /// may be 1 or 2.`
`  34:      /// `
`  35:      /// If w is complex, X and B are represented as NA x 2 matrices,`
`  36:      /// the first column of each being the real part and the second`
`  37:      /// being the imaginary part.`
`  38:      /// `
`  39:      /// "s" is a scaling factor (.LE. 1), computed by DLALN2, which is`
`  40:      /// so chosen that X can be computed without overflow.  X is further`
`  41:      /// scaled if necessary to assure that norm(ca A - w D)*norm(X) is less`
`  42:      /// than overflow.`
`  43:      /// `
`  44:      /// If both singular values of (ca A - w D) are less than SMIN,`
`  45:      /// SMIN*identity will be used instead of (ca A - w D).  If only one`
`  46:      /// singular value is less than SMIN, one element of (ca A - w D) will be`
`  47:      /// perturbed enough to make the smallest singular value roughly SMIN.`
`  48:      /// If both singular values are at least SMIN, (ca A - w D) will not be`
`  49:      /// perturbed.  In any case, the perturbation will be at most some small`
`  50:      /// multiple of max( SMIN, ulp*norm(ca A - w D) ).  The singular values`
`  51:      /// are computed by infinity-norm approximations, and thus will only be`
`  52:      /// correct to a factor of 2 or so.`
`  53:      /// `
`  54:      /// Note: all input quantities are assumed to be smaller than overflow`
`  55:      /// by a reasonable factor.  (See BIGNUM.)`
`  56:      /// `
`  57:      ///</summary>`
`  58:      public class DLALN2`
`  59:      {`
`  60:      `
`  61:   `
`  62:          #region Dependencies`
`  63:          `
`  64:          DLAMCH _dlamch; DLADIV _dladiv; `
`  65:   `
`  66:          #endregion`
`  67:   `
`  68:   `
`  69:          #region Fields`
`  70:          `
`  71:          const double ZERO = 0.0E0; const double ONE = 1.0E0; const double TWO = 2.0E0; int ICMAX = 0; int J = 0; double BBND = 0; `
`  72:          double BI1 = 0;double BI2 = 0; double BIGNUM = 0; double BNORM = 0; double BR1 = 0; double BR2 = 0; double CI21 = 0; `
`  73:          double CI22 = 0;double CMAX = 0; double CNORM = 0; double CR21 = 0; double CR22 = 0; double CSI = 0; double CSR = 0; `
`  74:          double LI21 = 0;double LR21 = 0; double SMINI = 0; double SMLNUM = 0; double TEMP = 0; double U22ABS = 0; double UI11 = 0; `
`  75:          double UI11R = 0;double UI12 = 0; double UI12S = 0; double UI22 = 0; double UR11 = 0; double UR11R = 0; double UR12 = 0; `
`  76:          double UR12S = 0;double UR22 = 0; double XI1 = 0; double XI2 = 0; double XR1 = 0; double XR2 = 0; `
`  77:          bool[] RSWAP = new bool[4]; int o_rswap = -1;bool[] ZSWAP = new bool[4];  int o_zswap = -1; `
`  78:          int[] IPIVOT = new int[4 * 4];  int o_ipivot = -5;double[] CIV = new double[4]; int o_civ = -1; `
`  79:          double[] CRV = new double[4]; int o_crv = -1;`
`  80:   `
`  81:          #endregion`
`  82:   `
`  83:          public DLALN2(DLAMCH dlamch, DLADIV dladiv)`
`  84:          {`
`  85:      `
`  86:   `
`  87:              #region Set Dependencies`
`  88:              `
`  89:              this._dlamch = dlamch; this._dladiv = dladiv; `
`  90:   `
`  91:              #endregion`
`  92:   `
`  93:   `
`  94:              #region Data Inicializacion`
`  95:              `
`  96:              //ZSWAP/.FALSE.,.FALSE.,.TRUE.,.TRUE.`
`  97:              ZSWAP[1 + o_zswap] = false;`
`  98:              ZSWAP[2 + o_zswap] = false;`
`  99:              ZSWAP[3 + o_zswap] = true;`
` 100:              ZSWAP[4 + o_zswap] = true;`
` 101:              //RSWAP/.FALSE.,.TRUE.,.FALSE.,.TRUE.`
` 102:              RSWAP[1 + o_rswap] = false;`
` 103:              RSWAP[2 + o_rswap] = true;`
` 104:              RSWAP[3 + o_rswap] = false;`
` 105:              RSWAP[4 + o_rswap] = true;`
` 106:              //IPIVOT/1,2,3,4,2,1,4,3,3,4,1,2,4,3,2,1`
` 107:              IPIVOT[0] = 1;`
` 108:              IPIVOT[1] = 2;`
` 109:              IPIVOT[2] = 3;`
` 110:              IPIVOT[3] = 4;`
` 111:              IPIVOT[4] = 2;`
` 112:              IPIVOT[5] = 1;`
` 113:              IPIVOT[6] = 4;`
` 114:              IPIVOT[7] = 3;`
` 115:              IPIVOT[8] = 3;`
` 116:              IPIVOT[9] = 4;`
` 117:              IPIVOT[10] = 1;`
` 118:              IPIVOT[11] = 2;`
` 119:              IPIVOT[12] = 4;`
` 120:              IPIVOT[12] = 3;`
` 121:              IPIVOT[14] = 2;`
` 122:              IPIVOT[15] = 1;`
` 123:   `
` 124:              #endregion`
` 125:   `
` 126:          }`
` 127:      `
` 128:          public DLALN2()`
` 129:          {`
` 130:      `
` 131:   `
` 132:              #region Dependencies (Initialization)`
` 133:              `
` 134:              LSAME lsame = new LSAME();`
` 135:              DLAMC3 dlamc3 = new DLAMC3();`
` 136:              DLADIV dladiv = new DLADIV();`
` 137:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);`
` 138:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);`
` 139:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);`
` 140:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);`
` 141:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);`
` 142:   `
` 143:              #endregion`
` 144:   `
` 145:   `
` 146:              #region Set Dependencies`
` 147:              `
` 148:              this._dlamch = dlamch; this._dladiv = dladiv; `
` 149:   `
` 150:              #endregion`
` 151:   `
` 152:   `
` 153:              #region Data Inicializacion`
` 154:              `
` 155:              //ZSWAP/.FALSE.,.FALSE.,.TRUE.,.TRUE.`
` 156:              ZSWAP[1 + o_zswap] = false;`
` 157:              ZSWAP[2 + o_zswap] = false;`
` 158:              ZSWAP[3 + o_zswap] = true;`
` 159:              ZSWAP[4 + o_zswap] = true;`
` 160:              //RSWAP/.FALSE.,.TRUE.,.FALSE.,.TRUE.`
` 161:              RSWAP[1 + o_rswap] = false;`
` 162:              RSWAP[2 + o_rswap] = true;`
` 163:              RSWAP[3 + o_rswap] = false;`
` 164:              RSWAP[4 + o_rswap] = true;`
` 165:              //IPIVOT/1,2,3,4,2,1,4,3,3,4,1,2,4,3,2,1`
` 166:              IPIVOT[0] = 1;`
` 167:              IPIVOT[1] = 2;`
` 168:              IPIVOT[2] = 3;`
` 169:              IPIVOT[3] = 4;`
` 170:              IPIVOT[4] = 2;`
` 171:              IPIVOT[5] = 1;`
` 172:              IPIVOT[6] = 4;`
` 173:              IPIVOT[7] = 3;`
` 174:              IPIVOT[8] = 3;`
` 175:              IPIVOT[9] = 4;`
` 176:              IPIVOT[10] = 1;`
` 177:              IPIVOT[11] = 2;`
` 178:              IPIVOT[12] = 4;`
` 179:              IPIVOT[12] = 3;`
` 180:              IPIVOT[14] = 2;`
` 181:              IPIVOT[15] = 1;`
` 182:   `
` 183:              #endregion`
` 184:   `
` 185:          }`
` 186:          /// <summary>`
` 187:          /// Purpose`
` 188:          /// =======`
` 189:          /// `
` 190:          /// DLALN2 solves a system of the form  (ca A - w D ) X = s B`
` 191:          /// or (ca A' - w D) X = s B   with possible scaling ("s") and`
` 192:          /// perturbation of A.  (A' means A-transpose.)`
` 193:          /// `
` 194:          /// A is an NA x NA real matrix, ca is a real scalar, D is an NA x NA`
` 195:          /// real diagonal matrix, w is a real or complex value, and X and B are`
` 196:          /// NA x 1 matrices -- real if w is real, complex if w is complex.  NA`
` 197:          /// may be 1 or 2.`
` 198:          /// `
` 199:          /// If w is complex, X and B are represented as NA x 2 matrices,`
` 200:          /// the first column of each being the real part and the second`
` 201:          /// being the imaginary part.`
` 202:          /// `
` 203:          /// "s" is a scaling factor (.LE. 1), computed by DLALN2, which is`
` 204:          /// so chosen that X can be computed without overflow.  X is further`
` 205:          /// scaled if necessary to assure that norm(ca A - w D)*norm(X) is less`
` 206:          /// than overflow.`
` 207:          /// `
` 208:          /// If both singular values of (ca A - w D) are less than SMIN,`
` 209:          /// SMIN*identity will be used instead of (ca A - w D).  If only one`
` 210:          /// singular value is less than SMIN, one element of (ca A - w D) will be`
` 211:          /// perturbed enough to make the smallest singular value roughly SMIN.`
` 212:          /// If both singular values are at least SMIN, (ca A - w D) will not be`
` 213:          /// perturbed.  In any case, the perturbation will be at most some small`
` 214:          /// multiple of max( SMIN, ulp*norm(ca A - w D) ).  The singular values`
` 215:          /// are computed by infinity-norm approximations, and thus will only be`
` 216:          /// correct to a factor of 2 or so.`
` 217:          /// `
` 218:          /// Note: all input quantities are assumed to be smaller than overflow`
` 219:          /// by a reasonable factor.  (See BIGNUM.)`
` 220:          /// `
` 221:          ///</summary>`
` 222:          /// <param name="LTRANS">`
` 223:          /// (input) LOGICAL`
` 224:          /// =.TRUE.:  A-transpose will be used.`
` 225:          /// =.FALSE.: A will be used (not transposed.)`
` 226:          ///</param>`
` 227:          /// <param name="NA">`
` 228:          /// (input) INTEGER`
` 229:          /// The size of the matrix A.  It may (only) be 1 or 2.`
` 230:          ///</param>`
` 231:          /// <param name="NW">`
` 232:          /// (input) INTEGER`
` 233:          /// 1 if "w" is real, 2 if "w" is complex.  It may only be 1`
` 234:          /// or 2.`
` 235:          ///</param>`
` 236:          /// <param name="SMIN">`
` 237:          /// (input) DOUBLE PRECISION`
` 238:          /// The desired lower bound on the singular values of A.  This`
` 239:          /// should be a safe distance away from underflow or overflow,`
` 240:          /// say, between (underflow/machine precision) and  (machine`
` 241:          /// precision * overflow ).  (See BIGNUM and ULP.)`
` 242:          ///</param>`
` 243:          /// <param name="CA">`
` 244:          /// (input) DOUBLE PRECISION`
` 245:          /// The coefficient c, which A is multiplied by.`
` 246:          ///</param>`
` 247:          /// <param name="A">`
` 248:          /// is an NA x NA real matrix, ca is a real scalar, D is an NA x NA`
` 249:          ///</param>`
` 250:          /// <param name="LDA">`
` 251:          /// (input) INTEGER`
` 252:          /// The leading dimension of A.  It must be at least NA.`
` 253:          ///</param>`
` 254:          /// <param name="D1">`
` 255:          /// (input) DOUBLE PRECISION`
` 256:          /// The 1,1 element in the diagonal matrix D.`
` 257:          ///</param>`
` 258:          /// <param name="D2">`
` 259:          /// (input) DOUBLE PRECISION`
` 260:          /// The 2,2 element in the diagonal matrix D.  Not used if NW=1.`
` 261:          ///</param>`
` 262:          /// <param name="B">`
` 263:          /// (input) DOUBLE PRECISION array, dimension (LDB,NW)`
` 264:          /// The NA x NW matrix B (right-hand side).  If NW=2 ("w" is`
` 265:          /// complex), column 1 contains the real part of B and column 2`
` 266:          /// contains the imaginary part.`
` 267:          ///</param>`
` 268:          /// <param name="LDB">`
` 269:          /// (input) INTEGER`
` 270:          /// The leading dimension of B.  It must be at least NA.`
` 271:          ///</param>`
` 272:          /// <param name="WR">`
` 273:          /// (input) DOUBLE PRECISION`
` 274:          /// The real part of the scalar "w".`
` 275:          ///</param>`
` 276:          /// <param name="WI">`
` 277:          /// (input) DOUBLE PRECISION`
` 278:          /// The imaginary part of the scalar "w".  Not used if NW=1.`
` 279:          ///</param>`
` 280:          /// <param name="X">`
` 281:          /// (output) DOUBLE PRECISION array, dimension (LDX,NW)`
` 282:          /// The NA x NW matrix X (unknowns), as computed by DLALN2.`
` 283:          /// If NW=2 ("w" is complex), on exit, column 1 will contain`
` 284:          /// the real part of X and column 2 will contain the imaginary`
` 285:          /// part.`
` 286:          ///</param>`
` 287:          /// <param name="LDX">`
` 288:          /// (input) INTEGER`
` 289:          /// The leading dimension of X.  It must be at least NA.`
` 290:          ///</param>`
` 291:          /// <param name="SCALE">`
` 292:          /// (output) DOUBLE PRECISION`
` 293:          /// The scale factor that B must be multiplied by to insure`
` 294:          /// that overflow does not occur when computing X.  Thus,`
` 295:          /// (ca A - w D) X  will be SCALE*B, not B (ignoring`
` 296:          /// perturbations of A.)  It will be at most 1.`
` 297:          ///</param>`
` 298:          /// <param name="XNORM">`
` 299:          /// (output) DOUBLE PRECISION`
` 300:          /// The infinity-norm of X, when X is regarded as an NA x NW`
` 301:          /// real matrix.`
` 302:          ///</param>`
` 303:          /// <param name="INFO">`
` 304:          /// (output) INTEGER`
` 305:          /// An error flag.  It will be set to zero if no error occurs,`
` 306:          /// a negative number if an argument is in error, or a positive`
` 307:          /// number if  ca A - w D  had to be perturbed.`
` 308:          /// The possible values are:`
` 309:          /// = 0: No error occurred, and (ca A - w D) did not have to be`
` 310:          /// perturbed.`
` 311:          /// = 1: (ca A - w D) had to be perturbed to make its smallest`
` 312:          /// (or only) singular value greater than SMIN.`
` 313:          /// NOTE: In the interests of speed, this routine does not`
` 314:          /// check the inputs for errors.`
` 315:          ///</param>`
` 316:          public void Run(bool LTRANS, int NA, int NW, double SMIN, double CA, double[] A, int offset_a`
` 317:                           , int LDA, double D1, double D2, double[] B, int offset_b, int LDB, double WR`
` 318:                           , double WI, ref double[] X, int offset_x, int LDX, ref double SCALE, ref double XNORM, ref int INFO)`
` 319:          {`
` 320:   `
` 321:              #region Array Index Correction`
` 322:              `
` 323:               int o_a = -1 - LDA + offset_a;  int o_b = -1 - LDB + offset_b;  int o_x = -1 - LDX + offset_x; `
` 324:   `
` 325:              #endregion`
` 326:   `
` 327:   `
` 328:              #region Prolog`
` 329:              `
` 330:              `
` 331:              `
` 332:              // *`
` 333:              // *    SUBRUTINA MODIFICADA, SE ELIMINO LA EQUIVALENCIA `
` 334:              // *`
` 335:              // *      EQUIVALENCE        ( CI( 1, 1 ), CIV( 1 ) ),`
` 336:              // *     \$                   ( CR( 1, 1 ), CRV( 1 ) )`
` 337:              // *`
` 338:              // *    DEBIDO A QUE EL TRADUCTOR NO LA SOPORTA`
` 339:              `
` 340:              // *`
` 341:              // *  -- LAPACK auxiliary routine (version 3.1) --`
` 342:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..`
` 343:              // *     November 2006`
` 344:              // *`
` 345:              // *     .. Scalar Arguments ..`
` 346:              // *     ..`
` 347:              // *     .. Array Arguments ..`
` 348:              // *     ..`
` 349:              // *`
` 350:              // *  Purpose`
` 351:              // *  =======`
` 352:              // *`
` 353:              // *  DLALN2 solves a system of the form  (ca A - w D ) X = s B`
` 354:              // *  or (ca A' - w D) X = s B   with possible scaling ("s") and`
` 355:              // *  perturbation of A.  (A' means A-transpose.)`
` 356:              // *`
` 357:              // *  A is an NA x NA real matrix, ca is a real scalar, D is an NA x NA`
` 358:              // *  real diagonal matrix, w is a real or complex value, and X and B are`
` 359:              // *  NA x 1 matrices -- real if w is real, complex if w is complex.  NA`
` 360:              // *  may be 1 or 2.`
` 361:              // *`
` 362:              // *  If w is complex, X and B are represented as NA x 2 matrices,`
` 363:              // *  the first column of each being the real part and the second`
` 364:              // *  being the imaginary part.`
` 365:              // *`
` 366:              // *  "s" is a scaling factor (.LE. 1), computed by DLALN2, which is`
` 367:              // *  so chosen that X can be computed without overflow.  X is further`
` 368:              // *  scaled if necessary to assure that norm(ca A - w D)*norm(X) is less`
` 369:              // *  than overflow.`
` 370:              // *`
` 371:              // *  If both singular values of (ca A - w D) are less than SMIN,`
` 372:              // *  SMIN*identity will be used instead of (ca A - w D).  If only one`
` 373:              // *  singular value is less than SMIN, one element of (ca A - w D) will be`
` 374:              // *  perturbed enough to make the smallest singular value roughly SMIN.`
` 375:              // *  If both singular values are at least SMIN, (ca A - w D) will not be`
` 376:              // *  perturbed.  In any case, the perturbation will be at most some small`
` 377:              // *  multiple of max( SMIN, ulp*norm(ca A - w D) ).  The singular values`
` 378:              // *  are computed by infinity-norm approximations, and thus will only be`
` 379:              // *  correct to a factor of 2 or so.`
` 380:              // *`
` 381:              // *  Note: all input quantities are assumed to be smaller than overflow`
` 382:              // *  by a reasonable factor.  (See BIGNUM.)`
` 383:              // *`
` 384:              // *  Arguments`
` 385:              // *  ==========`
` 386:              // *`
` 387:              // *  LTRANS  (input) LOGICAL`
` 388:              // *          =.TRUE.:  A-transpose will be used.`
` 389:              // *          =.FALSE.: A will be used (not transposed.)`
` 390:              // *`
` 391:              // *  NA      (input) INTEGER`
` 392:              // *          The size of the matrix A.  It may (only) be 1 or 2.`
` 393:              // *`
` 394:              // *  NW      (input) INTEGER`
` 395:              // *          1 if "w" is real, 2 if "w" is complex.  It may only be 1`
` 396:              // *          or 2.`
` 397:              // *`
` 398:              // *  SMIN    (input) DOUBLE PRECISION`
` 399:              // *          The desired lower bound on the singular values of A.  This`
` 400:              // *          should be a safe distance away from underflow or overflow,`
` 401:              // *          say, between (underflow/machine precision) and  (machine`
` 402:              // *          precision * overflow ).  (See BIGNUM and ULP.)`
` 403:              // *`
` 404:              // *  CA      (input) DOUBLE PRECISION`
` 405:              // *          The coefficient c, which A is multiplied by.`
` 406:              // *`
` 407:              // *  A       (input) DOUBLE PRECISION array, dimension (LDA,NA)`
` 408:              // *          The NA x NA matrix A.`
` 409:              // *`
` 410:              // *  LDA     (input) INTEGER`
` 411:              // *          The leading dimension of A.  It must be at least NA.`
` 412:              // *`
` 413:              // *  D1      (input) DOUBLE PRECISION`
` 414:              // *          The 1,1 element in the diagonal matrix D.`
` 415:              // *`
` 416:              // *  D2      (input) DOUBLE PRECISION`
` 417:              // *          The 2,2 element in the diagonal matrix D.  Not used if NW=1.`
` 418:              // *`
` 419:              // *  B       (input) DOUBLE PRECISION array, dimension (LDB,NW)`
` 420:              // *          The NA x NW matrix B (right-hand side).  If NW=2 ("w" is`
` 421:              // *          complex), column 1 contains the real part of B and column 2`
` 422:              // *          contains the imaginary part.`
` 423:              // *`
` 424:              // *  LDB     (input) INTEGER`
` 425:              // *          The leading dimension of B.  It must be at least NA.`
` 426:              // *`
` 427:              // *  WR      (input) DOUBLE PRECISION`
` 428:              // *          The real part of the scalar "w".`
` 429:              // *`
` 430:              // *  WI      (input) DOUBLE PRECISION`
` 431:              // *          The imaginary part of the scalar "w".  Not used if NW=1.`
` 432:              // *`
` 433:              // *  X       (output) DOUBLE PRECISION array, dimension (LDX,NW)`
` 434:              // *          The NA x NW matrix X (unknowns), as computed by DLALN2.`
` 435:              // *          If NW=2 ("w" is complex), on exit, column 1 will contain`
` 436:              // *          the real part of X and column 2 will contain the imaginary`
` 437:              // *          part.`
` 438:              // *`
` 439:              // *  LDX     (input) INTEGER`
` 440:              // *          The leading dimension of X.  It must be at least NA.`
` 441:              // *`
` 442:              // *  SCALE   (output) DOUBLE PRECISION`
` 443:              // *          The scale factor that B must be multiplied by to insure`
` 444:              // *          that overflow does not occur when computing X.  Thus,`
` 445:              // *          (ca A - w D) X  will be SCALE*B, not B (ignoring`
` 446:              // *          perturbations of A.)  It will be at most 1.`
` 447:              // *`
` 448:              // *  XNORM   (output) DOUBLE PRECISION`
` 449:              // *          The infinity-norm of X, when X is regarded as an NA x NW`
` 450:              // *          real matrix.`
` 451:              // *`
` 452:              // *  INFO    (output) INTEGER`
` 453:              // *          An error flag.  It will be set to zero if no error occurs,`
` 454:              // *          a negative number if an argument is in error, or a positive`
` 455:              // *          number if  ca A - w D  had to be perturbed.`
` 456:              // *          The possible values are:`
` 457:              // *          = 0: No error occurred, and (ca A - w D) did not have to be`
` 458:              // *                 perturbed.`
` 459:              // *          = 1: (ca A - w D) had to be perturbed to make its smallest`
` 460:              // *               (or only) singular value greater than SMIN.`
` 461:              // *          NOTE: In the interests of speed, this routine does not`
` 462:              // *                check the inputs for errors.`
` 463:              // *`
` 464:              // * =====================================================================`
` 465:              // *`
` 466:              // *     .. Parameters ..`
` 467:              // *     ..`
` 468:              // *     .. Local Scalars ..`
` 469:              // *     ..`
` 470:              // *     .. Local Arrays ..`
` 471:              `
` 472:              `
` 473:              `
` 474:              // *-----------------------------INICIA MI MODIFICACION--------------------------`
` 475:              `
` 476:              // *      DOUBLE PRECISION   CI( 2, 2 ), CIV( 4 ), CR( 2, 2 ), CRV( 4 )  `
` 477:              // *-----------------------------TERMINA MI MODIFICACION-------------------------`
` 478:              `
` 479:              `
` 480:              `
` 481:              // *     ..`
` 482:              // *     .. External Functions ..`
` 483:              // *     ..`
` 484:              // *     .. External Subroutines ..`
` 485:              // *     ..`
` 486:              // *     .. Intrinsic Functions ..`
` 487:              //      INTRINSIC          ABS, MAX;`
` 488:              // *     ..`
` 489:              // *     .. Equivalences ..`
` 490:              `
` 491:              `
` 492:              // *-----------------------------INICIA MI MODIFICACION--------------------------`
` 493:              // *      EQUIVALENCE        ( CI( 1, 1 ), CIV( 1 ) ),       `
` 494:              // *     \$                   ( CR( 1, 1 ), CRV( 1 ) )`
` 495:              // *-----------------------------TERMINA MI MODIFICACION-------------------------`
` 496:              `
` 497:              `
` 498:              `
` 499:              // *     .. Data statements ..`
` 500:              // *     ..`
` 501:              // *     .. Executable Statements ..`
` 502:              // *`
` 503:              // *     Compute BIGNUM`
` 504:              // *`
` 505:   `
` 506:              #endregion`
` 507:   `
` 508:   `
` 509:              #region Body`
` 510:              `
` 511:              SMLNUM = TWO * this._dlamch.Run("Safe minimum");`
` 512:              BIGNUM = ONE / SMLNUM;`
` 513:              SMINI = Math.Max(SMIN, SMLNUM);`
` 514:              // *`
` 515:              // *     Don't check for input errors`
` 516:              // *`
` 517:              INFO = 0;`
` 518:              // *`
` 519:              // *     Standard Initializations`
` 520:              // *`
` 521:              SCALE = ONE;`
` 522:              // *`
` 523:              if (NA == 1)`
` 524:              {`
` 525:                  // *`
` 526:                  // *        1 x 1  (i.e., scalar) system   C X = B`
` 527:                  // *`
` 528:                  if (NW == 1)`
` 529:                  {`
` 530:                      // *`
` 531:                      // *           Real 1x1 system.`
` 532:                      // *`
` 533:                      // *           C = ca A - w D`
` 534:                      // *`
` 535:                      CSR = CA * A[1+1 * LDA + o_a] - WR * D1;`
` 536:                      CNORM = Math.Abs(CSR);`
` 537:                      // *`
` 538:                      // *           If | C | < SMINI, use C = SMINI`
` 539:                      // *`
` 540:                      if (CNORM < SMINI)`
` 541:                      {`
` 542:                          CSR = SMINI;`
` 543:                          CNORM = SMINI;`
` 544:                          INFO = 1;`
` 545:                      }`
` 546:                      // *`
` 547:                      // *           Check scaling for  X = B / C`
` 548:                      // *`
` 549:                      BNORM = Math.Abs(B[1+1 * LDB + o_b]);`
` 550:                      if (CNORM < ONE && BNORM > ONE)`
` 551:                      {`
` 552:                          if (BNORM > BIGNUM * CNORM) SCALE = ONE / BNORM;`
` 553:                      }`
` 554:                      // *`
` 555:                      // *           Compute X`
` 556:                      // *`
` 557:                      X[1+1 * LDX + o_x] = (B[1+1 * LDB + o_b] * SCALE) / CSR;`
` 558:                      XNORM = Math.Abs(X[1+1 * LDX + o_x]);`
` 559:                  }`
` 560:                  else`
` 561:                  {`
` 562:                      // *`
` 563:                      // *           Complex 1x1 system (w is complex)`
` 564:                      // *`
` 565:                      // *           C = ca A - w D`
` 566:                      // *`
` 567:                      CSR = CA * A[1+1 * LDA + o_a] - WR * D1;`
` 568:                      CSI =  - WI * D1;`
` 569:                      CNORM = Math.Abs(CSR) + Math.Abs(CSI);`
` 570:                      // *`
` 571:                      // *           If | C | < SMINI, use C = SMINI`
` 572:                      // *`
` 573:                      if (CNORM < SMINI)`
` 574:                      {`
` 575:                          CSR = SMINI;`
` 576:                          CSI = ZERO;`
` 577:                          CNORM = SMINI;`
` 578:                          INFO = 1;`
` 579:                      }`
` 580:                      // *`
` 581:                      // *           Check scaling for  X = B / C`
` 582:                      // *`
` 583:                      BNORM = Math.Abs(B[1+1 * LDB + o_b]) + Math.Abs(B[1+2 * LDB + o_b]);`
` 584:                      if (CNORM < ONE && BNORM > ONE)`
` 585:                      {`
` 586:                          if (BNORM > BIGNUM * CNORM) SCALE = ONE / BNORM;`
` 587:                      }`
` 588:                      // *`
` 589:                      // *           Compute X`
` 590:                      // *`
` 591:                      this._dladiv.Run(SCALE * B[1+1 * LDB + o_b], SCALE * B[1+2 * LDB + o_b], CSR, CSI, ref X[1+1 * LDX + o_x], ref X[1+2 * LDX + o_x]);`
` 592:                      XNORM = Math.Abs(X[1+1 * LDX + o_x]) + Math.Abs(X[1+2 * LDX + o_x]);`
` 593:                  }`
` 594:                  // *`
` 595:              }`
` 596:              else`
` 597:              {`
` 598:                  // *`
` 599:                  // *        2x2 System`
` 600:                  // *`
` 601:                  // *        Compute the real part of  C = ca A - w D  (or  ca A' - w D )`
` 602:                  // *`
` 603:                  `
` 604:                  `
` 605:                  `
` 606:                  `
` 607:                  // *-----------------------------INICIA MI MODIFICACION--------------------------`
` 608:                  `
` 609:                  // *         CR( 1, 1 ) = CA*A( 1, 1 ) - WR*D1`
` 610:                  // *         CR( 2, 2 ) = CA*A( 2, 2 ) - WR*D2`
` 611:                  `
` 612:                  CRV[1 + o_crv] = CA * A[1+1 * LDA + o_a] - WR * D1;`
` 613:                  CRV[4 + o_crv] = CA * A[2+2 * LDA + o_a] - WR * D2;`
` 614:                  `
` 615:                  if (LTRANS)`
` 616:                  {`
` 617:                      // *            CR( 1, 2 ) = CA*A( 2, 1 )`
` 618:                      // *            CR( 2, 1 ) = CA*A( 1, 2 )`
` 619:                      `
` 620:                      CRV[3 + o_crv] = CA * A[2+1 * LDA + o_a];`
` 621:                      CRV[2 + o_crv] = CA * A[1+2 * LDA + o_a];`
` 622:                      `
` 623:                  }`
` 624:                  else`
` 625:                  {`
` 626:                      // *            CR( 2, 1 ) = CA*A( 2, 1 )`
` 627:                      // *            CR( 1, 2 ) = CA*A( 1, 2 )`
` 628:                      `
` 629:                      `
` 630:                      CRV[2 + o_crv] = CA * A[2+1 * LDA + o_a];`
` 631:                      CRV[3 + o_crv] = CA * A[1+2 * LDA + o_a];`
` 632:                      `
` 633:                  }`
` 634:                  `
` 635:                  // *-----------------------------TERMINA MI MODIFICACION-------------------------`
` 636:                  `
` 637:                  `
` 638:                  `
` 639:                  // *`
` 640:                  if (NW == 1)`
` 641:                  {`
` 642:                      // *`
` 643:                      // *           Real 2x2 system  (w is real)`
` 644:                      // *`
` 645:                      // *           Find the largest element in C`
` 646:                      // *`
` 647:                      CMAX = ZERO;`
` 648:                      ICMAX = 0;`
` 649:                      // *`
` 650:                      for (J = 1; J <= 4; J++)`
` 651:                      {`
` 652:                          if (Math.Abs(CRV[J + o_crv]) > CMAX)`
` 653:                          {`
` 654:                              CMAX = Math.Abs(CRV[J + o_crv]);`
` 655:                              ICMAX = J;`
` 656:                          }`
` 657:                      }`
` 658:                      // *`
` 659:                      // *           If norm(C) < SMINI, use SMINI*identity.`
` 660:                      // *`
` 661:                      if (CMAX < SMINI)`
` 662:                      {`
` 663:                          BNORM = Math.Max(Math.Abs(B[1+1 * LDB + o_b]), Math.Abs(B[2+1 * LDB + o_b]));`
` 664:                          if (SMINI < ONE && BNORM > ONE)`
` 665:                          {`
` 666:                              if (BNORM > BIGNUM * SMINI) SCALE = ONE / BNORM;`
` 667:                          }`
` 668:                          TEMP = SCALE / SMINI;`
` 669:                          X[1+1 * LDX + o_x] = TEMP * B[1+1 * LDB + o_b];`
` 670:                          X[2+1 * LDX + o_x] = TEMP * B[2+1 * LDB + o_b];`
` 671:                          XNORM = TEMP * BNORM;`
` 672:                          INFO = 1;`
` 673:                          return;`
` 674:                      }`
` 675:                      // *`
` 676:                      // *           Gaussian elimination with complete pivoting.`
` 677:                      // *`
` 678:                      UR11 = CRV[ICMAX + o_crv];`
` 679:                      CR21 = CRV[IPIVOT[2+ICMAX * 4 + o_ipivot] + o_crv];`
` 680:                      UR12 = CRV[IPIVOT[3+ICMAX * 4 + o_ipivot] + o_crv];`
` 681:                      CR22 = CRV[IPIVOT[4+ICMAX * 4 + o_ipivot] + o_crv];`
` 682:                      UR11R = ONE / UR11;`
` 683:                      LR21 = UR11R * CR21;`
` 684:                      UR22 = CR22 - UR12 * LR21;`
` 685:                      // *`
` 686:                      // *           If smaller pivot < SMINI, use SMINI`
` 687:                      // *`
` 688:                      if (Math.Abs(UR22) < SMINI)`
` 689:                      {`
` 690:                          UR22 = SMINI;`
` 691:                          INFO = 1;`
` 692:                      }`
` 693:                      if (RSWAP[ICMAX + o_rswap])`
` 694:                      {`
` 695:                          BR1 = B[2+1 * LDB + o_b];`
` 696:                          BR2 = B[1+1 * LDB + o_b];`
` 697:                      }`
` 698:                      else`
` 699:                      {`
` 700:                          BR1 = B[1+1 * LDB + o_b];`
` 701:                          BR2 = B[2+1 * LDB + o_b];`
` 702:                      }`
` 703:                      BR2 = BR2 - LR21 * BR1;`
` 704:                      BBND = Math.Max(Math.Abs(BR1 * (UR22 * UR11R)), Math.Abs(BR2));`
` 705:                      if (BBND > ONE && Math.Abs(UR22) < ONE)`
` 706:                      {`
` 707:                          if (BBND >= BIGNUM * Math.Abs(UR22)) SCALE = ONE / BBND;`
` 708:                      }`
` 709:                      // *`
` 710:                      XR2 = (BR2 * SCALE) / UR22;`
` 711:                      XR1 = (SCALE * BR1) * UR11R - XR2 * (UR11R * UR12);`
` 712:                      if (ZSWAP[ICMAX + o_zswap])`
` 713:                      {`
` 714:                          X[1+1 * LDX + o_x] = XR2;`
` 715:                          X[2+1 * LDX + o_x] = XR1;`
` 716:                      }`
` 717:                      else`
` 718:                      {`
` 719:                          X[1+1 * LDX + o_x] = XR1;`
` 720:                          X[2+1 * LDX + o_x] = XR2;`
` 721:                      }`
` 722:                      XNORM = Math.Max(Math.Abs(XR1), Math.Abs(XR2));`
` 723:                      // *`
` 724:                      // *           Further scaling if  norm(A) norm(X) > overflow`
` 725:                      // *`
` 726:                      if (XNORM > ONE && CMAX > ONE)`
` 727:                      {`
` 728:                          if (XNORM > BIGNUM / CMAX)`
` 729:                          {`
` 730:                              TEMP = CMAX / BIGNUM;`
` 731:                              X[1+1 * LDX + o_x] = TEMP * X[1+1 * LDX + o_x];`
` 732:                              X[2+1 * LDX + o_x] = TEMP * X[2+1 * LDX + o_x];`
` 733:                              XNORM = TEMP * XNORM;`
` 734:                              SCALE = TEMP * SCALE;`
` 735:                          }`
` 736:                      }`
` 737:                  }`
` 738:                  else`
` 739:                  {`
` 740:                      // *`
` 741:                      // *           Complex 2x2 system  (w is complex)`
` 742:                      // *`
` 743:                      // *           Find the largest element in C`
` 744:                      // *`
` 745:                      `
` 746:                      `
` 747:                      // *-----------------------------INICIA MI MODIFICACION--------------------------`
` 748:                      // *            CI( 1, 1 ) = -WI*D1`
` 749:                      // *            CI( 2, 1 ) = ZERO`
` 750:                      // *            CI( 1, 2 ) = ZERO`
` 751:                      // *            CI( 2, 2 ) = -WI*D2`
` 752:                      CIV[1 + o_civ] =  - WI * D1;`
` 753:                      CIV[2 + o_civ] = ZERO;`
` 754:                      CIV[3 + o_civ] = ZERO;`
` 755:                      CIV[4 + o_civ] =  - WI * D2;`
` 756:                      // *-----------------------------TERMINA MI MODIFICACION-------------------------`
` 757:                      `
` 758:                      `
` 759:                      `
` 760:                      `
` 761:                      `
` 762:                      CMAX = ZERO;`
` 763:                      ICMAX = 0;`
` 764:                      // *`
` 765:                      for (J = 1; J <= 4; J++)`
` 766:                      {`
` 767:                          if (Math.Abs(CRV[J + o_crv]) + Math.Abs(CIV[J + o_civ]) > CMAX)`
` 768:                          {`
` 769:                              CMAX = Math.Abs(CRV[J + o_crv]) + Math.Abs(CIV[J + o_civ]);`
` 770:                              ICMAX = J;`
` 771:                          }`
` 772:                      }`
` 773:                      // *`
` 774:                      // *           If norm(C) < SMINI, use SMINI*identity.`
` 775:                      // *`
` 776:                      if (CMAX < SMINI)`
` 777:                      {`
` 778:                          BNORM = Math.Max(Math.Abs(B[1+1 * LDB + o_b]) + Math.Abs(B[1+2 * LDB + o_b]), Math.Abs(B[2+1 * LDB + o_b]) + Math.Abs(B[2+2 * LDB + o_b]));`
` 779:                          if (SMINI < ONE && BNORM > ONE)`
` 780:                          {`
` 781:                              if (BNORM > BIGNUM * SMINI) SCALE = ONE / BNORM;`
` 782:                          }`
` 783:                          TEMP = SCALE / SMINI;`
` 784:                          X[1+1 * LDX + o_x] = TEMP * B[1+1 * LDB + o_b];`
` 785:                          X[2+1 * LDX + o_x] = TEMP * B[2+1 * LDB + o_b];`
` 786:                          X[1+2 * LDX + o_x] = TEMP * B[1+2 * LDB + o_b];`
` 787:                          X[2+2 * LDX + o_x] = TEMP * B[2+2 * LDB + o_b];`
` 788:                          XNORM = TEMP * BNORM;`
` 789:                          INFO = 1;`
` 790:                          return;`
` 791:                      }`
` 792:                      // *`
` 793:                      // *           Gaussian elimination with complete pivoting.`
` 794:                      // *`
` 795:                      UR11 = CRV[ICMAX + o_crv];`
` 796:                      UI11 = CIV[ICMAX + o_civ];`
` 797:                      CR21 = CRV[IPIVOT[2+ICMAX * 4 + o_ipivot] + o_crv];`
` 798:                      CI21 = CIV[IPIVOT[2+ICMAX * 4 + o_ipivot] + o_civ];`
` 799:                      UR12 = CRV[IPIVOT[3+ICMAX * 4 + o_ipivot] + o_crv];`
` 800:                      UI12 = CIV[IPIVOT[3+ICMAX * 4 + o_ipivot] + o_civ];`
` 801:                      CR22 = CRV[IPIVOT[4+ICMAX * 4 + o_ipivot] + o_crv];`
` 802:                      CI22 = CIV[IPIVOT[4+ICMAX * 4 + o_ipivot] + o_civ];`
` 803:                      if (ICMAX == 1 || ICMAX == 4)`
` 804:                      {`
` 805:                          // *`
` 806:                          // *              Code when off-diagonals of pivoted C are real`
` 807:                          // *`
` 808:                          if (Math.Abs(UR11) > Math.Abs(UI11))`
` 809:                          {`
` 810:                              TEMP = UI11 / UR11;`
` 811:                              UR11R = ONE / (UR11 * (ONE + Math.Pow(TEMP,2)));`
` 812:                              UI11R =  - TEMP * UR11R;`
` 813:                          }`
` 814:                          else`
` 815:                          {`
` 816:                              TEMP = UR11 / UI11;`
` 817:                              UI11R =  - ONE / (UI11 * (ONE + Math.Pow(TEMP,2)));`
` 818:                              UR11R =  - TEMP * UI11R;`
` 819:                          }`
` 820:                          LR21 = CR21 * UR11R;`
` 821:                          LI21 = CR21 * UI11R;`
` 822:                          UR12S = UR12 * UR11R;`
` 823:                          UI12S = UR12 * UI11R;`
` 824:                          UR22 = CR22 - UR12 * LR21;`
` 825:                          UI22 = CI22 - UR12 * LI21;`
` 826:                      }`
` 827:                      else`
` 828:                      {`
` 829:                          // *`
` 830:                          // *              Code when diagonals of pivoted C are real`
` 831:                          // *`
` 832:                          UR11R = ONE / UR11;`
` 833:                          UI11R = ZERO;`
` 834:                          LR21 = CR21 * UR11R;`
` 835:                          LI21 = CI21 * UR11R;`
` 836:                          UR12S = UR12 * UR11R;`
` 837:                          UI12S = UI12 * UR11R;`
` 838:                          UR22 = CR22 - UR12 * LR21 + UI12 * LI21;`
` 839:                          UI22 =  - UR12 * LI21 - UI12 * LR21;`
` 840:                      }`
` 841:                      U22ABS = Math.Abs(UR22) + Math.Abs(UI22);`
` 842:                      // *`
` 843:                      // *           If smaller pivot < SMINI, use SMINI`
` 844:                      // *`
` 845:                      if (U22ABS < SMINI)`
` 846:                      {`
` 847:                          UR22 = SMINI;`
` 848:                          UI22 = ZERO;`
` 849:                          INFO = 1;`
` 850:                      }`
` 851:                      if (RSWAP[ICMAX + o_rswap])`
` 852:                      {`
` 853:                          BR2 = B[1+1 * LDB + o_b];`
` 854:                          BR1 = B[2+1 * LDB + o_b];`
` 855:                          BI2 = B[1+2 * LDB + o_b];`
` 856:                          BI1 = B[2+2 * LDB + o_b];`
` 857:                      }`
` 858:                      else`
` 859:                      {`
` 860:                          BR1 = B[1+1 * LDB + o_b];`
` 861:                          BR2 = B[2+1 * LDB + o_b];`
` 862:                          BI1 = B[1+2 * LDB + o_b];`
` 863:                          BI2 = B[2+2 * LDB + o_b];`
` 864:                      }`
` 865:                      BR2 = BR2 - LR21 * BR1 + LI21 * BI1;`
` 866:                      BI2 = BI2 - LI21 * BR1 - LR21 * BI1;`
` 867:                      BBND = Math.Max((Math.Abs(BR1) + Math.Abs(BI1)) * (U22ABS * (Math.Abs(UR11R) + Math.Abs(UI11R))), Math.Abs(BR2) + Math.Abs(BI2));`
` 868:                      if (BBND > ONE && U22ABS < ONE)`
` 869:                      {`
` 870:                          if (BBND >= BIGNUM * U22ABS)`
` 871:                          {`
` 872:                              SCALE = ONE / BBND;`
` 873:                              BR1 = SCALE * BR1;`
` 874:                              BI1 = SCALE * BI1;`
` 875:                              BR2 = SCALE * BR2;`
` 876:                              BI2 = SCALE * BI2;`
` 877:                          }`
` 878:                      }`
` 879:                      // *`
` 880:                      this._dladiv.Run(BR2, BI2, UR22, UI22, ref XR2, ref XI2);`
` 881:                      XR1 = UR11R * BR1 - UI11R * BI1 - UR12S * XR2 + UI12S * XI2;`
` 882:                      XI1 = UI11R * BR1 + UR11R * BI1 - UI12S * XR2 - UR12S * XI2;`
` 883:                      if (ZSWAP[ICMAX + o_zswap])`
` 884:                      {`
` 885:                          X[1+1 * LDX + o_x] = XR2;`
` 886:                          X[2+1 * LDX + o_x] = XR1;`
` 887:                          X[1+2 * LDX + o_x] = XI2;`
` 888:                          X[2+2 * LDX + o_x] = XI1;`
` 889:                      }`
` 890:                      else`
` 891:                      {`
` 892:                          X[1+1 * LDX + o_x] = XR1;`
` 893:                          X[2+1 * LDX + o_x] = XR2;`
` 894:                          X[1+2 * LDX + o_x] = XI1;`
` 895:                          X[2+2 * LDX + o_x] = XI2;`
` 896:                      }`
` 897:                      XNORM = Math.Max(Math.Abs(XR1) + Math.Abs(XI1), Math.Abs(XR2) + Math.Abs(XI2));`
` 898:                      // *`
` 899:                      // *           Further scaling if  norm(A) norm(X) > overflow`
` 900:                      // *`
` 901:                      if (XNORM > ONE && CMAX > ONE)`
` 902:                      {`
` 903:                          if (XNORM > BIGNUM / CMAX)`
` 904:                          {`
` 905:                              TEMP = CMAX / BIGNUM;`
` 906:                              X[1+1 * LDX + o_x] = TEMP * X[1+1 * LDX + o_x];`
` 907:                              X[2+1 * LDX + o_x] = TEMP * X[2+1 * LDX + o_x];`
` 908:                              X[1+2 * LDX + o_x] = TEMP * X[1+2 * LDX + o_x];`
` 909:                              X[2+2 * LDX + o_x] = TEMP * X[2+2 * LDX + o_x];`
` 910:                              XNORM = TEMP * XNORM;`
` 911:                              SCALE = TEMP * SCALE;`
` 912:                          }`
` 913:                      }`
` 914:                  }`
` 915:              }`
` 916:              // *`
` 917:              return;`
` 918:              // *`
` 919:              // *     End of DLALN2`
` 920:              // *`
` 921:   `
` 922:              #endregion`
` 923:   `
` 924:          }`
` 925:      }`
` 926:  }`