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