148
149
150
151
152
153
154 INTEGER LDA, LDB, LWORK, M, N, NRHS, RANK, RKSEL, SCALE
155 REAL NORMA, NORMB
156
157
158 INTEGER ISEED( 4 )
159 REAL A( LDA, * ), B( LDB, * ), S( * ), WORK( LWORK )
160
161
162
163
164
165 REAL ZERO, ONE, TWO, SVMIN
166 parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0,
167 $ svmin = 0.1e0 )
168
169
170 INTEGER INFO, J, MN
171 REAL BIGNUM, EPS, SMLNUM, TEMP
172
173
174 REAL DUMMY( 1 )
175
176
177 REAL SASUM, SLAMCH, SLANGE, SLARND, SNRM2
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(
'SQRT15', 16 )
192 RETURN
193 END IF
194
195 smlnum =
slamch(
'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(
'SQRT15', 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 slaord(
'Decreasing', rank, s, 1 )
229
230
231
232 CALL slarnv( 2, iseed, m, work )
233 CALL sscal( m, one /
snrm2( m, work, 1 ), work, 1 )
234 CALL slaset(
'Full', m, rank, zero, one, a, lda )
235 CALL slarf(
'Left', m, rank, work, 1, two, a, lda,
236 $ work( m+1 ) )
237
238
239
240
241
242 CALL slarnv( 2, iseed, rank*nrhs, work )
243 CALL sgemm(
'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 sscal( m, s( j ), a( 1, j ), 1 )
252 40 CONTINUE
253 IF( rank.LT.n )
254 $
CALL slaset(
'Full', m, n-rank, zero, zero, a( 1, rank+1 ),
255 $ lda )
256 CALL slaror(
'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 slaset(
'Full', m, n, zero, zero, a, lda )
269 CALL slaset(
'Full', m, nrhs, zero, zero, b, ldb )
270
271 END IF
272
273
274
275 IF( scale.NE.1 ) THEN
276 norma =
slange(
'Max', m, n, a, lda, dummy )
277 IF( norma.NE.zero ) THEN
278 IF( scale.EQ.2 ) THEN
279
280
281
282 CALL slascl(
'General', 0, 0, norma, bignum, m, n, a,
283 $ lda, info )
284 CALL slascl(
'General', 0, 0, norma, bignum, mn, 1, s,
285 $ mn, info )
286 CALL slascl(
'General', 0, 0, norma, bignum, m, nrhs, b,
287 $ ldb, info )
288 ELSE IF( scale.EQ.3 ) THEN
289
290
291
292 CALL slascl(
'General', 0, 0, norma, smlnum, m, n, a,
293 $ lda, info )
294 CALL slascl(
'General', 0, 0, norma, smlnum, mn, 1, s,
295 $ mn, info )
296 CALL slascl(
'General', 0, 0, norma, smlnum, m, nrhs, b,
297 $ ldb, info )
298 ELSE
299 CALL xerbla(
'SQRT15', 1 )
300 RETURN
301 END IF
302 END IF
303 END IF
304
305 norma =
sasum( mn, s, 1 )
306 normb =
slange(
'One-norm', m, nrhs, b, ldb, dummy )
307
308 RETURN
309
310
311
subroutine xerbla(srname, info)
real function sasum(n, sx, incx)
SASUM
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
real function slamch(cmach)
SLAMCH
real function slange(norm, m, n, a, lda, work)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine slarf(side, m, n, v, incv, tau, c, ldc, work)
SLARF applies an elementary reflector to a general rectangular matrix.
subroutine slarnv(idist, iseed, n, x)
SLARNV returns a vector of random numbers from a uniform or normal distribution.
subroutine slascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine slaset(uplo, m, n, alpha, beta, a, lda)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
real(wp) function snrm2(n, x, incx)
SNRM2
subroutine sscal(n, sa, sx, incx)
SSCAL
subroutine slaord(job, n, x, incx)
SLAORD
real function slarnd(idist, iseed)
SLARND
subroutine slaror(side, init, m, n, a, lda, iseed, x, info)
SLAROR