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