149
150
151
152
153
154
155 INTEGER LDA, LDB, LWORK, M, N, NRHS, RANK, RKSEL, SCALE
156 DOUBLE PRECISION NORMA, NORMB
157
158
159 INTEGER ISEED( 4 )
160 DOUBLE PRECISION S( * )
161 COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( LWORK )
162
163
164
165
166
167 DOUBLE PRECISION ZERO, ONE, TWO, SVMIN
168 parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0,
169 $ svmin = 0.1d+0 )
170 COMPLEX*16 CZERO, CONE
171 parameter( czero = ( 0.0d+0, 0.0d+0 ),
172 $ cone = ( 1.0d+0, 0.0d+0 ) )
173
174
175 INTEGER INFO, J, MN
176 DOUBLE PRECISION BIGNUM, EPS, SMLNUM, TEMP
177
178
179 DOUBLE PRECISION DUMMY( 1 )
180
181
182 DOUBLE PRECISION DASUM, DLAMCH, DLARND, DZNRM2, ZLANGE
184
185
188
189
190 INTRINSIC abs, dcmplx, max, min
191
192
193
194 mn = min( m, n )
195 IF( lwork.LT.max( m+mn, mn*nrhs, 2*n+m ) ) THEN
196 CALL xerbla(
'ZQRT15', 16 )
197 RETURN
198 END IF
199
200 smlnum =
dlamch(
'Safe minimum' )
201 bignum = one / smlnum
203 smlnum = ( smlnum / eps ) / eps
204 bignum = one / smlnum
205
206
207
208 IF( rksel.EQ.1 ) THEN
209 rank = mn
210 ELSE IF( rksel.EQ.2 ) THEN
211 rank = ( 3*mn ) / 4
212 DO 10 j = rank + 1, mn
213 s( j ) = zero
214 10 CONTINUE
215 ELSE
216 CALL xerbla(
'ZQRT15', 2 )
217 END IF
218
219 IF( rank.GT.0 ) THEN
220
221
222
223 s( 1 ) = one
224 DO 30 j = 2, rank
225 20 CONTINUE
227 IF( temp.GT.svmin ) THEN
228 s( j ) = abs( temp )
229 ELSE
230 GO TO 20
231 END IF
232 30 CONTINUE
233 CALL dlaord(
'Decreasing', rank, s, 1 )
234
235
236
237 CALL zlarnv( 2, iseed, m, work )
239 CALL zlaset(
'Full', m, rank, czero, cone, a, lda )
240 CALL zlarf(
'Left', m, rank, work, 1, dcmplx( two ), a, lda,
241 $ work( m+1 ) )
242
243
244
245
246
247 CALL zlarnv( 2, iseed, rank*nrhs, work )
248 CALL zgemm(
'No transpose',
'No transpose', m, nrhs, rank,
249 $ cone, a, lda, work, rank, czero, b, ldb )
250
251
252
253
254
255 DO 40 j = 1, rank
256 CALL zdscal( m, s( j ), a( 1, j ), 1 )
257 40 CONTINUE
258 IF( rank.LT.n )
259 $
CALL zlaset(
'Full', m, n-rank, czero, czero,
260 $ a( 1, rank+1 ), lda )
261 CALL zlaror(
'Right',
'No initialization', m, n, a, lda, iseed,
262 $ work, info )
263
264 ELSE
265
266
267
268
269
270 DO 50 j = 1, mn
271 s( j ) = zero
272 50 CONTINUE
273 CALL zlaset(
'Full', m, n, czero, czero, a, lda )
274 CALL zlaset(
'Full', m, nrhs, czero, czero, b, ldb )
275
276 END IF
277
278
279
280 IF( scale.NE.1 ) THEN
281 norma =
zlange(
'Max', m, n, a, lda, dummy )
282 IF( norma.NE.zero ) THEN
283 IF( scale.EQ.2 ) THEN
284
285
286
287 CALL zlascl(
'General', 0, 0, norma, bignum, m, n, a,
288 $ lda, info )
289 CALL dlascl(
'General', 0, 0, norma, bignum, mn, 1, s,
290 $ mn, info )
291 CALL zlascl(
'General', 0, 0, norma, bignum, m, nrhs, b,
292 $ ldb, info )
293 ELSE IF( scale.EQ.3 ) THEN
294
295
296
297 CALL zlascl(
'General', 0, 0, norma, smlnum, m, n, a,
298 $ lda, info )
299 CALL dlascl(
'General', 0, 0, norma, smlnum, mn, 1, s,
300 $ mn, info )
301 CALL zlascl(
'General', 0, 0, norma, smlnum, m, nrhs, b,
302 $ ldb, info )
303 ELSE
304 CALL xerbla(
'ZQRT15', 1 )
305 RETURN
306 END IF
307 END IF
308 END IF
309
310 norma =
dasum( mn, s, 1 )
311 normb =
zlange(
'One-norm', m, nrhs, b, ldb, dummy )
312
313 RETURN
314
315
316
subroutine xerbla(srname, info)
subroutine dlaord(job, n, x, incx)
DLAORD
double precision function dlarnd(idist, iseed)
DLARND
double precision function dasum(n, dx, incx)
DASUM
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
double precision function dlamch(cmach)
DLAMCH
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 zlarf(side, m, n, v, incv, tau, c, ldc, work)
ZLARF applies an elementary reflector to a general rectangular matrix.
subroutine zlarnv(idist, iseed, n, x)
ZLARNV returns a vector of random numbers from a uniform or normal distribution.
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 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 zlaset(uplo, m, n, alpha, beta, a, lda)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
real(wp) function dznrm2(n, x, incx)
DZNRM2
subroutine zdscal(n, da, zx, incx)
ZDSCAL
subroutine zlaror(side, init, m, n, a, lda, iseed, x, info)
ZLAROR