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