1: SUBROUTINE DGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK, LWORK,
2: $ INFO )
3: *
4: * -- LAPACK driver routine (version 3.2) --
5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7: * November 2006
8: *
9: * .. Scalar Arguments ..
10: INTEGER INFO, LDA, LDB, LWORK, M, N, P
11: * ..
12: * .. Array Arguments ..
13: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), D( * ), WORK( * ),
14: $ X( * ), Y( * )
15: * ..
16: *
17: * Purpose
18: * =======
19: *
20: * DGGGLM solves a general Gauss-Markov linear model (GLM) problem:
21: *
22: * minimize || y ||_2 subject to d = A*x + B*y
23: * x
24: *
25: * where A is an N-by-M matrix, B is an N-by-P matrix, and d is a
26: * given N-vector. It is assumed that M <= N <= M+P, and
27: *
28: * rank(A) = M and rank( A B ) = N.
29: *
30: * Under these assumptions, the constrained equation is always
31: * consistent, and there is a unique solution x and a minimal 2-norm
32: * solution y, which is obtained using a generalized QR factorization
33: * of the matrices (A, B) given by
34: *
35: * A = Q*(R), B = Q*T*Z.
36: * (0)
37: *
38: * In particular, if matrix B is square nonsingular, then the problem
39: * GLM is equivalent to the following weighted linear least squares
40: * problem
41: *
42: * minimize || inv(B)*(d-A*x) ||_2
43: * x
44: *
45: * where inv(B) denotes the inverse of B.
46: *
47: * Arguments
48: * =========
49: *
50: * N (input) INTEGER
51: * The number of rows of the matrices A and B. N >= 0.
52: *
53: * M (input) INTEGER
54: * The number of columns of the matrix A. 0 <= M <= N.
55: *
56: * P (input) INTEGER
57: * The number of columns of the matrix B. P >= N-M.
58: *
59: * A (input/output) DOUBLE PRECISION array, dimension (LDA,M)
60: * On entry, the N-by-M matrix A.
61: * On exit, the upper triangular part of the array A contains
62: * the M-by-M upper triangular matrix R.
63: *
64: * LDA (input) INTEGER
65: * The leading dimension of the array A. LDA >= max(1,N).
66: *
67: * B (input/output) DOUBLE PRECISION array, dimension (LDB,P)
68: * On entry, the N-by-P matrix B.
69: * On exit, if N <= P, the upper triangle of the subarray
70: * B(1:N,P-N+1:P) contains the N-by-N upper triangular matrix T;
71: * if N > P, the elements on and above the (N-P)th subdiagonal
72: * contain the N-by-P upper trapezoidal matrix T.
73: *
74: * LDB (input) INTEGER
75: * The leading dimension of the array B. LDB >= max(1,N).
76: *
77: * D (input/output) DOUBLE PRECISION array, dimension (N)
78: * On entry, D is the left hand side of the GLM equation.
79: * On exit, D is destroyed.
80: *
81: * X (output) DOUBLE PRECISION array, dimension (M)
82: * Y (output) DOUBLE PRECISION array, dimension (P)
83: * On exit, X and Y are the solutions of the GLM problem.
84: *
85: * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
86: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
87: *
88: * LWORK (input) INTEGER
89: * The dimension of the array WORK. LWORK >= max(1,N+M+P).
90: * For optimum performance, LWORK >= M+min(N,P)+max(N,P)*NB,
91: * where NB is an upper bound for the optimal blocksizes for
92: * DGEQRF, SGERQF, DORMQR and SORMRQ.
93: *
94: * If LWORK = -1, then a workspace query is assumed; the routine
95: * only calculates the optimal size of the WORK array, returns
96: * this value as the first entry of the WORK array, and no error
97: * message related to LWORK is issued by XERBLA.
98: *
99: * INFO (output) INTEGER
100: * = 0: successful exit.
101: * < 0: if INFO = -i, the i-th argument had an illegal value.
102: * = 1: the upper triangular factor R associated with A in the
103: * generalized QR factorization of the pair (A, B) is
104: * singular, so that rank(A) < M; the least squares
105: * solution could not be computed.
106: * = 2: the bottom (N-M) by (N-M) part of the upper trapezoidal
107: * factor T associated with B in the generalized QR
108: * factorization of the pair (A, B) is singular, so that
109: * rank( A B ) < N; the least squares solution could not
110: * be computed.
111: *
112: * ===================================================================
113: *
114: * .. Parameters ..
115: DOUBLE PRECISION ZERO, ONE
116: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
117: * ..
118: * .. Local Scalars ..
119: LOGICAL LQUERY
120: INTEGER I, LOPT, LWKMIN, LWKOPT, NB, NB1, NB2, NB3,
121: $ NB4, NP
122: * ..
123: * .. External Subroutines ..
124: EXTERNAL DCOPY, DGEMV, DGGQRF, DORMQR, DORMRQ, DTRTRS,
125: $ XERBLA
126: * ..
127: * .. External Functions ..
128: INTEGER ILAENV
129: EXTERNAL ILAENV
130: * ..
131: * .. Intrinsic Functions ..
132: INTRINSIC INT, MAX, MIN
133: * ..
134: * .. Executable Statements ..
135: *
136: * Test the input parameters
137: *
138: INFO = 0
139: NP = MIN( N, P )
140: LQUERY = ( LWORK.EQ.-1 )
141: IF( N.LT.0 ) THEN
142: INFO = -1
143: ELSE IF( M.LT.0 .OR. M.GT.N ) THEN
144: INFO = -2
145: ELSE IF( P.LT.0 .OR. P.LT.N-M ) THEN
146: INFO = -3
147: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
148: INFO = -5
149: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
150: INFO = -7
151: END IF
152: *
153: * Calculate workspace
154: *
155: IF( INFO.EQ.0) THEN
156: IF( N.EQ.0 ) THEN
157: LWKMIN = 1
158: LWKOPT = 1
159: ELSE
160: NB1 = ILAENV( 1, 'DGEQRF', ' ', N, M, -1, -1 )
161: NB2 = ILAENV( 1, 'DGERQF', ' ', N, M, -1, -1 )
162: NB3 = ILAENV( 1, 'DORMQR', ' ', N, M, P, -1 )
163: NB4 = ILAENV( 1, 'DORMRQ', ' ', N, M, P, -1 )
164: NB = MAX( NB1, NB2, NB3, NB4 )
165: LWKMIN = M + N + P
166: LWKOPT = M + NP + MAX( N, P )*NB
167: END IF
168: WORK( 1 ) = LWKOPT
169: *
170: IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
171: INFO = -12
172: END IF
173: END IF
174: *
175: IF( INFO.NE.0 ) THEN
176: CALL XERBLA( 'DGGGLM', -INFO )
177: RETURN
178: ELSE IF( LQUERY ) THEN
179: RETURN
180: END IF
181: *
182: * Quick return if possible
183: *
184: IF( N.EQ.0 )
185: $ RETURN
186: *
187: * Compute the GQR factorization of matrices A and B:
188: *
189: * Q'*A = ( R11 ) M, Q'*B*Z' = ( T11 T12 ) M
190: * ( 0 ) N-M ( 0 T22 ) N-M
191: * M M+P-N N-M
192: *
193: * where R11 and T22 are upper triangular, and Q and Z are
194: * orthogonal.
195: *
196: CALL DGGQRF( N, M, P, A, LDA, WORK, B, LDB, WORK( M+1 ),
197: $ WORK( M+NP+1 ), LWORK-M-NP, INFO )
198: LOPT = WORK( M+NP+1 )
199: *
200: * Update left-hand-side vector d = Q'*d = ( d1 ) M
201: * ( d2 ) N-M
202: *
203: CALL DORMQR( 'Left', 'Transpose', N, 1, M, A, LDA, WORK, D,
204: $ MAX( 1, N ), WORK( M+NP+1 ), LWORK-M-NP, INFO )
205: LOPT = MAX( LOPT, INT( WORK( M+NP+1 ) ) )
206: *
207: * Solve T22*y2 = d2 for y2
208: *
209: IF( N.GT.M ) THEN
210: CALL DTRTRS( 'Upper', 'No transpose', 'Non unit', N-M, 1,
211: $ B( M+1, M+P-N+1 ), LDB, D( M+1 ), N-M, INFO )
212: *
213: IF( INFO.GT.0 ) THEN
214: INFO = 1
215: RETURN
216: END IF
217: *
218: CALL DCOPY( N-M, D( M+1 ), 1, Y( M+P-N+1 ), 1 )
219: END IF
220: *
221: * Set y1 = 0
222: *
223: DO 10 I = 1, M + P - N
224: Y( I ) = ZERO
225: 10 CONTINUE
226: *
227: * Update d1 = d1 - T12*y2
228: *
229: CALL DGEMV( 'No transpose', M, N-M, -ONE, B( 1, M+P-N+1 ), LDB,
230: $ Y( M+P-N+1 ), 1, ONE, D, 1 )
231: *
232: * Solve triangular system: R11*x = d1
233: *
234: IF( M.GT.0 ) THEN
235: CALL DTRTRS( 'Upper', 'No Transpose', 'Non unit', M, 1, A, LDA,
236: $ D, M, INFO )
237: *
238: IF( INFO.GT.0 ) THEN
239: INFO = 2
240: RETURN
241: END IF
242: *
243: * Copy D to X
244: *
245: CALL DCOPY( M, D, 1, X, 1 )
246: END IF
247: *
248: * Backward transformation y = Z'*y
249: *
250: CALL DORMRQ( 'Left', 'Transpose', P, 1, NP,
251: $ B( MAX( 1, N-P+1 ), 1 ), LDB, WORK( M+1 ), Y,
252: $ MAX( 1, P ), WORK( M+NP+1 ), LWORK-M-NP, INFO )
253: WORK( 1 ) = M + NP + MAX( LOPT, INT( WORK( M+NP+1 ) ) )
254: *
255: RETURN
256: *
257: * End of DGGGLM
258: *
259: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dggglm.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack