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();
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;
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:  }