Annotation of rpl/lapack/lapack/zggglm.f, revision 1.5
1.1 bertrand 1: SUBROUTINE ZGGGLM( 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: COMPLEX*16 A( LDA, * ), B( LDB, * ), D( * ), WORK( * ),
14: $ X( * ), Y( * )
15: * ..
16: *
17: * Purpose
18: * =======
19: *
20: * ZGGGLM 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) COMPLEX*16 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) COMPLEX*16 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) COMPLEX*16 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) COMPLEX*16 array, dimension (M)
82: * Y (output) COMPLEX*16 array, dimension (P)
83: * On exit, X and Y are the solutions of the GLM problem.
84: *
85: * WORK (workspace/output) COMPLEX*16 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: * ZGEQRF, ZGERQF, ZUNMQR and ZUNMRQ.
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: COMPLEX*16 CZERO, CONE
116: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
117: $ CONE = ( 1.0D+0, 0.0D+0 ) )
118: * ..
119: * .. Local Scalars ..
120: LOGICAL LQUERY
121: INTEGER I, LOPT, LWKMIN, LWKOPT, NB, NB1, NB2, NB3,
122: $ NB4, NP
123: * ..
124: * .. External Subroutines ..
125: EXTERNAL XERBLA, ZCOPY, ZGEMV, ZGGQRF, ZTRTRS, ZUNMQR,
126: $ ZUNMRQ
127: * ..
128: * .. External Functions ..
129: INTEGER ILAENV
130: EXTERNAL ILAENV
131: * ..
132: * .. Intrinsic Functions ..
133: INTRINSIC INT, MAX, MIN
134: * ..
135: * .. Executable Statements ..
136: *
137: * Test the input parameters
138: *
139: INFO = 0
140: NP = MIN( N, P )
141: LQUERY = ( LWORK.EQ.-1 )
142: IF( N.LT.0 ) THEN
143: INFO = -1
144: ELSE IF( M.LT.0 .OR. M.GT.N ) THEN
145: INFO = -2
146: ELSE IF( P.LT.0 .OR. P.LT.N-M ) THEN
147: INFO = -3
148: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
149: INFO = -5
150: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
151: INFO = -7
152: END IF
153: *
154: * Calculate workspace
155: *
156: IF( INFO.EQ.0) THEN
157: IF( N.EQ.0 ) THEN
158: LWKMIN = 1
159: LWKOPT = 1
160: ELSE
161: NB1 = ILAENV( 1, 'ZGEQRF', ' ', N, M, -1, -1 )
162: NB2 = ILAENV( 1, 'ZGERQF', ' ', N, M, -1, -1 )
163: NB3 = ILAENV( 1, 'ZUNMQR', ' ', N, M, P, -1 )
164: NB4 = ILAENV( 1, 'ZUNMRQ', ' ', N, M, P, -1 )
165: NB = MAX( NB1, NB2, NB3, NB4 )
166: LWKMIN = M + N + P
167: LWKOPT = M + NP + MAX( N, P )*NB
168: END IF
169: WORK( 1 ) = LWKOPT
170: *
171: IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
172: INFO = -12
173: END IF
174: END IF
175: *
176: IF( INFO.NE.0 ) THEN
177: CALL XERBLA( 'ZGGGLM', -INFO )
178: RETURN
179: ELSE IF( LQUERY ) THEN
180: RETURN
181: END IF
182: *
183: * Quick return if possible
184: *
185: IF( N.EQ.0 )
186: $ RETURN
187: *
188: * Compute the GQR factorization of matrices A and B:
189: *
190: * Q'*A = ( R11 ) M, Q'*B*Z' = ( T11 T12 ) M
191: * ( 0 ) N-M ( 0 T22 ) N-M
192: * M M+P-N N-M
193: *
194: * where R11 and T22 are upper triangular, and Q and Z are
195: * unitary.
196: *
197: CALL ZGGQRF( N, M, P, A, LDA, WORK, B, LDB, WORK( M+1 ),
198: $ WORK( M+NP+1 ), LWORK-M-NP, INFO )
199: LOPT = WORK( M+NP+1 )
200: *
201: * Update left-hand-side vector d = Q'*d = ( d1 ) M
202: * ( d2 ) N-M
203: *
204: CALL ZUNMQR( 'Left', 'Conjugate transpose', N, 1, M, A, LDA, WORK,
205: $ D, MAX( 1, N ), WORK( M+NP+1 ), LWORK-M-NP, INFO )
206: LOPT = MAX( LOPT, INT( WORK( M+NP+1 ) ) )
207: *
208: * Solve T22*y2 = d2 for y2
209: *
210: IF( N.GT.M ) THEN
211: CALL ZTRTRS( 'Upper', 'No transpose', 'Non unit', N-M, 1,
212: $ B( M+1, M+P-N+1 ), LDB, D( M+1 ), N-M, INFO )
213: *
214: IF( INFO.GT.0 ) THEN
215: INFO = 1
216: RETURN
217: END IF
218: *
219: CALL ZCOPY( N-M, D( M+1 ), 1, Y( M+P-N+1 ), 1 )
220: END IF
221: *
222: * Set y1 = 0
223: *
224: DO 10 I = 1, M + P - N
225: Y( I ) = CZERO
226: 10 CONTINUE
227: *
228: * Update d1 = d1 - T12*y2
229: *
230: CALL ZGEMV( 'No transpose', M, N-M, -CONE, B( 1, M+P-N+1 ), LDB,
231: $ Y( M+P-N+1 ), 1, CONE, D, 1 )
232: *
233: * Solve triangular system: R11*x = d1
234: *
235: IF( M.GT.0 ) THEN
236: CALL ZTRTRS( 'Upper', 'No Transpose', 'Non unit', M, 1, A, LDA,
237: $ D, M, INFO )
238: *
239: IF( INFO.GT.0 ) THEN
240: INFO = 2
241: RETURN
242: END IF
243: *
244: * Copy D to X
245: *
246: CALL ZCOPY( M, D, 1, X, 1 )
247: END IF
248: *
249: * Backward transformation y = Z'*y
250: *
251: CALL ZUNMRQ( 'Left', 'Conjugate transpose', P, 1, NP,
252: $ B( MAX( 1, N-P+1 ), 1 ), LDB, WORK( M+1 ), Y,
253: $ MAX( 1, P ), WORK( M+NP+1 ), LWORK-M-NP, INFO )
254: WORK( 1 ) = M + NP + MAX( LOPT, INT( WORK( M+NP+1 ) ) )
255: *
256: RETURN
257: *
258: * End of ZGGGLM
259: *
260: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zggglm.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack