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: }