`   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 driver routine (version 3.1) --`
`  21:      /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..`
`  22:      /// November 2006`
`  23:      /// Purpose`
`  24:      /// =======`
`  25:      /// `
`  26:      /// DGELSY computes the minimum-norm solution to a real linear least`
`  27:      /// squares problem:`
`  28:      /// minimize || A * X - B ||`
`  29:      /// using a complete orthogonal factorization of A.  A is an M-by-N`
`  30:      /// matrix which may be rank-deficient.`
`  31:      /// `
`  32:      /// Several right hand side vectors b and solution vectors x can be`
`  33:      /// handled in a single call; they are stored as the columns of the`
`  34:      /// M-by-NRHS right hand side matrix B and the N-by-NRHS solution`
`  35:      /// matrix X.`
`  36:      /// `
`  37:      /// The routine first computes a QR factorization with column pivoting:`
`  38:      /// A * P = Q * [ R11 R12 ]`
`  39:      /// [  0  R22 ]`
`  40:      /// with R11 defined as the largest leading submatrix whose estimated`
`  41:      /// condition number is less than 1/RCOND.  The order of R11, RANK,`
`  42:      /// is the effective rank of A.`
`  43:      /// `
`  44:      /// Then, R22 is considered to be negligible, and R12 is annihilated`
`  45:      /// by orthogonal transformations from the right, arriving at the`
`  46:      /// complete orthogonal factorization:`
`  47:      /// A * P = Q * [ T11 0 ] * Z`
`  48:      /// [  0  0 ]`
`  49:      /// The minimum-norm solution is then`
`  50:      /// X = P * Z' [ inv(T11)*Q1'*B ]`
`  51:      /// [        0       ]`
`  52:      /// where Q1 consists of the first RANK columns of Q.`
`  53:      /// `
`  54:      /// This routine is basically identical to the original xGELSX except`
`  55:      /// three differences:`
`  56:      /// o The call to the subroutine xGEQPF has been substituted by the`
`  57:      /// the call to the subroutine xGEQP3. This subroutine is a Blas-3`
`  58:      /// version of the QR factorization with column pivoting.`
`  59:      /// o Matrix B (the right hand side) is updated with Blas-3.`
`  60:      /// o The permutation of matrix B (the right hand side) is faster and`
`  61:      /// more simple.`
`  62:      /// `
`  63:      ///</summary>`
`  64:      public class DGELSY`
`  65:      {`
`  66:      `
`  67:   `
`  68:          #region Dependencies`
`  69:          `
`  70:          ILAENV _ilaenv; DLAMCH _dlamch; DLANGE _dlange; DCOPY _dcopy; DGEQP3 _dgeqp3; DLABAD _dlabad; DLAIC1 _dlaic1; `
`  71:          DLASCL _dlascl;DLASET _dlaset; DORMQR _dormqr; DORMRZ _dormrz; DTRSM _dtrsm; DTZRZF _dtzrzf; XERBLA _xerbla; `
`  72:   `
`  73:          #endregion`
`  74:   `
`  75:   `
`  76:          #region Fields`
`  77:          `
`  78:          const int IMAX = 1; const int IMIN = 2; const double ZERO = 0.0E+0; const double ONE = 1.0E+0; bool LQUERY = false; `
`  79:          int I = 0;int IASCL = 0; int IBSCL = 0; int ISMAX = 0; int ISMIN = 0; int J = 0; int LWKMIN = 0; int LWKOPT = 0; `
`  80:          int MN = 0;int NB = 0; int NB1 = 0; int NB2 = 0; int NB3 = 0; int NB4 = 0; double ANRM = 0; double BIGNUM = 0; `
`  81:          double BNRM = 0;double C1 = 0; double C2 = 0; double S1 = 0; double S2 = 0; double SMAX = 0; double SMAXPR = 0; `
`  82:          double SMIN = 0;double SMINPR = 0; double SMLNUM = 0; double WSIZE = 0; `
`  83:   `
`  84:          #endregion`
`  85:   `
`  86:          public DGELSY(ILAENV ilaenv, DLAMCH dlamch, DLANGE dlange, DCOPY dcopy, DGEQP3 dgeqp3, DLABAD dlabad, DLAIC1 dlaic1, DLASCL dlascl, DLASET dlaset, DORMQR dormqr`
`  87:                        , DORMRZ dormrz, DTRSM dtrsm, DTZRZF dtzrzf, XERBLA xerbla)`
`  88:          {`
`  89:      `
`  90:   `
`  91:              #region Set Dependencies`
`  92:              `
`  93:              this._ilaenv = ilaenv; this._dlamch = dlamch; this._dlange = dlange; this._dcopy = dcopy; this._dgeqp3 = dgeqp3; `
`  94:              this._dlabad = dlabad;this._dlaic1 = dlaic1; this._dlascl = dlascl; this._dlaset = dlaset; this._dormqr = dormqr; `
`  95:              this._dormrz = dormrz;this._dtrsm = dtrsm; this._dtzrzf = dtzrzf; this._xerbla = xerbla; `
`  96:   `
`  97:              #endregion`
`  98:   `
`  99:          }`
` 100:      `
` 101:          public DGELSY()`
` 102:          {`
` 103:      `
` 104:   `
` 105:              #region Dependencies (Initialization)`
` 106:              `
` 107:              IEEECK ieeeck = new IEEECK();`
` 108:              IPARMQ iparmq = new IPARMQ();`
` 109:              LSAME lsame = new LSAME();`
` 110:              DLAMC3 dlamc3 = new DLAMC3();`
` 111:              DLASSQ dlassq = new DLASSQ();`
` 112:              DCOPY dcopy = new DCOPY();`
` 113:              XERBLA xerbla = new XERBLA();`
` 114:              DLAPY2 dlapy2 = new DLAPY2();`
` 115:              DNRM2 dnrm2 = new DNRM2();`
` 116:              DSCAL dscal = new DSCAL();`
` 117:              DSWAP dswap = new DSWAP();`
` 118:              IDAMAX idamax = new IDAMAX();`
` 119:              DLABAD dlabad = new DLABAD();`
` 120:              DDOT ddot = new DDOT();`
` 121:              DAXPY daxpy = new DAXPY();`
` 122:              ILAENV ilaenv = new ILAENV(ieeeck, iparmq);`
` 123:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);`
` 124:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);`
` 125:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);`
` 126:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);`
` 127:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);`
` 128:              DLANGE dlange = new DLANGE(dlassq, lsame);`
` 129:              DGEMV dgemv = new DGEMV(lsame, xerbla);`
` 130:              DGER dger = new DGER(xerbla);`
` 131:              DLARF dlarf = new DLARF(dgemv, dger, lsame);`
` 132:              DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);`
` 133:              DGEQR2 dgeqr2 = new DGEQR2(dlarf, dlarfg, xerbla);`
` 134:              DGEMM dgemm = new DGEMM(lsame, xerbla);`
` 135:              DTRMM dtrmm = new DTRMM(lsame, xerbla);`
` 136:              DLARFB dlarfb = new DLARFB(lsame, dcopy, dgemm, dtrmm);`
` 137:              DTRMV dtrmv = new DTRMV(lsame, xerbla);`
` 138:              DLARFT dlarft = new DLARFT(dgemv, dtrmv, lsame);`
` 139:              DGEQRF dgeqrf = new DGEQRF(dgeqr2, dlarfb, dlarft, xerbla, ilaenv);`
` 140:              DLAQP2 dlaqp2 = new DLAQP2(dlarf, dlarfg, dswap, idamax, dlamch, dnrm2);`
` 141:              DLAQPS dlaqps = new DLAQPS(dgemm, dgemv, dlarfg, dswap, idamax, dlamch, dnrm2);`
` 142:              DORM2R dorm2r = new DORM2R(lsame, dlarf, xerbla);`
` 143:              DORMQR dormqr = new DORMQR(lsame, ilaenv, dlarfb, dlarft, dorm2r, xerbla);`
` 144:              DGEQP3 dgeqp3 = new DGEQP3(dgeqrf, dlaqp2, dlaqps, dormqr, dswap, xerbla, ilaenv, dnrm2);`
` 145:              DLAIC1 dlaic1 = new DLAIC1(ddot, dlamch);`
` 146:              DLASCL dlascl = new DLASCL(lsame, dlamch, xerbla);`
` 147:              DLASET dlaset = new DLASET(lsame);`
` 148:              DLARZB dlarzb = new DLARZB(lsame, dcopy, dgemm, dtrmm, xerbla);`
` 149:              DLARZT dlarzt = new DLARZT(dgemv, dtrmv, xerbla, lsame);`
` 150:              DLARZ dlarz = new DLARZ(daxpy, dcopy, dgemv, dger, lsame);`
` 151:              DORMR3 dormr3 = new DORMR3(lsame, dlarz, xerbla);`
` 152:              DORMRZ dormrz = new DORMRZ(lsame, ilaenv, dlarzb, dlarzt, dormr3, xerbla);`
` 153:              DTRSM dtrsm = new DTRSM(lsame, xerbla);`
` 154:              DLATRZ dlatrz = new DLATRZ(dlarfg, dlarz);`
` 155:              DTZRZF dtzrzf = new DTZRZF(dlarzb, dlarzt, dlatrz, xerbla, ilaenv);`
` 156:   `
` 157:              #endregion`
` 158:   `
` 159:   `
` 160:              #region Set Dependencies`
` 161:              `
` 162:              this._ilaenv = ilaenv; this._dlamch = dlamch; this._dlange = dlange; this._dcopy = dcopy; this._dgeqp3 = dgeqp3; `
` 163:              this._dlabad = dlabad;this._dlaic1 = dlaic1; this._dlascl = dlascl; this._dlaset = dlaset; this._dormqr = dormqr; `
` 164:              this._dormrz = dormrz;this._dtrsm = dtrsm; this._dtzrzf = dtzrzf; this._xerbla = xerbla; `
` 165:   `
` 166:              #endregion`
` 167:   `
` 168:          }`
` 169:          /// <summary>`
` 170:          /// Purpose`
` 171:          /// =======`
` 172:          /// `
` 173:          /// DGELSY computes the minimum-norm solution to a real linear least`
` 174:          /// squares problem:`
` 175:          /// minimize || A * X - B ||`
` 176:          /// using a complete orthogonal factorization of A.  A is an M-by-N`
` 177:          /// matrix which may be rank-deficient.`
` 178:          /// `
` 179:          /// Several right hand side vectors b and solution vectors x can be`
` 180:          /// handled in a single call; they are stored as the columns of the`
` 181:          /// M-by-NRHS right hand side matrix B and the N-by-NRHS solution`
` 182:          /// matrix X.`
` 183:          /// `
` 184:          /// The routine first computes a QR factorization with column pivoting:`
` 185:          /// A * P = Q * [ R11 R12 ]`
` 186:          /// [  0  R22 ]`
` 187:          /// with R11 defined as the largest leading submatrix whose estimated`
` 188:          /// condition number is less than 1/RCOND.  The order of R11, RANK,`
` 189:          /// is the effective rank of A.`
` 190:          /// `
` 191:          /// Then, R22 is considered to be negligible, and R12 is annihilated`
` 192:          /// by orthogonal transformations from the right, arriving at the`
` 193:          /// complete orthogonal factorization:`
` 194:          /// A * P = Q * [ T11 0 ] * Z`
` 195:          /// [  0  0 ]`
` 196:          /// The minimum-norm solution is then`
` 197:          /// X = P * Z' [ inv(T11)*Q1'*B ]`
` 198:          /// [        0       ]`
` 199:          /// where Q1 consists of the first RANK columns of Q.`
` 200:          /// `
` 201:          /// This routine is basically identical to the original xGELSX except`
` 202:          /// three differences:`
` 203:          /// o The call to the subroutine xGEQPF has been substituted by the`
` 204:          /// the call to the subroutine xGEQP3. This subroutine is a Blas-3`
` 205:          /// version of the QR factorization with column pivoting.`
` 206:          /// o Matrix B (the right hand side) is updated with Blas-3.`
` 207:          /// o The permutation of matrix B (the right hand side) is faster and`
` 208:          /// more simple.`
` 209:          /// `
` 210:          ///</summary>`
` 211:          /// <param name="M">`
` 212:          /// (input) INTEGER`
` 213:          /// The number of rows of the matrix A.  M .GE. 0.`
` 214:          ///</param>`
` 215:          /// <param name="N">`
` 216:          /// (input) INTEGER`
` 217:          /// The number of columns of the matrix A.  N .GE. 0.`
` 218:          ///</param>`
` 219:          /// <param name="NRHS">`
` 220:          /// (input) INTEGER`
` 221:          /// The number of right hand sides, i.e., the number of`
` 222:          /// columns of matrices B and X. NRHS .GE. 0.`
` 223:          ///</param>`
` 224:          /// <param name="A">`
` 225:          /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)`
` 226:          /// On entry, the M-by-N matrix A.`
` 227:          /// On exit, A has been overwritten by details of its`
` 228:          /// complete orthogonal factorization.`
` 229:          ///</param>`
` 230:          /// <param name="LDA">`
` 231:          /// (input) INTEGER`
` 232:          /// The leading dimension of the array A.  LDA .GE. max(1,M).`
` 233:          ///</param>`
` 234:          /// <param name="B">`
` 235:          /// (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)`
` 236:          /// On entry, the M-by-NRHS right hand side matrix B.`
` 237:          /// On exit, the N-by-NRHS solution matrix X.`
` 238:          ///</param>`
` 239:          /// <param name="LDB">`
` 240:          /// (input) INTEGER`
` 241:          /// The leading dimension of the array B. LDB .GE. max(1,M,N).`
` 242:          ///</param>`
` 243:          /// <param name="JPVT">`
` 244:          /// (input/output) INTEGER array, dimension (N)`
` 245:          /// On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted`
` 246:          /// to the front of AP, otherwise column i is a free column.`
` 247:          /// On exit, if JPVT(i) = k, then the i-th column of AP`
` 248:          /// was the k-th column of A.`
` 249:          ///</param>`
` 250:          /// <param name="RCOND">`
` 251:          /// (input) DOUBLE PRECISION`
` 252:          /// RCOND is used to determine the effective rank of A, which`
` 253:          /// is defined as the order of the largest leading triangular`
` 254:          /// submatrix R11 in the QR factorization with pivoting of A,`
` 255:          /// whose estimated condition number .LT. 1/RCOND.`
` 256:          ///</param>`
` 257:          /// <param name="RANK">`
` 258:          /// (output) INTEGER`
` 259:          /// The effective rank of A, i.e., the order of the submatrix`
` 260:          /// R11.  This is the same as the order of the submatrix T11`
` 261:          /// in the complete orthogonal factorization of A.`
` 262:          ///</param>`
` 263:          /// <param name="WORK">`
` 264:          /// (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))`
` 265:          /// On exit, if INFO = 0, WORK(1) returns the optimal LWORK.`
` 266:          ///</param>`
` 267:          /// <param name="LWORK">`
` 268:          /// (input) INTEGER`
` 269:          /// The dimension of the array WORK.`
` 270:          /// The unblocked strategy requires that:`
` 271:          /// LWORK .GE. MAX( MN+3*N+1, 2*MN+NRHS ),`
` 272:          /// where MN = min( M, N ).`
` 273:          /// The block algorithm requires that:`
` 274:          /// LWORK .GE. MAX( MN+2*N+NB*(N+1), 2*MN+NB*NRHS ),`
` 275:          /// where NB is an upper bound on the blocksize returned`
` 276:          /// by ILAENV for the routines DGEQP3, DTZRZF, STZRQF, DORMQR,`
` 277:          /// and DORMRZ.`
` 278:          /// `
` 279:          /// If LWORK = -1, then a workspace query is assumed; the routine`
` 280:          /// only calculates the optimal size of the WORK array, returns`
` 281:          /// this value as the first entry of the WORK array, and no error`
` 282:          /// message related to LWORK is issued by XERBLA.`
` 283:          ///</param>`
` 284:          /// <param name="INFO">`
` 285:          /// (output) INTEGER`
` 286:          /// = 0: successful exit`
` 287:          /// .LT. 0: If INFO = -i, the i-th argument had an illegal value.`
` 288:          ///</param>`
` 289:          public void Run(int M, int N, int NRHS, ref double[] A, int offset_a, int LDA, ref double[] B, int offset_b`
` 290:                           , int LDB, ref int[] JPVT, int offset_jpvt, double RCOND, ref int RANK, ref double[] WORK, int offset_work, int LWORK`
` 291:                           , ref int INFO)`
` 292:          {`
` 293:   `
` 294:              #region Array Index Correction`
` 295:              `
` 296:               int o_a = -1 - LDA + offset_a;  int o_b = -1 - LDB + offset_b;  int o_jpvt = -1 + offset_jpvt; `
` 297:               int o_work = -1 + offset_work;`
` 298:   `
` 299:              #endregion`
` 300:   `
` 301:   `
` 302:              #region Prolog`
` 303:              `
` 304:              // *`
` 305:              // *  -- LAPACK driver routine (version 3.1) --`
` 306:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..`
` 307:              // *     November 2006`
` 308:              // *`
` 309:              // *     .. Scalar Arguments ..`
` 310:              // *     ..`
` 311:              // *     .. Array Arguments ..`
` 312:              // *     ..`
` 313:              // *`
` 314:              // *  Purpose`
` 315:              // *  =======`
` 316:              // *`
` 317:              // *  DGELSY computes the minimum-norm solution to a real linear least`
` 318:              // *  squares problem:`
` 319:              // *      minimize || A * X - B ||`
` 320:              // *  using a complete orthogonal factorization of A.  A is an M-by-N`
` 321:              // *  matrix which may be rank-deficient.`
` 322:              // *`
` 323:              // *  Several right hand side vectors b and solution vectors x can be`
` 324:              // *  handled in a single call; they are stored as the columns of the`
` 325:              // *  M-by-NRHS right hand side matrix B and the N-by-NRHS solution`
` 326:              // *  matrix X.`
` 327:              // *`
` 328:              // *  The routine first computes a QR factorization with column pivoting:`
` 329:              // *      A * P = Q * [ R11 R12 ]`
` 330:              // *                  [  0  R22 ]`
` 331:              // *  with R11 defined as the largest leading submatrix whose estimated`
` 332:              // *  condition number is less than 1/RCOND.  The order of R11, RANK,`
` 333:              // *  is the effective rank of A.`
` 334:              // *`
` 335:              // *  Then, R22 is considered to be negligible, and R12 is annihilated`
` 336:              // *  by orthogonal transformations from the right, arriving at the`
` 337:              // *  complete orthogonal factorization:`
` 338:              // *     A * P = Q * [ T11 0 ] * Z`
` 339:              // *                 [  0  0 ]`
` 340:              // *  The minimum-norm solution is then`
` 341:              // *     X = P * Z' [ inv(T11)*Q1'*B ]`
` 342:              // *                [        0       ]`
` 343:              // *  where Q1 consists of the first RANK columns of Q.`
` 344:              // *`
` 345:              // *  This routine is basically identical to the original xGELSX except`
` 346:              // *  three differences:`
` 347:              // *    o The call to the subroutine xGEQPF has been substituted by the`
` 348:              // *      the call to the subroutine xGEQP3. This subroutine is a Blas-3`
` 349:              // *      version of the QR factorization with column pivoting.`
` 350:              // *    o Matrix B (the right hand side) is updated with Blas-3.`
` 351:              // *    o The permutation of matrix B (the right hand side) is faster and`
` 352:              // *      more simple.`
` 353:              // *`
` 354:              // *  Arguments`
` 355:              // *  =========`
` 356:              // *`
` 357:              // *  M       (input) INTEGER`
` 358:              // *          The number of rows of the matrix A.  M >= 0.`
` 359:              // *`
` 360:              // *  N       (input) INTEGER`
` 361:              // *          The number of columns of the matrix A.  N >= 0.`
` 362:              // *`
` 363:              // *  NRHS    (input) INTEGER`
` 364:              // *          The number of right hand sides, i.e., the number of`
` 365:              // *          columns of matrices B and X. NRHS >= 0.`
` 366:              // *`
` 367:              // *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)`
` 368:              // *          On entry, the M-by-N matrix A.`
` 369:              // *          On exit, A has been overwritten by details of its`
` 370:              // *          complete orthogonal factorization.`
` 371:              // *`
` 372:              // *  LDA     (input) INTEGER`
` 373:              // *          The leading dimension of the array A.  LDA >= max(1,M).`
` 374:              // *`
` 375:              // *  B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)`
` 376:              // *          On entry, the M-by-NRHS right hand side matrix B.`
` 377:              // *          On exit, the N-by-NRHS solution matrix X.`
` 378:              // *`
` 379:              // *  LDB     (input) INTEGER`
` 380:              // *          The leading dimension of the array B. LDB >= max(1,M,N).`
` 381:              // *`
` 382:              // *  JPVT    (input/output) INTEGER array, dimension (N)`
` 383:              // *          On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted`
` 384:              // *          to the front of AP, otherwise column i is a free column.`
` 385:              // *          On exit, if JPVT(i) = k, then the i-th column of AP`
` 386:              // *          was the k-th column of A.`
` 387:              // *`
` 388:              // *  RCOND   (input) DOUBLE PRECISION`
` 389:              // *          RCOND is used to determine the effective rank of A, which`
` 390:              // *          is defined as the order of the largest leading triangular`
` 391:              // *          submatrix R11 in the QR factorization with pivoting of A,`
` 392:              // *          whose estimated condition number < 1/RCOND.`
` 393:              // *`
` 394:              // *  RANK    (output) INTEGER`
` 395:              // *          The effective rank of A, i.e., the order of the submatrix`
` 396:              // *          R11.  This is the same as the order of the submatrix T11`
` 397:              // *          in the complete orthogonal factorization of A.`
` 398:              // *`
` 399:              // *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))`
` 400:              // *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.`
` 401:              // *`
` 402:              // *  LWORK   (input) INTEGER`
` 403:              // *          The dimension of the array WORK.`
` 404:              // *          The unblocked strategy requires that:`
` 405:              // *             LWORK >= MAX( MN+3*N+1, 2*MN+NRHS ),`
` 406:              // *          where MN = min( M, N ).`
` 407:              // *          The block algorithm requires that:`
` 408:              // *             LWORK >= MAX( MN+2*N+NB*(N+1), 2*MN+NB*NRHS ),`
` 409:              // *          where NB is an upper bound on the blocksize returned`
` 410:              // *          by ILAENV for the routines DGEQP3, DTZRZF, STZRQF, DORMQR,`
` 411:              // *          and DORMRZ.`
` 412:              // *`
` 413:              // *          If LWORK = -1, then a workspace query is assumed; the routine`
` 414:              // *          only calculates the optimal size of the WORK array, returns`
` 415:              // *          this value as the first entry of the WORK array, and no error`
` 416:              // *          message related to LWORK is issued by XERBLA.`
` 417:              // *`
` 418:              // *  INFO    (output) INTEGER`
` 419:              // *          = 0: successful exit`
` 420:              // *          < 0: If INFO = -i, the i-th argument had an illegal value.`
` 421:              // *`
` 422:              // *  Further Details`
` 423:              // *  ===============`
` 424:              // *`
` 425:              // *  Based on contributions by`
` 426:              // *    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA`
` 427:              // *    E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain`
` 428:              // *    G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain`
` 429:              // *`
` 430:              // *  =====================================================================`
` 431:              // *`
` 432:              // *     .. Parameters ..`
` 433:              // *     ..`
` 434:              // *     .. Local Scalars ..`
` 435:              // *     ..`
` 436:              // *     .. External Functions ..`
` 437:              // *     ..`
` 438:              // *     .. External Subroutines ..`
` 439:              // *     ..`
` 440:              // *     .. Intrinsic Functions ..`
` 441:              //      INTRINSIC          ABS, MAX, MIN;`
` 442:              // *     ..`
` 443:              // *     .. Executable Statements ..`
` 444:              // *`
` 445:   `
` 446:              #endregion`
` 447:   `
` 448:   `
` 449:              #region Body`
` 450:              `
` 451:              MN = Math.Min(M, N);`
` 452:              ISMIN = MN + 1;`
` 453:              ISMAX = 2 * MN + 1;`
` 454:              // *`
` 455:              // *     Test the input arguments.`
` 456:              // *`
` 457:              INFO = 0;`
` 458:              LQUERY = (LWORK ==  - 1);`
` 459:              if (M < 0)`
` 460:              {`
` 461:                  INFO =  - 1;`
` 462:              }`
` 463:              else`
` 464:              {`
` 465:                  if (N < 0)`
` 466:                  {`
` 467:                      INFO =  - 2;`
` 468:                  }`
` 469:                  else`
` 470:                  {`
` 471:                      if (NRHS < 0)`
` 472:                      {`
` 473:                          INFO =  - 3;`
` 474:                      }`
` 475:                      else`
` 476:                      {`
` 477:                          if (LDA < Math.Max(1, M))`
` 478:                          {`
` 479:                              INFO =  - 5;`
` 480:                          }`
` 481:                          else`
` 482:                          {`
` 483:                              if (LDB < Math.Max(1, Math.Max(M, N)))`
` 484:                              {`
` 485:                                  INFO =  - 7;`
` 486:                              }`
` 487:                          }`
` 488:                      }`
` 489:                  }`
` 490:              }`
` 491:              // *`
` 492:              // *     Figure out optimal block size`
` 493:              // *`
` 494:              if (INFO == 0)`
` 495:              {`
` 496:                  if (MN == 0 || NRHS == 0)`
` 497:                  {`
` 498:                      LWKMIN = 1;`
` 499:                      LWKOPT = 1;`
` 500:                  }`
` 501:                  else`
` 502:                  {`
` 503:                      NB1 = this._ilaenv.Run(1, "DGEQRF", " ", M, N,  - 1,  - 1);`
` 504:                      NB2 = this._ilaenv.Run(1, "DGERQF", " ", M, N,  - 1,  - 1);`
` 505:                      NB3 = this._ilaenv.Run(1, "DORMQR", " ", M, N, NRHS,  - 1);`
` 506:                      NB4 = this._ilaenv.Run(1, "DORMRQ", " ", M, N, NRHS,  - 1);`
` 507:                      NB = Math.Max(NB1, Math.Max(NB2, Math.Max(NB3, NB4)));`
` 508:                      LWKMIN = MN + Math.Max(2 * MN, Math.Max(N + 1, MN + NRHS));`
` 509:                      LWKOPT = Math.Max(LWKMIN, Math.Max(MN + 2 * N + NB * (N + 1), 2 * MN + NB * NRHS));`
` 510:                  }`
` 511:                  WORK[1 + o_work] = LWKOPT;`
` 512:                  // *`
` 513:                  if (LWORK < LWKMIN && !LQUERY)`
` 514:                  {`
` 515:                      INFO =  - 12;`
` 516:                  }`
` 517:              }`
` 518:              // *`
` 519:              if (INFO != 0)`
` 520:              {`
` 521:                  this._xerbla.Run("DGELSY",  - INFO);`
` 522:                  return;`
` 523:              }`
` 524:              else`
` 525:              {`
` 526:                  if (LQUERY)`
` 527:                  {`
` 528:                      return;`
` 529:                  }`
` 530:              }`
` 531:              // *`
` 532:              // *     Quick return if possible`
` 533:              // *`
` 534:              if (MN == 0 || NRHS == 0)`
` 535:              {`
` 536:                  RANK = 0;`
` 537:                  return;`
` 538:              }`
` 539:              // *`
` 540:              // *     Get machine parameters`
` 541:              // *`
` 542:              SMLNUM = this._dlamch.Run("S") / this._dlamch.Run("P");`
` 543:              BIGNUM = ONE / SMLNUM;`
` 544:              this._dlabad.Run(ref SMLNUM, ref BIGNUM);`
` 545:              // *`
` 546:              // *     Scale A, B if max entries outside range [SMLNUM,BIGNUM]`
` 547:              // *`
` 548:              ANRM = this._dlange.Run("M", M, N, A, offset_a, LDA, ref WORK, offset_work);`
` 549:              IASCL = 0;`
` 550:              if (ANRM > ZERO && ANRM < SMLNUM)`
` 551:              {`
` 552:                  // *`
` 553:                  // *        Scale matrix norm up to SMLNUM`
` 554:                  // *`
` 555:                  this._dlascl.Run("G", 0, 0, ANRM, SMLNUM, M`
` 556:                                   , N, ref A, offset_a, LDA, ref INFO);`
` 557:                  IASCL = 1;`
` 558:              }`
` 559:              else`
` 560:              {`
` 561:                  if (ANRM > BIGNUM)`
` 562:                  {`
` 563:                      // *`
` 564:                      // *        Scale matrix norm down to BIGNUM`
` 565:                      // *`
` 566:                      this._dlascl.Run("G", 0, 0, ANRM, BIGNUM, M`
` 567:                                       , N, ref A, offset_a, LDA, ref INFO);`
` 568:                      IASCL = 2;`
` 569:                  }`
` 570:                  else`
` 571:                  {`
` 572:                      if (ANRM == ZERO)`
` 573:                      {`
` 574:                          // *`
` 575:                          // *        Matrix all zero. Return zero solution.`
` 576:                          // *`
` 577:                          this._dlaset.Run("F", Math.Max(M, N), NRHS, ZERO, ZERO, ref B, offset_b`
` 578:                                           , LDB);`
` 579:                          RANK = 0;`
` 580:                          goto LABEL70;`
` 581:                      }`
` 582:                  }`
` 583:              }`
` 584:              // *`
` 585:              BNRM = this._dlange.Run("M", M, NRHS, B, offset_b, LDB, ref WORK, offset_work);`
` 586:              IBSCL = 0;`
` 587:              if (BNRM > ZERO && BNRM < SMLNUM)`
` 588:              {`
` 589:                  // *`
` 590:                  // *        Scale matrix norm up to SMLNUM`
` 591:                  // *`
` 592:                  this._dlascl.Run("G", 0, 0, BNRM, SMLNUM, M`
` 593:                                   , NRHS, ref B, offset_b, LDB, ref INFO);`
` 594:                  IBSCL = 1;`
` 595:              }`
` 596:              else`
` 597:              {`
` 598:                  if (BNRM > BIGNUM)`
` 599:                  {`
` 600:                      // *`
` 601:                      // *        Scale matrix norm down to BIGNUM`
` 602:                      // *`
` 603:                      this._dlascl.Run("G", 0, 0, BNRM, BIGNUM, M`
` 604:                                       , NRHS, ref B, offset_b, LDB, ref INFO);`
` 605:                      IBSCL = 2;`
` 606:                  }`
` 607:              }`
` 608:              // *`
` 609:              // *     Compute QR factorization with column pivoting of A:`
` 610:              // *        A * P = Q * R`
` 611:              // *`
` 612:              this._dgeqp3.Run(M, N, ref A, offset_a, LDA, ref JPVT, offset_jpvt, ref WORK, 1 + o_work`
` 613:                               , ref WORK, MN + 1 + o_work, LWORK - MN, ref INFO);`
` 614:              WSIZE = MN + WORK[MN + 1 + o_work];`
` 615:              // *`
` 616:              // *     workspace: MN+2*N+NB*(N+1).`
` 617:              // *     Details of Householder rotations stored in WORK(1:MN).`
` 618:              // *`
` 619:              // *     Determine RANK using incremental condition estimation`
` 620:              // *`
` 621:              WORK[ISMIN + o_work] = ONE;`
` 622:              WORK[ISMAX + o_work] = ONE;`
` 623:              SMAX = Math.Abs(A[1+1 * LDA + o_a]);`
` 624:              SMIN = SMAX;`
` 625:              if (Math.Abs(A[1+1 * LDA + o_a]) == ZERO)`
` 626:              {`
` 627:                  RANK = 0;`
` 628:                  this._dlaset.Run("F", Math.Max(M, N), NRHS, ZERO, ZERO, ref B, offset_b`
` 629:                                   , LDB);`
` 630:                  goto LABEL70;`
` 631:              }`
` 632:              else`
` 633:              {`
` 634:                  RANK = 1;`
` 635:              }`
` 636:              // *`
` 637:          LABEL10:;`
` 638:              if (RANK < MN)`
` 639:              {`
` 640:                  I = RANK + 1;`
` 641:                  this._dlaic1.Run(IMIN, RANK, WORK, ISMIN + o_work, SMIN, A, 1+I * LDA + o_a, A[I+I * LDA + o_a]`
` 642:                                   , ref SMINPR, ref S1, ref C1);`
` 643:                  this._dlaic1.Run(IMAX, RANK, WORK, ISMAX + o_work, SMAX, A, 1+I * LDA + o_a, A[I+I * LDA + o_a]`
` 644:                                   , ref SMAXPR, ref S2, ref C2);`
` 645:                  // *`
` 646:                  if (SMAXPR * RCOND <= SMINPR)`
` 647:                  {`
` 648:                      for (I = 1; I <= RANK; I++)`
` 649:                      {`
` 650:                          WORK[ISMIN + I - 1 + o_work] = S1 * WORK[ISMIN + I - 1 + o_work];`
` 651:                          WORK[ISMAX + I - 1 + o_work] = S2 * WORK[ISMAX + I - 1 + o_work];`
` 652:                      }`
` 653:                      WORK[ISMIN + RANK + o_work] = C1;`
` 654:                      WORK[ISMAX + RANK + o_work] = C2;`
` 655:                      SMIN = SMINPR;`
` 656:                      SMAX = SMAXPR;`
` 657:                      RANK = RANK + 1;`
` 658:                      goto LABEL10;`
` 659:                  }`
` 660:              }`
` 661:              // *`
` 662:              // *     workspace: 3*MN.`
` 663:              // *`
` 664:              // *     Logically partition R = [ R11 R12 ]`
` 665:              // *                             [  0  R22 ]`
` 666:              // *     where R11 = R(1:RANK,1:RANK)`
` 667:              // *`
` 668:              // *     [R11,R12] = [ T11, 0 ] * Y`
` 669:              // *`
` 670:              if (RANK < N)`
` 671:              {`
` 672:                  this._dtzrzf.Run(RANK, N, ref A, offset_a, LDA, ref WORK, MN + 1 + o_work, ref WORK, 2 * MN + 1 + o_work`
` 673:                                   , LWORK - 2 * MN, ref INFO);`
` 674:              }`
` 675:              // *`
` 676:              // *     workspace: 2*MN.`
` 677:              // *     Details of Householder rotations stored in WORK(MN+1:2*MN)`
` 678:              // *`
` 679:              // *     B(1:M,1:NRHS) := Q' * B(1:M,1:NRHS)`
` 680:              // *`
` 681:              this._dormqr.Run("Left", "Transpose", M, NRHS, MN, ref A, offset_a`
` 682:                               , LDA, WORK, 1 + o_work, ref B, offset_b, LDB, ref WORK, 2 * MN + 1 + o_work, LWORK - 2 * MN`
` 683:                               , ref INFO);`
` 684:              WSIZE = Math.Max(WSIZE, 2 * MN + WORK[2 * MN + 1 + o_work]);`
` 685:              // *`
` 686:              // *     workspace: 2*MN+NB*NRHS.`
` 687:              // *`
` 688:              // *     B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS)`
` 689:              // *`
` 690:              this._dtrsm.Run("Left", "Upper", "No transpose", "Non-unit", RANK, NRHS`
` 691:                              , ONE, A, offset_a, LDA, ref B, offset_b, LDB);`
` 692:              // *`
` 693:              for (J = 1; J <= NRHS; J++)`
` 694:              {`
` 695:                  for (I = RANK + 1; I <= N; I++)`
` 696:                  {`
` 697:                      B[I+J * LDB + o_b] = ZERO;`
` 698:                  }`
` 699:              }`
` 700:              // *`
` 701:              // *     B(1:N,1:NRHS) := Y' * B(1:N,1:NRHS)`
` 702:              // *`
` 703:              if (RANK < N)`
` 704:              {`
` 705:                  this._dormrz.Run("Left", "Transpose", N, NRHS, RANK, N - RANK`
` 706:                                   , A, offset_a, LDA, WORK, MN + 1 + o_work, ref B, offset_b, LDB, ref WORK, 2 * MN + 1 + o_work`
` 707:                                   , LWORK - 2 * MN, ref INFO);`
` 708:              }`
` 709:              // *`
` 710:              // *     workspace: 2*MN+NRHS.`
` 711:              // *`
` 712:              // *     B(1:N,1:NRHS) := P * B(1:N,1:NRHS)`
` 713:              // *`
` 714:              for (J = 1; J <= NRHS; J++)`
` 715:              {`
` 716:                  for (I = 1; I <= N; I++)`
` 717:                  {`
` 718:                      WORK[JPVT[I + o_jpvt] + o_work] = B[I+J * LDB + o_b];`
` 719:                  }`
` 720:                  this._dcopy.Run(N, WORK, 1 + o_work, 1, ref B, 1+J * LDB + o_b, 1);`
` 721:              }`
` 722:              // *`
` 723:              // *     workspace: N.`
` 724:              // *`
` 725:              // *     Undo scaling`
` 726:              // *`
` 727:              if (IASCL == 1)`
` 728:              {`
` 729:                  this._dlascl.Run("G", 0, 0, ANRM, SMLNUM, N`
` 730:                                   , NRHS, ref B, offset_b, LDB, ref INFO);`
` 731:                  this._dlascl.Run("U", 0, 0, SMLNUM, ANRM, RANK`
` 732:                                   , RANK, ref A, offset_a, LDA, ref INFO);`
` 733:              }`
` 734:              else`
` 735:              {`
` 736:                  if (IASCL == 2)`
` 737:                  {`
` 738:                      this._dlascl.Run("G", 0, 0, ANRM, BIGNUM, N`
` 739:                                       , NRHS, ref B, offset_b, LDB, ref INFO);`
` 740:                      this._dlascl.Run("U", 0, 0, BIGNUM, ANRM, RANK`
` 741:                                       , RANK, ref A, offset_a, LDA, ref INFO);`
` 742:                  }`
` 743:              }`
` 744:              if (IBSCL == 1)`
` 745:              {`
` 746:                  this._dlascl.Run("G", 0, 0, SMLNUM, BNRM, N`
` 747:                                   , NRHS, ref B, offset_b, LDB, ref INFO);`
` 748:              }`
` 749:              else`
` 750:              {`
` 751:                  if (IBSCL == 2)`
` 752:                  {`
` 753:                      this._dlascl.Run("G", 0, 0, BIGNUM, BNRM, N`
` 754:                                       , NRHS, ref B, offset_b, LDB, ref INFO);`
` 755:                  }`
` 756:              }`
` 757:              // *`
` 758:          LABEL70:;`
` 759:              WORK[1 + o_work] = LWKOPT;`
` 760:              // *`
` 761:              return;`
` 762:              // *`
` 763:              // *     End of DGELSY`
` 764:              // *`
` 765:   `
` 766:              #endregion`
` 767:   `
` 768:          }`
` 769:      }`
` 770:  }`