`   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:      /// DGGSVD computes the generalized singular value decomposition (GSVD)`
`  27:      /// of an M-by-N real matrix A and P-by-N real matrix B:`
`  28:      /// `
`  29:      /// U'*A*Q = D1*( 0 R ),    V'*B*Q = D2*( 0 R )`
`  30:      /// `
`  31:      /// where U, V and Q are orthogonal matrices, and Z' is the transpose`
`  32:      /// of Z.  Let K+L = the effective numerical rank of the matrix (A',B')',`
`  33:      /// then R is a K+L-by-K+L nonsingular upper triangular matrix, D1 and`
`  34:      /// D2 are M-by-(K+L) and P-by-(K+L) "diagonal" matrices and of the`
`  35:      /// following structures, respectively:`
`  36:      /// `
`  37:      /// If M-K-L .GE. 0,`
`  38:      /// `
`  39:      /// K  L`
`  40:      /// D1 =     K ( I  0 )`
`  41:      /// L ( 0  C )`
`  42:      /// M-K-L ( 0  0 )`
`  43:      /// `
`  44:      /// K  L`
`  45:      /// D2 =   L ( 0  S )`
`  46:      /// P-L ( 0  0 )`
`  47:      /// `
`  48:      /// N-K-L  K    L`
`  49:      /// ( 0 R ) = K (  0   R11  R12 )`
`  50:      /// L (  0    0   R22 )`
`  51:      /// `
`  52:      /// where`
`  53:      /// `
`  54:      /// C = diag( ALPHA(K+1), ... , ALPHA(K+L) ),`
`  55:      /// S = diag( BETA(K+1),  ... , BETA(K+L) ),`
`  56:      /// C**2 + S**2 = I.`
`  57:      /// `
`  58:      /// R is stored in A(1:K+L,N-K-L+1:N) on exit.`
`  59:      /// `
`  60:      /// If M-K-L .LT. 0,`
`  61:      /// `
`  62:      /// K M-K K+L-M`
`  63:      /// D1 =   K ( I  0    0   )`
`  64:      /// M-K ( 0  C    0   )`
`  65:      /// `
`  66:      /// K M-K K+L-M`
`  67:      /// D2 =   M-K ( 0  S    0  )`
`  68:      /// K+L-M ( 0  0    I  )`
`  69:      /// P-L ( 0  0    0  )`
`  70:      /// `
`  71:      /// N-K-L  K   M-K  K+L-M`
`  72:      /// ( 0 R ) =     K ( 0    R11  R12  R13  )`
`  73:      /// M-K ( 0     0   R22  R23  )`
`  74:      /// K+L-M ( 0     0    0   R33  )`
`  75:      /// `
`  76:      /// where`
`  77:      /// `
`  78:      /// C = diag( ALPHA(K+1), ... , ALPHA(M) ),`
`  79:      /// S = diag( BETA(K+1),  ... , BETA(M) ),`
`  80:      /// C**2 + S**2 = I.`
`  81:      /// `
`  82:      /// (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N), and R33 is stored`
`  83:      /// ( 0  R22 R23 )`
`  84:      /// in B(M-K+1:L,N+M-K-L+1:N) on exit.`
`  85:      /// `
`  86:      /// The routine computes C, S, R, and optionally the orthogonal`
`  87:      /// transformation matrices U, V and Q.`
`  88:      /// `
`  89:      /// In particular, if B is an N-by-N nonsingular matrix, then the GSVD of`
`  90:      /// A and B implicitly gives the SVD of A*inv(B):`
`  91:      /// A*inv(B) = U*(D1*inv(D2))*V'.`
`  92:      /// If ( A',B')' has orthonormal columns, then the GSVD of A and B is`
`  93:      /// also equal to the CS decomposition of A and B. Furthermore, the GSVD`
`  94:      /// can be used to derive the solution of the eigenvalue problem:`
`  95:      /// A'*A x = lambda* B'*B x.`
`  96:      /// In some literature, the GSVD of A and B is presented in the form`
`  97:      /// U'*A*X = ( 0 D1 ),   V'*B*X = ( 0 D2 )`
`  98:      /// where U and V are orthogonal and X is nonsingular, D1 and D2 are`
`  99:      /// ``diagonal''.  The former GSVD form can be converted to the latter`
` 100:      /// form by taking the nonsingular matrix X as`
` 101:      /// `
` 102:      /// X = Q*( I   0    )`
` 103:      /// ( 0 inv(R) ).`
` 104:      /// `
` 105:      ///</summary>`
` 106:      public class DGGSVD`
` 107:      {`
` 108:      `
` 109:   `
` 110:          #region Dependencies`
` 111:          `
` 112:          LSAME _lsame; DLAMCH _dlamch; DLANGE _dlange; DCOPY _dcopy; DGGSVP _dggsvp; DTGSJA _dtgsja; XERBLA _xerbla; `
` 113:   `
` 114:          #endregion`
` 115:   `
` 116:   `
` 117:          #region Fields`
` 118:          `
` 119:          bool WANTQ = false; bool WANTU = false; bool WANTV = false; int I = 0; int IBND = 0; int ISUB = 0; int J = 0; `
` 120:          int NCYCLE = 0;double ANORM = 0; double BNORM = 0; double SMAX = 0; double TEMP = 0; double TOLA = 0; double TOLB = 0; `
` 121:          double ULP = 0;double UNFL = 0; `
` 122:   `
` 123:          #endregion`
` 124:   `
` 125:          public DGGSVD(LSAME lsame, DLAMCH dlamch, DLANGE dlange, DCOPY dcopy, DGGSVP dggsvp, DTGSJA dtgsja, XERBLA xerbla)`
` 126:          {`
` 127:      `
` 128:   `
` 129:              #region Set Dependencies`
` 130:              `
` 131:              this._lsame = lsame; this._dlamch = dlamch; this._dlange = dlange; this._dcopy = dcopy; this._dggsvp = dggsvp; `
` 132:              this._dtgsja = dtgsja;this._xerbla = xerbla; `
` 133:   `
` 134:              #endregion`
` 135:   `
` 136:          }`
` 137:      `
` 138:          public DGGSVD()`
` 139:          {`
` 140:      `
` 141:   `
` 142:              #region Dependencies (Initialization)`
` 143:              `
` 144:              LSAME lsame = new LSAME();`
` 145:              DLAMC3 dlamc3 = new DLAMC3();`
` 146:              DLASSQ dlassq = new DLASSQ();`
` 147:              DCOPY dcopy = new DCOPY();`
` 148:              XERBLA xerbla = new XERBLA();`
` 149:              DLAPY2 dlapy2 = new DLAPY2();`
` 150:              DNRM2 dnrm2 = new DNRM2();`
` 151:              DSCAL dscal = new DSCAL();`
` 152:              DSWAP dswap = new DSWAP();`
` 153:              IDAMAX idamax = new IDAMAX();`
` 154:              DLAPMT dlapmt = new DLAPMT();`
` 155:              DDOT ddot = new DDOT();`
` 156:              DAXPY daxpy = new DAXPY();`
` 157:              DLAS2 dlas2 = new DLAS2();`
` 158:              DROT drot = new DROT();`
` 159:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);`
` 160:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);`
` 161:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);`
` 162:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);`
` 163:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);`
` 164:              DLANGE dlange = new DLANGE(dlassq, lsame);`
` 165:              DGEMV dgemv = new DGEMV(lsame, xerbla);`
` 166:              DGER dger = new DGER(xerbla);`
` 167:              DLARF dlarf = new DLARF(dgemv, dger, lsame);`
` 168:              DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);`
` 169:              DGEQR2 dgeqr2 = new DGEQR2(dlarf, dlarfg, xerbla);`
` 170:              DORM2R dorm2r = new DORM2R(lsame, dlarf, xerbla);`
` 171:              DGEQPF dgeqpf = new DGEQPF(dgeqr2, dlarf, dlarfg, dorm2r, dswap, xerbla, idamax, dlamch, dnrm2);`
` 172:              DGERQ2 dgerq2 = new DGERQ2(dlarf, dlarfg, xerbla);`
` 173:              DLACPY dlacpy = new DLACPY(lsame);`
` 174:              DLASET dlaset = new DLASET(lsame);`
` 175:              DORG2R dorg2r = new DORG2R(dlarf, dscal, xerbla);`
` 176:              DORMR2 dormr2 = new DORMR2(lsame, dlarf, xerbla);`
` 177:              DGGSVP dggsvp = new DGGSVP(lsame, dgeqpf, dgeqr2, dgerq2, dlacpy, dlapmt, dlaset, dorg2r, dorm2r, dormr2`
` 178:                                         , xerbla);`
` 179:              DLARTG dlartg = new DLARTG(dlamch);`
` 180:              DLASV2 dlasv2 = new DLASV2(dlamch);`
` 181:              DLAGS2 dlags2 = new DLAGS2(dlartg, dlasv2);`
` 182:              DLAPLL dlapll = new DLAPLL(ddot, daxpy, dlarfg, dlas2);`
` 183:              DTGSJA dtgsja = new DTGSJA(lsame, dcopy, dlags2, dlapll, dlartg, dlaset, drot, dscal, xerbla);`
` 184:   `
` 185:              #endregion`
` 186:   `
` 187:   `
` 188:              #region Set Dependencies`
` 189:              `
` 190:              this._lsame = lsame; this._dlamch = dlamch; this._dlange = dlange; this._dcopy = dcopy; this._dggsvp = dggsvp; `
` 191:              this._dtgsja = dtgsja;this._xerbla = xerbla; `
` 192:   `
` 193:              #endregion`
` 194:   `
` 195:          }`
` 196:          /// <summary>`
` 197:          /// Purpose`
` 198:          /// =======`
` 199:          /// `
` 200:          /// DGGSVD computes the generalized singular value decomposition (GSVD)`
` 201:          /// of an M-by-N real matrix A and P-by-N real matrix B:`
` 202:          /// `
` 203:          /// U'*A*Q = D1*( 0 R ),    V'*B*Q = D2*( 0 R )`
` 204:          /// `
` 205:          /// where U, V and Q are orthogonal matrices, and Z' is the transpose`
` 206:          /// of Z.  Let K+L = the effective numerical rank of the matrix (A',B')',`
` 207:          /// then R is a K+L-by-K+L nonsingular upper triangular matrix, D1 and`
` 208:          /// D2 are M-by-(K+L) and P-by-(K+L) "diagonal" matrices and of the`
` 209:          /// following structures, respectively:`
` 210:          /// `
` 211:          /// If M-K-L .GE. 0,`
` 212:          /// `
` 213:          /// K  L`
` 214:          /// D1 =     K ( I  0 )`
` 215:          /// L ( 0  C )`
` 216:          /// M-K-L ( 0  0 )`
` 217:          /// `
` 218:          /// K  L`
` 219:          /// D2 =   L ( 0  S )`
` 220:          /// P-L ( 0  0 )`
` 221:          /// `
` 222:          /// N-K-L  K    L`
` 223:          /// ( 0 R ) = K (  0   R11  R12 )`
` 224:          /// L (  0    0   R22 )`
` 225:          /// `
` 226:          /// where`
` 227:          /// `
` 228:          /// C = diag( ALPHA(K+1), ... , ALPHA(K+L) ),`
` 229:          /// S = diag( BETA(K+1),  ... , BETA(K+L) ),`
` 230:          /// C**2 + S**2 = I.`
` 231:          /// `
` 232:          /// R is stored in A(1:K+L,N-K-L+1:N) on exit.`
` 233:          /// `
` 234:          /// If M-K-L .LT. 0,`
` 235:          /// `
` 236:          /// K M-K K+L-M`
` 237:          /// D1 =   K ( I  0    0   )`
` 238:          /// M-K ( 0  C    0   )`
` 239:          /// `
` 240:          /// K M-K K+L-M`
` 241:          /// D2 =   M-K ( 0  S    0  )`
` 242:          /// K+L-M ( 0  0    I  )`
` 243:          /// P-L ( 0  0    0  )`
` 244:          /// `
` 245:          /// N-K-L  K   M-K  K+L-M`
` 246:          /// ( 0 R ) =     K ( 0    R11  R12  R13  )`
` 247:          /// M-K ( 0     0   R22  R23  )`
` 248:          /// K+L-M ( 0     0    0   R33  )`
` 249:          /// `
` 250:          /// where`
` 251:          /// `
` 252:          /// C = diag( ALPHA(K+1), ... , ALPHA(M) ),`
` 253:          /// S = diag( BETA(K+1),  ... , BETA(M) ),`
` 254:          /// C**2 + S**2 = I.`
` 255:          /// `
` 256:          /// (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N), and R33 is stored`
` 257:          /// ( 0  R22 R23 )`
` 258:          /// in B(M-K+1:L,N+M-K-L+1:N) on exit.`
` 259:          /// `
` 260:          /// The routine computes C, S, R, and optionally the orthogonal`
` 261:          /// transformation matrices U, V and Q.`
` 262:          /// `
` 263:          /// In particular, if B is an N-by-N nonsingular matrix, then the GSVD of`
` 264:          /// A and B implicitly gives the SVD of A*inv(B):`
` 265:          /// A*inv(B) = U*(D1*inv(D2))*V'.`
` 266:          /// If ( A',B')' has orthonormal columns, then the GSVD of A and B is`
` 267:          /// also equal to the CS decomposition of A and B. Furthermore, the GSVD`
` 268:          /// can be used to derive the solution of the eigenvalue problem:`
` 269:          /// A'*A x = lambda* B'*B x.`
` 270:          /// In some literature, the GSVD of A and B is presented in the form`
` 271:          /// U'*A*X = ( 0 D1 ),   V'*B*X = ( 0 D2 )`
` 272:          /// where U and V are orthogonal and X is nonsingular, D1 and D2 are`
` 273:          /// ``diagonal''.  The former GSVD form can be converted to the latter`
` 274:          /// form by taking the nonsingular matrix X as`
` 275:          /// `
` 276:          /// X = Q*( I   0    )`
` 277:          /// ( 0 inv(R) ).`
` 278:          /// `
` 279:          ///</summary>`
` 280:          /// <param name="JOBU">`
` 281:          /// (input) CHARACTER*1`
` 282:          /// = 'U':  Orthogonal matrix U is computed;`
` 283:          /// = 'N':  U is not computed.`
` 284:          ///</param>`
` 285:          /// <param name="JOBV">`
` 286:          /// (input) CHARACTER*1`
` 287:          /// = 'V':  Orthogonal matrix V is computed;`
` 288:          /// = 'N':  V is not computed.`
` 289:          ///</param>`
` 290:          /// <param name="JOBQ">`
` 291:          /// (input) CHARACTER*1`
` 292:          /// = 'Q':  Orthogonal matrix Q is computed;`
` 293:          /// = 'N':  Q is not computed.`
` 294:          ///</param>`
` 295:          /// <param name="M">`
` 296:          /// (input) INTEGER`
` 297:          /// The number of rows of the matrix A.  M .GE. 0.`
` 298:          ///</param>`
` 299:          /// <param name="N">`
` 300:          /// (input) INTEGER`
` 301:          /// The number of columns of the matrices A and B.  N .GE. 0.`
` 302:          ///</param>`
` 303:          /// <param name="P">`
` 304:          /// (input) INTEGER`
` 305:          /// The number of rows of the matrix B.  P .GE. 0.`
` 306:          ///</param>`
` 307:          /// <param name="K">`
` 308:          /// L`
` 309:          ///</param>`
` 310:          /// <param name="L">`
` 311:          /// ( 0  C )`
` 312:          ///</param>`
` 313:          /// <param name="A">`
` 314:          /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)`
` 315:          /// On entry, the M-by-N matrix A.`
` 316:          /// On exit, A contains the triangular matrix R, or part of R.`
` 317:          /// See Purpose for details.`
` 318:          ///</param>`
` 319:          /// <param name="LDA">`
` 320:          /// (input) INTEGER`
` 321:          /// The leading dimension of the array A. LDA .GE. max(1,M).`
` 322:          ///</param>`
` 323:          /// <param name="B">`
` 324:          /// (input/output) DOUBLE PRECISION array, dimension (LDB,N)`
` 325:          /// On entry, the P-by-N matrix B.`
` 326:          /// On exit, B contains the triangular matrix R if M-K-L .LT. 0.`
` 327:          /// See Purpose for details.`
` 328:          ///</param>`
` 329:          /// <param name="LDB">`
` 330:          /// (input) INTEGER`
` 331:          /// The leading dimension of the array B. LDB .GE. max(1,P).`
` 332:          ///</param>`
` 333:          /// <param name="ALPHA">`
` 334:          /// (output) DOUBLE PRECISION array, dimension (N)`
` 335:          ///</param>`
` 336:          /// <param name="BETA">`
` 337:          /// (output) DOUBLE PRECISION array, dimension (N)`
` 338:          /// On exit, ALPHA and BETA contain the generalized singular`
` 339:          /// value pairs of A and B;`
` 340:          /// ALPHA(1:K) = 1,`
` 341:          /// BETA(1:K)  = 0,`
` 342:          /// and if M-K-L .GE. 0,`
` 343:          /// ALPHA(K+1:K+L) = C,`
` 344:          /// BETA(K+1:K+L)  = S,`
` 345:          /// or if M-K-L .LT. 0,`
` 346:          /// ALPHA(K+1:M)=C, ALPHA(M+1:K+L)=0`
` 347:          /// BETA(K+1:M) =S, BETA(M+1:K+L) =1`
` 348:          /// and`
` 349:          /// ALPHA(K+L+1:N) = 0`
` 350:          /// BETA(K+L+1:N)  = 0`
` 351:          ///</param>`
` 352:          /// <param name="U">`
` 353:          /// (output) DOUBLE PRECISION array, dimension (LDU,M)`
` 354:          /// If JOBU = 'U', U contains the M-by-M orthogonal matrix U.`
` 355:          /// If JOBU = 'N', U is not referenced.`
` 356:          ///</param>`
` 357:          /// <param name="LDU">`
` 358:          /// (input) INTEGER`
` 359:          /// The leading dimension of the array U. LDU .GE. max(1,M) if`
` 360:          /// JOBU = 'U'; LDU .GE. 1 otherwise.`
` 361:          ///</param>`
` 362:          /// <param name="V">`
` 363:          /// (output) DOUBLE PRECISION array, dimension (LDV,P)`
` 364:          /// If JOBV = 'V', V contains the P-by-P orthogonal matrix V.`
` 365:          /// If JOBV = 'N', V is not referenced.`
` 366:          ///</param>`
` 367:          /// <param name="LDV">`
` 368:          /// (input) INTEGER`
` 369:          /// The leading dimension of the array V. LDV .GE. max(1,P) if`
` 370:          /// JOBV = 'V'; LDV .GE. 1 otherwise.`
` 371:          ///</param>`
` 372:          /// <param name="Q">`
` 373:          /// (output) DOUBLE PRECISION array, dimension (LDQ,N)`
` 374:          /// If JOBQ = 'Q', Q contains the N-by-N orthogonal matrix Q.`
` 375:          /// If JOBQ = 'N', Q is not referenced.`
` 376:          ///</param>`
` 377:          /// <param name="LDQ">`
` 378:          /// (input) INTEGER`
` 379:          /// The leading dimension of the array Q. LDQ .GE. max(1,N) if`
` 380:          /// JOBQ = 'Q'; LDQ .GE. 1 otherwise.`
` 381:          ///</param>`
` 382:          /// <param name="WORK">`
` 383:          /// (workspace) DOUBLE PRECISION array,`
` 384:          /// dimension (max(3*N,M,P)+N)`
` 385:          ///</param>`
` 386:          /// <param name="IWORK">`
` 387:          /// (workspace/output) INTEGER array, dimension (N)`
` 388:          /// On exit, IWORK stores the sorting information. More`
` 389:          /// precisely, the following loop will sort ALPHA`
` 390:          /// for I = K+1, min(M,K+L)`
` 391:          /// swap ALPHA(I) and ALPHA(IWORK(I))`
` 392:          /// endfor`
` 393:          /// such that ALPHA(1) .GE. ALPHA(2) .GE. ... .GE. ALPHA(N).`
` 394:          ///</param>`
` 395:          /// <param name="INFO">`
` 396:          /// (output) INTEGER`
` 397:          /// = 0:  successful exit`
` 398:          /// .LT. 0:  if INFO = -i, the i-th argument had an illegal value.`
` 399:          /// .GT. 0:  if INFO = 1, the Jacobi-type procedure failed to`
` 400:          /// converge.  For further details, see subroutine DTGSJA.`
` 401:          ///</param>`
` 402:          public void Run(string JOBU, string JOBV, string JOBQ, int M, int N, int P`
` 403:                           , ref int K, ref int L, ref double[] A, int offset_a, int LDA, ref double[] B, int offset_b, int LDB`
` 404:                           , ref double[] ALPHA, int offset_alpha, ref double[] BETA, int offset_beta, ref double[] U, int offset_u, int LDU, ref double[] V, int offset_v, int LDV`
` 405:                           , ref double[] Q, int offset_q, int LDQ, ref double[] WORK, int offset_work, ref int[] IWORK, int offset_iwork, ref int INFO)`
` 406:          {`
` 407:   `
` 408:              #region Array Index Correction`
` 409:              `
` 410:               int o_a = -1 - LDA + offset_a;  int o_b = -1 - LDB + offset_b;  int o_alpha = -1 + offset_alpha; `
` 411:               int o_beta = -1 + offset_beta; int o_u = -1 - LDU + offset_u;  int o_v = -1 - LDV + offset_v; `
` 412:               int o_q = -1 - LDQ + offset_q; int o_work = -1 + offset_work;  int o_iwork = -1 + offset_iwork; `
` 413:   `
` 414:              #endregion`
` 415:   `
` 416:   `
` 417:              #region Strings`
` 418:              `
` 419:              JOBU = JOBU.Substring(0, 1);  JOBV = JOBV.Substring(0, 1);  JOBQ = JOBQ.Substring(0, 1);  `
` 420:   `
` 421:              #endregion`
` 422:   `
` 423:   `
` 424:              #region Prolog`
` 425:              `
` 426:              // *`
` 427:              // *  -- LAPACK driver routine (version 3.1) --`
` 428:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..`
` 429:              // *     November 2006`
` 430:              // *`
` 431:              // *     .. Scalar Arguments ..`
` 432:              // *     ..`
` 433:              // *     .. Array Arguments ..`
` 434:              // *     ..`
` 435:              // *`
` 436:              // *  Purpose`
` 437:              // *  =======`
` 438:              // *`
` 439:              // *  DGGSVD computes the generalized singular value decomposition (GSVD)`
` 440:              // *  of an M-by-N real matrix A and P-by-N real matrix B:`
` 441:              // *`
` 442:              // *      U'*A*Q = D1*( 0 R ),    V'*B*Q = D2*( 0 R )`
` 443:              // *`
` 444:              // *  where U, V and Q are orthogonal matrices, and Z' is the transpose`
` 445:              // *  of Z.  Let K+L = the effective numerical rank of the matrix (A',B')',`
` 446:              // *  then R is a K+L-by-K+L nonsingular upper triangular matrix, D1 and`
` 447:              // *  D2 are M-by-(K+L) and P-by-(K+L) "diagonal" matrices and of the`
` 448:              // *  following structures, respectively:`
` 449:              // *`
` 450:              // *  If M-K-L >= 0,`
` 451:              // *`
` 452:              // *                      K  L`
` 453:              // *         D1 =     K ( I  0 )`
` 454:              // *                  L ( 0  C )`
` 455:              // *              M-K-L ( 0  0 )`
` 456:              // *`
` 457:              // *                    K  L`
` 458:              // *         D2 =   L ( 0  S )`
` 459:              // *              P-L ( 0  0 )`
` 460:              // *`
` 461:              // *                  N-K-L  K    L`
` 462:              // *    ( 0 R ) = K (  0   R11  R12 )`
` 463:              // *              L (  0    0   R22 )`
` 464:              // *`
` 465:              // *  where`
` 466:              // *`
` 467:              // *    C = diag( ALPHA(K+1), ... , ALPHA(K+L) ),`
` 468:              // *    S = diag( BETA(K+1),  ... , BETA(K+L) ),`
` 469:              // *    C**2 + S**2 = I.`
` 470:              // *`
` 471:              // *    R is stored in A(1:K+L,N-K-L+1:N) on exit.`
` 472:              // *`
` 473:              // *  If M-K-L < 0,`
` 474:              // *`
` 475:              // *                    K M-K K+L-M`
` 476:              // *         D1 =   K ( I  0    0   )`
` 477:              // *              M-K ( 0  C    0   )`
` 478:              // *`
` 479:              // *                      K M-K K+L-M`
` 480:              // *         D2 =   M-K ( 0  S    0  )`
` 481:              // *              K+L-M ( 0  0    I  )`
` 482:              // *                P-L ( 0  0    0  )`
` 483:              // *`
` 484:              // *                     N-K-L  K   M-K  K+L-M`
` 485:              // *    ( 0 R ) =     K ( 0    R11  R12  R13  )`
` 486:              // *                M-K ( 0     0   R22  R23  )`
` 487:              // *              K+L-M ( 0     0    0   R33  )`
` 488:              // *`
` 489:              // *  where`
` 490:              // *`
` 491:              // *    C = diag( ALPHA(K+1), ... , ALPHA(M) ),`
` 492:              // *    S = diag( BETA(K+1),  ... , BETA(M) ),`
` 493:              // *    C**2 + S**2 = I.`
` 494:              // *`
` 495:              // *    (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N), and R33 is stored`
` 496:              // *    ( 0  R22 R23 )`
` 497:              // *    in B(M-K+1:L,N+M-K-L+1:N) on exit.`
` 498:              // *`
` 499:              // *  The routine computes C, S, R, and optionally the orthogonal`
` 500:              // *  transformation matrices U, V and Q.`
` 501:              // *`
` 502:              // *  In particular, if B is an N-by-N nonsingular matrix, then the GSVD of`
` 503:              // *  A and B implicitly gives the SVD of A*inv(B):`
` 504:              // *                       A*inv(B) = U*(D1*inv(D2))*V'.`
` 505:              // *  If ( A',B')' has orthonormal columns, then the GSVD of A and B is`
` 506:              // *  also equal to the CS decomposition of A and B. Furthermore, the GSVD`
` 507:              // *  can be used to derive the solution of the eigenvalue problem:`
` 508:              // *                       A'*A x = lambda* B'*B x.`
` 509:              // *  In some literature, the GSVD of A and B is presented in the form`
` 510:              // *                   U'*A*X = ( 0 D1 ),   V'*B*X = ( 0 D2 )`
` 511:              // *  where U and V are orthogonal and X is nonsingular, D1 and D2 are`
` 512:              // *  ``diagonal''.  The former GSVD form can be converted to the latter`
` 513:              // *  form by taking the nonsingular matrix X as`
` 514:              // *`
` 515:              // *                       X = Q*( I   0    )`
` 516:              // *                             ( 0 inv(R) ).`
` 517:              // *`
` 518:              // *  Arguments`
` 519:              // *  =========`
` 520:              // *`
` 521:              // *  JOBU    (input) CHARACTER*1`
` 522:              // *          = 'U':  Orthogonal matrix U is computed;`
` 523:              // *          = 'N':  U is not computed.`
` 524:              // *`
` 525:              // *  JOBV    (input) CHARACTER*1`
` 526:              // *          = 'V':  Orthogonal matrix V is computed;`
` 527:              // *          = 'N':  V is not computed.`
` 528:              // *`
` 529:              // *  JOBQ    (input) CHARACTER*1`
` 530:              // *          = 'Q':  Orthogonal matrix Q is computed;`
` 531:              // *          = 'N':  Q is not computed.`
` 532:              // *`
` 533:              // *  M       (input) INTEGER`
` 534:              // *          The number of rows of the matrix A.  M >= 0.`
` 535:              // *`
` 536:              // *  N       (input) INTEGER`
` 537:              // *          The number of columns of the matrices A and B.  N >= 0.`
` 538:              // *`
` 539:              // *  P       (input) INTEGER`
` 540:              // *          The number of rows of the matrix B.  P >= 0.`
` 541:              // *`
` 542:              // *  K       (output) INTEGER`
` 543:              // *  L       (output) INTEGER`
` 544:              // *          On exit, K and L specify the dimension of the subblocks`
` 545:              // *          described in the Purpose section.`
` 546:              // *          K + L = effective numerical rank of (A',B')'.`
` 547:              // *`
` 548:              // *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)`
` 549:              // *          On entry, the M-by-N matrix A.`
` 550:              // *          On exit, A contains the triangular matrix R, or part of R.`
` 551:              // *          See Purpose for details.`
` 552:              // *`
` 553:              // *  LDA     (input) INTEGER`
` 554:              // *          The leading dimension of the array A. LDA >= max(1,M).`
` 555:              // *`
` 556:              // *  B       (input/output) DOUBLE PRECISION array, dimension (LDB,N)`
` 557:              // *          On entry, the P-by-N matrix B.`
` 558:              // *          On exit, B contains the triangular matrix R if M-K-L < 0.`
` 559:              // *          See Purpose for details.`
` 560:              // *`
` 561:              // *  LDB     (input) INTEGER`
` 562:              // *          The leading dimension of the array B. LDB >= max(1,P).`
` 563:              // *`
` 564:              // *  ALPHA   (output) DOUBLE PRECISION array, dimension (N)`
` 565:              // *  BETA    (output) DOUBLE PRECISION array, dimension (N)`
` 566:              // *          On exit, ALPHA and BETA contain the generalized singular`
` 567:              // *          value pairs of A and B;`
` 568:              // *            ALPHA(1:K) = 1,`
` 569:              // *            BETA(1:K)  = 0,`
` 570:              // *          and if M-K-L >= 0,`
` 571:              // *            ALPHA(K+1:K+L) = C,`
` 572:              // *            BETA(K+1:K+L)  = S,`
` 573:              // *          or if M-K-L < 0,`
` 574:              // *            ALPHA(K+1:M)=C, ALPHA(M+1:K+L)=0`
` 575:              // *            BETA(K+1:M) =S, BETA(M+1:K+L) =1`
` 576:              // *          and`
` 577:              // *            ALPHA(K+L+1:N) = 0`
` 578:              // *            BETA(K+L+1:N)  = 0`
` 579:              // *`
` 580:              // *  U       (output) DOUBLE PRECISION array, dimension (LDU,M)`
` 581:              // *          If JOBU = 'U', U contains the M-by-M orthogonal matrix U.`
` 582:              // *          If JOBU = 'N', U is not referenced.`
` 583:              // *`
` 584:              // *  LDU     (input) INTEGER`
` 585:              // *          The leading dimension of the array U. LDU >= max(1,M) if`
` 586:              // *          JOBU = 'U'; LDU >= 1 otherwise.`
` 587:              // *`
` 588:              // *  V       (output) DOUBLE PRECISION array, dimension (LDV,P)`
` 589:              // *          If JOBV = 'V', V contains the P-by-P orthogonal matrix V.`
` 590:              // *          If JOBV = 'N', V is not referenced.`
` 591:              // *`
` 592:              // *  LDV     (input) INTEGER`
` 593:              // *          The leading dimension of the array V. LDV >= max(1,P) if`
` 594:              // *          JOBV = 'V'; LDV >= 1 otherwise.`
` 595:              // *`
` 596:              // *  Q       (output) DOUBLE PRECISION array, dimension (LDQ,N)`
` 597:              // *          If JOBQ = 'Q', Q contains the N-by-N orthogonal matrix Q.`
` 598:              // *          If JOBQ = 'N', Q is not referenced.`
` 599:              // *`
` 600:              // *  LDQ     (input) INTEGER`
` 601:              // *          The leading dimension of the array Q. LDQ >= max(1,N) if`
` 602:              // *          JOBQ = 'Q'; LDQ >= 1 otherwise.`
` 603:              // *`
` 604:              // *  WORK    (workspace) DOUBLE PRECISION array,`
` 605:              // *                      dimension (max(3*N,M,P)+N)`
` 606:              // *`
` 607:              // *  IWORK   (workspace/output) INTEGER array, dimension (N)`
` 608:              // *          On exit, IWORK stores the sorting information. More`
` 609:              // *          precisely, the following loop will sort ALPHA`
` 610:              // *             for I = K+1, min(M,K+L)`
` 611:              // *                 swap ALPHA(I) and ALPHA(IWORK(I))`
` 612:              // *             endfor`
` 613:              // *          such that ALPHA(1) >= ALPHA(2) >= ... >= ALPHA(N).`
` 614:              // *`
` 615:              // *  INFO    (output) INTEGER`
` 616:              // *          = 0:  successful exit`
` 617:              // *          < 0:  if INFO = -i, the i-th argument had an illegal value.`
` 618:              // *          > 0:  if INFO = 1, the Jacobi-type procedure failed to`
` 619:              // *                converge.  For further details, see subroutine DTGSJA.`
` 620:              // *`
` 621:              // *  Internal Parameters`
` 622:              // *  ===================`
` 623:              // *`
` 624:              // *  TOLA    DOUBLE PRECISION`
` 625:              // *  TOLB    DOUBLE PRECISION`
` 626:              // *          TOLA and TOLB are the thresholds to determine the effective`
` 627:              // *          rank of (A',B')'. Generally, they are set to`
` 628:              // *                   TOLA = MAX(M,N)*norm(A)*MAZHEPS,`
` 629:              // *                   TOLB = MAX(P,N)*norm(B)*MAZHEPS.`
` 630:              // *          The size of TOLA and TOLB may affect the size of backward`
` 631:              // *          errors of the decomposition.`
` 632:              // *`
` 633:              // *  Further Details`
` 634:              // *  ===============`
` 635:              // *`
` 636:              // *  2-96 Based on modifications by`
` 637:              // *     Ming Gu and Huan Ren, Computer Science Division, University of`
` 638:              // *     California at Berkeley, USA`
` 639:              // *`
` 640:              // *  =====================================================================`
` 641:              // *`
` 642:              // *     .. Local Scalars ..`
` 643:              // *     ..`
` 644:              // *     .. External Functions ..`
` 645:              // *     ..`
` 646:              // *     .. External Subroutines ..`
` 647:              // *     ..`
` 648:              // *     .. Intrinsic Functions ..`
` 649:              //      INTRINSIC          MAX, MIN;`
` 650:              // *     ..`
` 651:              // *     .. Executable Statements ..`
` 652:              // *`
` 653:              // *     Test the input parameters`
` 654:              // *`
` 655:   `
` 656:              #endregion`
` 657:   `
` 658:   `
` 659:              #region Body`
` 660:              `
` 661:              WANTU = this._lsame.Run(JOBU, "U");`
` 662:              WANTV = this._lsame.Run(JOBV, "V");`
` 663:              WANTQ = this._lsame.Run(JOBQ, "Q");`
` 664:              // *`
` 665:              INFO = 0;`
` 666:              if (!(WANTU || this._lsame.Run(JOBU, "N")))`
` 667:              {`
` 668:                  INFO =  - 1;`
` 669:              }`
` 670:              else`
` 671:              {`
` 672:                  if (!(WANTV || this._lsame.Run(JOBV, "N")))`
` 673:                  {`
` 674:                      INFO =  - 2;`
` 675:                  }`
` 676:                  else`
` 677:                  {`
` 678:                      if (!(WANTQ || this._lsame.Run(JOBQ, "N")))`
` 679:                      {`
` 680:                          INFO =  - 3;`
` 681:                      }`
` 682:                      else`
` 683:                      {`
` 684:                          if (M < 0)`
` 685:                          {`
` 686:                              INFO =  - 4;`
` 687:                          }`
` 688:                          else`
` 689:                          {`
` 690:                              if (N < 0)`
` 691:                              {`
` 692:                                  INFO =  - 5;`
` 693:                              }`
` 694:                              else`
` 695:                              {`
` 696:                                  if (P < 0)`
` 697:                                  {`
` 698:                                      INFO =  - 6;`
` 699:                                  }`
` 700:                                  else`
` 701:                                  {`
` 702:                                      if (LDA < Math.Max(1, M))`
` 703:                                      {`
` 704:                                          INFO =  - 10;`
` 705:                                      }`
` 706:                                      else`
` 707:                                      {`
` 708:                                          if (LDB < Math.Max(1, P))`
` 709:                                          {`
` 710:                                              INFO =  - 12;`
` 711:                                          }`
` 712:                                          else`
` 713:                                          {`
` 714:                                              if (LDU < 1 || (WANTU && LDU < M))`
` 715:                                              {`
` 716:                                                  INFO =  - 16;`
` 717:                                              }`
` 718:                                              else`
` 719:                                              {`
` 720:                                                  if (LDV < 1 || (WANTV && LDV < P))`
` 721:                                                  {`
` 722:                                                      INFO =  - 18;`
` 723:                                                  }`
` 724:                                                  else`
` 725:                                                  {`
` 726:                                                      if (LDQ < 1 || (WANTQ && LDQ < N))`
` 727:                                                      {`
` 728:                                                          INFO =  - 20;`
` 729:                                                      }`
` 730:                                                  }`
` 731:                                              }`
` 732:                                          }`
` 733:                                      }`
` 734:                                  }`
` 735:                              }`
` 736:                          }`
` 737:                      }`
` 738:                  }`
` 739:              }`
` 740:              if (INFO != 0)`
` 741:              {`
` 742:                  this._xerbla.Run("DGGSVD",  - INFO);`
` 743:                  return;`
` 744:              }`
` 745:              // *`
` 746:              // *     Compute the Frobenius norm of matrices A and B`
` 747:              // *`
` 748:              ANORM = this._dlange.Run("1", M, N, A, offset_a, LDA, ref WORK, offset_work);`
` 749:              BNORM = this._dlange.Run("1", P, N, B, offset_b, LDB, ref WORK, offset_work);`
` 750:              // *`
` 751:              // *     Get machine precision and set up threshold for determining`
` 752:              // *     the effective numerical rank of the matrices A and B.`
` 753:              // *`
` 754:              ULP = this._dlamch.Run("Precision");`
` 755:              UNFL = this._dlamch.Run("Safe Minimum");`
` 756:              TOLA = Math.Max(M, N) * Math.Max(ANORM, UNFL) * ULP;`
` 757:              TOLB = Math.Max(P, N) * Math.Max(BNORM, UNFL) * ULP;`
` 758:              // *`
` 759:              // *     Preprocessing`
` 760:              // *`
` 761:              this._dggsvp.Run(JOBU, JOBV, JOBQ, M, P, N`
` 762:                               , ref A, offset_a, LDA, ref B, offset_b, LDB, TOLA, TOLB`
` 763:                               , ref K, ref L, ref U, offset_u, LDU, ref V, offset_v, LDV`
` 764:                               , ref Q, offset_q, LDQ, ref IWORK, offset_iwork, ref WORK, offset_work, ref WORK, N + 1 + o_work, ref INFO);`
` 765:              // *`
` 766:              // *     Compute the GSVD of two upper "triangular" matrices`
` 767:              // *`
` 768:              this._dtgsja.Run(JOBU, JOBV, JOBQ, M, P, N`
` 769:                               , K, L, ref A, offset_a, LDA, ref B, offset_b, LDB`
` 770:                               , TOLA, TOLB, ref ALPHA, offset_alpha, ref BETA, offset_beta, ref U, offset_u, LDU`
` 771:                               , ref V, offset_v, LDV, ref Q, offset_q, LDQ, ref WORK, offset_work, ref NCYCLE`
` 772:                               , ref INFO);`
` 773:              // *`
` 774:              // *     Sort the singular values and store the pivot indices in IWORK`
` 775:              // *     Copy ALPHA to WORK, then sort ALPHA in WORK`
` 776:              // *`
` 777:              this._dcopy.Run(N, ALPHA, offset_alpha, 1, ref WORK, offset_work, 1);`
` 778:              IBND = Math.Min(L, M - K);`
` 779:              for (I = 1; I <= IBND; I++)`
` 780:              {`
` 781:                  // *`
` 782:                  // *        Scan for largest ALPHA(K+I)`
` 783:                  // *`
` 784:                  ISUB = I;`
` 785:                  SMAX = WORK[K + I + o_work];`
` 786:                  for (J = I + 1; J <= IBND; J++)`
` 787:                  {`
` 788:                      TEMP = WORK[K + J + o_work];`
` 789:                      if (TEMP > SMAX)`
` 790:                      {`
` 791:                          ISUB = J;`
` 792:                          SMAX = TEMP;`
` 793:                      }`
` 794:                  }`
` 795:                  if (ISUB != I)`
` 796:                  {`
` 797:                      WORK[K + ISUB + o_work] = WORK[K + I + o_work];`
` 798:                      WORK[K + I + o_work] = SMAX;`
` 799:                      IWORK[K + I + o_iwork] = K + ISUB;`
` 800:                  }`
` 801:                  else`
` 802:                  {`
` 803:                      IWORK[K + I + o_iwork] = K + I;`
` 804:                  }`
` 805:              }`
` 806:              // *`
` 807:              return;`
` 808:              // *`
` 809:              // *     End of DGGSVD`
` 810:              // *`
` 811:   `
` 812:              #endregion`
` 813:   `
` 814:          }`
` 815:      }`
` 816:  }`