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: /// DGGLSE solves the linear equality-constrained least squares (LSE)
27: /// problem:
28: ///
29: /// minimize || c - A*x ||_2 subject to B*x = d
30: ///
31: /// where A is an M-by-N matrix, B is a P-by-N matrix, c is a given
32: /// M-vector, and d is a given P-vector. It is assumed that
33: /// P .LE. N .LE. M+P, and
34: ///
35: /// rank(B) = P and rank( (A) ) = N.
36: /// ( (B) )
37: ///
38: /// These conditions ensure that the LSE problem has a unique solution,
39: /// which is obtained using a generalized RQ factorization of the
40: /// matrices (B, A) given by
41: ///
42: /// B = (0 R)*Q, A = Z*T*Q.
43: ///
44: ///</summary>
45: public class DGGLSE
46: {
47:
48:
49: #region Dependencies
50:
51: DAXPY _daxpy; DCOPY _dcopy; DGEMV _dgemv; DGGRQF _dggrqf; DORMQR _dormqr; DORMRQ _dormrq; DTRMV _dtrmv; DTRTRS _dtrtrs;
52: XERBLA _xerbla;ILAENV _ilaenv;
53:
54: #endregion
55:
56:
57: #region Fields
58:
59: const double ONE = 1.0E+0; bool LQUERY = false; int LOPT = 0; int LWKMIN = 0; int LWKOPT = 0; int MN = 0; int NB = 0;
60: int NB1 = 0;int NB2 = 0; int NB3 = 0; int NB4 = 0; int NR = 0;
61:
62: #endregion
63:
64: public DGGLSE(DAXPY daxpy, DCOPY dcopy, DGEMV dgemv, DGGRQF dggrqf, DORMQR dormqr, DORMRQ dormrq, DTRMV dtrmv, DTRTRS dtrtrs, XERBLA xerbla, ILAENV ilaenv)
65: {
66:
67:
68: #region Set Dependencies
69:
70: this._daxpy = daxpy; this._dcopy = dcopy; this._dgemv = dgemv; this._dggrqf = dggrqf; this._dormqr = dormqr;
71: this._dormrq = dormrq;this._dtrmv = dtrmv; this._dtrtrs = dtrtrs; this._xerbla = xerbla; this._ilaenv = ilaenv;
72:
73: #endregion
74:
75: }
76:
77: public DGGLSE()
78: {
79:
80:
81: #region Dependencies (Initialization)
82:
83: DAXPY daxpy = new DAXPY();
84: DCOPY dcopy = new DCOPY();
85: LSAME lsame = new LSAME();
86: XERBLA xerbla = new XERBLA();
87: DLAMC3 dlamc3 = new DLAMC3();
88: DLAPY2 dlapy2 = new DLAPY2();
89: DNRM2 dnrm2 = new DNRM2();
90: DSCAL dscal = new DSCAL();
91: IEEECK ieeeck = new IEEECK();
92: IPARMQ iparmq = new IPARMQ();
93: DGEMV dgemv = new DGEMV(lsame, xerbla);
94: DGER dger = new DGER(xerbla);
95: DLARF dlarf = new DLARF(dgemv, dger, lsame);
96: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
97: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
98: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
99: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
100: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
101: DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
102: DGEQR2 dgeqr2 = new DGEQR2(dlarf, dlarfg, xerbla);
103: DGEMM dgemm = new DGEMM(lsame, xerbla);
104: DTRMM dtrmm = new DTRMM(lsame, xerbla);
105: DLARFB dlarfb = new DLARFB(lsame, dcopy, dgemm, dtrmm);
106: DTRMV dtrmv = new DTRMV(lsame, xerbla);
107: DLARFT dlarft = new DLARFT(dgemv, dtrmv, lsame);
108: ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
109: DGEQRF dgeqrf = new DGEQRF(dgeqr2, dlarfb, dlarft, xerbla, ilaenv);
110: DGERQ2 dgerq2 = new DGERQ2(dlarf, dlarfg, xerbla);
111: DGERQF dgerqf = new DGERQF(dgerq2, dlarfb, dlarft, xerbla, ilaenv);
112: DORMR2 dormr2 = new DORMR2(lsame, dlarf, xerbla);
113: DORMRQ dormrq = new DORMRQ(lsame, ilaenv, dlarfb, dlarft, dormr2, xerbla);
114: DGGRQF dggrqf = new DGGRQF(dgeqrf, dgerqf, dormrq, xerbla, ilaenv);
115: DORM2R dorm2r = new DORM2R(lsame, dlarf, xerbla);
116: DORMQR dormqr = new DORMQR(lsame, ilaenv, dlarfb, dlarft, dorm2r, xerbla);
117: DTRSM dtrsm = new DTRSM(lsame, xerbla);
118: DTRTRS dtrtrs = new DTRTRS(lsame, dtrsm, xerbla);
119:
120: #endregion
121:
122:
123: #region Set Dependencies
124:
125: this._daxpy = daxpy; this._dcopy = dcopy; this._dgemv = dgemv; this._dggrqf = dggrqf; this._dormqr = dormqr;
126: this._dormrq = dormrq;this._dtrmv = dtrmv; this._dtrtrs = dtrtrs; this._xerbla = xerbla; this._ilaenv = ilaenv;
127:
128: #endregion
129:
130: }
131: /// <summary>
132: /// Purpose
133: /// =======
134: ///
135: /// DGGLSE solves the linear equality-constrained least squares (LSE)
136: /// problem:
137: ///
138: /// minimize || c - A*x ||_2 subject to B*x = d
139: ///
140: /// where A is an M-by-N matrix, B is a P-by-N matrix, c is a given
141: /// M-vector, and d is a given P-vector. It is assumed that
142: /// P .LE. N .LE. M+P, and
143: ///
144: /// rank(B) = P and rank( (A) ) = N.
145: /// ( (B) )
146: ///
147: /// These conditions ensure that the LSE problem has a unique solution,
148: /// which is obtained using a generalized RQ factorization of the
149: /// matrices (B, A) given by
150: ///
151: /// B = (0 R)*Q, A = Z*T*Q.
152: ///
153: ///</summary>
154: /// <param name="M">
155: /// (input) INTEGER
156: /// The number of rows of the matrix A. M .GE. 0.
157: ///</param>
158: /// <param name="N">
159: /// (input) INTEGER
160: /// The number of columns of the matrices A and B. N .GE. 0.
161: ///</param>
162: /// <param name="P">
163: /// (input) INTEGER
164: /// The number of rows of the matrix B. 0 .LE. P .LE. N .LE. M+P.
165: ///</param>
166: /// <param name="A">
167: /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)
168: /// On entry, the M-by-N matrix A.
169: /// On exit, the elements on and above the diagonal of the array
170: /// contain the min(M,N)-by-N upper trapezoidal matrix T.
171: ///</param>
172: /// <param name="LDA">
173: /// (input) INTEGER
174: /// The leading dimension of the array A. LDA .GE. max(1,M).
175: ///</param>
176: /// <param name="B">
177: /// = (0 R)*Q, A = Z*T*Q.
178: ///</param>
179: /// <param name="LDB">
180: /// (input) INTEGER
181: /// The leading dimension of the array B. LDB .GE. max(1,P).
182: ///</param>
183: /// <param name="C">
184: /// (input/output) DOUBLE PRECISION array, dimension (M)
185: /// On entry, C contains the right hand side vector for the
186: /// least squares part of the LSE problem.
187: /// On exit, the residual sum of squares for the solution
188: /// is given by the sum of squares of elements N-P+1 to M of
189: /// vector C.
190: ///</param>
191: /// <param name="D">
192: /// (input/output) DOUBLE PRECISION array, dimension (P)
193: /// On entry, D contains the right hand side vector for the
194: /// constrained equation.
195: /// On exit, D is destroyed.
196: ///</param>
197: /// <param name="X">
198: /// (output) DOUBLE PRECISION array, dimension (N)
199: /// On exit, X is the solution of the LSE problem.
200: ///</param>
201: /// <param name="WORK">
202: /// (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
203: /// On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
204: ///</param>
205: /// <param name="LWORK">
206: /// (input) INTEGER
207: /// The dimension of the array WORK. LWORK .GE. max(1,M+N+P).
208: /// For optimum performance LWORK .GE. P+min(M,N)+max(M,N)*NB,
209: /// where NB is an upper bound for the optimal blocksizes for
210: /// DGEQRF, SGERQF, DORMQR and SORMRQ.
211: ///
212: /// If LWORK = -1, then a workspace query is assumed; the routine
213: /// only calculates the optimal size of the WORK array, returns
214: /// this value as the first entry of the WORK array, and no error
215: /// message related to LWORK is issued by XERBLA.
216: ///</param>
217: /// <param name="INFO">
218: /// (output) INTEGER
219: /// = 0: successful exit.
220: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value.
221: /// = 1: the upper triangular factor R associated with B in the
222: /// generalized RQ factorization of the pair (B, A) is
223: /// singular, so that rank(B) .LT. P; the least squares
224: /// solution could not be computed.
225: /// = 2: the (N-P) by (N-P) part of the upper trapezoidal factor
226: /// T associated with A in the generalized RQ factorization
227: /// of the pair (B, A) is singular, so that
228: /// rank( (A) ) .LT. N; the least squares solution could not
229: /// ( (B) )
230: /// be computed.
231: ///</param>
232: public void Run(int M, int N, int P, ref double[] A, int offset_a, int LDA, ref double[] B, int offset_b
233: , int LDB, ref double[] C, int offset_c, ref double[] D, int offset_d, ref double[] X, int offset_x, ref double[] WORK, int offset_work, int LWORK
234: , ref int INFO)
235: {
236:
237: #region Array Index Correction
238:
239: int o_a = -1 - LDA + offset_a; int o_b = -1 - LDB + offset_b; int o_c = -1 + offset_c; int o_d = -1 + offset_d;
240: int o_x = -1 + offset_x; int o_work = -1 + offset_work;
241:
242: #endregion
243:
244:
245: #region Prolog
246:
247: // *
248: // * -- LAPACK driver routine (version 3.1) --
249: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
250: // * November 2006
251: // *
252: // * .. Scalar Arguments ..
253: // * ..
254: // * .. Array Arguments ..
255: // * ..
256: // *
257: // * Purpose
258: // * =======
259: // *
260: // * DGGLSE solves the linear equality-constrained least squares (LSE)
261: // * problem:
262: // *
263: // * minimize || c - A*x ||_2 subject to B*x = d
264: // *
265: // * where A is an M-by-N matrix, B is a P-by-N matrix, c is a given
266: // * M-vector, and d is a given P-vector. It is assumed that
267: // * P <= N <= M+P, and
268: // *
269: // * rank(B) = P and rank( (A) ) = N.
270: // * ( (B) )
271: // *
272: // * These conditions ensure that the LSE problem has a unique solution,
273: // * which is obtained using a generalized RQ factorization of the
274: // * matrices (B, A) given by
275: // *
276: // * B = (0 R)*Q, A = Z*T*Q.
277: // *
278: // * Arguments
279: // * =========
280: // *
281: // * M (input) INTEGER
282: // * The number of rows of the matrix A. M >= 0.
283: // *
284: // * N (input) INTEGER
285: // * The number of columns of the matrices A and B. N >= 0.
286: // *
287: // * P (input) INTEGER
288: // * The number of rows of the matrix B. 0 <= P <= N <= M+P.
289: // *
290: // * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
291: // * On entry, the M-by-N matrix A.
292: // * On exit, the elements on and above the diagonal of the array
293: // * contain the min(M,N)-by-N upper trapezoidal matrix T.
294: // *
295: // * LDA (input) INTEGER
296: // * The leading dimension of the array A. LDA >= max(1,M).
297: // *
298: // * B (input/output) DOUBLE PRECISION array, dimension (LDB,N)
299: // * On entry, the P-by-N matrix B.
300: // * On exit, the upper triangle of the subarray B(1:P,N-P+1:N)
301: // * contains the P-by-P upper triangular matrix R.
302: // *
303: // * LDB (input) INTEGER
304: // * The leading dimension of the array B. LDB >= max(1,P).
305: // *
306: // * C (input/output) DOUBLE PRECISION array, dimension (M)
307: // * On entry, C contains the right hand side vector for the
308: // * least squares part of the LSE problem.
309: // * On exit, the residual sum of squares for the solution
310: // * is given by the sum of squares of elements N-P+1 to M of
311: // * vector C.
312: // *
313: // * D (input/output) DOUBLE PRECISION array, dimension (P)
314: // * On entry, D contains the right hand side vector for the
315: // * constrained equation.
316: // * On exit, D is destroyed.
317: // *
318: // * X (output) DOUBLE PRECISION array, dimension (N)
319: // * On exit, X is the solution of the LSE problem.
320: // *
321: // * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
322: // * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
323: // *
324: // * LWORK (input) INTEGER
325: // * The dimension of the array WORK. LWORK >= max(1,M+N+P).
326: // * For optimum performance LWORK >= P+min(M,N)+max(M,N)*NB,
327: // * where NB is an upper bound for the optimal blocksizes for
328: // * DGEQRF, SGERQF, DORMQR and SORMRQ.
329: // *
330: // * If LWORK = -1, then a workspace query is assumed; the routine
331: // * only calculates the optimal size of the WORK array, returns
332: // * this value as the first entry of the WORK array, and no error
333: // * message related to LWORK is issued by XERBLA.
334: // *
335: // * INFO (output) INTEGER
336: // * = 0: successful exit.
337: // * < 0: if INFO = -i, the i-th argument had an illegal value.
338: // * = 1: the upper triangular factor R associated with B in the
339: // * generalized RQ factorization of the pair (B, A) is
340: // * singular, so that rank(B) < P; the least squares
341: // * solution could not be computed.
342: // * = 2: the (N-P) by (N-P) part of the upper trapezoidal factor
343: // * T associated with A in the generalized RQ factorization
344: // * of the pair (B, A) is singular, so that
345: // * rank( (A) ) < N; the least squares solution could not
346: // * ( (B) )
347: // * be computed.
348: // *
349: // * =====================================================================
350: // *
351: // * .. Parameters ..
352: // * ..
353: // * .. Local Scalars ..
354: // * ..
355: // * .. External Subroutines ..
356: // * ..
357: // * .. External Functions ..
358: // * ..
359: // * .. Intrinsic Functions ..
360: // INTRINSIC INT, MAX, MIN;
361: // * ..
362: // * .. Executable Statements ..
363: // *
364: // * Test the input parameters
365: // *
366:
367: #endregion
368:
369:
370: #region Body
371:
372: INFO = 0;
373: MN = Math.Min(M, N);
374: LQUERY = (LWORK == - 1);
375: if (M < 0)
376: {
377: INFO = - 1;
378: }
379: else
380: {
381: if (N < 0)
382: {
383: INFO = - 2;
384: }
385: else
386: {
387: if (P < 0 || P > N || P < N - M)
388: {
389: INFO = - 3;
390: }
391: else
392: {
393: if (LDA < Math.Max(1, M))
394: {
395: INFO = - 5;
396: }
397: else
398: {
399: if (LDB < Math.Max(1, P))
400: {
401: INFO = - 7;
402: }
403: }
404: }
405: }
406: }
407: // *
408: // * Calculate workspace
409: // *
410: if (INFO == 0)
411: {
412: if (N == 0)
413: {
414: LWKMIN = 1;
415: LWKOPT = 1;
416: }
417: else
418: {
419: NB1 = this._ilaenv.Run(1, "DGEQRF", " ", M, N, - 1, - 1);
420: NB2 = this._ilaenv.Run(1, "DGERQF", " ", M, N, - 1, - 1);
421: NB3 = this._ilaenv.Run(1, "DORMQR", " ", M, N, P, - 1);
422: NB4 = this._ilaenv.Run(1, "DORMRQ", " ", M, N, P, - 1);
423: NB = Math.Max(NB1, Math.Max(NB2, Math.Max(NB3, NB4)));
424: LWKMIN = M + N + P;
425: LWKOPT = P + MN + Math.Max(M, N) * NB;
426: }
427: WORK[1 + o_work] = LWKOPT;
428: // *
429: if (LWORK < LWKMIN && !LQUERY)
430: {
431: INFO = - 12;
432: }
433: }
434: // *
435: if (INFO != 0)
436: {
437: this._xerbla.Run("DGGLSE", - INFO);
438: return;
439: }
440: else
441: {
442: if (LQUERY)
443: {
444: return;
445: }
446: }
447: // *
448: // * Quick return if possible
449: // *
450: if (N == 0) return;
451: // *
452: // * Compute the GRQ factorization of matrices B and A:
453: // *
454: // * B*Q' = ( 0 T12 ) P Z'*A*Q' = ( R11 R12 ) N-P
455: // * N-P P ( 0 R22 ) M+P-N
456: // * N-P P
457: // *
458: // * where T12 and R11 are upper triangular, and Q and Z are
459: // * orthogonal.
460: // *
461: this._dggrqf.Run(P, M, N, ref B, offset_b, LDB, ref WORK, offset_work
462: , ref A, offset_a, LDA, ref WORK, P + 1 + o_work, ref WORK, P + MN + 1 + o_work, LWORK - P - MN, ref INFO);
463: LOPT = (int)WORK[P + MN + 1 + o_work];
464: // *
465: // * Update c = Z'*c = ( c1 ) N-P
466: // * ( c2 ) M+P-N
467: // *
468: this._dormqr.Run("Left", "Transpose", M, 1, MN, ref A, offset_a
469: , LDA, WORK, P + 1 + o_work, ref C, offset_c, Math.Max(1, M), ref WORK, P + MN + 1 + o_work, LWORK - P - MN
470: , ref INFO);
471: LOPT = Math.Max(LOPT, Convert.ToInt32(Math.Truncate(WORK[P + MN + 1 + o_work])));
472: // *
473: // * Solve T12*x2 = d for x2
474: // *
475: if (P > 0)
476: {
477: this._dtrtrs.Run("Upper", "No transpose", "Non-unit", P, 1, B, 1+(N - P + 1) * LDB + o_b
478: , LDB, ref D, offset_d, P, ref INFO);
479: // *
480: if (INFO > 0)
481: {
482: INFO = 1;
483: return;
484: }
485: // *
486: // * Put the solution in X
487: // *
488: this._dcopy.Run(P, D, offset_d, 1, ref X, N - P + 1 + o_x, 1);
489: // *
490: // * Update c1
491: // *
492: this._dgemv.Run("No transpose", N - P, P, - ONE, A, 1+(N - P + 1) * LDA + o_a, LDA
493: , D, offset_d, 1, ONE, ref C, offset_c, 1);
494: }
495: // *
496: // * Solve R11*x1 = c1 for x1
497: // *
498: if (N > P)
499: {
500: this._dtrtrs.Run("Upper", "No transpose", "Non-unit", N - P, 1, A, offset_a
501: , LDA, ref C, offset_c, N - P, ref INFO);
502: // *
503: if (INFO > 0)
504: {
505: INFO = 2;
506: return;
507: }
508: // *
509: // * Put the solutions in X
510: // *
511: this._dcopy.Run(N - P, C, offset_c, 1, ref X, offset_x, 1);
512: }
513: // *
514: // * Compute the residual vector:
515: // *
516: if (M < N)
517: {
518: NR = M + P - N;
519: if (NR > 0)
520: {
521: this._dgemv.Run("No transpose", NR, N - M, - ONE, A, N - P + 1+(M + 1) * LDA + o_a, LDA
522: , D, NR + 1 + o_d, 1, ONE, ref C, N - P + 1 + o_c, 1);
523: }
524: }
525: else
526: {
527: NR = P;
528: }
529: if (NR > 0)
530: {
531: this._dtrmv.Run("Upper", "No transpose", "Non unit", NR, A, N - P + 1+(N - P + 1) * LDA + o_a, LDA
532: , ref D, offset_d, 1);
533: this._daxpy.Run(NR, - ONE, D, offset_d, 1, ref C, N - P + 1 + o_c, 1);
534: }
535: // *
536: // * Backward transformation x = Q'*x
537: // *
538: this._dormrq.Run("Left", "Transpose", N, 1, P, ref B, offset_b
539: , LDB, WORK, 1 + o_work, ref X, offset_x, N, ref WORK, P + MN + 1 + o_work, LWORK - P - MN
540: , ref INFO);
541: WORK[1 + o_work] = P + MN + Math.Max(LOPT, Convert.ToInt32(Math.Truncate(WORK[P + MN + 1 + o_work])));
542: // *
543: return;
544: // *
545: // * End of DGGLSE
546: // *
547:
548: #endregion
549:
550: }
551: }
552: }