116
117
118
119
120
121
122 CHARACTER TRANS
123 INTEGER LDA, LDX, LWORK, M, N, NRHS
124
125
126 COMPLEX*16 A( LDA, * ), WORK( LWORK ), X( LDX, * )
127
128
129
130
131
132 DOUBLE PRECISION ZERO, ONE
133 parameter( zero = 0.0d0, one = 1.0d0 )
134
135
136 LOGICAL TPSD
137 INTEGER I, INFO, J, LDWORK
138 DOUBLE PRECISION ANRM, ERR, XNRM
139
140
141 DOUBLE PRECISION RWORK( 1 )
142
143
144 LOGICAL LSAME
145 DOUBLE PRECISION DLAMCH, ZLANGE
147
148
150
151
152 INTRINSIC abs, dble, dconjg, max, min
153
154
155
157 IF(
lsame( trans,
'N' ) )
THEN
158 ldwork = m + nrhs
159 tpsd = .false.
160 IF( lwork.LT.( m+nrhs )*( n+2 ) ) THEN
161 CALL xerbla(
'ZQRT14', 10 )
162 RETURN
163 ELSE IF( n.LE.0 .OR. nrhs.LE.0 ) THEN
164 RETURN
165 END IF
166 ELSE IF(
lsame( trans,
'C' ) )
THEN
167 ldwork = m
168 tpsd = .true.
169 IF( lwork.LT.( n+nrhs )*( m+2 ) ) THEN
170 CALL xerbla(
'ZQRT14', 10 )
171 RETURN
172 ELSE IF( m.LE.0 .OR. nrhs.LE.0 ) THEN
173 RETURN
174 END IF
175 ELSE
176 CALL xerbla(
'ZQRT14', 1 )
177 RETURN
178 END IF
179
180
181
182 CALL zlacpy(
'All', m, n, a, lda, work, ldwork )
183 anrm =
zlange(
'M', m, n, work, ldwork, rwork )
184 IF( anrm.NE.zero )
185 $
CALL zlascl(
'G', 0, 0, anrm, one, m, n, work, ldwork, info )
186
187
188
189 IF( tpsd ) THEN
190
191
192
193 CALL zlacpy(
'All', m, nrhs, x, ldx, work( n*ldwork+1 ),
194 $ ldwork )
195 xnrm =
zlange(
'M', m, nrhs, work( n*ldwork+1 ), ldwork,
196 $ rwork )
197 IF( xnrm.NE.zero )
198 $
CALL zlascl(
'G', 0, 0, xnrm, one, m, nrhs,
199 $ work( n*ldwork+1 ), ldwork, info )
200
201
202
203 CALL zgeqr2( m, n+nrhs, work, ldwork,
204 $ work( ldwork*( n+nrhs )+1 ),
205 $ work( ldwork*( n+nrhs )+min( m, n+nrhs )+1 ),
206 $ info )
207
208
209
210
211 err = zero
212 DO 20 j = n + 1, n + nrhs
213 DO 10 i = n + 1, min( m, j )
214 err = max( err, abs( work( i+( j-1 )*m ) ) )
215 10 CONTINUE
216 20 CONTINUE
217
218 ELSE
219
220
221
222 DO 40 i = 1, n
223 DO 30 j = 1, nrhs
224 work( m+j+( i-1 )*ldwork ) = dconjg( x( i, j ) )
225 30 CONTINUE
226 40 CONTINUE
227
228 xnrm =
zlange(
'M', nrhs, n, work( m+1 ), ldwork, rwork )
229 IF( xnrm.NE.zero )
230 $
CALL zlascl(
'G', 0, 0, xnrm, one, nrhs, n, work( m+1 ),
231 $ ldwork, info )
232
233
234
235 CALL zgelq2( ldwork, n, work, ldwork, work( ldwork*n+1 ),
236 $ work( ldwork*( n+1 )+1 ), info )
237
238
239
240
241 err = zero
242 DO 60 j = m + 1, n
243 DO 50 i = j, ldwork
244 err = max( err, abs( work( i+( j-1 )*ldwork ) ) )
245 50 CONTINUE
246 60 CONTINUE
247
248 END IF
249
250 zqrt14 = err / ( dble( max( m, n, nrhs ) )*
dlamch(
'Epsilon' ) )
251
252 RETURN
253
254
255
double precision function dlamch(CMACH)
DLAMCH
subroutine xerbla(SRNAME, INFO)
XERBLA
logical function lsame(CA, CB)
LSAME
double precision function zqrt14(TRANS, M, N, NRHS, A, LDA, X, LDX, WORK, LWORK)
ZQRT14
double precision function zlange(NORM, M, N, A, LDA, WORK)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine zgelq2(M, N, A, LDA, TAU, WORK, INFO)
ZGELQ2 computes the LQ factorization of a general rectangular matrix using an unblocked algorithm.
subroutine zgeqr2(M, N, A, LDA, TAU, WORK, INFO)
ZGEQR2 computes the QR factorization of a general rectangular matrix using an unblocked algorithm.
subroutine zlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.