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:      /// DGELSY computes the minimum-norm solution to a real linear least
  27:      /// squares problem:
  28:      /// minimize || A * X - B ||
  29:      /// using a complete orthogonal factorization of A.  A is an M-by-N
  30:      /// matrix which may be rank-deficient.
  31:      /// 
  32:      /// Several right hand side vectors b and solution vectors x can be
  33:      /// handled in a single call; they are stored as the columns of the
  34:      /// M-by-NRHS right hand side matrix B and the N-by-NRHS solution
  35:      /// matrix X.
  36:      /// 
  37:      /// The routine first computes a QR factorization with column pivoting:
  38:      /// A * P = Q * [ R11 R12 ]
  39:      /// [  0  R22 ]
  40:      /// with R11 defined as the largest leading submatrix whose estimated
  41:      /// condition number is less than 1/RCOND.  The order of R11, RANK,
  42:      /// is the effective rank of A.
  43:      /// 
  44:      /// Then, R22 is considered to be negligible, and R12 is annihilated
  45:      /// by orthogonal transformations from the right, arriving at the
  46:      /// complete orthogonal factorization:
  47:      /// A * P = Q * [ T11 0 ] * Z
  48:      /// [  0  0 ]
  49:      /// The minimum-norm solution is then
  50:      /// X = P * Z' [ inv(T11)*Q1'*B ]
  51:      /// [        0       ]
  52:      /// where Q1 consists of the first RANK columns of Q.
  53:      /// 
  54:      /// This routine is basically identical to the original xGELSX except
  55:      /// three differences:
  56:      /// o The call to the subroutine xGEQPF has been substituted by the
  57:      /// the call to the subroutine xGEQP3. This subroutine is a Blas-3
  58:      /// version of the QR factorization with column pivoting.
  59:      /// o Matrix B (the right hand side) is updated with Blas-3.
  60:      /// o The permutation of matrix B (the right hand side) is faster and
  61:      /// more simple.
  62:      /// 
  63:      ///</summary>
  64:      public class DGELSY
  65:      {
  66:      
  67:   
  68:          #region Dependencies
  69:          
  70:          ILAENV _ilaenv; DLAMCH _dlamch; DLANGE _dlange; DCOPY _dcopy; DGEQP3 _dgeqp3; DLABAD _dlabad; DLAIC1 _dlaic1; 
  71:          DLASCL _dlascl;DLASET _dlaset; DORMQR _dormqr; DORMRZ _dormrz; DTRSM _dtrsm; DTZRZF _dtzrzf; XERBLA _xerbla; 
  72:   
  73:          #endregion
  74:   
  75:   
  76:          #region Fields
  77:          
  78:          const int IMAX = 1; const int IMIN = 2; const double ZERO = 0.0E+0; const double ONE = 1.0E+0; bool LQUERY = false; 
  79:          int I = 0;int IASCL = 0; int IBSCL = 0; int ISMAX = 0; int ISMIN = 0; int J = 0; int LWKMIN = 0; int LWKOPT = 0; 
  80:          int MN = 0;int NB = 0; int NB1 = 0; int NB2 = 0; int NB3 = 0; int NB4 = 0; double ANRM = 0; double BIGNUM = 0; 
  81:          double BNRM = 0;double C1 = 0; double C2 = 0; double S1 = 0; double S2 = 0; double SMAX = 0; double SMAXPR = 0; 
  82:          double SMIN = 0;double SMINPR = 0; double SMLNUM = 0; double WSIZE = 0; 
  83:   
  84:          #endregion
  85:   
  86:          public DGELSY(ILAENV ilaenv, DLAMCH dlamch, DLANGE dlange, DCOPY dcopy, DGEQP3 dgeqp3, DLABAD dlabad, DLAIC1 dlaic1, DLASCL dlascl, DLASET dlaset, DORMQR dormqr
  87:                        , DORMRZ dormrz, DTRSM dtrsm, DTZRZF dtzrzf, XERBLA xerbla)
  88:          {
  89:      
  90:   
  91:              #region Set Dependencies
  92:              
  93:              this._ilaenv = ilaenv; this._dlamch = dlamch; this._dlange = dlange; this._dcopy = dcopy; this._dgeqp3 = dgeqp3; 
  94:              this._dlabad = dlabad;this._dlaic1 = dlaic1; this._dlascl = dlascl; this._dlaset = dlaset; this._dormqr = dormqr; 
  95:              this._dormrz = dormrz;this._dtrsm = dtrsm; this._dtzrzf = dtzrzf; this._xerbla = xerbla; 
  96:   
  97:              #endregion
  98:   
  99:          }
 100:      
 101:          public DGELSY()
 102:          {
 103:      
 104:   
 105:              #region Dependencies (Initialization)
 106:              
 107:              IEEECK ieeeck = new IEEECK();
 108:              IPARMQ iparmq = new IPARMQ();
 109:              LSAME lsame = new LSAME();
 110:              DLAMC3 dlamc3 = new DLAMC3();
 111:              DLASSQ dlassq = new DLASSQ();
 112:              DCOPY dcopy = new DCOPY();
 113:              XERBLA xerbla = new XERBLA();
 114:              DLAPY2 dlapy2 = new DLAPY2();
 115:              DNRM2 dnrm2 = new DNRM2();
 116:              DSCAL dscal = new DSCAL();
 117:              DSWAP dswap = new DSWAP();
 118:              IDAMAX idamax = new IDAMAX();
 119:              DLABAD dlabad = new DLABAD();
 120:              DDOT ddot = new DDOT();
 121:              DAXPY daxpy = new DAXPY();
 122:              ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
 123:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);
 124:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);
 125:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);
 126:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
 127:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
 128:              DLANGE dlange = new DLANGE(dlassq, lsame);
 129:              DGEMV dgemv = new DGEMV(lsame, xerbla);
 130:              DGER dger = new DGER(xerbla);
 131:              DLARF dlarf = new DLARF(dgemv, dger, lsame);
 132:              DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
 133:              DGEQR2 dgeqr2 = new DGEQR2(dlarf, dlarfg, xerbla);
 134:              DGEMM dgemm = new DGEMM(lsame, xerbla);
 135:              DTRMM dtrmm = new DTRMM(lsame, xerbla);
 136:              DLARFB dlarfb = new DLARFB(lsame, dcopy, dgemm, dtrmm);
 137:              DTRMV dtrmv = new DTRMV(lsame, xerbla);
 138:              DLARFT dlarft = new DLARFT(dgemv, dtrmv, lsame);
 139:              DGEQRF dgeqrf = new DGEQRF(dgeqr2, dlarfb, dlarft, xerbla, ilaenv);
 140:              DLAQP2 dlaqp2 = new DLAQP2(dlarf, dlarfg, dswap, idamax, dlamch, dnrm2);
 141:              DLAQPS dlaqps = new DLAQPS(dgemm, dgemv, dlarfg, dswap, idamax, dlamch, dnrm2);
 142:              DORM2R dorm2r = new DORM2R(lsame, dlarf, xerbla);
 143:              DORMQR dormqr = new DORMQR(lsame, ilaenv, dlarfb, dlarft, dorm2r, xerbla);
 144:              DGEQP3 dgeqp3 = new DGEQP3(dgeqrf, dlaqp2, dlaqps, dormqr, dswap, xerbla, ilaenv, dnrm2);
 145:              DLAIC1 dlaic1 = new DLAIC1(ddot, dlamch);
 146:              DLASCL dlascl = new DLASCL(lsame, dlamch, xerbla);
 147:              DLASET dlaset = new DLASET(lsame);
 148:              DLARZB dlarzb = new DLARZB(lsame, dcopy, dgemm, dtrmm, xerbla);
 149:              DLARZT dlarzt = new DLARZT(dgemv, dtrmv, xerbla, lsame);
 150:              DLARZ dlarz = new DLARZ(daxpy, dcopy, dgemv, dger, lsame);
 151:              DORMR3 dormr3 = new DORMR3(lsame, dlarz, xerbla);
 152:              DORMRZ dormrz = new DORMRZ(lsame, ilaenv, dlarzb, dlarzt, dormr3, xerbla);
 153:              DTRSM dtrsm = new DTRSM(lsame, xerbla);
 154:              DLATRZ dlatrz = new DLATRZ(dlarfg, dlarz);
 155:              DTZRZF dtzrzf = new DTZRZF(dlarzb, dlarzt, dlatrz, xerbla, ilaenv);
 156:   
 157:              #endregion
 158:   
 159:   
 160:              #region Set Dependencies
 161:              
 162:              this._ilaenv = ilaenv; this._dlamch = dlamch; this._dlange = dlange; this._dcopy = dcopy; this._dgeqp3 = dgeqp3; 
 163:              this._dlabad = dlabad;this._dlaic1 = dlaic1; this._dlascl = dlascl; this._dlaset = dlaset; this._dormqr = dormqr; 
 164:              this._dormrz = dormrz;this._dtrsm = dtrsm; this._dtzrzf = dtzrzf; this._xerbla = xerbla; 
 165:   
 166:              #endregion
 167:   
 168:          }
 169:          /// <summary>
 170:          /// Purpose
 171:          /// =======
 172:          /// 
 173:          /// DGELSY computes the minimum-norm solution to a real linear least
 174:          /// squares problem:
 175:          /// minimize || A * X - B ||
 176:          /// using a complete orthogonal factorization of A.  A is an M-by-N
 177:          /// matrix which may be rank-deficient.
 178:          /// 
 179:          /// Several right hand side vectors b and solution vectors x can be
 180:          /// handled in a single call; they are stored as the columns of the
 181:          /// M-by-NRHS right hand side matrix B and the N-by-NRHS solution
 182:          /// matrix X.
 183:          /// 
 184:          /// The routine first computes a QR factorization with column pivoting:
 185:          /// A * P = Q * [ R11 R12 ]
 186:          /// [  0  R22 ]
 187:          /// with R11 defined as the largest leading submatrix whose estimated
 188:          /// condition number is less than 1/RCOND.  The order of R11, RANK,
 189:          /// is the effective rank of A.
 190:          /// 
 191:          /// Then, R22 is considered to be negligible, and R12 is annihilated
 192:          /// by orthogonal transformations from the right, arriving at the
 193:          /// complete orthogonal factorization:
 194:          /// A * P = Q * [ T11 0 ] * Z
 195:          /// [  0  0 ]
 196:          /// The minimum-norm solution is then
 197:          /// X = P * Z' [ inv(T11)*Q1'*B ]
 198:          /// [        0       ]
 199:          /// where Q1 consists of the first RANK columns of Q.
 200:          /// 
 201:          /// This routine is basically identical to the original xGELSX except
 202:          /// three differences:
 203:          /// o The call to the subroutine xGEQPF has been substituted by the
 204:          /// the call to the subroutine xGEQP3. This subroutine is a Blas-3
 205:          /// version of the QR factorization with column pivoting.
 206:          /// o Matrix B (the right hand side) is updated with Blas-3.
 207:          /// o The permutation of matrix B (the right hand side) is faster and
 208:          /// more simple.
 209:          /// 
 210:          ///</summary>
 211:          /// <param name="M">
 212:          /// (input) INTEGER
 213:          /// The number of rows of the matrix A.  M .GE. 0.
 214:          ///</param>
 215:          /// <param name="N">
 216:          /// (input) INTEGER
 217:          /// The number of columns of the matrix A.  N .GE. 0.
 218:          ///</param>
 219:          /// <param name="NRHS">
 220:          /// (input) INTEGER
 221:          /// The number of right hand sides, i.e., the number of
 222:          /// columns of matrices B and X. NRHS .GE. 0.
 223:          ///</param>
 224:          /// <param name="A">
 225:          /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)
 226:          /// On entry, the M-by-N matrix A.
 227:          /// On exit, A has been overwritten by details of its
 228:          /// complete orthogonal factorization.
 229:          ///</param>
 230:          /// <param name="LDA">
 231:          /// (input) INTEGER
 232:          /// The leading dimension of the array A.  LDA .GE. max(1,M).
 233:          ///</param>
 234:          /// <param name="B">
 235:          /// (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
 236:          /// On entry, the M-by-NRHS right hand side matrix B.
 237:          /// On exit, the N-by-NRHS solution matrix X.
 238:          ///</param>
 239:          /// <param name="LDB">
 240:          /// (input) INTEGER
 241:          /// The leading dimension of the array B. LDB .GE. max(1,M,N).
 242:          ///</param>
 243:          /// <param name="JPVT">
 244:          /// (input/output) INTEGER array, dimension (N)
 245:          /// On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
 246:          /// to the front of AP, otherwise column i is a free column.
 247:          /// On exit, if JPVT(i) = k, then the i-th column of AP
 248:          /// was the k-th column of A.
 249:          ///</param>
 250:          /// <param name="RCOND">
 251:          /// (input) DOUBLE PRECISION
 252:          /// RCOND is used to determine the effective rank of A, which
 253:          /// is defined as the order of the largest leading triangular
 254:          /// submatrix R11 in the QR factorization with pivoting of A,
 255:          /// whose estimated condition number .LT. 1/RCOND.
 256:          ///</param>
 257:          /// <param name="RANK">
 258:          /// (output) INTEGER
 259:          /// The effective rank of A, i.e., the order of the submatrix
 260:          /// R11.  This is the same as the order of the submatrix T11
 261:          /// in the complete orthogonal factorization of A.
 262:          ///</param>
 263:          /// <param name="WORK">
 264:          /// (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 265:          /// On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 266:          ///</param>
 267:          /// <param name="LWORK">
 268:          /// (input) INTEGER
 269:          /// The dimension of the array WORK.
 270:          /// The unblocked strategy requires that:
 271:          /// LWORK .GE. MAX( MN+3*N+1, 2*MN+NRHS ),
 272:          /// where MN = min( M, N ).
 273:          /// The block algorithm requires that:
 274:          /// LWORK .GE. MAX( MN+2*N+NB*(N+1), 2*MN+NB*NRHS ),
 275:          /// where NB is an upper bound on the blocksize returned
 276:          /// by ILAENV for the routines DGEQP3, DTZRZF, STZRQF, DORMQR,
 277:          /// and DORMRZ.
 278:          /// 
 279:          /// If LWORK = -1, then a workspace query is assumed; the routine
 280:          /// only calculates the optimal size of the WORK array, returns
 281:          /// this value as the first entry of the WORK array, and no error
 282:          /// message related to LWORK is issued by XERBLA.
 283:          ///</param>
 284:          /// <param name="INFO">
 285:          /// (output) INTEGER
 286:          /// = 0: successful exit
 287:          /// .LT. 0: If INFO = -i, the i-th argument had an illegal value.
 288:          ///</param>
 289:          public void Run(int M, int N, int NRHS, ref double[] A, int offset_a, int LDA, ref double[] B, int offset_b
 290:                           , int LDB, ref int[] JPVT, int offset_jpvt, double RCOND, ref int RANK, ref double[] WORK, int offset_work, int LWORK
 291:                           , ref int INFO)
 292:          {
 293:   
 294:              #region Array Index Correction
 295:              
 296:               int o_a = -1 - LDA + offset_a;  int o_b = -1 - LDB + offset_b;  int o_jpvt = -1 + offset_jpvt; 
 297:               int o_work = -1 + offset_work;
 298:   
 299:              #endregion
 300:   
 301:   
 302:              #region Prolog
 303:              
 304:              // *
 305:              // *  -- LAPACK driver routine (version 3.1) --
 306:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 307:              // *     November 2006
 308:              // *
 309:              // *     .. Scalar Arguments ..
 310:              // *     ..
 311:              // *     .. Array Arguments ..
 312:              // *     ..
 313:              // *
 314:              // *  Purpose
 315:              // *  =======
 316:              // *
 317:              // *  DGELSY computes the minimum-norm solution to a real linear least
 318:              // *  squares problem:
 319:              // *      minimize || A * X - B ||
 320:              // *  using a complete orthogonal factorization of A.  A is an M-by-N
 321:              // *  matrix which may be rank-deficient.
 322:              // *
 323:              // *  Several right hand side vectors b and solution vectors x can be
 324:              // *  handled in a single call; they are stored as the columns of the
 325:              // *  M-by-NRHS right hand side matrix B and the N-by-NRHS solution
 326:              // *  matrix X.
 327:              // *
 328:              // *  The routine first computes a QR factorization with column pivoting:
 329:              // *      A * P = Q * [ R11 R12 ]
 330:              // *                  [  0  R22 ]
 331:              // *  with R11 defined as the largest leading submatrix whose estimated
 332:              // *  condition number is less than 1/RCOND.  The order of R11, RANK,
 333:              // *  is the effective rank of A.
 334:              // *
 335:              // *  Then, R22 is considered to be negligible, and R12 is annihilated
 336:              // *  by orthogonal transformations from the right, arriving at the
 337:              // *  complete orthogonal factorization:
 338:              // *     A * P = Q * [ T11 0 ] * Z
 339:              // *                 [  0  0 ]
 340:              // *  The minimum-norm solution is then
 341:              // *     X = P * Z' [ inv(T11)*Q1'*B ]
 342:              // *                [        0       ]
 343:              // *  where Q1 consists of the first RANK columns of Q.
 344:              // *
 345:              // *  This routine is basically identical to the original xGELSX except
 346:              // *  three differences:
 347:              // *    o The call to the subroutine xGEQPF has been substituted by the
 348:              // *      the call to the subroutine xGEQP3. This subroutine is a Blas-3
 349:              // *      version of the QR factorization with column pivoting.
 350:              // *    o Matrix B (the right hand side) is updated with Blas-3.
 351:              // *    o The permutation of matrix B (the right hand side) is faster and
 352:              // *      more simple.
 353:              // *
 354:              // *  Arguments
 355:              // *  =========
 356:              // *
 357:              // *  M       (input) INTEGER
 358:              // *          The number of rows of the matrix A.  M >= 0.
 359:              // *
 360:              // *  N       (input) INTEGER
 361:              // *          The number of columns of the matrix A.  N >= 0.
 362:              // *
 363:              // *  NRHS    (input) INTEGER
 364:              // *          The number of right hand sides, i.e., the number of
 365:              // *          columns of matrices B and X. NRHS >= 0.
 366:              // *
 367:              // *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
 368:              // *          On entry, the M-by-N matrix A.
 369:              // *          On exit, A has been overwritten by details of its
 370:              // *          complete orthogonal factorization.
 371:              // *
 372:              // *  LDA     (input) INTEGER
 373:              // *          The leading dimension of the array A.  LDA >= max(1,M).
 374:              // *
 375:              // *  B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
 376:              // *          On entry, the M-by-NRHS right hand side matrix B.
 377:              // *          On exit, the N-by-NRHS solution matrix X.
 378:              // *
 379:              // *  LDB     (input) INTEGER
 380:              // *          The leading dimension of the array B. LDB >= max(1,M,N).
 381:              // *
 382:              // *  JPVT    (input/output) INTEGER array, dimension (N)
 383:              // *          On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
 384:              // *          to the front of AP, otherwise column i is a free column.
 385:              // *          On exit, if JPVT(i) = k, then the i-th column of AP
 386:              // *          was the k-th column of A.
 387:              // *
 388:              // *  RCOND   (input) DOUBLE PRECISION
 389:              // *          RCOND is used to determine the effective rank of A, which
 390:              // *          is defined as the order of the largest leading triangular
 391:              // *          submatrix R11 in the QR factorization with pivoting of A,
 392:              // *          whose estimated condition number < 1/RCOND.
 393:              // *
 394:              // *  RANK    (output) INTEGER
 395:              // *          The effective rank of A, i.e., the order of the submatrix
 396:              // *          R11.  This is the same as the order of the submatrix T11
 397:              // *          in the complete orthogonal factorization of A.
 398:              // *
 399:              // *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 400:              // *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 401:              // *
 402:              // *  LWORK   (input) INTEGER
 403:              // *          The dimension of the array WORK.
 404:              // *          The unblocked strategy requires that:
 405:              // *             LWORK >= MAX( MN+3*N+1, 2*MN+NRHS ),
 406:              // *          where MN = min( M, N ).
 407:              // *          The block algorithm requires that:
 408:              // *             LWORK >= MAX( MN+2*N+NB*(N+1), 2*MN+NB*NRHS ),
 409:              // *          where NB is an upper bound on the blocksize returned
 410:              // *          by ILAENV for the routines DGEQP3, DTZRZF, STZRQF, DORMQR,
 411:              // *          and DORMRZ.
 412:              // *
 413:              // *          If LWORK = -1, then a workspace query is assumed; the routine
 414:              // *          only calculates the optimal size of the WORK array, returns
 415:              // *          this value as the first entry of the WORK array, and no error
 416:              // *          message related to LWORK is issued by XERBLA.
 417:              // *
 418:              // *  INFO    (output) INTEGER
 419:              // *          = 0: successful exit
 420:              // *          < 0: If INFO = -i, the i-th argument had an illegal value.
 421:              // *
 422:              // *  Further Details
 423:              // *  ===============
 424:              // *
 425:              // *  Based on contributions by
 426:              // *    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
 427:              // *    E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
 428:              // *    G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
 429:              // *
 430:              // *  =====================================================================
 431:              // *
 432:              // *     .. Parameters ..
 433:              // *     ..
 434:              // *     .. Local Scalars ..
 435:              // *     ..
 436:              // *     .. External Functions ..
 437:              // *     ..
 438:              // *     .. External Subroutines ..
 439:              // *     ..
 440:              // *     .. Intrinsic Functions ..
 441:              //      INTRINSIC          ABS, MAX, MIN;
 442:              // *     ..
 443:              // *     .. Executable Statements ..
 444:              // *
 445:   
 446:              #endregion
 447:   
 448:   
 449:              #region Body
 450:              
 451:              MN = Math.Min(M, N);
 452:              ISMIN = MN + 1;
 453:              ISMAX = 2 * MN + 1;
 454:              // *
 455:              // *     Test the input arguments.
 456:              // *
 457:              INFO = 0;
 458:              LQUERY = (LWORK ==  - 1);
 459:              if (M < 0)
 460:              {
 461:                  INFO =  - 1;
 462:              }
 463:              else
 464:              {
 465:                  if (N < 0)
 466:                  {
 467:                      INFO =  - 2;
 468:                  }
 469:                  else
 470:                  {
 471:                      if (NRHS < 0)
 472:                      {
 473:                          INFO =  - 3;
 474:                      }
 475:                      else
 476:                      {
 477:                          if (LDA < Math.Max(1, M))
 478:                          {
 479:                              INFO =  - 5;
 480:                          }
 481:                          else
 482:                          {
 483:                              if (LDB < Math.Max(1, Math.Max(M, N)))
 484:                              {
 485:                                  INFO =  - 7;
 486:                              }
 487:                          }
 488:                      }
 489:                  }
 490:              }
 491:              // *
 492:              // *     Figure out optimal block size
 493:              // *
 494:              if (INFO == 0)
 495:              {
 496:                  if (MN == 0 || NRHS == 0)
 497:                  {
 498:                      LWKMIN = 1;
 499:                      LWKOPT = 1;
 500:                  }
 501:                  else
 502:                  {
 503:                      NB1 = this._ilaenv.Run(1, "DGEQRF", " ", M, N,  - 1,  - 1);
 504:                      NB2 = this._ilaenv.Run(1, "DGERQF", " ", M, N,  - 1,  - 1);
 505:                      NB3 = this._ilaenv.Run(1, "DORMQR", " ", M, N, NRHS,  - 1);
 506:                      NB4 = this._ilaenv.Run(1, "DORMRQ", " ", M, N, NRHS,  - 1);
 507:                      NB = Math.Max(NB1, Math.Max(NB2, Math.Max(NB3, NB4)));
 508:                      LWKMIN = MN + Math.Max(2 * MN, Math.Max(N + 1, MN + NRHS));
 509:                      LWKOPT = Math.Max(LWKMIN, Math.Max(MN + 2 * N + NB * (N + 1), 2 * MN + NB * NRHS));
 510:                  }
 511:                  WORK[1 + o_work] = LWKOPT;
 512:                  // *
 513:                  if (LWORK < LWKMIN && !LQUERY)
 514:                  {
 515:                      INFO =  - 12;
 516:                  }
 517:              }
 518:              // *
 519:              if (INFO != 0)
 520:              {
 521:                  this._xerbla.Run("DGELSY",  - INFO);
 522:                  return;
 523:              }
 524:              else
 525:              {
 526:                  if (LQUERY)
 527:                  {
 528:                      return;
 529:                  }
 530:              }
 531:              // *
 532:              // *     Quick return if possible
 533:              // *
 534:              if (MN == 0 || NRHS == 0)
 535:              {
 536:                  RANK = 0;
 537:                  return;
 538:              }
 539:              // *
 540:              // *     Get machine parameters
 541:              // *
 542:              SMLNUM = this._dlamch.Run("S") / this._dlamch.Run("P");
 543:              BIGNUM = ONE / SMLNUM;
 544:              this._dlabad.Run(ref SMLNUM, ref BIGNUM);
 545:              // *
 546:              // *     Scale A, B if max entries outside range [SMLNUM,BIGNUM]
 547:              // *
 548:              ANRM = this._dlange.Run("M", M, N, A, offset_a, LDA, ref WORK, offset_work);
 549:              IASCL = 0;
 550:              if (ANRM > ZERO && ANRM < SMLNUM)
 551:              {
 552:                  // *
 553:                  // *        Scale matrix norm up to SMLNUM
 554:                  // *
 555:                  this._dlascl.Run("G", 0, 0, ANRM, SMLNUM, M
 556:                                   , N, ref A, offset_a, LDA, ref INFO);
 557:                  IASCL = 1;
 558:              }
 559:              else
 560:              {
 561:                  if (ANRM > BIGNUM)
 562:                  {
 563:                      // *
 564:                      // *        Scale matrix norm down to BIGNUM
 565:                      // *
 566:                      this._dlascl.Run("G", 0, 0, ANRM, BIGNUM, M
 567:                                       , N, ref A, offset_a, LDA, ref INFO);
 568:                      IASCL = 2;
 569:                  }
 570:                  else
 571:                  {
 572:                      if (ANRM == ZERO)
 573:                      {
 574:                          // *
 575:                          // *        Matrix all zero. Return zero solution.
 576:                          // *
 577:                          this._dlaset.Run("F", Math.Max(M, N), NRHS, ZERO, ZERO, ref B, offset_b
 578:                                           , LDB);
 579:                          RANK = 0;
 580:                          goto LABEL70;
 581:                      }
 582:                  }
 583:              }
 584:              // *
 585:              BNRM = this._dlange.Run("M", M, NRHS, B, offset_b, LDB, ref WORK, offset_work);
 586:              IBSCL = 0;
 587:              if (BNRM > ZERO && BNRM < SMLNUM)
 588:              {
 589:                  // *
 590:                  // *        Scale matrix norm up to SMLNUM
 591:                  // *
 592:                  this._dlascl.Run("G", 0, 0, BNRM, SMLNUM, M
 593:                                   , NRHS, ref B, offset_b, LDB, ref INFO);
 594:                  IBSCL = 1;
 595:              }
 596:              else
 597:              {
 598:                  if (BNRM > BIGNUM)
 599:                  {
 600:                      // *
 601:                      // *        Scale matrix norm down to BIGNUM
 602:                      // *
 603:                      this._dlascl.Run("G", 0, 0, BNRM, BIGNUM, M
 604:                                       , NRHS, ref B, offset_b, LDB, ref INFO);
 605:                      IBSCL = 2;
 606:                  }
 607:              }
 608:              // *
 609:              // *     Compute QR factorization with column pivoting of A:
 610:              // *        A * P = Q * R
 611:              // *
 612:              this._dgeqp3.Run(M, N, ref A, offset_a, LDA, ref JPVT, offset_jpvt, ref WORK, 1 + o_work
 613:                               , ref WORK, MN + 1 + o_work, LWORK - MN, ref INFO);
 614:              WSIZE = MN + WORK[MN + 1 + o_work];
 615:              // *
 616:              // *     workspace: MN+2*N+NB*(N+1).
 617:              // *     Details of Householder rotations stored in WORK(1:MN).
 618:              // *
 619:              // *     Determine RANK using incremental condition estimation
 620:              // *
 621:              WORK[ISMIN + o_work] = ONE;
 622:              WORK[ISMAX + o_work] = ONE;
 623:              SMAX = Math.Abs(A[1+1 * LDA + o_a]);
 624:              SMIN = SMAX;
 625:              if (Math.Abs(A[1+1 * LDA + o_a]) == ZERO)
 626:              {
 627:                  RANK = 0;
 628:                  this._dlaset.Run("F", Math.Max(M, N), NRHS, ZERO, ZERO, ref B, offset_b
 629:                                   , LDB);
 630:                  goto LABEL70;
 631:              }
 632:              else
 633:              {
 634:                  RANK = 1;
 635:              }
 636:              // *
 637:          LABEL10:;
 638:              if (RANK < MN)
 639:              {
 640:                  I = RANK + 1;
 641:                  this._dlaic1.Run(IMIN, RANK, WORK, ISMIN + o_work, SMIN, A, 1+I * LDA + o_a, A[I+I * LDA + o_a]
 642:                                   , ref SMINPR, ref S1, ref C1);
 643:                  this._dlaic1.Run(IMAX, RANK, WORK, ISMAX + o_work, SMAX, A, 1+I * LDA + o_a, A[I+I * LDA + o_a]
 644:                                   , ref SMAXPR, ref S2, ref C2);
 645:                  // *
 646:                  if (SMAXPR * RCOND <= SMINPR)
 647:                  {
 648:                      for (I = 1; I <= RANK; I++)
 649:                      {
 650:                          WORK[ISMIN + I - 1 + o_work] = S1 * WORK[ISMIN + I - 1 + o_work];
 651:                          WORK[ISMAX + I - 1 + o_work] = S2 * WORK[ISMAX + I - 1 + o_work];
 652:                      }
 653:                      WORK[ISMIN + RANK + o_work] = C1;
 654:                      WORK[ISMAX + RANK + o_work] = C2;
 655:                      SMIN = SMINPR;
 656:                      SMAX = SMAXPR;
 657:                      RANK = RANK + 1;
 658:                      goto LABEL10;
 659:                  }
 660:              }
 661:              // *
 662:              // *     workspace: 3*MN.
 663:              // *
 664:              // *     Logically partition R = [ R11 R12 ]
 665:              // *                             [  0  R22 ]
 666:              // *     where R11 = R(1:RANK,1:RANK)
 667:              // *
 668:              // *     [R11,R12] = [ T11, 0 ] * Y
 669:              // *
 670:              if (RANK < N)
 671:              {
 672:                  this._dtzrzf.Run(RANK, N, ref A, offset_a, LDA, ref WORK, MN + 1 + o_work, ref WORK, 2 * MN + 1 + o_work
 673:                                   , LWORK - 2 * MN, ref INFO);
 674:              }
 675:              // *
 676:              // *     workspace: 2*MN.
 677:              // *     Details of Householder rotations stored in WORK(MN+1:2*MN)
 678:              // *
 679:              // *     B(1:M,1:NRHS) := Q' * B(1:M,1:NRHS)
 680:              // *
 681:              this._dormqr.Run("Left", "Transpose", M, NRHS, MN, ref A, offset_a
 682:                               , LDA, WORK, 1 + o_work, ref B, offset_b, LDB, ref WORK, 2 * MN + 1 + o_work, LWORK - 2 * MN
 683:                               , ref INFO);
 684:              WSIZE = Math.Max(WSIZE, 2 * MN + WORK[2 * MN + 1 + o_work]);
 685:              // *
 686:              // *     workspace: 2*MN+NB*NRHS.
 687:              // *
 688:              // *     B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS)
 689:              // *
 690:              this._dtrsm.Run("Left", "Upper", "No transpose", "Non-unit", RANK, NRHS
 691:                              , ONE, A, offset_a, LDA, ref B, offset_b, LDB);
 692:              // *
 693:              for (J = 1; J <= NRHS; J++)
 694:              {
 695:                  for (I = RANK + 1; I <= N; I++)
 696:                  {
 697:                      B[I+J * LDB + o_b] = ZERO;
 698:                  }
 699:              }
 700:              // *
 701:              // *     B(1:N,1:NRHS) := Y' * B(1:N,1:NRHS)
 702:              // *
 703:              if (RANK < N)
 704:              {
 705:                  this._dormrz.Run("Left", "Transpose", N, NRHS, RANK, N - RANK
 706:                                   , A, offset_a, LDA, WORK, MN + 1 + o_work, ref B, offset_b, LDB, ref WORK, 2 * MN + 1 + o_work
 707:                                   , LWORK - 2 * MN, ref INFO);
 708:              }
 709:              // *
 710:              // *     workspace: 2*MN+NRHS.
 711:              // *
 712:              // *     B(1:N,1:NRHS) := P * B(1:N,1:NRHS)
 713:              // *
 714:              for (J = 1; J <= NRHS; J++)
 715:              {
 716:                  for (I = 1; I <= N; I++)
 717:                  {
 718:                      WORK[JPVT[I + o_jpvt] + o_work] = B[I+J * LDB + o_b];
 719:                  }
 720:                  this._dcopy.Run(N, WORK, 1 + o_work, 1, ref B, 1+J * LDB + o_b, 1);
 721:              }
 722:              // *
 723:              // *     workspace: N.
 724:              // *
 725:              // *     Undo scaling
 726:              // *
 727:              if (IASCL == 1)
 728:              {
 729:                  this._dlascl.Run("G", 0, 0, ANRM, SMLNUM, N
 730:                                   , NRHS, ref B, offset_b, LDB, ref INFO);
 731:                  this._dlascl.Run("U", 0, 0, SMLNUM, ANRM, RANK
 732:                                   , RANK, ref A, offset_a, LDA, ref INFO);
 733:              }
 734:              else
 735:              {
 736:                  if (IASCL == 2)
 737:                  {
 738:                      this._dlascl.Run("G", 0, 0, ANRM, BIGNUM, N
 739:                                       , NRHS, ref B, offset_b, LDB, ref INFO);
 740:                      this._dlascl.Run("U", 0, 0, BIGNUM, ANRM, RANK
 741:                                       , RANK, ref A, offset_a, LDA, ref INFO);
 742:                  }
 743:              }
 744:              if (IBSCL == 1)
 745:              {
 746:                  this._dlascl.Run("G", 0, 0, SMLNUM, BNRM, N
 747:                                   , NRHS, ref B, offset_b, LDB, ref INFO);
 748:              }
 749:              else
 750:              {
 751:                  if (IBSCL == 2)
 752:                  {
 753:                      this._dlascl.Run("G", 0, 0, BIGNUM, BNRM, N
 754:                                       , NRHS, ref B, offset_b, LDB, ref INFO);
 755:                  }
 756:              }
 757:              // *
 758:          LABEL70:;
 759:              WORK[1 + o_work] = LWKOPT;
 760:              // *
 761:              return;
 762:              // *
 763:              // *     End of DGELSY
 764:              // *
 765:   
 766:              #endregion
 767:   
 768:          }
 769:      }
 770:  }