116
117
118
119
120
121
122 CHARACTER TRANS
123 INTEGER LDA, LDX, LWORK, M, N, NRHS
124
125
126 DOUBLE PRECISION 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, DLANGE
147
148
150
151
152 INTRINSIC abs, dble, 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(
'DQRT14', 10 )
162 RETURN
163 ELSE IF( n.LE.0 .OR. nrhs.LE.0 ) THEN
164 RETURN
165 END IF
166 ELSE IF(
lsame( trans,
'T' ) )
THEN
167 ldwork = m
168 tpsd = .true.
169 IF( lwork.LT.( n+nrhs )*( m+2 ) ) THEN
170 CALL xerbla(
'DQRT14', 10 )
171 RETURN
172 ELSE IF( m.LE.0 .OR. nrhs.LE.0 ) THEN
173 RETURN
174 END IF
175 ELSE
176 CALL xerbla(
'DQRT14', 1 )
177 RETURN
178 END IF
179
180
181
182 CALL dlacpy(
'All', m, n, a, lda, work, ldwork )
183 anrm =
dlange(
'M', m, n, work, ldwork, rwork )
184 IF( anrm.NE.zero )
185 $
CALL dlascl(
'G', 0, 0, anrm, one, m, n, work, ldwork, info )
186
187
188
189 IF( tpsd ) THEN
190
191
192
193 CALL dlacpy(
'All', m, nrhs, x, ldx, work( n*ldwork+1 ),
194 $ ldwork )
195 xnrm =
dlange(
'M', m, nrhs, work( n*ldwork+1 ), ldwork,
196 $ rwork )
197 IF( xnrm.NE.zero )
198 $
CALL dlascl(
'G', 0, 0, xnrm, one, m, nrhs,
199 $ work( n*ldwork+1 ), ldwork, info )
200
201
202
203 CALL dgeqr2( 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 ) = x( i, j )
225 30 CONTINUE
226 40 CONTINUE
227
228 xnrm =
dlange(
'M', nrhs, n, work( m+1 ), ldwork, rwork )
229 IF( xnrm.NE.zero )
230 $
CALL dlascl(
'G', 0, 0, xnrm, one, nrhs, n, work( m+1 ),
231 $ ldwork, info )
232
233
234
235 CALL dgelq2( 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 dqrt14 = err / ( dble( max( m, n, nrhs ) )*
dlamch(
'Epsilon' ) )
251
252 RETURN
253
254
255
double precision function dlamch(CMACH)
DLAMCH
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
subroutine xerbla(SRNAME, INFO)
XERBLA
logical function lsame(CA, CB)
LSAME
double precision function dqrt14(TRANS, M, N, NRHS, A, LDA, X, LDX, WORK, LWORK)
DQRT14
double precision function dlange(NORM, M, N, A, LDA, WORK)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine dgelq2(M, N, A, LDA, TAU, WORK, INFO)
DGELQ2 computes the LQ factorization of a general rectangular matrix using an unblocked algorithm.
subroutine dgeqr2(M, N, A, LDA, TAU, WORK, INFO)
DGEQR2 computes the QR factorization of a general rectangular matrix using an unblocked algorithm.