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:      /// DGELS solves overdetermined or underdetermined real linear systems
  27:      /// involving an M-by-N matrix A, or its transpose, using a QR or LQ
  28:      /// factorization of A.  It is assumed that A has full rank.
  29:      /// 
  30:      /// The following options are provided:
  31:      /// 
  32:      /// 1. If TRANS = 'N' and m .GE. n:  find the least squares solution of
  33:      /// an overdetermined system, i.e., solve the least squares problem
  34:      /// minimize || B - A*X ||.
  35:      /// 
  36:      /// 2. If TRANS = 'N' and m .LT. n:  find the minimum norm solution of
  37:      /// an underdetermined system A * X = B.
  38:      /// 
  39:      /// 3. If TRANS = 'T' and m .GE. n:  find the minimum norm solution of
  40:      /// an undetermined system A**T * X = B.
  41:      /// 
  42:      /// 4. If TRANS = 'T' and m .LT. n:  find the least squares solution of
  43:      /// an overdetermined system, i.e., solve the least squares problem
  44:      /// minimize || B - A**T * X ||.
  45:      /// 
  46:      /// Several right hand side vectors b and solution vectors x can be
  47:      /// handled in a single call; they are stored as the columns of the
  48:      /// M-by-NRHS right hand side matrix B and the N-by-NRHS solution
  49:      /// matrix X.
  50:      /// 
  51:      ///</summary>
  52:      public class DGELS
  53:      {
  54:      
  55:   
  56:          #region Dependencies
  57:          
  58:          LSAME _lsame; ILAENV _ilaenv; DLABAD _dlabad; DLAMCH _dlamch; DLANGE _dlange; DGELQF _dgelqf; DGEQRF _dgeqrf; 
  59:          DLASCL _dlascl;DLASET _dlaset; DORMLQ _dormlq; DORMQR _dormqr; DTRTRS _dtrtrs; XERBLA _xerbla; 
  60:   
  61:          #endregion
  62:   
  63:   
  64:          #region Fields
  65:          
  66:          const double ZERO = 0.0E0; const double ONE = 1.0E0; bool LQUERY = false; bool TPSD = false; int BROW = 0; int I = 0; 
  67:          int IASCL = 0;int IBSCL = 0; int J = 0; int MN = 0; int NB = 0; int SCLLEN = 0; int WSIZE = 0; double ANRM = 0; 
  68:          double BIGNUM = 0;double BNRM = 0; double SMLNUM = 0; double[] RWORK = new double[1]; int offset_rwork = 0; 
  69:   
  70:          #endregion
  71:   
  72:          public DGELS(LSAME lsame, ILAENV ilaenv, DLABAD dlabad, DLAMCH dlamch, DLANGE dlange, DGELQF dgelqf, DGEQRF dgeqrf, DLASCL dlascl, DLASET dlaset, DORMLQ dormlq
  73:                       , DORMQR dormqr, DTRTRS dtrtrs, XERBLA xerbla)
  74:          {
  75:      
  76:   
  77:              #region Set Dependencies
  78:              
  79:              this._lsame = lsame; this._ilaenv = ilaenv; this._dlabad = dlabad; this._dlamch = dlamch; this._dlange = dlange; 
  80:              this._dgelqf = dgelqf;this._dgeqrf = dgeqrf; this._dlascl = dlascl; this._dlaset = dlaset; this._dormlq = dormlq; 
  81:              this._dormqr = dormqr;this._dtrtrs = dtrtrs; this._xerbla = xerbla; 
  82:   
  83:              #endregion
  84:   
  85:          }
  86:      
  87:          public DGELS()
  88:          {
  89:      
  90:   
  91:              #region Dependencies (Initialization)
  92:              
  93:              LSAME lsame = new LSAME();
  94:              IEEECK ieeeck = new IEEECK();
  95:              IPARMQ iparmq = new IPARMQ();
  96:              DLABAD dlabad = new DLABAD();
  97:              DLAMC3 dlamc3 = new DLAMC3();
  98:              DLASSQ dlassq = new DLASSQ();
  99:              XERBLA xerbla = new XERBLA();
 100:              DLAPY2 dlapy2 = new DLAPY2();
 101:              DNRM2 dnrm2 = new DNRM2();
 102:              DSCAL dscal = new DSCAL();
 103:              DCOPY dcopy = new DCOPY();
 104:              ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
 105:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);
 106:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);
 107:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);
 108:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
 109:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
 110:              DLANGE dlange = new DLANGE(dlassq, lsame);
 111:              DGEMV dgemv = new DGEMV(lsame, xerbla);
 112:              DGER dger = new DGER(xerbla);
 113:              DLARF dlarf = new DLARF(dgemv, dger, lsame);
 114:              DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
 115:              DGELQ2 dgelq2 = new DGELQ2(dlarf, dlarfg, xerbla);
 116:              DGEMM dgemm = new DGEMM(lsame, xerbla);
 117:              DTRMM dtrmm = new DTRMM(lsame, xerbla);
 118:              DLARFB dlarfb = new DLARFB(lsame, dcopy, dgemm, dtrmm);
 119:              DTRMV dtrmv = new DTRMV(lsame, xerbla);
 120:              DLARFT dlarft = new DLARFT(dgemv, dtrmv, lsame);
 121:              DGELQF dgelqf = new DGELQF(dgelq2, dlarfb, dlarft, xerbla, ilaenv);
 122:              DGEQR2 dgeqr2 = new DGEQR2(dlarf, dlarfg, xerbla);
 123:              DGEQRF dgeqrf = new DGEQRF(dgeqr2, dlarfb, dlarft, xerbla, ilaenv);
 124:              DLASCL dlascl = new DLASCL(lsame, dlamch, xerbla);
 125:              DLASET dlaset = new DLASET(lsame);
 126:              DORML2 dorml2 = new DORML2(lsame, dlarf, xerbla);
 127:              DORMLQ dormlq = new DORMLQ(lsame, ilaenv, dlarfb, dlarft, dorml2, xerbla);
 128:              DORM2R dorm2r = new DORM2R(lsame, dlarf, xerbla);
 129:              DORMQR dormqr = new DORMQR(lsame, ilaenv, dlarfb, dlarft, dorm2r, xerbla);
 130:              DTRSM dtrsm = new DTRSM(lsame, xerbla);
 131:              DTRTRS dtrtrs = new DTRTRS(lsame, dtrsm, xerbla);
 132:   
 133:              #endregion
 134:   
 135:   
 136:              #region Set Dependencies
 137:              
 138:              this._lsame = lsame; this._ilaenv = ilaenv; this._dlabad = dlabad; this._dlamch = dlamch; this._dlange = dlange; 
 139:              this._dgelqf = dgelqf;this._dgeqrf = dgeqrf; this._dlascl = dlascl; this._dlaset = dlaset; this._dormlq = dormlq; 
 140:              this._dormqr = dormqr;this._dtrtrs = dtrtrs; this._xerbla = xerbla; 
 141:   
 142:              #endregion
 143:   
 144:          }
 145:          /// <summary>
 146:          /// Purpose
 147:          /// =======
 148:          /// 
 149:          /// DGELS solves overdetermined or underdetermined real linear systems
 150:          /// involving an M-by-N matrix A, or its transpose, using a QR or LQ
 151:          /// factorization of A.  It is assumed that A has full rank.
 152:          /// 
 153:          /// The following options are provided:
 154:          /// 
 155:          /// 1. If TRANS = 'N' and m .GE. n:  find the least squares solution of
 156:          /// an overdetermined system, i.e., solve the least squares problem
 157:          /// minimize || B - A*X ||.
 158:          /// 
 159:          /// 2. If TRANS = 'N' and m .LT. n:  find the minimum norm solution of
 160:          /// an underdetermined system A * X = B.
 161:          /// 
 162:          /// 3. If TRANS = 'T' and m .GE. n:  find the minimum norm solution of
 163:          /// an undetermined system A**T * X = B.
 164:          /// 
 165:          /// 4. If TRANS = 'T' and m .LT. n:  find the least squares solution of
 166:          /// an overdetermined system, i.e., solve the least squares problem
 167:          /// minimize || B - A**T * X ||.
 168:          /// 
 169:          /// Several right hand side vectors b and solution vectors x can be
 170:          /// handled in a single call; they are stored as the columns of the
 171:          /// M-by-NRHS right hand side matrix B and the N-by-NRHS solution
 172:          /// matrix X.
 173:          /// 
 174:          ///</summary>
 175:          /// <param name="TRANS">
 176:          /// (input) CHARACTER*1
 177:          /// = 'N': the linear system involves A;
 178:          /// = 'T': the linear system involves A**T.
 179:          ///</param>
 180:          /// <param name="M">
 181:          /// (input) INTEGER
 182:          /// The number of rows of the matrix A.  M .GE. 0.
 183:          ///</param>
 184:          /// <param name="N">
 185:          /// (input) INTEGER
 186:          /// The number of columns of the matrix A.  N .GE. 0.
 187:          ///</param>
 188:          /// <param name="NRHS">
 189:          /// (input) INTEGER
 190:          /// The number of right hand sides, i.e., the number of
 191:          /// columns of the matrices B and X. NRHS .GE.0.
 192:          ///</param>
 193:          /// <param name="A">
 194:          /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)
 195:          /// On entry, the M-by-N matrix A.
 196:          /// On exit,
 197:          /// if M .GE. N, A is overwritten by details of its QR
 198:          /// factorization as returned by DGEQRF;
 199:          /// if M .LT.  N, A is overwritten by details of its LQ
 200:          /// factorization as returned by DGELQF.
 201:          ///</param>
 202:          /// <param name="LDA">
 203:          /// (input) INTEGER
 204:          /// The leading dimension of the array A.  LDA .GE. max(1,M).
 205:          ///</param>
 206:          /// <param name="B">
 207:          /// (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
 208:          /// On entry, the matrix B of right hand side vectors, stored
 209:          /// columnwise; B is M-by-NRHS if TRANS = 'N', or N-by-NRHS
 210:          /// if TRANS = 'T'.
 211:          /// On exit, if INFO = 0, B is overwritten by the solution
 212:          /// vectors, stored columnwise:
 213:          /// if TRANS = 'N' and m .GE. n, rows 1 to n of B contain the least
 214:          /// squares solution vectors; the residual sum of squares for the
 215:          /// solution in each column is given by the sum of squares of
 216:          /// elements N+1 to M in that column;
 217:          /// if TRANS = 'N' and m .LT. n, rows 1 to N of B contain the
 218:          /// minimum norm solution vectors;
 219:          /// if TRANS = 'T' and m .GE. n, rows 1 to M of B contain the
 220:          /// minimum norm solution vectors;
 221:          /// if TRANS = 'T' and m .LT. n, rows 1 to M of B contain the
 222:          /// least squares solution vectors; the residual sum of squares
 223:          /// for the solution in each column is given by the sum of
 224:          /// squares of elements M+1 to N in that column.
 225:          ///</param>
 226:          /// <param name="LDB">
 227:          /// (input) INTEGER
 228:          /// The leading dimension of the array B. LDB .GE. MAX(1,M,N).
 229:          ///</param>
 230:          /// <param name="WORK">
 231:          /// (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 232:          /// On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 233:          ///</param>
 234:          /// <param name="LWORK">
 235:          /// (input) INTEGER
 236:          /// The dimension of the array WORK.
 237:          /// LWORK .GE. max( 1, MN + max( MN, NRHS ) ).
 238:          /// For optimal performance,
 239:          /// LWORK .GE. max( 1, MN + max( MN, NRHS )*NB ).
 240:          /// where MN = min(M,N) and NB is the optimum block size.
 241:          /// 
 242:          /// If LWORK = -1, then a workspace query is assumed; the routine
 243:          /// only calculates the optimal size of the WORK array, returns
 244:          /// this value as the first entry of the WORK array, and no error
 245:          /// message related to LWORK is issued by XERBLA.
 246:          ///</param>
 247:          /// <param name="INFO">
 248:          /// (output) INTEGER
 249:          /// = 0:  successful exit
 250:          /// .LT. 0:  if INFO = -i, the i-th argument had an illegal value
 251:          /// .GT. 0:  if INFO =  i, the i-th diagonal element of the
 252:          /// triangular factor of A is zero, so that A does not have
 253:          /// full rank; the least squares solution could not be
 254:          /// computed.
 255:          ///</param>
 256:          public void Run(string TRANS, int M, int N, int NRHS, ref double[] A, int offset_a, int LDA
 257:                           , ref double[] B, int offset_b, int LDB, ref double[] WORK, int offset_work, int LWORK, ref int INFO)
 258:          {
 259:   
 260:              #region Array Index Correction
 261:              
 262:               int o_a = -1 - LDA + offset_a;  int o_b = -1 - LDB + offset_b;  int o_work = -1 + offset_work; 
 263:   
 264:              #endregion
 265:   
 266:   
 267:              #region Strings
 268:              
 269:              TRANS = TRANS.Substring(0, 1);  
 270:   
 271:              #endregion
 272:   
 273:   
 274:              #region Prolog
 275:              
 276:              // *
 277:              // *  -- LAPACK driver routine (version 3.1) --
 278:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 279:              // *     November 2006
 280:              // *
 281:              // *     .. Scalar Arguments ..
 282:              // *     ..
 283:              // *     .. Array Arguments ..
 284:              // *     ..
 285:              // *
 286:              // *  Purpose
 287:              // *  =======
 288:              // *
 289:              // *  DGELS solves overdetermined or underdetermined real linear systems
 290:              // *  involving an M-by-N matrix A, or its transpose, using a QR or LQ
 291:              // *  factorization of A.  It is assumed that A has full rank.
 292:              // *
 293:              // *  The following options are provided:
 294:              // *
 295:              // *  1. If TRANS = 'N' and m >= n:  find the least squares solution of
 296:              // *     an overdetermined system, i.e., solve the least squares problem
 297:              // *                  minimize || B - A*X ||.
 298:              // *
 299:              // *  2. If TRANS = 'N' and m < n:  find the minimum norm solution of
 300:              // *     an underdetermined system A * X = B.
 301:              // *
 302:              // *  3. If TRANS = 'T' and m >= n:  find the minimum norm solution of
 303:              // *     an undetermined system A**T * X = B.
 304:              // *
 305:              // *  4. If TRANS = 'T' and m < n:  find the least squares solution of
 306:              // *     an overdetermined system, i.e., solve the least squares problem
 307:              // *                  minimize || B - A**T * X ||.
 308:              // *
 309:              // *  Several right hand side vectors b and solution vectors x can be
 310:              // *  handled in a single call; they are stored as the columns of the
 311:              // *  M-by-NRHS right hand side matrix B and the N-by-NRHS solution
 312:              // *  matrix X.
 313:              // *
 314:              // *  Arguments
 315:              // *  =========
 316:              // *
 317:              // *  TRANS   (input) CHARACTER*1
 318:              // *          = 'N': the linear system involves A;
 319:              // *          = 'T': the linear system involves A**T.
 320:              // *
 321:              // *  M       (input) INTEGER
 322:              // *          The number of rows of the matrix A.  M >= 0.
 323:              // *
 324:              // *  N       (input) INTEGER
 325:              // *          The number of columns of the matrix A.  N >= 0.
 326:              // *
 327:              // *  NRHS    (input) INTEGER
 328:              // *          The number of right hand sides, i.e., the number of
 329:              // *          columns of the matrices B and X. NRHS >=0.
 330:              // *
 331:              // *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
 332:              // *          On entry, the M-by-N matrix A.
 333:              // *          On exit,
 334:              // *            if M >= N, A is overwritten by details of its QR
 335:              // *                       factorization as returned by DGEQRF;
 336:              // *            if M <  N, A is overwritten by details of its LQ
 337:              // *                       factorization as returned by DGELQF.
 338:              // *
 339:              // *  LDA     (input) INTEGER
 340:              // *          The leading dimension of the array A.  LDA >= max(1,M).
 341:              // *
 342:              // *  B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
 343:              // *          On entry, the matrix B of right hand side vectors, stored
 344:              // *          columnwise; B is M-by-NRHS if TRANS = 'N', or N-by-NRHS
 345:              // *          if TRANS = 'T'.
 346:              // *          On exit, if INFO = 0, B is overwritten by the solution
 347:              // *          vectors, stored columnwise:
 348:              // *          if TRANS = 'N' and m >= n, rows 1 to n of B contain the least
 349:              // *          squares solution vectors; the residual sum of squares for the
 350:              // *          solution in each column is given by the sum of squares of
 351:              // *          elements N+1 to M in that column;
 352:              // *          if TRANS = 'N' and m < n, rows 1 to N of B contain the
 353:              // *          minimum norm solution vectors;
 354:              // *          if TRANS = 'T' and m >= n, rows 1 to M of B contain the
 355:              // *          minimum norm solution vectors;
 356:              // *          if TRANS = 'T' and m < n, rows 1 to M of B contain the
 357:              // *          least squares solution vectors; the residual sum of squares
 358:              // *          for the solution in each column is given by the sum of
 359:              // *          squares of elements M+1 to N in that column.
 360:              // *
 361:              // *  LDB     (input) INTEGER
 362:              // *          The leading dimension of the array B. LDB >= MAX(1,M,N).
 363:              // *
 364:              // *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 365:              // *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 366:              // *
 367:              // *  LWORK   (input) INTEGER
 368:              // *          The dimension of the array WORK.
 369:              // *          LWORK >= max( 1, MN + max( MN, NRHS ) ).
 370:              // *          For optimal performance,
 371:              // *          LWORK >= max( 1, MN + max( MN, NRHS )*NB ).
 372:              // *          where MN = min(M,N) and NB is the optimum block size.
 373:              // *
 374:              // *          If LWORK = -1, then a workspace query is assumed; the routine
 375:              // *          only calculates the optimal size of the WORK array, returns
 376:              // *          this value as the first entry of the WORK array, and no error
 377:              // *          message related to LWORK is issued by XERBLA.
 378:              // *
 379:              // *  INFO    (output) INTEGER
 380:              // *          = 0:  successful exit
 381:              // *          < 0:  if INFO = -i, the i-th argument had an illegal value
 382:              // *          > 0:  if INFO =  i, the i-th diagonal element of the
 383:              // *                triangular factor of A is zero, so that A does not have
 384:              // *                full rank; the least squares solution could not be
 385:              // *                computed.
 386:              // *
 387:              // *  =====================================================================
 388:              // *
 389:              // *     .. Parameters ..
 390:              // *     ..
 391:              // *     .. Local Scalars ..
 392:              // *     ..
 393:              // *     .. Local Arrays ..
 394:              // *     ..
 395:              // *     .. External Functions ..
 396:              // *     ..
 397:              // *     .. External Subroutines ..
 398:              // *     ..
 399:              // *     .. Intrinsic Functions ..
 400:              //      INTRINSIC          DBLE, MAX, MIN;
 401:              // *     ..
 402:              // *     .. Executable Statements ..
 403:              // *
 404:              // *     Test the input arguments.
 405:              // *
 406:   
 407:              #endregion
 408:   
 409:   
 410:              #region Body
 411:              
 412:              INFO = 0;
 413:              MN = Math.Min(M, N);
 414:              LQUERY = (LWORK ==  - 1);
 415:              if (!(this._lsame.Run(TRANS, "N") || this._lsame.Run(TRANS, "T")))
 416:              {
 417:                  INFO =  - 1;
 418:              }
 419:              else
 420:              {
 421:                  if (M < 0)
 422:                  {
 423:                      INFO =  - 2;
 424:                  }
 425:                  else
 426:                  {
 427:                      if (N < 0)
 428:                      {
 429:                          INFO =  - 3;
 430:                      }
 431:                      else
 432:                      {
 433:                          if (NRHS < 0)
 434:                          {
 435:                              INFO =  - 4;
 436:                          }
 437:                          else
 438:                          {
 439:                              if (LDA < Math.Max(1, M))
 440:                              {
 441:                                  INFO =  - 6;
 442:                              }
 443:                              else
 444:                              {
 445:                                  if (LDB < Math.Max(1, Math.Max(M, N)))
 446:                                  {
 447:                                      INFO =  - 8;
 448:                                  }
 449:                                  else
 450:                                  {
 451:                                      if (LWORK < Math.Max(1, MN + Math.Max(MN, NRHS)) && !LQUERY)
 452:                                      {
 453:                                          INFO =  - 10;
 454:                                      }
 455:                                  }
 456:                              }
 457:                          }
 458:                      }
 459:                  }
 460:              }
 461:              // *
 462:              // *     Figure out optimal block size
 463:              // *
 464:              if (INFO == 0 || INFO ==  - 10)
 465:              {
 466:                  // *
 467:                  TPSD = true;
 468:                  if (this._lsame.Run(TRANS, "N")) TPSD = false;
 469:                  // *
 470:                  if (M >= N)
 471:                  {
 472:                      NB = this._ilaenv.Run(1, "DGEQRF", " ", M, N,  - 1,  - 1);
 473:                      if (TPSD)
 474:                      {
 475:                          NB = Math.Max(NB, this._ilaenv.Run(1, "DORMQR", "LN", M, NRHS, N,  - 1));
 476:                      }
 477:                      else
 478:                      {
 479:                          NB = Math.Max(NB, this._ilaenv.Run(1, "DORMQR", "LT", M, NRHS, N,  - 1));
 480:                      }
 481:                  }
 482:                  else
 483:                  {
 484:                      NB = this._ilaenv.Run(1, "DGELQF", " ", M, N,  - 1,  - 1);
 485:                      if (TPSD)
 486:                      {
 487:                          NB = Math.Max(NB, this._ilaenv.Run(1, "DORMLQ", "LT", N, NRHS, M,  - 1));
 488:                      }
 489:                      else
 490:                      {
 491:                          NB = Math.Max(NB, this._ilaenv.Run(1, "DORMLQ", "LN", N, NRHS, M,  - 1));
 492:                      }
 493:                  }
 494:                  // *
 495:                  WSIZE = Math.Max(1, MN + Math.Max(MN, NRHS) * NB);
 496:                  WORK[1 + o_work] = Convert.ToDouble(WSIZE);
 497:                  // *
 498:              }
 499:              // *
 500:              if (INFO != 0)
 501:              {
 502:                  this._xerbla.Run("DGELS ",  - INFO);
 503:                  return;
 504:              }
 505:              else
 506:              {
 507:                  if (LQUERY)
 508:                  {
 509:                      return;
 510:                  }
 511:              }
 512:              // *
 513:              // *     Quick return if possible
 514:              // *
 515:              if (Math.Min(M, Math.Min(N, NRHS)) == 0)
 516:              {
 517:                  this._dlaset.Run("Full", Math.Max(M, N), NRHS, ZERO, ZERO, ref B, offset_b
 518:                                   , LDB);
 519:                  return;
 520:              }
 521:              // *
 522:              // *     Get machine parameters
 523:              // *
 524:              SMLNUM = this._dlamch.Run("S") / this._dlamch.Run("P");
 525:              BIGNUM = ONE / SMLNUM;
 526:              this._dlabad.Run(ref SMLNUM, ref BIGNUM);
 527:              // *
 528:              // *     Scale A, B if max element outside range [SMLNUM,BIGNUM]
 529:              // *
 530:              ANRM = this._dlange.Run("M", M, N, A, offset_a, LDA, ref RWORK, offset_rwork);
 531:              IASCL = 0;
 532:              if (ANRM > ZERO && ANRM < SMLNUM)
 533:              {
 534:                  // *
 535:                  // *        Scale matrix norm up to SMLNUM
 536:                  // *
 537:                  this._dlascl.Run("G", 0, 0, ANRM, SMLNUM, M
 538:                                   , N, ref A, offset_a, LDA, ref INFO);
 539:                  IASCL = 1;
 540:              }
 541:              else
 542:              {
 543:                  if (ANRM > BIGNUM)
 544:                  {
 545:                      // *
 546:                      // *        Scale matrix norm down to BIGNUM
 547:                      // *
 548:                      this._dlascl.Run("G", 0, 0, ANRM, BIGNUM, M
 549:                                       , N, ref A, offset_a, LDA, ref INFO);
 550:                      IASCL = 2;
 551:                  }
 552:                  else
 553:                  {
 554:                      if (ANRM == ZERO)
 555:                      {
 556:                          // *
 557:                          // *        Matrix all zero. Return zero solution.
 558:                          // *
 559:                          this._dlaset.Run("F", Math.Max(M, N), NRHS, ZERO, ZERO, ref B, offset_b
 560:                                           , LDB);
 561:                          goto LABEL50;
 562:                      }
 563:                  }
 564:              }
 565:              // *
 566:              BROW = M;
 567:              if (TPSD) BROW = N;
 568:              BNRM = this._dlange.Run("M", BROW, NRHS, B, offset_b, LDB, ref RWORK, offset_rwork);
 569:              IBSCL = 0;
 570:              if (BNRM > ZERO && BNRM < SMLNUM)
 571:              {
 572:                  // *
 573:                  // *        Scale matrix norm up to SMLNUM
 574:                  // *
 575:                  this._dlascl.Run("G", 0, 0, BNRM, SMLNUM, BROW
 576:                                   , NRHS, ref B, offset_b, LDB, ref INFO);
 577:                  IBSCL = 1;
 578:              }
 579:              else
 580:              {
 581:                  if (BNRM > BIGNUM)
 582:                  {
 583:                      // *
 584:                      // *        Scale matrix norm down to BIGNUM
 585:                      // *
 586:                      this._dlascl.Run("G", 0, 0, BNRM, BIGNUM, BROW
 587:                                       , NRHS, ref B, offset_b, LDB, ref INFO);
 588:                      IBSCL = 2;
 589:                  }
 590:              }
 591:              // *
 592:              if (M >= N)
 593:              {
 594:                  // *
 595:                  // *        compute QR factorization of A
 596:                  // *
 597:                  this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, 1 + o_work, ref WORK, MN + 1 + o_work
 598:                                   , LWORK - MN, ref INFO);
 599:                  // *
 600:                  // *        workspace at least N, optimally N*NB
 601:                  // *
 602:                  if (!TPSD)
 603:                  {
 604:                      // *
 605:                      // *           Least-Squares Problem min || A * X - B ||
 606:                      // *
 607:                      // *           B(1:M,1:NRHS) := Q' * B(1:M,1:NRHS)
 608:                      // *
 609:                      this._dormqr.Run("Left", "Transpose", M, NRHS, N, ref A, offset_a
 610:                                       , LDA, WORK, 1 + o_work, ref B, offset_b, LDB, ref WORK, MN + 1 + o_work, LWORK - MN
 611:                                       , ref INFO);
 612:                      // *
 613:                      // *           workspace at least NRHS, optimally NRHS*NB
 614:                      // *
 615:                      // *           B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
 616:                      // *
 617:                      this._dtrtrs.Run("Upper", "No transpose", "Non-unit", N, NRHS, A, offset_a
 618:                                       , LDA, ref B, offset_b, LDB, ref INFO);
 619:                      // *
 620:                      if (INFO > 0)
 621:                      {
 622:                          return;
 623:                      }
 624:                      // *
 625:                      SCLLEN = N;
 626:                      // *
 627:                  }
 628:                  else
 629:                  {
 630:                      // *
 631:                      // *           Overdetermined system of equations A' * X = B
 632:                      // *
 633:                      // *           B(1:N,1:NRHS) := inv(R') * B(1:N,1:NRHS)
 634:                      // *
 635:                      this._dtrtrs.Run("Upper", "Transpose", "Non-unit", N, NRHS, A, offset_a
 636:                                       , LDA, ref B, offset_b, LDB, ref INFO);
 637:                      // *
 638:                      if (INFO > 0)
 639:                      {
 640:                          return;
 641:                      }
 642:                      // *
 643:                      // *           B(N+1:M,1:NRHS) = ZERO
 644:                      // *
 645:                      for (J = 1; J <= NRHS; J++)
 646:                      {
 647:                          for (I = N + 1; I <= M; I++)
 648:                          {
 649:                              B[I+J * LDB + o_b] = ZERO;
 650:                          }
 651:                      }
 652:                      // *
 653:                      // *           B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS)
 654:                      // *
 655:                      this._dormqr.Run("Left", "No transpose", M, NRHS, N, ref A, offset_a
 656:                                       , LDA, WORK, 1 + o_work, ref B, offset_b, LDB, ref WORK, MN + 1 + o_work, LWORK - MN
 657:                                       , ref INFO);
 658:                      // *
 659:                      // *           workspace at least NRHS, optimally NRHS*NB
 660:                      // *
 661:                      SCLLEN = M;
 662:                      // *
 663:                  }
 664:                  // *
 665:              }
 666:              else
 667:              {
 668:                  // *
 669:                  // *        Compute LQ factorization of A
 670:                  // *
 671:                  this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, 1 + o_work, ref WORK, MN + 1 + o_work
 672:                                   , LWORK - MN, ref INFO);
 673:                  // *
 674:                  // *        workspace at least M, optimally M*NB.
 675:                  // *
 676:                  if (!TPSD)
 677:                  {
 678:                      // *
 679:                      // *           underdetermined system of equations A * X = B
 680:                      // *
 681:                      // *           B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS)
 682:                      // *
 683:                      this._dtrtrs.Run("Lower", "No transpose", "Non-unit", M, NRHS, A, offset_a
 684:                                       , LDA, ref B, offset_b, LDB, ref INFO);
 685:                      // *
 686:                      if (INFO > 0)
 687:                      {
 688:                          return;
 689:                      }
 690:                      // *
 691:                      // *           B(M+1:N,1:NRHS) = 0
 692:                      // *
 693:                      for (J = 1; J <= NRHS; J++)
 694:                      {
 695:                          for (I = M + 1; I <= N; I++)
 696:                          {
 697:                              B[I+J * LDB + o_b] = ZERO;
 698:                          }
 699:                      }
 700:                      // *
 701:                      // *           B(1:N,1:NRHS) := Q(1:N,:)' * B(1:M,1:NRHS)
 702:                      // *
 703:                      this._dormlq.Run("Left", "Transpose", N, NRHS, M, ref A, offset_a
 704:                                       , LDA, WORK, 1 + o_work, ref B, offset_b, LDB, ref WORK, MN + 1 + o_work, LWORK - MN
 705:                                       , ref INFO);
 706:                      // *
 707:                      // *           workspace at least NRHS, optimally NRHS*NB
 708:                      // *
 709:                      SCLLEN = N;
 710:                      // *
 711:                  }
 712:                  else
 713:                  {
 714:                      // *
 715:                      // *           overdetermined system min || A' * X - B ||
 716:                      // *
 717:                      // *           B(1:N,1:NRHS) := Q * B(1:N,1:NRHS)
 718:                      // *
 719:                      this._dormlq.Run("Left", "No transpose", N, NRHS, M, ref A, offset_a
 720:                                       , LDA, WORK, 1 + o_work, ref B, offset_b, LDB, ref WORK, MN + 1 + o_work, LWORK - MN
 721:                                       , ref INFO);
 722:                      // *
 723:                      // *           workspace at least NRHS, optimally NRHS*NB
 724:                      // *
 725:                      // *           B(1:M,1:NRHS) := inv(L') * B(1:M,1:NRHS)
 726:                      // *
 727:                      this._dtrtrs.Run("Lower", "Transpose", "Non-unit", M, NRHS, A, offset_a
 728:                                       , LDA, ref B, offset_b, LDB, ref INFO);
 729:                      // *
 730:                      if (INFO > 0)
 731:                      {
 732:                          return;
 733:                      }
 734:                      // *
 735:                      SCLLEN = M;
 736:                      // *
 737:                  }
 738:                  // *
 739:              }
 740:              // *
 741:              // *     Undo scaling
 742:              // *
 743:              if (IASCL == 1)
 744:              {
 745:                  this._dlascl.Run("G", 0, 0, ANRM, SMLNUM, SCLLEN
 746:                                   , NRHS, ref B, offset_b, LDB, ref INFO);
 747:              }
 748:              else
 749:              {
 750:                  if (IASCL == 2)
 751:                  {
 752:                      this._dlascl.Run("G", 0, 0, ANRM, BIGNUM, SCLLEN
 753:                                       , NRHS, ref B, offset_b, LDB, ref INFO);
 754:                  }
 755:              }
 756:              if (IBSCL == 1)
 757:              {
 758:                  this._dlascl.Run("G", 0, 0, SMLNUM, BNRM, SCLLEN
 759:                                   , NRHS, ref B, offset_b, LDB, ref INFO);
 760:              }
 761:              else
 762:              {
 763:                  if (IBSCL == 2)
 764:                  {
 765:                      this._dlascl.Run("G", 0, 0, BIGNUM, BNRM, SCLLEN
 766:                                       , NRHS, ref B, offset_b, LDB, ref INFO);
 767:                  }
 768:              }
 769:              // *
 770:          LABEL50:;
 771:              WORK[1 + o_work] = Convert.ToDouble(WSIZE);
 772:              // *
 773:              return;
 774:              // *
 775:              // *     End of DGELS
 776:              // *
 777:   
 778:              #endregion
 779:   
 780:          }
 781:      }
 782:  }