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