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 routine (version 3.1) --
  21:      /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
  22:      /// November 2006
  23:      /// Purpose
  24:      /// =======
  25:      /// 
  26:      /// DGGSVP computes orthogonal matrices U, V and Q such that
  27:      /// 
  28:      /// N-K-L  K    L
  29:      /// U'*A*Q =     K ( 0    A12  A13 )  if M-K-L .GE. 0;
  30:      /// L ( 0     0   A23 )
  31:      /// M-K-L ( 0     0    0  )
  32:      /// 
  33:      /// N-K-L  K    L
  34:      /// =     K ( 0    A12  A13 )  if M-K-L .LT. 0;
  35:      /// M-K ( 0     0   A23 )
  36:      /// 
  37:      /// N-K-L  K    L
  38:      /// V'*B*Q =   L ( 0     0   B13 )
  39:      /// P-L ( 0     0    0  )
  40:      /// 
  41:      /// where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular
  42:      /// upper triangular; A23 is L-by-L upper triangular if M-K-L .GE. 0,
  43:      /// otherwise A23 is (M-K)-by-L upper trapezoidal.  K+L = the effective
  44:      /// numerical rank of the (M+P)-by-N matrix (A',B')'.  Z' denotes the
  45:      /// transpose of Z.
  46:      /// 
  47:      /// This decomposition is the preprocessing step for computing the
  48:      /// Generalized Singular Value Decomposition (GSVD), see subroutine
  49:      /// DGGSVD.
  50:      /// 
  51:      ///</summary>
  52:      public class DGGSVP
  53:      {
  54:      
  55:   
  56:          #region Dependencies
  57:          
  58:          LSAME _lsame; DGEQPF _dgeqpf; DGEQR2 _dgeqr2; DGERQ2 _dgerq2; DLACPY _dlacpy; DLAPMT _dlapmt; DLASET _dlaset; 
  59:          DORG2R _dorg2r;DORM2R _dorm2r; DORMR2 _dormr2; XERBLA _xerbla; 
  60:   
  61:          #endregion
  62:   
  63:   
  64:          #region Fields
  65:          
  66:          const double ZERO = 0.0E+0; const double ONE = 1.0E+0; bool FORWRD = false; bool WANTQ = false; bool WANTU = false; 
  67:          bool WANTV = false;int I = 0; int J = 0; 
  68:   
  69:          #endregion
  70:   
  71:          public DGGSVP(LSAME lsame, DGEQPF dgeqpf, DGEQR2 dgeqr2, DGERQ2 dgerq2, DLACPY dlacpy, DLAPMT dlapmt, DLASET dlaset, DORG2R dorg2r, DORM2R dorm2r, DORMR2 dormr2
  72:                        , XERBLA xerbla)
  73:          {
  74:      
  75:   
  76:              #region Set Dependencies
  77:              
  78:              this._lsame = lsame; this._dgeqpf = dgeqpf; this._dgeqr2 = dgeqr2; this._dgerq2 = dgerq2; this._dlacpy = dlacpy; 
  79:              this._dlapmt = dlapmt;this._dlaset = dlaset; this._dorg2r = dorg2r; this._dorm2r = dorm2r; this._dormr2 = dormr2; 
  80:              this._xerbla = xerbla;
  81:   
  82:              #endregion
  83:   
  84:          }
  85:      
  86:          public DGGSVP()
  87:          {
  88:      
  89:   
  90:              #region Dependencies (Initialization)
  91:              
  92:              LSAME lsame = new LSAME();
  93:              XERBLA xerbla = new XERBLA();
  94:              DLAMC3 dlamc3 = new DLAMC3();
  95:              DLAPY2 dlapy2 = new DLAPY2();
  96:              DNRM2 dnrm2 = new DNRM2();
  97:              DSCAL dscal = new DSCAL();
  98:              DSWAP dswap = new DSWAP();
  99:              IDAMAX idamax = new IDAMAX();
 100:              DLAPMT dlapmt = new DLAPMT();
 101:              DGEMV dgemv = new DGEMV(lsame, xerbla);
 102:              DGER dger = new DGER(xerbla);
 103:              DLARF dlarf = new DLARF(dgemv, dger, lsame);
 104:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);
 105:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);
 106:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);
 107:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
 108:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
 109:              DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
 110:              DGEQR2 dgeqr2 = new DGEQR2(dlarf, dlarfg, xerbla);
 111:              DORM2R dorm2r = new DORM2R(lsame, dlarf, xerbla);
 112:              DGEQPF dgeqpf = new DGEQPF(dgeqr2, dlarf, dlarfg, dorm2r, dswap, xerbla, idamax, dlamch, dnrm2);
 113:              DGERQ2 dgerq2 = new DGERQ2(dlarf, dlarfg, xerbla);
 114:              DLACPY dlacpy = new DLACPY(lsame);
 115:              DLASET dlaset = new DLASET(lsame);
 116:              DORG2R dorg2r = new DORG2R(dlarf, dscal, xerbla);
 117:              DORMR2 dormr2 = new DORMR2(lsame, dlarf, xerbla);
 118:   
 119:              #endregion
 120:   
 121:   
 122:              #region Set Dependencies
 123:              
 124:              this._lsame = lsame; this._dgeqpf = dgeqpf; this._dgeqr2 = dgeqr2; this._dgerq2 = dgerq2; this._dlacpy = dlacpy; 
 125:              this._dlapmt = dlapmt;this._dlaset = dlaset; this._dorg2r = dorg2r; this._dorm2r = dorm2r; this._dormr2 = dormr2; 
 126:              this._xerbla = xerbla;
 127:   
 128:              #endregion
 129:   
 130:          }
 131:          /// <summary>
 132:          /// Purpose
 133:          /// =======
 134:          /// 
 135:          /// DGGSVP computes orthogonal matrices U, V and Q such that
 136:          /// 
 137:          /// N-K-L  K    L
 138:          /// U'*A*Q =     K ( 0    A12  A13 )  if M-K-L .GE. 0;
 139:          /// L ( 0     0   A23 )
 140:          /// M-K-L ( 0     0    0  )
 141:          /// 
 142:          /// N-K-L  K    L
 143:          /// =     K ( 0    A12  A13 )  if M-K-L .LT. 0;
 144:          /// M-K ( 0     0   A23 )
 145:          /// 
 146:          /// N-K-L  K    L
 147:          /// V'*B*Q =   L ( 0     0   B13 )
 148:          /// P-L ( 0     0    0  )
 149:          /// 
 150:          /// where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular
 151:          /// upper triangular; A23 is L-by-L upper triangular if M-K-L .GE. 0,
 152:          /// otherwise A23 is (M-K)-by-L upper trapezoidal.  K+L = the effective
 153:          /// numerical rank of the (M+P)-by-N matrix (A',B')'.  Z' denotes the
 154:          /// transpose of Z.
 155:          /// 
 156:          /// This decomposition is the preprocessing step for computing the
 157:          /// Generalized Singular Value Decomposition (GSVD), see subroutine
 158:          /// DGGSVD.
 159:          /// 
 160:          ///</summary>
 161:          /// <param name="JOBU">
 162:          /// (input) CHARACTER*1
 163:          /// = 'U':  Orthogonal matrix U is computed;
 164:          /// = 'N':  U is not computed.
 165:          ///</param>
 166:          /// <param name="JOBV">
 167:          /// (input) CHARACTER*1
 168:          /// = 'V':  Orthogonal matrix V is computed;
 169:          /// = 'N':  V is not computed.
 170:          ///</param>
 171:          /// <param name="JOBQ">
 172:          /// (input) CHARACTER*1
 173:          /// = 'Q':  Orthogonal matrix Q is computed;
 174:          /// = 'N':  Q is not computed.
 175:          ///</param>
 176:          /// <param name="M">
 177:          /// (input) INTEGER
 178:          /// The number of rows of the matrix A.  M .GE. 0.
 179:          ///</param>
 180:          /// <param name="P">
 181:          /// (input) INTEGER
 182:          /// The number of rows of the matrix B.  P .GE. 0.
 183:          ///</param>
 184:          /// <param name="N">
 185:          /// (input) INTEGER
 186:          /// The number of columns of the matrices A and B.  N .GE. 0.
 187:          ///</param>
 188:          /// <param name="A">
 189:          /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)
 190:          /// On entry, the M-by-N matrix A.
 191:          /// On exit, A contains the triangular (or trapezoidal) matrix
 192:          /// described in the Purpose section.
 193:          ///</param>
 194:          /// <param name="LDA">
 195:          /// (input) INTEGER
 196:          /// The leading dimension of the array A. LDA .GE. max(1,M).
 197:          ///</param>
 198:          /// <param name="B">
 199:          /// (input/output) DOUBLE PRECISION array, dimension (LDB,N)
 200:          /// On entry, the P-by-N matrix B.
 201:          /// On exit, B contains the triangular matrix described in
 202:          /// the Purpose section.
 203:          ///</param>
 204:          /// <param name="LDB">
 205:          /// (input) INTEGER
 206:          /// The leading dimension of the array B. LDB .GE. max(1,P).
 207:          ///</param>
 208:          /// <param name="TOLA">
 209:          /// (input) DOUBLE PRECISION
 210:          ///</param>
 211:          /// <param name="TOLB">
 212:          /// (input) DOUBLE PRECISION
 213:          /// TOLA and TOLB are the thresholds to determine the effective
 214:          /// numerical rank of matrix B and a subblock of A. Generally,
 215:          /// they are set to
 216:          /// TOLA = MAX(M,N)*norm(A)*MAZHEPS,
 217:          /// TOLB = MAX(P,N)*norm(B)*MAZHEPS.
 218:          /// The size of TOLA and TOLB may affect the size of backward
 219:          /// errors of the decomposition.
 220:          ///</param>
 221:          /// <param name="K">
 222:          /// (output) INTEGER
 223:          ///</param>
 224:          /// <param name="L">
 225:          /// ( 0     0   A23 )
 226:          ///</param>
 227:          /// <param name="U">
 228:          /// (output) DOUBLE PRECISION array, dimension (LDU,M)
 229:          /// If JOBU = 'U', U contains the orthogonal matrix U.
 230:          /// If JOBU = 'N', U is not referenced.
 231:          ///</param>
 232:          /// <param name="LDU">
 233:          /// (input) INTEGER
 234:          /// The leading dimension of the array U. LDU .GE. max(1,M) if
 235:          /// JOBU = 'U'; LDU .GE. 1 otherwise.
 236:          ///</param>
 237:          /// <param name="V">
 238:          /// (output) DOUBLE PRECISION array, dimension (LDV,M)
 239:          /// If JOBV = 'V', V contains the orthogonal matrix V.
 240:          /// If JOBV = 'N', V is not referenced.
 241:          ///</param>
 242:          /// <param name="LDV">
 243:          /// (input) INTEGER
 244:          /// The leading dimension of the array V. LDV .GE. max(1,P) if
 245:          /// JOBV = 'V'; LDV .GE. 1 otherwise.
 246:          ///</param>
 247:          /// <param name="Q">
 248:          /// (output) DOUBLE PRECISION array, dimension (LDQ,N)
 249:          /// If JOBQ = 'Q', Q contains the orthogonal matrix Q.
 250:          /// If JOBQ = 'N', Q is not referenced.
 251:          ///</param>
 252:          /// <param name="LDQ">
 253:          /// (input) INTEGER
 254:          /// The leading dimension of the array Q. LDQ .GE. max(1,N) if
 255:          /// JOBQ = 'Q'; LDQ .GE. 1 otherwise.
 256:          ///</param>
 257:          /// <param name="IWORK">
 258:          /// (workspace) INTEGER array, dimension (N)
 259:          ///</param>
 260:          /// <param name="TAU">
 261:          /// (workspace) DOUBLE PRECISION array, dimension (N)
 262:          ///</param>
 263:          /// <param name="WORK">
 264:          /// (workspace) DOUBLE PRECISION array, dimension (max(3*N,M,P))
 265:          ///</param>
 266:          /// <param name="INFO">
 267:          /// (output) INTEGER
 268:          /// = 0:  successful exit
 269:          /// .LT. 0:  if INFO = -i, the i-th argument had an illegal value.
 270:          /// 
 271:          ///</param>
 272:          public void Run(string JOBU, string JOBV, string JOBQ, int M, int P, int N
 273:                           , ref double[] A, int offset_a, int LDA, ref double[] B, int offset_b, int LDB, double TOLA, double TOLB
 274:                           , ref int K, ref int L, ref double[] U, int offset_u, int LDU, ref double[] V, int offset_v, int LDV
 275:                           , ref double[] Q, int offset_q, int LDQ, ref int[] IWORK, int offset_iwork, ref double[] TAU, int offset_tau, ref double[] WORK, int offset_work, ref int INFO)
 276:          {
 277:   
 278:              #region Array Index Correction
 279:              
 280:               int o_a = -1 - LDA + offset_a;  int o_b = -1 - LDB + offset_b;  int o_u = -1 - LDU + offset_u; 
 281:               int o_v = -1 - LDV + offset_v; int o_q = -1 - LDQ + offset_q;  int o_iwork = -1 + offset_iwork; 
 282:               int o_tau = -1 + offset_tau; int o_work = -1 + offset_work; 
 283:   
 284:              #endregion
 285:   
 286:   
 287:              #region Strings
 288:              
 289:              JOBU = JOBU.Substring(0, 1);  JOBV = JOBV.Substring(0, 1);  JOBQ = JOBQ.Substring(0, 1);  
 290:   
 291:              #endregion
 292:   
 293:   
 294:              #region Prolog
 295:              
 296:              // *
 297:              // *  -- LAPACK routine (version 3.1) --
 298:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 299:              // *     November 2006
 300:              // *
 301:              // *     .. Scalar Arguments ..
 302:              // *     ..
 303:              // *     .. Array Arguments ..
 304:              // *     ..
 305:              // *
 306:              // *  Purpose
 307:              // *  =======
 308:              // *
 309:              // *  DGGSVP computes orthogonal matrices U, V and Q such that
 310:              // *
 311:              // *                   N-K-L  K    L
 312:              // *   U'*A*Q =     K ( 0    A12  A13 )  if M-K-L >= 0;
 313:              // *                L ( 0     0   A23 )
 314:              // *            M-K-L ( 0     0    0  )
 315:              // *
 316:              // *                   N-K-L  K    L
 317:              // *          =     K ( 0    A12  A13 )  if M-K-L < 0;
 318:              // *              M-K ( 0     0   A23 )
 319:              // *
 320:              // *                 N-K-L  K    L
 321:              // *   V'*B*Q =   L ( 0     0   B13 )
 322:              // *            P-L ( 0     0    0  )
 323:              // *
 324:              // *  where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular
 325:              // *  upper triangular; A23 is L-by-L upper triangular if M-K-L >= 0,
 326:              // *  otherwise A23 is (M-K)-by-L upper trapezoidal.  K+L = the effective
 327:              // *  numerical rank of the (M+P)-by-N matrix (A',B')'.  Z' denotes the
 328:              // *  transpose of Z.
 329:              // *
 330:              // *  This decomposition is the preprocessing step for computing the
 331:              // *  Generalized Singular Value Decomposition (GSVD), see subroutine
 332:              // *  DGGSVD.
 333:              // *
 334:              // *  Arguments
 335:              // *  =========
 336:              // *
 337:              // *  JOBU    (input) CHARACTER*1
 338:              // *          = 'U':  Orthogonal matrix U is computed;
 339:              // *          = 'N':  U is not computed.
 340:              // *
 341:              // *  JOBV    (input) CHARACTER*1
 342:              // *          = 'V':  Orthogonal matrix V is computed;
 343:              // *          = 'N':  V is not computed.
 344:              // *
 345:              // *  JOBQ    (input) CHARACTER*1
 346:              // *          = 'Q':  Orthogonal matrix Q is computed;
 347:              // *          = 'N':  Q is not computed.
 348:              // *
 349:              // *  M       (input) INTEGER
 350:              // *          The number of rows of the matrix A.  M >= 0.
 351:              // *
 352:              // *  P       (input) INTEGER
 353:              // *          The number of rows of the matrix B.  P >= 0.
 354:              // *
 355:              // *  N       (input) INTEGER
 356:              // *          The number of columns of the matrices A and B.  N >= 0.
 357:              // *
 358:              // *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
 359:              // *          On entry, the M-by-N matrix A.
 360:              // *          On exit, A contains the triangular (or trapezoidal) matrix
 361:              // *          described in the Purpose section.
 362:              // *
 363:              // *  LDA     (input) INTEGER
 364:              // *          The leading dimension of the array A. LDA >= max(1,M).
 365:              // *
 366:              // *  B       (input/output) DOUBLE PRECISION array, dimension (LDB,N)
 367:              // *          On entry, the P-by-N matrix B.
 368:              // *          On exit, B contains the triangular matrix described in
 369:              // *          the Purpose section.
 370:              // *
 371:              // *  LDB     (input) INTEGER
 372:              // *          The leading dimension of the array B. LDB >= max(1,P).
 373:              // *
 374:              // *  TOLA    (input) DOUBLE PRECISION
 375:              // *  TOLB    (input) DOUBLE PRECISION
 376:              // *          TOLA and TOLB are the thresholds to determine the effective
 377:              // *          numerical rank of matrix B and a subblock of A. Generally,
 378:              // *          they are set to
 379:              // *             TOLA = MAX(M,N)*norm(A)*MAZHEPS,
 380:              // *             TOLB = MAX(P,N)*norm(B)*MAZHEPS.
 381:              // *          The size of TOLA and TOLB may affect the size of backward
 382:              // *          errors of the decomposition.
 383:              // *
 384:              // *  K       (output) INTEGER
 385:              // *  L       (output) INTEGER
 386:              // *          On exit, K and L specify the dimension of the subblocks
 387:              // *          described in Purpose.
 388:              // *          K + L = effective numerical rank of (A',B')'.
 389:              // *
 390:              // *  U       (output) DOUBLE PRECISION array, dimension (LDU,M)
 391:              // *          If JOBU = 'U', U contains the orthogonal matrix U.
 392:              // *          If JOBU = 'N', U is not referenced.
 393:              // *
 394:              // *  LDU     (input) INTEGER
 395:              // *          The leading dimension of the array U. LDU >= max(1,M) if
 396:              // *          JOBU = 'U'; LDU >= 1 otherwise.
 397:              // *
 398:              // *  V       (output) DOUBLE PRECISION array, dimension (LDV,M)
 399:              // *          If JOBV = 'V', V contains the orthogonal matrix V.
 400:              // *          If JOBV = 'N', V is not referenced.
 401:              // *
 402:              // *  LDV     (input) INTEGER
 403:              // *          The leading dimension of the array V. LDV >= max(1,P) if
 404:              // *          JOBV = 'V'; LDV >= 1 otherwise.
 405:              // *
 406:              // *  Q       (output) DOUBLE PRECISION array, dimension (LDQ,N)
 407:              // *          If JOBQ = 'Q', Q contains the orthogonal matrix Q.
 408:              // *          If JOBQ = 'N', Q is not referenced.
 409:              // *
 410:              // *  LDQ     (input) INTEGER
 411:              // *          The leading dimension of the array Q. LDQ >= max(1,N) if
 412:              // *          JOBQ = 'Q'; LDQ >= 1 otherwise.
 413:              // *
 414:              // *  IWORK   (workspace) INTEGER array, dimension (N)
 415:              // *
 416:              // *  TAU     (workspace) DOUBLE PRECISION array, dimension (N)
 417:              // *
 418:              // *  WORK    (workspace) DOUBLE PRECISION array, dimension (max(3*N,M,P))
 419:              // *
 420:              // *  INFO    (output) INTEGER
 421:              // *          = 0:  successful exit
 422:              // *          < 0:  if INFO = -i, the i-th argument had an illegal value.
 423:              // *
 424:              // *
 425:              // *  Further Details
 426:              // *  ===============
 427:              // *
 428:              // *  The subroutine uses LAPACK subroutine DGEQPF for the QR factorization
 429:              // *  with column pivoting to detect the effective numerical rank of the
 430:              // *  a matrix. It may be replaced by a better rank determination strategy.
 431:              // *
 432:              // *  =====================================================================
 433:              // *
 434:              // *     .. Parameters ..
 435:              // *     ..
 436:              // *     .. Local Scalars ..
 437:              // *     ..
 438:              // *     .. External Functions ..
 439:              // *     ..
 440:              // *     .. External Subroutines ..
 441:              // *     ..
 442:              // *     .. Intrinsic Functions ..
 443:              //      INTRINSIC          ABS, MAX, MIN;
 444:              // *     ..
 445:              // *     .. Executable Statements ..
 446:              // *
 447:              // *     Test the input parameters
 448:              // *
 449:   
 450:              #endregion
 451:   
 452:   
 453:              #region Body
 454:              
 455:              WANTU = this._lsame.Run(JOBU, "U");
 456:              WANTV = this._lsame.Run(JOBV, "V");
 457:              WANTQ = this._lsame.Run(JOBQ, "Q");
 458:              FORWRD = true;
 459:              // *
 460:              INFO = 0;
 461:              if (!(WANTU || this._lsame.Run(JOBU, "N")))
 462:              {
 463:                  INFO =  - 1;
 464:              }
 465:              else
 466:              {
 467:                  if (!(WANTV || this._lsame.Run(JOBV, "N")))
 468:                  {
 469:                      INFO =  - 2;
 470:                  }
 471:                  else
 472:                  {
 473:                      if (!(WANTQ || this._lsame.Run(JOBQ, "N")))
 474:                      {
 475:                          INFO =  - 3;
 476:                      }
 477:                      else
 478:                      {
 479:                          if (M < 0)
 480:                          {
 481:                              INFO =  - 4;
 482:                          }
 483:                          else
 484:                          {
 485:                              if (P < 0)
 486:                              {
 487:                                  INFO =  - 5;
 488:                              }
 489:                              else
 490:                              {
 491:                                  if (N < 0)
 492:                                  {
 493:                                      INFO =  - 6;
 494:                                  }
 495:                                  else
 496:                                  {
 497:                                      if (LDA < Math.Max(1, M))
 498:                                      {
 499:                                          INFO =  - 8;
 500:                                      }
 501:                                      else
 502:                                      {
 503:                                          if (LDB < Math.Max(1, P))
 504:                                          {
 505:                                              INFO =  - 10;
 506:                                          }
 507:                                          else
 508:                                          {
 509:                                              if (LDU < 1 || (WANTU && LDU < M))
 510:                                              {
 511:                                                  INFO =  - 16;
 512:                                              }
 513:                                              else
 514:                                              {
 515:                                                  if (LDV < 1 || (WANTV && LDV < P))
 516:                                                  {
 517:                                                      INFO =  - 18;
 518:                                                  }
 519:                                                  else
 520:                                                  {
 521:                                                      if (LDQ < 1 || (WANTQ && LDQ < N))
 522:                                                      {
 523:                                                          INFO =  - 20;
 524:                                                      }
 525:                                                  }
 526:                                              }
 527:                                          }
 528:                                      }
 529:                                  }
 530:                              }
 531:                          }
 532:                      }
 533:                  }
 534:              }
 535:              if (INFO != 0)
 536:              {
 537:                  this._xerbla.Run("DGGSVP",  - INFO);
 538:                  return;
 539:              }
 540:              // *
 541:              // *     QR with column pivoting of B: B*P = V*( S11 S12 )
 542:              // *                                           (  0   0  )
 543:              // *
 544:              for (I = 1; I <= N; I++)
 545:              {
 546:                  IWORK[I + o_iwork] = 0;
 547:              }
 548:              this._dgeqpf.Run(P, N, ref B, offset_b, LDB, ref IWORK, offset_iwork, ref TAU, offset_tau
 549:                               , ref WORK, offset_work, ref INFO);
 550:              // *
 551:              // *     Update A := A*P
 552:              // *
 553:              this._dlapmt.Run(FORWRD, M, N, ref A, offset_a, LDA, ref IWORK, offset_iwork);
 554:              // *
 555:              // *     Determine the effective rank of matrix B.
 556:              // *
 557:              L = 0;
 558:              for (I = 1; I <= Math.Min(P, N); I++)
 559:              {
 560:                  if (Math.Abs(B[I+I * LDB + o_b]) > TOLB) L = L + 1;
 561:              }
 562:              // *
 563:              if (WANTV)
 564:              {
 565:                  // *
 566:                  // *        Copy the details of V, and form V.
 567:                  // *
 568:                  this._dlaset.Run("Full", P, P, ZERO, ZERO, ref V, offset_v
 569:                                   , LDV);
 570:                  if (P > 1)
 571:                  {
 572:                      this._dlacpy.Run("Lower", P - 1, N, B, 2+1 * LDB + o_b, LDB, ref V, 2+1 * LDV + o_v
 573:                                       , LDV);
 574:                  }
 575:                  this._dorg2r.Run(P, P, Math.Min(P, N), ref V, offset_v, LDV, TAU, offset_tau
 576:                                   , ref WORK, offset_work, ref INFO);
 577:              }
 578:              // *
 579:              // *     Clean up B
 580:              // *
 581:              for (J = 1; J <= L - 1; J++)
 582:              {
 583:                  for (I = J + 1; I <= L; I++)
 584:                  {
 585:                      B[I+J * LDB + o_b] = ZERO;
 586:                  }
 587:              }
 588:              if (P > L)
 589:              {
 590:                  this._dlaset.Run("Full", P - L, N, ZERO, ZERO, ref B, L + 1+1 * LDB + o_b
 591:                                   , LDB);
 592:              }
 593:              // *
 594:              if (WANTQ)
 595:              {
 596:                  // *
 597:                  // *        Set Q = I and Update Q := Q*P
 598:                  // *
 599:                  this._dlaset.Run("Full", N, N, ZERO, ONE, ref Q, offset_q
 600:                                   , LDQ);
 601:                  this._dlapmt.Run(FORWRD, N, N, ref Q, offset_q, LDQ, ref IWORK, offset_iwork);
 602:              }
 603:              // *
 604:              if (P >= L && N != L)
 605:              {
 606:                  // *
 607:                  // *        RQ factorization of (S11 S12): ( S11 S12 ) = ( 0 S12 )*Z
 608:                  // *
 609:                  this._dgerq2.Run(L, N, ref B, offset_b, LDB, ref TAU, offset_tau, ref WORK, offset_work
 610:                                   , ref INFO);
 611:                  // *
 612:                  // *        Update A := A*Z'
 613:                  // *
 614:                  this._dormr2.Run("Right", "Transpose", M, N, L, ref B, offset_b
 615:                                   , LDB, TAU, offset_tau, ref A, offset_a, LDA, ref WORK, offset_work, ref INFO);
 616:                  // *
 617:                  if (WANTQ)
 618:                  {
 619:                      // *
 620:                      // *           Update Q := Q*Z'
 621:                      // *
 622:                      this._dormr2.Run("Right", "Transpose", N, N, L, ref B, offset_b
 623:                                       , LDB, TAU, offset_tau, ref Q, offset_q, LDQ, ref WORK, offset_work, ref INFO);
 624:                  }
 625:                  // *
 626:                  // *        Clean up B
 627:                  // *
 628:                  this._dlaset.Run("Full", L, N - L, ZERO, ZERO, ref B, offset_b
 629:                                   , LDB);
 630:                  for (J = N - L + 1; J <= N; J++)
 631:                  {
 632:                      for (I = J - N + L + 1; I <= L; I++)
 633:                      {
 634:                          B[I+J * LDB + o_b] = ZERO;
 635:                      }
 636:                  }
 637:                  // *
 638:              }
 639:              // *
 640:              // *     Let              N-L     L
 641:              // *                A = ( A11    A12 ) M,
 642:              // *
 643:              // *     then the following does the complete QR decomposition of A11:
 644:              // *
 645:              // *              A11 = U*(  0  T12 )*P1'
 646:              // *                      (  0   0  )
 647:              // *
 648:              for (I = 1; I <= N - L; I++)
 649:              {
 650:                  IWORK[I + o_iwork] = 0;
 651:              }
 652:              this._dgeqpf.Run(M, N - L, ref A, offset_a, LDA, ref IWORK, offset_iwork, ref TAU, offset_tau
 653:                               , ref WORK, offset_work, ref INFO);
 654:              // *
 655:              // *     Determine the effective rank of A11
 656:              // *
 657:              K = 0;
 658:              for (I = 1; I <= Math.Min(M, N - L); I++)
 659:              {
 660:                  if (Math.Abs(A[I+I * LDA + o_a]) > TOLA) K = K + 1;
 661:              }
 662:              // *
 663:              // *     Update A12 := U'*A12, where A12 = A( 1:M, N-L+1:N )
 664:              // *
 665:              this._dorm2r.Run("Left", "Transpose", M, L, Math.Min(M, N - L), ref A, offset_a
 666:                               , LDA, TAU, offset_tau, ref A, 1+(N - L + 1) * LDA + o_a, LDA, ref WORK, offset_work, ref INFO);
 667:              // *
 668:              if (WANTU)
 669:              {
 670:                  // *
 671:                  // *        Copy the details of U, and form U
 672:                  // *
 673:                  this._dlaset.Run("Full", M, M, ZERO, ZERO, ref U, offset_u
 674:                                   , LDU);
 675:                  if (M > 1)
 676:                  {
 677:                      this._dlacpy.Run("Lower", M - 1, N - L, A, 2+1 * LDA + o_a, LDA, ref U, 2+1 * LDU + o_u
 678:                                       , LDU);
 679:                  }
 680:                  this._dorg2r.Run(M, M, Math.Min(M, N - L), ref U, offset_u, LDU, TAU, offset_tau
 681:                                   , ref WORK, offset_work, ref INFO);
 682:              }
 683:              // *
 684:              if (WANTQ)
 685:              {
 686:                  // *
 687:                  // *        Update Q( 1:N, 1:N-L )  = Q( 1:N, 1:N-L )*P1
 688:                  // *
 689:                  this._dlapmt.Run(FORWRD, N, N - L, ref Q, offset_q, LDQ, ref IWORK, offset_iwork);
 690:              }
 691:              // *
 692:              // *     Clean up A: set the strictly lower triangular part of
 693:              // *     A(1:K, 1:K) = 0, and A( K+1:M, 1:N-L ) = 0.
 694:              // *
 695:              for (J = 1; J <= K - 1; J++)
 696:              {
 697:                  for (I = J + 1; I <= K; I++)
 698:                  {
 699:                      A[I+J * LDA + o_a] = ZERO;
 700:                  }
 701:              }
 702:              if (M > K)
 703:              {
 704:                  this._dlaset.Run("Full", M - K, N - L, ZERO, ZERO, ref A, K + 1+1 * LDA + o_a
 705:                                   , LDA);
 706:              }
 707:              // *
 708:              if (N - L > K)
 709:              {
 710:                  // *
 711:                  // *        RQ factorization of ( T11 T12 ) = ( 0 T12 )*Z1
 712:                  // *
 713:                  this._dgerq2.Run(K, N - L, ref A, offset_a, LDA, ref TAU, offset_tau, ref WORK, offset_work
 714:                                   , ref INFO);
 715:                  // *
 716:                  if (WANTQ)
 717:                  {
 718:                      // *
 719:                      // *           Update Q( 1:N,1:N-L ) = Q( 1:N,1:N-L )*Z1'
 720:                      // *
 721:                      this._dormr2.Run("Right", "Transpose", N, N - L, K, ref A, offset_a
 722:                                       , LDA, TAU, offset_tau, ref Q, offset_q, LDQ, ref WORK, offset_work, ref INFO);
 723:                  }
 724:                  // *
 725:                  // *        Clean up A
 726:                  // *
 727:                  this._dlaset.Run("Full", K, N - L - K, ZERO, ZERO, ref A, offset_a
 728:                                   , LDA);
 729:                  for (J = N - L - K + 1; J <= N - L; J++)
 730:                  {
 731:                      for (I = J - N + L + K + 1; I <= K; I++)
 732:                      {
 733:                          A[I+J * LDA + o_a] = ZERO;
 734:                      }
 735:                  }
 736:                  // *
 737:              }
 738:              // *
 739:              if (M > K)
 740:              {
 741:                  // *
 742:                  // *        QR factorization of A( K+1:M,N-L+1:N )
 743:                  // *
 744:                  this._dgeqr2.Run(M - K, L, ref A, K + 1+(N - L + 1) * LDA + o_a, LDA, ref TAU, offset_tau, ref WORK, offset_work
 745:                                   , ref INFO);
 746:                  // *
 747:                  if (WANTU)
 748:                  {
 749:                      // *
 750:                      // *           Update U(:,K+1:M) := U(:,K+1:M)*U1
 751:                      // *
 752:                      this._dorm2r.Run("Right", "No transpose", M, M - K, Math.Min(M - K, L), ref A, K + 1+(N - L + 1) * LDA + o_a
 753:                                       , LDA, TAU, offset_tau, ref U, 1+(K + 1) * LDU + o_u, LDU, ref WORK, offset_work, ref INFO);
 754:                  }
 755:                  // *
 756:                  // *        Clean up
 757:                  // *
 758:                  for (J = N - L + 1; J <= N; J++)
 759:                  {
 760:                      for (I = J - N + K + L + 1; I <= M; I++)
 761:                      {
 762:                          A[I+J * LDA + o_a] = ZERO;
 763:                      }
 764:                  }
 765:                  // *
 766:              }
 767:              // *
 768:              return;
 769:              // *
 770:              // *     End of DGGSVP
 771:              // *
 772:   
 773:              #endregion
 774:   
 775:          }
 776:      }
 777:  }