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