Skip Navigation Links
Numerical Libraries
Linear Algebra
Differential Equations
Optimization
Samples
Skip Navigation Links
Linear Algebra
CSLapack
CSBlas
   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:  }