`   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:      /// DGGLSE solves the linear equality-constrained least squares (LSE)`
`  27:      /// problem:`
`  28:      /// `
`  29:      /// minimize || c - A*x ||_2   subject to   B*x = d`
`  30:      /// `
`  31:      /// where A is an M-by-N matrix, B is a P-by-N matrix, c is a given`
`  32:      /// M-vector, and d is a given P-vector. It is assumed that`
`  33:      /// P .LE. N .LE. M+P, and`
`  34:      /// `
`  35:      /// rank(B) = P and  rank( (A) ) = N.`
`  36:      /// ( (B) )`
`  37:      /// `
`  38:      /// These conditions ensure that the LSE problem has a unique solution,`
`  39:      /// which is obtained using a generalized RQ factorization of the`
`  40:      /// matrices (B, A) given by`
`  41:      /// `
`  42:      /// B = (0 R)*Q,   A = Z*T*Q.`
`  43:      /// `
`  44:      ///</summary>`
`  45:      public class DGGLSE`
`  46:      {`
`  47:      `
`  48:   `
`  49:          #region Dependencies`
`  50:          `
`  51:          DAXPY _daxpy; DCOPY _dcopy; DGEMV _dgemv; DGGRQF _dggrqf; DORMQR _dormqr; DORMRQ _dormrq; DTRMV _dtrmv; DTRTRS _dtrtrs; `
`  52:          XERBLA _xerbla;ILAENV _ilaenv; `
`  53:   `
`  54:          #endregion`
`  55:   `
`  56:   `
`  57:          #region Fields`
`  58:          `
`  59:          const double ONE = 1.0E+0; bool LQUERY = false; int LOPT = 0; int LWKMIN = 0; int LWKOPT = 0; int MN = 0; int NB = 0; `
`  60:          int NB1 = 0;int NB2 = 0; int NB3 = 0; int NB4 = 0; int NR = 0; `
`  61:   `
`  62:          #endregion`
`  63:   `
`  64:          public DGGLSE(DAXPY daxpy, DCOPY dcopy, DGEMV dgemv, DGGRQF dggrqf, DORMQR dormqr, DORMRQ dormrq, DTRMV dtrmv, DTRTRS dtrtrs, XERBLA xerbla, ILAENV ilaenv)`
`  65:          {`
`  66:      `
`  67:   `
`  68:              #region Set Dependencies`
`  69:              `
`  70:              this._daxpy = daxpy; this._dcopy = dcopy; this._dgemv = dgemv; this._dggrqf = dggrqf; this._dormqr = dormqr; `
`  71:              this._dormrq = dormrq;this._dtrmv = dtrmv; this._dtrtrs = dtrtrs; this._xerbla = xerbla; this._ilaenv = ilaenv; `
`  72:   `
`  73:              #endregion`
`  74:   `
`  75:          }`
`  76:      `
`  77:          public DGGLSE()`
`  78:          {`
`  79:      `
`  80:   `
`  81:              #region Dependencies (Initialization)`
`  82:              `
`  83:              DAXPY daxpy = new DAXPY();`
`  84:              DCOPY dcopy = new DCOPY();`
`  85:              LSAME lsame = new LSAME();`
`  86:              XERBLA xerbla = new XERBLA();`
`  87:              DLAMC3 dlamc3 = new DLAMC3();`
`  88:              DLAPY2 dlapy2 = new DLAPY2();`
`  89:              DNRM2 dnrm2 = new DNRM2();`
`  90:              DSCAL dscal = new DSCAL();`
`  91:              IEEECK ieeeck = new IEEECK();`
`  92:              IPARMQ iparmq = new IPARMQ();`
`  93:              DGEMV dgemv = new DGEMV(lsame, xerbla);`
`  94:              DGER dger = new DGER(xerbla);`
`  95:              DLARF dlarf = new DLARF(dgemv, dger, lsame);`
`  96:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);`
`  97:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);`
`  98:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);`
`  99:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);`
` 100:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);`
` 101:              DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);`
` 102:              DGEQR2 dgeqr2 = new DGEQR2(dlarf, dlarfg, xerbla);`
` 103:              DGEMM dgemm = new DGEMM(lsame, xerbla);`
` 104:              DTRMM dtrmm = new DTRMM(lsame, xerbla);`
` 105:              DLARFB dlarfb = new DLARFB(lsame, dcopy, dgemm, dtrmm);`
` 106:              DTRMV dtrmv = new DTRMV(lsame, xerbla);`
` 107:              DLARFT dlarft = new DLARFT(dgemv, dtrmv, lsame);`
` 108:              ILAENV ilaenv = new ILAENV(ieeeck, iparmq);`
` 109:              DGEQRF dgeqrf = new DGEQRF(dgeqr2, dlarfb, dlarft, xerbla, ilaenv);`
` 110:              DGERQ2 dgerq2 = new DGERQ2(dlarf, dlarfg, xerbla);`
` 111:              DGERQF dgerqf = new DGERQF(dgerq2, dlarfb, dlarft, xerbla, ilaenv);`
` 112:              DORMR2 dormr2 = new DORMR2(lsame, dlarf, xerbla);`
` 113:              DORMRQ dormrq = new DORMRQ(lsame, ilaenv, dlarfb, dlarft, dormr2, xerbla);`
` 114:              DGGRQF dggrqf = new DGGRQF(dgeqrf, dgerqf, dormrq, xerbla, ilaenv);`
` 115:              DORM2R dorm2r = new DORM2R(lsame, dlarf, xerbla);`
` 116:              DORMQR dormqr = new DORMQR(lsame, ilaenv, dlarfb, dlarft, dorm2r, xerbla);`
` 117:              DTRSM dtrsm = new DTRSM(lsame, xerbla);`
` 118:              DTRTRS dtrtrs = new DTRTRS(lsame, dtrsm, xerbla);`
` 119:   `
` 120:              #endregion`
` 121:   `
` 122:   `
` 123:              #region Set Dependencies`
` 124:              `
` 125:              this._daxpy = daxpy; this._dcopy = dcopy; this._dgemv = dgemv; this._dggrqf = dggrqf; this._dormqr = dormqr; `
` 126:              this._dormrq = dormrq;this._dtrmv = dtrmv; this._dtrtrs = dtrtrs; this._xerbla = xerbla; this._ilaenv = ilaenv; `
` 127:   `
` 128:              #endregion`
` 129:   `
` 130:          }`
` 131:          /// <summary>`
` 132:          /// Purpose`
` 133:          /// =======`
` 134:          /// `
` 135:          /// DGGLSE solves the linear equality-constrained least squares (LSE)`
` 136:          /// problem:`
` 137:          /// `
` 138:          /// minimize || c - A*x ||_2   subject to   B*x = d`
` 139:          /// `
` 140:          /// where A is an M-by-N matrix, B is a P-by-N matrix, c is a given`
` 141:          /// M-vector, and d is a given P-vector. It is assumed that`
` 142:          /// P .LE. N .LE. M+P, and`
` 143:          /// `
` 144:          /// rank(B) = P and  rank( (A) ) = N.`
` 145:          /// ( (B) )`
` 146:          /// `
` 147:          /// These conditions ensure that the LSE problem has a unique solution,`
` 148:          /// which is obtained using a generalized RQ factorization of the`
` 149:          /// matrices (B, A) given by`
` 150:          /// `
` 151:          /// B = (0 R)*Q,   A = Z*T*Q.`
` 152:          /// `
` 153:          ///</summary>`
` 154:          /// <param name="M">`
` 155:          /// (input) INTEGER`
` 156:          /// The number of rows of the matrix A.  M .GE. 0.`
` 157:          ///</param>`
` 158:          /// <param name="N">`
` 159:          /// (input) INTEGER`
` 160:          /// The number of columns of the matrices A and B. N .GE. 0.`
` 161:          ///</param>`
` 162:          /// <param name="P">`
` 163:          /// (input) INTEGER`
` 164:          /// The number of rows of the matrix B. 0 .LE. P .LE. N .LE. M+P.`
` 165:          ///</param>`
` 166:          /// <param name="A">`
` 167:          /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)`
` 168:          /// On entry, the M-by-N matrix A.`
` 169:          /// On exit, the elements on and above the diagonal of the array`
` 170:          /// contain the min(M,N)-by-N upper trapezoidal matrix T.`
` 171:          ///</param>`
` 172:          /// <param name="LDA">`
` 173:          /// (input) INTEGER`
` 174:          /// The leading dimension of the array A. LDA .GE. max(1,M).`
` 175:          ///</param>`
` 176:          /// <param name="B">`
` 177:          /// = (0 R)*Q,   A = Z*T*Q.`
` 178:          ///</param>`
` 179:          /// <param name="LDB">`
` 180:          /// (input) INTEGER`
` 181:          /// The leading dimension of the array B. LDB .GE. max(1,P).`
` 182:          ///</param>`
` 183:          /// <param name="C">`
` 184:          /// (input/output) DOUBLE PRECISION array, dimension (M)`
` 185:          /// On entry, C contains the right hand side vector for the`
` 186:          /// least squares part of the LSE problem.`
` 187:          /// On exit, the residual sum of squares for the solution`
` 188:          /// is given by the sum of squares of elements N-P+1 to M of`
` 189:          /// vector C.`
` 190:          ///</param>`
` 191:          /// <param name="D">`
` 192:          /// (input/output) DOUBLE PRECISION array, dimension (P)`
` 193:          /// On entry, D contains the right hand side vector for the`
` 194:          /// constrained equation.`
` 195:          /// On exit, D is destroyed.`
` 196:          ///</param>`
` 197:          /// <param name="X">`
` 198:          /// (output) DOUBLE PRECISION array, dimension (N)`
` 199:          /// On exit, X is the solution of the LSE problem.`
` 200:          ///</param>`
` 201:          /// <param name="WORK">`
` 202:          /// (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))`
` 203:          /// On exit, if INFO = 0, WORK(1) returns the optimal LWORK.`
` 204:          ///</param>`
` 205:          /// <param name="LWORK">`
` 206:          /// (input) INTEGER`
` 207:          /// The dimension of the array WORK. LWORK .GE. max(1,M+N+P).`
` 208:          /// For optimum performance LWORK .GE. P+min(M,N)+max(M,N)*NB,`
` 209:          /// where NB is an upper bound for the optimal blocksizes for`
` 210:          /// DGEQRF, SGERQF, DORMQR and SORMRQ.`
` 211:          /// `
` 212:          /// If LWORK = -1, then a workspace query is assumed; the routine`
` 213:          /// only calculates the optimal size of the WORK array, returns`
` 214:          /// this value as the first entry of the WORK array, and no error`
` 215:          /// message related to LWORK is issued by XERBLA.`
` 216:          ///</param>`
` 217:          /// <param name="INFO">`
` 218:          /// (output) INTEGER`
` 219:          /// = 0:  successful exit.`
` 220:          /// .LT. 0:  if INFO = -i, the i-th argument had an illegal value.`
` 221:          /// = 1:  the upper triangular factor R associated with B in the`
` 222:          /// generalized RQ factorization of the pair (B, A) is`
` 223:          /// singular, so that rank(B) .LT. P; the least squares`
` 224:          /// solution could not be computed.`
` 225:          /// = 2:  the (N-P) by (N-P) part of the upper trapezoidal factor`
` 226:          /// T associated with A in the generalized RQ factorization`
` 227:          /// of the pair (B, A) is singular, so that`
` 228:          /// rank( (A) ) .LT. N; the least squares solution could not`
` 229:          /// ( (B) )`
` 230:          /// be computed.`
` 231:          ///</param>`
` 232:          public void Run(int M, int N, int P, ref double[] A, int offset_a, int LDA, ref double[] B, int offset_b`
` 233:                           , int LDB, ref double[] C, int offset_c, ref double[] D, int offset_d, ref double[] X, int offset_x, ref double[] WORK, int offset_work, int LWORK`
` 234:                           , ref int INFO)`
` 235:          {`
` 236:   `
` 237:              #region Array Index Correction`
` 238:              `
` 239:               int o_a = -1 - LDA + offset_a;  int o_b = -1 - LDB + offset_b;  int o_c = -1 + offset_c;  int o_d = -1 + offset_d; `
` 240:               int o_x = -1 + offset_x; int o_work = -1 + offset_work; `
` 241:   `
` 242:              #endregion`
` 243:   `
` 244:   `
` 245:              #region Prolog`
` 246:              `
` 247:              // *`
` 248:              // *  -- LAPACK driver routine (version 3.1) --`
` 249:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..`
` 250:              // *     November 2006`
` 251:              // *`
` 252:              // *     .. Scalar Arguments ..`
` 253:              // *     ..`
` 254:              // *     .. Array Arguments ..`
` 255:              // *     ..`
` 256:              // *`
` 257:              // *  Purpose`
` 258:              // *  =======`
` 259:              // *`
` 260:              // *  DGGLSE solves the linear equality-constrained least squares (LSE)`
` 261:              // *  problem:`
` 262:              // *`
` 263:              // *          minimize || c - A*x ||_2   subject to   B*x = d`
` 264:              // *`
` 265:              // *  where A is an M-by-N matrix, B is a P-by-N matrix, c is a given`
` 266:              // *  M-vector, and d is a given P-vector. It is assumed that`
` 267:              // *  P <= N <= M+P, and`
` 268:              // *`
` 269:              // *           rank(B) = P and  rank( (A) ) = N.`
` 270:              // *                                ( (B) )`
` 271:              // *`
` 272:              // *  These conditions ensure that the LSE problem has a unique solution,`
` 273:              // *  which is obtained using a generalized RQ factorization of the`
` 274:              // *  matrices (B, A) given by`
` 275:              // *`
` 276:              // *     B = (0 R)*Q,   A = Z*T*Q.`
` 277:              // *`
` 278:              // *  Arguments`
` 279:              // *  =========`
` 280:              // *`
` 281:              // *  M       (input) INTEGER`
` 282:              // *          The number of rows of the matrix A.  M >= 0.`
` 283:              // *`
` 284:              // *  N       (input) INTEGER`
` 285:              // *          The number of columns of the matrices A and B. N >= 0.`
` 286:              // *`
` 287:              // *  P       (input) INTEGER`
` 288:              // *          The number of rows of the matrix B. 0 <= P <= N <= M+P.`
` 289:              // *`
` 290:              // *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)`
` 291:              // *          On entry, the M-by-N matrix A.`
` 292:              // *          On exit, the elements on and above the diagonal of the array`
` 293:              // *          contain the min(M,N)-by-N upper trapezoidal matrix T.`
` 294:              // *`
` 295:              // *  LDA     (input) INTEGER`
` 296:              // *          The leading dimension of the array A. LDA >= max(1,M).`
` 297:              // *`
` 298:              // *  B       (input/output) DOUBLE PRECISION array, dimension (LDB,N)`
` 299:              // *          On entry, the P-by-N matrix B.`
` 300:              // *          On exit, the upper triangle of the subarray B(1:P,N-P+1:N)`
` 301:              // *          contains the P-by-P upper triangular matrix R.`
` 302:              // *`
` 303:              // *  LDB     (input) INTEGER`
` 304:              // *          The leading dimension of the array B. LDB >= max(1,P).`
` 305:              // *`
` 306:              // *  C       (input/output) DOUBLE PRECISION array, dimension (M)`
` 307:              // *          On entry, C contains the right hand side vector for the`
` 308:              // *          least squares part of the LSE problem.`
` 309:              // *          On exit, the residual sum of squares for the solution`
` 310:              // *          is given by the sum of squares of elements N-P+1 to M of`
` 311:              // *          vector C.`
` 312:              // *`
` 313:              // *  D       (input/output) DOUBLE PRECISION array, dimension (P)`
` 314:              // *          On entry, D contains the right hand side vector for the`
` 315:              // *          constrained equation.`
` 316:              // *          On exit, D is destroyed.`
` 317:              // *`
` 318:              // *  X       (output) DOUBLE PRECISION array, dimension (N)`
` 319:              // *          On exit, X is the solution of the LSE problem.`
` 320:              // *`
` 321:              // *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))`
` 322:              // *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.`
` 323:              // *`
` 324:              // *  LWORK   (input) INTEGER`
` 325:              // *          The dimension of the array WORK. LWORK >= max(1,M+N+P).`
` 326:              // *          For optimum performance LWORK >= P+min(M,N)+max(M,N)*NB,`
` 327:              // *          where NB is an upper bound for the optimal blocksizes for`
` 328:              // *          DGEQRF, SGERQF, DORMQR and SORMRQ.`
` 329:              // *`
` 330:              // *          If LWORK = -1, then a workspace query is assumed; the routine`
` 331:              // *          only calculates the optimal size of the WORK array, returns`
` 332:              // *          this value as the first entry of the WORK array, and no error`
` 333:              // *          message related to LWORK is issued by XERBLA.`
` 334:              // *`
` 335:              // *  INFO    (output) INTEGER`
` 336:              // *          = 0:  successful exit.`
` 337:              // *          < 0:  if INFO = -i, the i-th argument had an illegal value.`
` 338:              // *          = 1:  the upper triangular factor R associated with B in the`
` 339:              // *                generalized RQ factorization of the pair (B, A) is`
` 340:              // *                singular, so that rank(B) < P; the least squares`
` 341:              // *                solution could not be computed.`
` 342:              // *          = 2:  the (N-P) by (N-P) part of the upper trapezoidal factor`
` 343:              // *                T associated with A in the generalized RQ factorization`
` 344:              // *                of the pair (B, A) is singular, so that`
` 345:              // *                rank( (A) ) < N; the least squares solution could not`
` 346:              // *                    ( (B) )`
` 347:              // *                be computed.`
` 348:              // *`
` 349:              // *  =====================================================================`
` 350:              // *`
` 351:              // *     .. Parameters ..`
` 352:              // *     ..`
` 353:              // *     .. Local Scalars ..`
` 354:              // *     ..`
` 355:              // *     .. External Subroutines ..`
` 356:              // *     ..`
` 357:              // *     .. External Functions ..`
` 358:              // *     ..`
` 359:              // *     .. Intrinsic Functions ..`
` 360:              //      INTRINSIC          INT, MAX, MIN;`
` 361:              // *     ..`
` 362:              // *     .. Executable Statements ..`
` 363:              // *`
` 364:              // *     Test the input parameters`
` 365:              // *`
` 366:   `
` 367:              #endregion`
` 368:   `
` 369:   `
` 370:              #region Body`
` 371:              `
` 372:              INFO = 0;`
` 373:              MN = Math.Min(M, N);`
` 374:              LQUERY = (LWORK ==  - 1);`
` 375:              if (M < 0)`
` 376:              {`
` 377:                  INFO =  - 1;`
` 378:              }`
` 379:              else`
` 380:              {`
` 381:                  if (N < 0)`
` 382:                  {`
` 383:                      INFO =  - 2;`
` 384:                  }`
` 385:                  else`
` 386:                  {`
` 387:                      if (P < 0 || P > N || P < N - M)`
` 388:                      {`
` 389:                          INFO =  - 3;`
` 390:                      }`
` 391:                      else`
` 392:                      {`
` 393:                          if (LDA < Math.Max(1, M))`
` 394:                          {`
` 395:                              INFO =  - 5;`
` 396:                          }`
` 397:                          else`
` 398:                          {`
` 399:                              if (LDB < Math.Max(1, P))`
` 400:                              {`
` 401:                                  INFO =  - 7;`
` 402:                              }`
` 403:                          }`
` 404:                      }`
` 405:                  }`
` 406:              }`
` 407:              // *`
` 408:              // *     Calculate workspace`
` 409:              // *`
` 410:              if (INFO == 0)`
` 411:              {`
` 412:                  if (N == 0)`
` 413:                  {`
` 414:                      LWKMIN = 1;`
` 415:                      LWKOPT = 1;`
` 416:                  }`
` 417:                  else`
` 418:                  {`
` 419:                      NB1 = this._ilaenv.Run(1, "DGEQRF", " ", M, N,  - 1,  - 1);`
` 420:                      NB2 = this._ilaenv.Run(1, "DGERQF", " ", M, N,  - 1,  - 1);`
` 421:                      NB3 = this._ilaenv.Run(1, "DORMQR", " ", M, N, P,  - 1);`
` 422:                      NB4 = this._ilaenv.Run(1, "DORMRQ", " ", M, N, P,  - 1);`
` 423:                      NB = Math.Max(NB1, Math.Max(NB2, Math.Max(NB3, NB4)));`
` 424:                      LWKMIN = M + N + P;`
` 425:                      LWKOPT = P + MN + Math.Max(M, N) * NB;`
` 426:                  }`
` 427:                  WORK[1 + o_work] = LWKOPT;`
` 428:                  // *`
` 429:                  if (LWORK < LWKMIN && !LQUERY)`
` 430:                  {`
` 431:                      INFO =  - 12;`
` 432:                  }`
` 433:              }`
` 434:              // *`
` 435:              if (INFO != 0)`
` 436:              {`
` 437:                  this._xerbla.Run("DGGLSE",  - INFO);`
` 438:                  return;`
` 439:              }`
` 440:              else`
` 441:              {`
` 442:                  if (LQUERY)`
` 443:                  {`
` 444:                      return;`
` 445:                  }`
` 446:              }`
` 447:              // *`
` 448:              // *     Quick return if possible`
` 449:              // *`
` 450:              if (N == 0) return;`
` 451:              // *`
` 452:              // *     Compute the GRQ factorization of matrices B and A:`
` 453:              // *`
` 454:              // *            B*Q' = (  0  T12 ) P   Z'*A*Q' = ( R11 R12 ) N-P`
` 455:              // *                     N-P  P                  (  0  R22 ) M+P-N`
` 456:              // *                                               N-P  P`
` 457:              // *`
` 458:              // *     where T12 and R11 are upper triangular, and Q and Z are`
` 459:              // *     orthogonal.`
` 460:              // *`
` 461:              this._dggrqf.Run(P, M, N, ref B, offset_b, LDB, ref WORK, offset_work`
` 462:                               , ref A, offset_a, LDA, ref WORK, P + 1 + o_work, ref WORK, P + MN + 1 + o_work, LWORK - P - MN, ref INFO);`
` 463:              LOPT = (int)WORK[P + MN + 1 + o_work];`
` 464:              // *`
` 465:              // *     Update c = Z'*c = ( c1 ) N-P`
` 466:              // *                       ( c2 ) M+P-N`
` 467:              // *`
` 468:              this._dormqr.Run("Left", "Transpose", M, 1, MN, ref A, offset_a`
` 469:                               , LDA, WORK, P + 1 + o_work, ref C, offset_c, Math.Max(1, M), ref WORK, P + MN + 1 + o_work, LWORK - P - MN`
` 470:                               , ref INFO);`
` 471:              LOPT = Math.Max(LOPT, Convert.ToInt32(Math.Truncate(WORK[P + MN + 1 + o_work])));`
` 472:              // *`
` 473:              // *     Solve T12*x2 = d for x2`
` 474:              // *`
` 475:              if (P > 0)`
` 476:              {`
` 477:                  this._dtrtrs.Run("Upper", "No transpose", "Non-unit", P, 1, B, 1+(N - P + 1) * LDB + o_b`
` 478:                                   , LDB, ref D, offset_d, P, ref INFO);`
` 479:                  // *`
` 480:                  if (INFO > 0)`
` 481:                  {`
` 482:                      INFO = 1;`
` 483:                      return;`
` 484:                  }`
` 485:                  // *`
` 486:                  // *        Put the solution in X`
` 487:                  // *`
` 488:                  this._dcopy.Run(P, D, offset_d, 1, ref X, N - P + 1 + o_x, 1);`
` 489:                  // *`
` 490:                  // *        Update c1`
` 491:                  // *`
` 492:                  this._dgemv.Run("No transpose", N - P, P,  - ONE, A, 1+(N - P + 1) * LDA + o_a, LDA`
` 493:                                  , D, offset_d, 1, ONE, ref C, offset_c, 1);`
` 494:              }`
` 495:              // *`
` 496:              // *     Solve R11*x1 = c1 for x1`
` 497:              // *`
` 498:              if (N > P)`
` 499:              {`
` 500:                  this._dtrtrs.Run("Upper", "No transpose", "Non-unit", N - P, 1, A, offset_a`
` 501:                                   , LDA, ref C, offset_c, N - P, ref INFO);`
` 502:                  // *`
` 503:                  if (INFO > 0)`
` 504:                  {`
` 505:                      INFO = 2;`
` 506:                      return;`
` 507:                  }`
` 508:                  // *`
` 509:                  // *        Put the solutions in X`
` 510:                  // *`
` 511:                  this._dcopy.Run(N - P, C, offset_c, 1, ref X, offset_x, 1);`
` 512:              }`
` 513:              // *`
` 514:              // *     Compute the residual vector:`
` 515:              // *`
` 516:              if (M < N)`
` 517:              {`
` 518:                  NR = M + P - N;`
` 519:                  if (NR > 0)`
` 520:                  {`
` 521:                      this._dgemv.Run("No transpose", NR, N - M,  - ONE, A, N - P + 1+(M + 1) * LDA + o_a, LDA`
` 522:                                      , D, NR + 1 + o_d, 1, ONE, ref C, N - P + 1 + o_c, 1);`
` 523:                  }`
` 524:              }`
` 525:              else`
` 526:              {`
` 527:                  NR = P;`
` 528:              }`
` 529:              if (NR > 0)`
` 530:              {`
` 531:                  this._dtrmv.Run("Upper", "No transpose", "Non unit", NR, A, N - P + 1+(N - P + 1) * LDA + o_a, LDA`
` 532:                                  , ref D, offset_d, 1);`
` 533:                  this._daxpy.Run(NR,  - ONE, D, offset_d, 1, ref C, N - P + 1 + o_c, 1);`
` 534:              }`
` 535:              // *`
` 536:              // *     Backward transformation x = Q'*x`
` 537:              // *`
` 538:              this._dormrq.Run("Left", "Transpose", N, 1, P, ref B, offset_b`
` 539:                               , LDB, WORK, 1 + o_work, ref X, offset_x, N, ref WORK, P + MN + 1 + o_work, LWORK - P - MN`
` 540:                               , ref INFO);`
` 541:              WORK[1 + o_work] = P + MN + Math.Max(LOPT, Convert.ToInt32(Math.Truncate(WORK[P + MN + 1 + o_work])));`
` 542:              // *`
` 543:              return;`
` 544:              // *`
` 545:              // *     End of DGGLSE`
` 546:              // *`
` 547:   `
` 548:              #endregion`
` 549:   `
` 550:          }`
` 551:      }`
` 552:  }`