`   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:      /// DLABRD reduces the first NB rows and columns of a real general`
`  27:      /// m by n matrix A to upper or lower bidiagonal form by an orthogonal`
`  28:      /// transformation Q' * A * P, and returns the matrices X and Y which`
`  29:      /// are needed to apply the transformation to the unreduced part of A.`
`  30:      /// `
`  31:      /// If m .GE. n, A is reduced to upper bidiagonal form; if m .LT. n, to lower`
`  32:      /// bidiagonal form.`
`  33:      /// `
`  34:      /// This is an auxiliary routine called by DGEBRD`
`  35:      /// `
`  36:      ///</summary>`
`  37:      public class DLABRD`
`  38:      {`
`  39:      `
`  40:   `
`  41:          #region Dependencies`
`  42:          `
`  43:          DGEMV _dgemv; DLARFG _dlarfg; DSCAL _dscal; `
`  44:   `
`  45:          #endregion`
`  46:   `
`  47:   `
`  48:          #region Fields`
`  49:          `
`  50:          const double ZERO = 0.0E0; const double ONE = 1.0E0; int I = 0; `
`  51:   `
`  52:          #endregion`
`  53:   `
`  54:          public DLABRD(DGEMV dgemv, DLARFG dlarfg, DSCAL dscal)`
`  55:          {`
`  56:      `
`  57:   `
`  58:              #region Set Dependencies`
`  59:              `
`  60:              this._dgemv = dgemv; this._dlarfg = dlarfg; this._dscal = dscal; `
`  61:   `
`  62:              #endregion`
`  63:   `
`  64:          }`
`  65:      `
`  66:          public DLABRD()`
`  67:          {`
`  68:      `
`  69:   `
`  70:              #region Dependencies (Initialization)`
`  71:              `
`  72:              LSAME lsame = new LSAME();`
`  73:              XERBLA xerbla = new XERBLA();`
`  74:              DLAMC3 dlamc3 = new DLAMC3();`
`  75:              DLAPY2 dlapy2 = new DLAPY2();`
`  76:              DNRM2 dnrm2 = new DNRM2();`
`  77:              DSCAL dscal = new DSCAL();`
`  78:              DGEMV dgemv = new DGEMV(lsame, xerbla);`
`  79:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);`
`  80:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);`
`  81:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);`
`  82:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);`
`  83:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);`
`  84:              DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);`
`  85:   `
`  86:              #endregion`
`  87:   `
`  88:   `
`  89:              #region Set Dependencies`
`  90:              `
`  91:              this._dgemv = dgemv; this._dlarfg = dlarfg; this._dscal = dscal; `
`  92:   `
`  93:              #endregion`
`  94:   `
`  95:          }`
`  96:          /// <summary>`
`  97:          /// Purpose`
`  98:          /// =======`
`  99:          /// `
` 100:          /// DLABRD reduces the first NB rows and columns of a real general`
` 101:          /// m by n matrix A to upper or lower bidiagonal form by an orthogonal`
` 102:          /// transformation Q' * A * P, and returns the matrices X and Y which`
` 103:          /// are needed to apply the transformation to the unreduced part of A.`
` 104:          /// `
` 105:          /// If m .GE. n, A is reduced to upper bidiagonal form; if m .LT. n, to lower`
` 106:          /// bidiagonal form.`
` 107:          /// `
` 108:          /// This is an auxiliary routine called by DGEBRD`
` 109:          /// `
` 110:          ///</summary>`
` 111:          /// <param name="M">`
` 112:          /// (input) INTEGER`
` 113:          /// The number of rows in the matrix A.`
` 114:          ///</param>`
` 115:          /// <param name="N">`
` 116:          /// (input) INTEGER`
` 117:          /// The number of columns in the matrix A.`
` 118:          ///</param>`
` 119:          /// <param name="NB">`
` 120:          /// (input) INTEGER`
` 121:          /// The number of leading rows and columns of A to be reduced.`
` 122:          ///</param>`
` 123:          /// <param name="A">`
` 124:          /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)`
` 125:          /// On entry, the m by n general matrix to be reduced.`
` 126:          /// On exit, the first NB rows and columns of the matrix are`
` 127:          /// overwritten; the rest of the array is unchanged.`
` 128:          /// If m .GE. n, elements on and below the diagonal in the first NB`
` 129:          /// columns, with the array TAUQ, represent the orthogonal`
` 130:          /// matrix Q as a product of elementary reflectors; and`
` 131:          /// elements above the diagonal in the first NB rows, with the`
` 132:          /// array TAUP, represent the orthogonal matrix P as a product`
` 133:          /// of elementary reflectors.`
` 134:          /// If m .LT. n, elements below the diagonal in the first NB`
` 135:          /// columns, with the array TAUQ, represent the orthogonal`
` 136:          /// matrix Q as a product of elementary reflectors, and`
` 137:          /// elements on and above the diagonal in the first NB rows,`
` 138:          /// with the array TAUP, represent the orthogonal matrix P as`
` 139:          /// a product of elementary reflectors.`
` 140:          /// See Further Details.`
` 141:          ///</param>`
` 142:          /// <param name="LDA">`
` 143:          /// (input) INTEGER`
` 144:          /// The leading dimension of the array A.  LDA .GE. max(1,M).`
` 145:          ///</param>`
` 146:          /// <param name="D">`
` 147:          /// (output) DOUBLE PRECISION array, dimension (NB)`
` 148:          /// The diagonal elements of the first NB rows and columns of`
` 149:          /// the reduced matrix.  D(i) = A(i,i).`
` 150:          ///</param>`
` 151:          /// <param name="E">`
` 152:          /// (output) DOUBLE PRECISION array, dimension (NB)`
` 153:          /// The off-diagonal elements of the first NB rows and columns of`
` 154:          /// the reduced matrix.`
` 155:          ///</param>`
` 156:          /// <param name="TAUQ">`
` 157:          /// (output) DOUBLE PRECISION array dimension (NB)`
` 158:          /// The scalar factors of the elementary reflectors which`
` 159:          /// represent the orthogonal matrix Q. See Further Details.`
` 160:          ///</param>`
` 161:          /// <param name="TAUP">`
` 162:          /// (output) DOUBLE PRECISION array, dimension (NB)`
` 163:          /// The scalar factors of the elementary reflectors which`
` 164:          /// represent the orthogonal matrix P. See Further Details.`
` 165:          ///</param>`
` 166:          /// <param name="X">`
` 167:          /// (output) DOUBLE PRECISION array, dimension (LDX,NB)`
` 168:          /// The m-by-nb matrix X required to update the unreduced part`
` 169:          /// of A.`
` 170:          ///</param>`
` 171:          /// <param name="LDX">`
` 172:          /// (input) INTEGER`
` 173:          /// The leading dimension of the array X. LDX .GE. M.`
` 174:          ///</param>`
` 175:          /// <param name="Y">`
` 176:          /// (output) DOUBLE PRECISION array, dimension (LDY,NB)`
` 177:          /// The n-by-nb matrix Y required to update the unreduced part`
` 178:          /// of A.`
` 179:          ///</param>`
` 180:          /// <param name="LDY">`
` 181:          /// (input) INTEGER`
` 182:          /// The leading dimension of the array Y. LDY .GE. N.`
` 183:          ///</param>`
` 184:          public void Run(int M, int N, int NB, ref double[] A, int offset_a, int LDA, ref double[] D, int offset_d`
` 185:                           , ref double[] E, int offset_e, ref double[] TAUQ, int offset_tauq, ref double[] TAUP, int offset_taup, ref double[] X, int offset_x, int LDX, ref double[] Y, int offset_y`
` 186:                           , int LDY)`
` 187:          {`
` 188:   `
` 189:              #region Array Index Correction`
` 190:              `
` 191:               int o_a = -1 - LDA + offset_a;  int o_d = -1 + offset_d;  int o_e = -1 + offset_e;  int o_tauq = -1 + offset_tauq; `
` 192:               int o_taup = -1 + offset_taup; int o_x = -1 - LDX + offset_x;  int o_y = -1 - LDY + offset_y; `
` 193:   `
` 194:              #endregion`
` 195:   `
` 196:   `
` 197:              #region Prolog`
` 198:              `
` 199:              // *`
` 200:              // *  -- LAPACK auxiliary routine (version 3.1) --`
` 201:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..`
` 202:              // *     November 2006`
` 203:              // *`
` 204:              // *     .. Scalar Arguments ..`
` 205:              // *     ..`
` 206:              // *     .. Array Arguments ..`
` 207:              // *     ..`
` 208:              // *`
` 209:              // *  Purpose`
` 210:              // *  =======`
` 211:              // *`
` 212:              // *  DLABRD reduces the first NB rows and columns of a real general`
` 213:              // *  m by n matrix A to upper or lower bidiagonal form by an orthogonal`
` 214:              // *  transformation Q' * A * P, and returns the matrices X and Y which`
` 215:              // *  are needed to apply the transformation to the unreduced part of A.`
` 216:              // *`
` 217:              // *  If m >= n, A is reduced to upper bidiagonal form; if m < n, to lower`
` 218:              // *  bidiagonal form.`
` 219:              // *`
` 220:              // *  This is an auxiliary routine called by DGEBRD`
` 221:              // *`
` 222:              // *  Arguments`
` 223:              // *  =========`
` 224:              // *`
` 225:              // *  M       (input) INTEGER`
` 226:              // *          The number of rows in the matrix A.`
` 227:              // *`
` 228:              // *  N       (input) INTEGER`
` 229:              // *          The number of columns in the matrix A.`
` 230:              // *`
` 231:              // *  NB      (input) INTEGER`
` 232:              // *          The number of leading rows and columns of A to be reduced.`
` 233:              // *`
` 234:              // *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)`
` 235:              // *          On entry, the m by n general matrix to be reduced.`
` 236:              // *          On exit, the first NB rows and columns of the matrix are`
` 237:              // *          overwritten; the rest of the array is unchanged.`
` 238:              // *          If m >= n, elements on and below the diagonal in the first NB`
` 239:              // *            columns, with the array TAUQ, represent the orthogonal`
` 240:              // *            matrix Q as a product of elementary reflectors; and`
` 241:              // *            elements above the diagonal in the first NB rows, with the`
` 242:              // *            array TAUP, represent the orthogonal matrix P as a product`
` 243:              // *            of elementary reflectors.`
` 244:              // *          If m < n, elements below the diagonal in the first NB`
` 245:              // *            columns, with the array TAUQ, represent the orthogonal`
` 246:              // *            matrix Q as a product of elementary reflectors, and`
` 247:              // *            elements on and above the diagonal in the first NB rows,`
` 248:              // *            with the array TAUP, represent the orthogonal matrix P as`
` 249:              // *            a product of elementary reflectors.`
` 250:              // *          See Further Details.`
` 251:              // *`
` 252:              // *  LDA     (input) INTEGER`
` 253:              // *          The leading dimension of the array A.  LDA >= max(1,M).`
` 254:              // *`
` 255:              // *  D       (output) DOUBLE PRECISION array, dimension (NB)`
` 256:              // *          The diagonal elements of the first NB rows and columns of`
` 257:              // *          the reduced matrix.  D(i) = A(i,i).`
` 258:              // *`
` 259:              // *  E       (output) DOUBLE PRECISION array, dimension (NB)`
` 260:              // *          The off-diagonal elements of the first NB rows and columns of`
` 261:              // *          the reduced matrix.`
` 262:              // *`
` 263:              // *  TAUQ    (output) DOUBLE PRECISION array dimension (NB)`
` 264:              // *          The scalar factors of the elementary reflectors which`
` 265:              // *          represent the orthogonal matrix Q. See Further Details.`
` 266:              // *`
` 267:              // *  TAUP    (output) DOUBLE PRECISION array, dimension (NB)`
` 268:              // *          The scalar factors of the elementary reflectors which`
` 269:              // *          represent the orthogonal matrix P. See Further Details.`
` 270:              // *`
` 271:              // *  X       (output) DOUBLE PRECISION array, dimension (LDX,NB)`
` 272:              // *          The m-by-nb matrix X required to update the unreduced part`
` 273:              // *          of A.`
` 274:              // *`
` 275:              // *  LDX     (input) INTEGER`
` 276:              // *          The leading dimension of the array X. LDX >= M.`
` 277:              // *`
` 278:              // *  Y       (output) DOUBLE PRECISION array, dimension (LDY,NB)`
` 279:              // *          The n-by-nb matrix Y required to update the unreduced part`
` 280:              // *          of A.`
` 281:              // *`
` 282:              // *  LDY     (input) INTEGER`
` 283:              // *          The leading dimension of the array Y. LDY >= N.`
` 284:              // *`
` 285:              // *  Further Details`
` 286:              // *  ===============`
` 287:              // *`
` 288:              // *  The matrices Q and P are represented as products of elementary`
` 289:              // *  reflectors:`
` 290:              // *`
` 291:              // *     Q = H(1) H(2) . . . H(nb)  and  P = G(1) G(2) . . . G(nb)`
` 292:              // *`
` 293:              // *  Each H(i) and G(i) has the form:`
` 294:              // *`
` 295:              // *     H(i) = I - tauq * v * v'  and G(i) = I - taup * u * u'`
` 296:              // *`
` 297:              // *  where tauq and taup are real scalars, and v and u are real vectors.`
` 298:              // *`
` 299:              // *  If m >= n, v(1:i-1) = 0, v(i) = 1, and v(i:m) is stored on exit in`
` 300:              // *  A(i:m,i); u(1:i) = 0, u(i+1) = 1, and u(i+1:n) is stored on exit in`
` 301:              // *  A(i,i+1:n); tauq is stored in TAUQ(i) and taup in TAUP(i).`
` 302:              // *`
` 303:              // *  If m < n, v(1:i) = 0, v(i+1) = 1, and v(i+1:m) is stored on exit in`
` 304:              // *  A(i+2:m,i); u(1:i-1) = 0, u(i) = 1, and u(i:n) is stored on exit in`
` 305:              // *  A(i,i+1:n); tauq is stored in TAUQ(i) and taup in TAUP(i).`
` 306:              // *`
` 307:              // *  The elements of the vectors v and u together form the m-by-nb matrix`
` 308:              // *  V and the nb-by-n matrix U' which are needed, with X and Y, to apply`
` 309:              // *  the transformation to the unreduced part of the matrix, using a block`
` 310:              // *  update of the form:  A := A - V*Y' - X*U'.`
` 311:              // *`
` 312:              // *  The contents of A on exit are illustrated by the following examples`
` 313:              // *  with nb = 2:`
` 314:              // *`
` 315:              // *  m = 6 and n = 5 (m > n):          m = 5 and n = 6 (m < n):`
` 316:              // *`
` 317:              // *    (  1   1   u1  u1  u1 )           (  1   u1  u1  u1  u1  u1 )`
` 318:              // *    (  v1  1   1   u2  u2 )           (  1   1   u2  u2  u2  u2 )`
` 319:              // *    (  v1  v2  a   a   a  )           (  v1  1   a   a   a   a  )`
` 320:              // *    (  v1  v2  a   a   a  )           (  v1  v2  a   a   a   a  )`
` 321:              // *    (  v1  v2  a   a   a  )           (  v1  v2  a   a   a   a  )`
` 322:              // *    (  v1  v2  a   a   a  )`
` 323:              // *`
` 324:              // *  where a denotes an element of the original matrix which is unchanged,`
` 325:              // *  vi denotes an element of the vector defining H(i), and ui an element`
` 326:              // *  of the vector defining G(i).`
` 327:              // *`
` 328:              // *  =====================================================================`
` 329:              // *`
` 330:              // *     .. Parameters ..`
` 331:              // *     ..`
` 332:              // *     .. Local Scalars ..`
` 333:              // *     ..`
` 334:              // *     .. External Subroutines ..`
` 335:              // *     ..`
` 336:              // *     .. Intrinsic Functions ..`
` 337:              //      INTRINSIC          MIN;`
` 338:              // *     ..`
` 339:              // *     .. Executable Statements ..`
` 340:              // *`
` 341:              // *     Quick return if possible`
` 342:              // *`
` 343:   `
` 344:              #endregion`
` 345:   `
` 346:   `
` 347:              #region Body`
` 348:              `
` 349:              if (M <= 0 || N <= 0) return;`
` 350:              // *`
` 351:              if (M >= N)`
` 352:              {`
` 353:                  // *`
` 354:                  // *        Reduce to upper bidiagonal form`
` 355:                  // *`
` 356:                  for (I = 1; I <= NB; I++)`
` 357:                  {`
` 358:                      // *`
` 359:                      // *           Update A(i:m,i)`
` 360:                      // *`
` 361:                      this._dgemv.Run("No transpose", M - I + 1, I - 1,  - ONE, A, I+1 * LDA + o_a, LDA`
` 362:                                      , Y, I+1 * LDY + o_y, LDY, ONE, ref A, I+I * LDA + o_a, 1);`
` 363:                      this._dgemv.Run("No transpose", M - I + 1, I - 1,  - ONE, X, I+1 * LDX + o_x, LDX`
` 364:                                      , A, 1+I * LDA + o_a, 1, ONE, ref A, I+I * LDA + o_a, 1);`
` 365:                      // *`
` 366:                      // *           Generate reflection Q(i) to annihilate A(i+1:m,i)`
` 367:                      // *`
` 368:                      this._dlarfg.Run(M - I + 1, ref A[I+I * LDA + o_a], ref A, Math.Min(I + 1, M)+I * LDA + o_a, 1, ref TAUQ[I + o_tauq]);`
` 369:                      D[I + o_d] = A[I+I * LDA + o_a];`
` 370:                      if (I < N)`
` 371:                      {`
` 372:                          A[I+I * LDA + o_a] = ONE;`
` 373:                          // *`
` 374:                          // *              Compute Y(i+1:n,i)`
` 375:                          // *`
` 376:                          this._dgemv.Run("Transpose", M - I + 1, N - I, ONE, A, I+(I + 1) * LDA + o_a, LDA`
` 377:                                          , A, I+I * LDA + o_a, 1, ZERO, ref Y, I + 1+I * LDY + o_y, 1);`
` 378:                          this._dgemv.Run("Transpose", M - I + 1, I - 1, ONE, A, I+1 * LDA + o_a, LDA`
` 379:                                          , A, I+I * LDA + o_a, 1, ZERO, ref Y, 1+I * LDY + o_y, 1);`
` 380:                          this._dgemv.Run("No transpose", N - I, I - 1,  - ONE, Y, I + 1+1 * LDY + o_y, LDY`
` 381:                                          , Y, 1+I * LDY + o_y, 1, ONE, ref Y, I + 1+I * LDY + o_y, 1);`
` 382:                          this._dgemv.Run("Transpose", M - I + 1, I - 1, ONE, X, I+1 * LDX + o_x, LDX`
` 383:                                          , A, I+I * LDA + o_a, 1, ZERO, ref Y, 1+I * LDY + o_y, 1);`
` 384:                          this._dgemv.Run("Transpose", I - 1, N - I,  - ONE, A, 1+(I + 1) * LDA + o_a, LDA`
` 385:                                          , Y, 1+I * LDY + o_y, 1, ONE, ref Y, I + 1+I * LDY + o_y, 1);`
` 386:                          this._dscal.Run(N - I, TAUQ[I + o_tauq], ref Y, I + 1+I * LDY + o_y, 1);`
` 387:                          // *`
` 388:                          // *              Update A(i,i+1:n)`
` 389:                          // *`
` 390:                          this._dgemv.Run("No transpose", N - I, I,  - ONE, Y, I + 1+1 * LDY + o_y, LDY`
` 391:                                          , A, I+1 * LDA + o_a, LDA, ONE, ref A, I+(I + 1) * LDA + o_a, LDA);`
` 392:                          this._dgemv.Run("Transpose", I - 1, N - I,  - ONE, A, 1+(I + 1) * LDA + o_a, LDA`
` 393:                                          , X, I+1 * LDX + o_x, LDX, ONE, ref A, I+(I + 1) * LDA + o_a, LDA);`
` 394:                          // *`
` 395:                          // *              Generate reflection P(i) to annihilate A(i,i+2:n)`
` 396:                          // *`
` 397:                          this._dlarfg.Run(N - I, ref A[I+(I + 1) * LDA + o_a], ref A, I+Math.Min(I + 2, N) * LDA + o_a, LDA, ref TAUP[I + o_taup]);`
` 398:                          E[I + o_e] = A[I+(I + 1) * LDA + o_a];`
` 399:                          A[I+(I + 1) * LDA + o_a] = ONE;`
` 400:                          // *`
` 401:                          // *              Compute X(i+1:m,i)`
` 402:                          // *`
` 403:                          this._dgemv.Run("No transpose", M - I, N - I, ONE, A, I + 1+(I + 1) * LDA + o_a, LDA`
` 404:                                          , A, I+(I + 1) * LDA + o_a, LDA, ZERO, ref X, I + 1+I * LDX + o_x, 1);`
` 405:                          this._dgemv.Run("Transpose", N - I, I, ONE, Y, I + 1+1 * LDY + o_y, LDY`
` 406:                                          , A, I+(I + 1) * LDA + o_a, LDA, ZERO, ref X, 1+I * LDX + o_x, 1);`
` 407:                          this._dgemv.Run("No transpose", M - I, I,  - ONE, A, I + 1+1 * LDA + o_a, LDA`
` 408:                                          , X, 1+I * LDX + o_x, 1, ONE, ref X, I + 1+I * LDX + o_x, 1);`
` 409:                          this._dgemv.Run("No transpose", I - 1, N - I, ONE, A, 1+(I + 1) * LDA + o_a, LDA`
` 410:                                          , A, I+(I + 1) * LDA + o_a, LDA, ZERO, ref X, 1+I * LDX + o_x, 1);`
` 411:                          this._dgemv.Run("No transpose", M - I, I - 1,  - ONE, X, I + 1+1 * LDX + o_x, LDX`
` 412:                                          , X, 1+I * LDX + o_x, 1, ONE, ref X, I + 1+I * LDX + o_x, 1);`
` 413:                          this._dscal.Run(M - I, TAUP[I + o_taup], ref X, I + 1+I * LDX + o_x, 1);`
` 414:                      }`
` 415:                  }`
` 416:              }`
` 417:              else`
` 418:              {`
` 419:                  // *`
` 420:                  // *        Reduce to lower bidiagonal form`
` 421:                  // *`
` 422:                  for (I = 1; I <= NB; I++)`
` 423:                  {`
` 424:                      // *`
` 425:                      // *           Update A(i,i:n)`
` 426:                      // *`
` 427:                      this._dgemv.Run("No transpose", N - I + 1, I - 1,  - ONE, Y, I+1 * LDY + o_y, LDY`
` 428:                                      , A, I+1 * LDA + o_a, LDA, ONE, ref A, I+I * LDA + o_a, LDA);`
` 429:                      this._dgemv.Run("Transpose", I - 1, N - I + 1,  - ONE, A, 1+I * LDA + o_a, LDA`
` 430:                                      , X, I+1 * LDX + o_x, LDX, ONE, ref A, I+I * LDA + o_a, LDA);`
` 431:                      // *`
` 432:                      // *           Generate reflection P(i) to annihilate A(i,i+1:n)`
` 433:                      // *`
` 434:                      this._dlarfg.Run(N - I + 1, ref A[I+I * LDA + o_a], ref A, I+Math.Min(I + 1, N) * LDA + o_a, LDA, ref TAUP[I + o_taup]);`
` 435:                      D[I + o_d] = A[I+I * LDA + o_a];`
` 436:                      if (I < M)`
` 437:                      {`
` 438:                          A[I+I * LDA + o_a] = ONE;`
` 439:                          // *`
` 440:                          // *              Compute X(i+1:m,i)`
` 441:                          // *`
` 442:                          this._dgemv.Run("No transpose", M - I, N - I + 1, ONE, A, I + 1+I * LDA + o_a, LDA`
` 443:                                          , A, I+I * LDA + o_a, LDA, ZERO, ref X, I + 1+I * LDX + o_x, 1);`
` 444:                          this._dgemv.Run("Transpose", N - I + 1, I - 1, ONE, Y, I+1 * LDY + o_y, LDY`
` 445:                                          , A, I+I * LDA + o_a, LDA, ZERO, ref X, 1+I * LDX + o_x, 1);`
` 446:                          this._dgemv.Run("No transpose", M - I, I - 1,  - ONE, A, I + 1+1 * LDA + o_a, LDA`
` 447:                                          , X, 1+I * LDX + o_x, 1, ONE, ref X, I + 1+I * LDX + o_x, 1);`
` 448:                          this._dgemv.Run("No transpose", I - 1, N - I + 1, ONE, A, 1+I * LDA + o_a, LDA`
` 449:                                          , A, I+I * LDA + o_a, LDA, ZERO, ref X, 1+I * LDX + o_x, 1);`
` 450:                          this._dgemv.Run("No transpose", M - I, I - 1,  - ONE, X, I + 1+1 * LDX + o_x, LDX`
` 451:                                          , X, 1+I * LDX + o_x, 1, ONE, ref X, I + 1+I * LDX + o_x, 1);`
` 452:                          this._dscal.Run(M - I, TAUP[I + o_taup], ref X, I + 1+I * LDX + o_x, 1);`
` 453:                          // *`
` 454:                          // *              Update A(i+1:m,i)`
` 455:                          // *`
` 456:                          this._dgemv.Run("No transpose", M - I, I - 1,  - ONE, A, I + 1+1 * LDA + o_a, LDA`
` 457:                                          , Y, I+1 * LDY + o_y, LDY, ONE, ref A, I + 1+I * LDA + o_a, 1);`
` 458:                          this._dgemv.Run("No transpose", M - I, I,  - ONE, X, I + 1+1 * LDX + o_x, LDX`
` 459:                                          , A, 1+I * LDA + o_a, 1, ONE, ref A, I + 1+I * LDA + o_a, 1);`
` 460:                          // *`
` 461:                          // *              Generate reflection Q(i) to annihilate A(i+2:m,i)`
` 462:                          // *`
` 463:                          this._dlarfg.Run(M - I, ref A[I + 1+I * LDA + o_a], ref A, Math.Min(I + 2, M)+I * LDA + o_a, 1, ref TAUQ[I + o_tauq]);`
` 464:                          E[I + o_e] = A[I + 1+I * LDA + o_a];`
` 465:                          A[I + 1+I * LDA + o_a] = ONE;`
` 466:                          // *`
` 467:                          // *              Compute Y(i+1:n,i)`
` 468:                          // *`
` 469:                          this._dgemv.Run("Transpose", M - I, N - I, ONE, A, I + 1+(I + 1) * LDA + o_a, LDA`
` 470:                                          , A, I + 1+I * LDA + o_a, 1, ZERO, ref Y, I + 1+I * LDY + o_y, 1);`
` 471:                          this._dgemv.Run("Transpose", M - I, I - 1, ONE, A, I + 1+1 * LDA + o_a, LDA`
` 472:                                          , A, I + 1+I * LDA + o_a, 1, ZERO, ref Y, 1+I * LDY + o_y, 1);`
` 473:                          this._dgemv.Run("No transpose", N - I, I - 1,  - ONE, Y, I + 1+1 * LDY + o_y, LDY`
` 474:                                          , Y, 1+I * LDY + o_y, 1, ONE, ref Y, I + 1+I * LDY + o_y, 1);`
` 475:                          this._dgemv.Run("Transpose", M - I, I, ONE, X, I + 1+1 * LDX + o_x, LDX`
` 476:                                          , A, I + 1+I * LDA + o_a, 1, ZERO, ref Y, 1+I * LDY + o_y, 1);`
` 477:                          this._dgemv.Run("Transpose", I, N - I,  - ONE, A, 1+(I + 1) * LDA + o_a, LDA`
` 478:                                          , Y, 1+I * LDY + o_y, 1, ONE, ref Y, I + 1+I * LDY + o_y, 1);`
` 479:                          this._dscal.Run(N - I, TAUQ[I + o_tauq], ref Y, I + 1+I * LDY + o_y, 1);`
` 480:                      }`
` 481:                  }`
` 482:              }`
` 483:              return;`
` 484:              // *`
` 485:              // *     End of DLABRD`
` 486:              // *`
` 487:   `
` 488:              #endregion`
` 489:   `
` 490:          }`
` 491:      }`
` 492:  }`