Annotation of rpl/lapack/lapack/zggrqf.f, revision 1.2
1.1 bertrand 1: SUBROUTINE ZGGRQF( M, P, N, A, LDA, TAUA, B, LDB, TAUB, WORK,
2: $ LWORK, INFO )
3: *
4: * -- LAPACK 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, * ), TAUA( * ), TAUB( * ),
14: $ WORK( * )
15: * ..
16: *
17: * Purpose
18: * =======
19: *
20: * ZGGRQF computes a generalized RQ factorization of an M-by-N matrix A
21: * and a P-by-N matrix B:
22: *
23: * A = R*Q, B = Z*T*Q,
24: *
25: * where Q is an N-by-N unitary matrix, Z is a P-by-P unitary
26: * matrix, and R and T assume one of the forms:
27: *
28: * if M <= N, R = ( 0 R12 ) M, or if M > N, R = ( R11 ) M-N,
29: * N-M M ( R21 ) N
30: * N
31: *
32: * where R12 or R21 is upper triangular, and
33: *
34: * if P >= N, T = ( T11 ) N , or if P < N, T = ( T11 T12 ) P,
35: * ( 0 ) P-N P N-P
36: * N
37: *
38: * where T11 is upper triangular.
39: *
40: * In particular, if B is square and nonsingular, the GRQ factorization
41: * of A and B implicitly gives the RQ factorization of A*inv(B):
42: *
43: * A*inv(B) = (R*inv(T))*Z'
44: *
45: * where inv(B) denotes the inverse of the matrix B, and Z' denotes the
46: * conjugate transpose of the matrix Z.
47: *
48: * Arguments
49: * =========
50: *
51: * M (input) INTEGER
52: * The number of rows of the matrix A. M >= 0.
53: *
54: * P (input) INTEGER
55: * The number of rows of the matrix B. P >= 0.
56: *
57: * N (input) INTEGER
58: * The number of columns of the matrices A and B. N >= 0.
59: *
60: * A (input/output) COMPLEX*16 array, dimension (LDA,N)
61: * On entry, the M-by-N matrix A.
62: * On exit, if M <= N, the upper triangle of the subarray
63: * A(1:M,N-M+1:N) contains the M-by-M upper triangular matrix R;
64: * if M > N, the elements on and above the (M-N)-th subdiagonal
65: * contain the M-by-N upper trapezoidal matrix R; the remaining
66: * elements, with the array TAUA, represent the unitary
67: * matrix Q as a product of elementary reflectors (see Further
68: * Details).
69: *
70: * LDA (input) INTEGER
71: * The leading dimension of the array A. LDA >= max(1,M).
72: *
73: * TAUA (output) COMPLEX*16 array, dimension (min(M,N))
74: * The scalar factors of the elementary reflectors which
75: * represent the unitary matrix Q (see Further Details).
76: *
77: * B (input/output) COMPLEX*16 array, dimension (LDB,N)
78: * On entry, the P-by-N matrix B.
79: * On exit, the elements on and above the diagonal of the array
80: * contain the min(P,N)-by-N upper trapezoidal matrix T (T is
81: * upper triangular if P >= N); the elements below the diagonal,
82: * with the array TAUB, represent the unitary matrix Z as a
83: * product of elementary reflectors (see Further Details).
84: *
85: * LDB (input) INTEGER
86: * The leading dimension of the array B. LDB >= max(1,P).
87: *
88: * TAUB (output) COMPLEX*16 array, dimension (min(P,N))
89: * The scalar factors of the elementary reflectors which
90: * represent the unitary matrix Z (see Further Details).
91: *
92: * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
93: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
94: *
95: * LWORK (input) INTEGER
96: * The dimension of the array WORK. LWORK >= max(1,N,M,P).
97: * For optimum performance LWORK >= max(N,M,P)*max(NB1,NB2,NB3),
98: * where NB1 is the optimal blocksize for the RQ factorization
99: * of an M-by-N matrix, NB2 is the optimal blocksize for the
100: * QR factorization of a P-by-N matrix, and NB3 is the optimal
101: * blocksize for a call of ZUNMRQ.
102: *
103: * If LWORK = -1, then a workspace query is assumed; the routine
104: * only calculates the optimal size of the WORK array, returns
105: * this value as the first entry of the WORK array, and no error
106: * message related to LWORK is issued by XERBLA.
107: *
108: * INFO (output) INTEGER
109: * = 0: successful exit
110: * < 0: if INFO=-i, the i-th argument had an illegal value.
111: *
112: * Further Details
113: * ===============
114: *
115: * The matrix Q is represented as a product of elementary reflectors
116: *
117: * Q = H(1) H(2) . . . H(k), where k = min(m,n).
118: *
119: * Each H(i) has the form
120: *
121: * H(i) = I - taua * v * v'
122: *
123: * where taua is a complex scalar, and v is a complex vector with
124: * v(n-k+i+1:n) = 0 and v(n-k+i) = 1; v(1:n-k+i-1) is stored on exit in
125: * A(m-k+i,1:n-k+i-1), and taua in TAUA(i).
126: * To form Q explicitly, use LAPACK subroutine ZUNGRQ.
127: * To use Q to update another matrix, use LAPACK subroutine ZUNMRQ.
128: *
129: * The matrix Z is represented as a product of elementary reflectors
130: *
131: * Z = H(1) H(2) . . . H(k), where k = min(p,n).
132: *
133: * Each H(i) has the form
134: *
135: * H(i) = I - taub * v * v'
136: *
137: * where taub is a complex scalar, and v is a complex vector with
138: * v(1:i-1) = 0 and v(i) = 1; v(i+1:p) is stored on exit in B(i+1:p,i),
139: * and taub in TAUB(i).
140: * To form Z explicitly, use LAPACK subroutine ZUNGQR.
141: * To use Z to update another matrix, use LAPACK subroutine ZUNMQR.
142: *
143: * =====================================================================
144: *
145: * .. Local Scalars ..
146: LOGICAL LQUERY
147: INTEGER LOPT, LWKOPT, NB, NB1, NB2, NB3
148: * ..
149: * .. External Subroutines ..
150: EXTERNAL XERBLA, ZGEQRF, ZGERQF, ZUNMRQ
151: * ..
152: * .. External Functions ..
153: INTEGER ILAENV
154: EXTERNAL ILAENV
155: * ..
156: * .. Intrinsic Functions ..
157: INTRINSIC INT, MAX, MIN
158: * ..
159: * .. Executable Statements ..
160: *
161: * Test the input parameters
162: *
163: INFO = 0
164: NB1 = ILAENV( 1, 'ZGERQF', ' ', M, N, -1, -1 )
165: NB2 = ILAENV( 1, 'ZGEQRF', ' ', P, N, -1, -1 )
166: NB3 = ILAENV( 1, 'ZUNMRQ', ' ', M, N, P, -1 )
167: NB = MAX( NB1, NB2, NB3 )
168: LWKOPT = MAX( N, M, P )*NB
169: WORK( 1 ) = LWKOPT
170: LQUERY = ( LWORK.EQ.-1 )
171: IF( M.LT.0 ) THEN
172: INFO = -1
173: ELSE IF( P.LT.0 ) THEN
174: INFO = -2
175: ELSE IF( N.LT.0 ) THEN
176: INFO = -3
177: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
178: INFO = -5
179: ELSE IF( LDB.LT.MAX( 1, P ) ) THEN
180: INFO = -8
181: ELSE IF( LWORK.LT.MAX( 1, M, P, N ) .AND. .NOT.LQUERY ) THEN
182: INFO = -11
183: END IF
184: IF( INFO.NE.0 ) THEN
185: CALL XERBLA( 'ZGGRQF', -INFO )
186: RETURN
187: ELSE IF( LQUERY ) THEN
188: RETURN
189: END IF
190: *
191: * RQ factorization of M-by-N matrix A: A = R*Q
192: *
193: CALL ZGERQF( M, N, A, LDA, TAUA, WORK, LWORK, INFO )
194: LOPT = WORK( 1 )
195: *
196: * Update B := B*Q'
197: *
198: CALL ZUNMRQ( 'Right', 'Conjugate Transpose', P, N, MIN( M, N ),
199: $ A( MAX( 1, M-N+1 ), 1 ), LDA, TAUA, B, LDB, WORK,
200: $ LWORK, INFO )
201: LOPT = MAX( LOPT, INT( WORK( 1 ) ) )
202: *
203: * QR factorization of P-by-N matrix B: B = Z*T
204: *
205: CALL ZGEQRF( P, N, B, LDB, TAUB, WORK, LWORK, INFO )
206: WORK( 1 ) = MAX( LOPT, INT( WORK( 1 ) ) )
207: *
208: RETURN
209: *
210: * End of ZGGRQF
211: *
212: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zggrqf.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack