1       SUBROUTINE CQRT15( SCALE, RKSEL, M, N, NRHS, A, LDA, B, LDB, S,
  2      $                   RANK, NORMA, NORMB, ISEED, WORK, LWORK )
  3 *
  4 *  -- LAPACK test routine (version 3.1) --
  5 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
  6 *     November 2006
  7 *
  8 *     .. Scalar Arguments ..
  9       INTEGER            LDA, LDB, LWORK, M, N, NRHS, RANK, RKSEL, SCALE
 10       REAL               NORMA, NORMB
 11 *     ..
 12 *     .. Array Arguments ..
 13       INTEGER            ISEED( 4 )
 14       REAL               S( * )
 15       COMPLEX            A( LDA, * ), B( LDB, * ), WORK( LWORK )
 16 *     ..
 17 *
 18 *  Purpose
 19 *  =======
 20 *
 21 *  CQRT15 generates a matrix with full or deficient rank and of various
 22 *  norms.
 23 *
 24 *  Arguments
 25 *  =========
 26 *
 27 *  SCALE   (input) INTEGER
 28 *          SCALE = 1: normally scaled matrix
 29 *          SCALE = 2: matrix scaled up
 30 *          SCALE = 3: matrix scaled down
 31 *
 32 *  RKSEL   (input) INTEGER
 33 *          RKSEL = 1: full rank matrix
 34 *          RKSEL = 2: rank-deficient matrix
 35 *
 36 *  M       (input) INTEGER
 37 *          The number of rows of the matrix A.
 38 *
 39 *  N       (input) INTEGER
 40 *          The number of columns of A.
 41 *
 42 *  NRHS    (input) INTEGER
 43 *          The number of columns of B.
 44 *
 45 *  A       (output) COMPLEX array, dimension (LDA,N)
 46 *          The M-by-N matrix A.
 47 *
 48 *  LDA     (input) INTEGER
 49 *          The leading dimension of the array A.
 50 *
 51 *  B       (output) COMPLEX array, dimension (LDB, NRHS)
 52 *          A matrix that is in the range space of matrix A.
 53 *
 54 *  LDB     (input) INTEGER
 55 *          The leading dimension of the array B.
 56 *
 57 *  S       (output) REAL array, dimension MIN(M,N)
 58 *          Singular values of A.
 59 *
 60 *  RANK    (output) INTEGER
 61 *          number of nonzero singular values of A.
 62 *
 63 *  NORMA   (output) REAL
 64 *          one-norm norm of A.
 65 *
 66 *  NORMB   (output) REAL
 67 *          one-norm norm of B.
 68 *
 69 *  ISEED   (input/output) integer array, dimension (4)
 70 *          seed for random number generator.
 71 *
 72 *  WORK    (workspace) COMPLEX array, dimension (LWORK)
 73 *
 74 *  LWORK   (input) INTEGER
 75 *          length of work space required.
 76 *          LWORK >= MAX(M+MIN(M,N),NRHS*MIN(M,N),2*N+M)
 77 *
 78 *  =====================================================================
 79 *
 80 *     .. Parameters ..
 81       REAL               ZERO, ONE, TWO, SVMIN
 82       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0, TWO = 2.0E+0,
 83      $                   SVMIN = 0.1E+0 )
 84       COMPLEX            CZERO, CONE
 85       PARAMETER          ( CZERO = ( 0.0E+00.0E+0 ),
 86      $                   CONE = ( 1.0E+00.0E+0 ) )
 87 *     ..
 88 *     .. Local Scalars ..
 89       INTEGER            INFO, J, MN
 90       REAL               BIGNUM, EPS, SMLNUM, TEMP
 91 *     ..
 92 *     .. Local Arrays ..
 93       REAL               DUMMY( 1 )
 94 *     ..
 95 *     .. External Functions ..
 96       REAL               CLANGE, SASUM, SCNRM2, SLAMCH, SLARND
 97       EXTERNAL           CLANGE, SASUM, SCNRM2, SLAMCH, SLARND
 98 *     ..
 99 *     .. External Subroutines ..
100       EXTERNAL           CGEMM, CLARF, CLARNV, CLAROR, CLASCL, CLASET,
101      $                   CSSCAL, SLABAD, SLAORD, SLASCL, XERBLA
102 *     ..
103 *     .. Intrinsic Functions ..
104       INTRINSIC          ABSCMPLXMAXMIN
105 *     ..
106 *     .. Executable Statements ..
107 *
108       MN = MIN( M, N )
109       IF( LWORK.LT.MAX( M+MN, MN*NRHS, 2*N+M ) ) THEN
110          CALL XERBLA( 'CQRT15'16 )
111          RETURN
112       END IF
113 *
114       SMLNUM = SLAMCH( 'Safe minimum' )
115       BIGNUM = ONE / SMLNUM
116       CALL SLABAD( SMLNUM, BIGNUM )
117       EPS = SLAMCH( 'Epsilon' )
118       SMLNUM = ( SMLNUM / EPS ) / EPS
119       BIGNUM = ONE / SMLNUM
120 *
121 *     Determine rank and (unscaled) singular values
122 *
123       IF( RKSEL.EQ.1 ) THEN
124          RANK = MN
125       ELSE IF( RKSEL.EQ.2 ) THEN
126          RANK = ( 3*MN ) / 4
127          DO 10 J = RANK + 1, MN
128             S( J ) = ZERO
129    10    CONTINUE
130       ELSE
131          CALL XERBLA( 'CQRT15'2 )
132       END IF
133 *
134       IF( RANK.GT.0 ) THEN
135 *
136 *        Nontrivial case
137 *
138          S( 1 ) = ONE
139          DO 30 J = 2, RANK
140    20       CONTINUE
141             TEMP = SLARND( 1, ISEED )
142             IF( TEMP.GT.SVMIN ) THEN
143                S( J ) = ABS( TEMP )
144             ELSE
145                GO TO 20
146             END IF
147    30    CONTINUE
148          CALL SLAORD( 'Decreasing', RANK, S, 1 )
149 *
150 *        Generate 'rank' columns of a random orthogonal matrix in A
151 *
152          CALL CLARNV( 2, ISEED, M, WORK )
153          CALL CSSCAL( M, ONE / SCNRM2( M, WORK, 1 ), WORK, 1 )
154          CALL CLASET( 'Full', M, RANK, CZERO, CONE, A, LDA )
155          CALL CLARF( 'Left', M, RANK, WORK, 1CMPLX( TWO ), A, LDA,
156      $               WORK( M+1 ) )
157 *
158 *        workspace used: m+mn
159 *
160 *        Generate consistent rhs in the range space of A
161 *
162          CALL CLARNV( 2, ISEED, RANK*NRHS, WORK )
163          CALL CGEMM( 'No transpose''No transpose', M, NRHS, RANK,
164      $               CONE, A, LDA, WORK, RANK, CZERO, B, LDB )
165 *
166 *        work space used: <= mn *nrhs
167 *
168 *        generate (unscaled) matrix A
169 *
170          DO 40 J = 1, RANK
171             CALL CSSCAL( M, S( J ), A( 1, J ), 1 )
172    40    CONTINUE
173          IF( RANK.LT.N )
174      $      CALL CLASET( 'Full', M, N-RANK, CZERO, CZERO,
175      $                   A( 1, RANK+1 ), LDA )
176          CALL CLAROR( 'Right''No initialization', M, N, A, LDA, ISEED,
177      $                WORK, INFO )
178 *
179       ELSE
180 *
181 *        work space used 2*n+m
182 *
183 *        Generate null matrix and rhs
184 *
185          DO 50 J = 1, MN
186             S( J ) = ZERO
187    50    CONTINUE
188          CALL CLASET( 'Full', M, N, CZERO, CZERO, A, LDA )
189          CALL CLASET( 'Full', M, NRHS, CZERO, CZERO, B, LDB )
190 *
191       END IF
192 *
193 *     Scale the matrix
194 *
195       IFSCALE.NE.1 ) THEN
196          NORMA = CLANGE( 'Max', M, N, A, LDA, DUMMY )
197          IF( NORMA.NE.ZERO ) THEN
198             IFSCALE.EQ.2 ) THEN
199 *
200 *              matrix scaled up
201 *
202                CALL CLASCL( 'General'00, NORMA, BIGNUM, M, N, A,
203      $                      LDA, INFO )
204                CALL SLASCL( 'General'00, NORMA, BIGNUM, MN, 1, S,
205      $                      MN, INFO )
206                CALL CLASCL( 'General'00, NORMA, BIGNUM, M, NRHS, B,
207      $                      LDB, INFO )
208             ELSE IFSCALE.EQ.3 ) THEN
209 *
210 *              matrix scaled down
211 *
212                CALL CLASCL( 'General'00, NORMA, SMLNUM, M, N, A,
213      $                      LDA, INFO )
214                CALL SLASCL( 'General'00, NORMA, SMLNUM, MN, 1, S,
215      $                      MN, INFO )
216                CALL CLASCL( 'General'00, NORMA, SMLNUM, M, NRHS, B,
217      $                      LDB, INFO )
218             ELSE
219                CALL XERBLA( 'CQRT15'1 )
220                RETURN
221             END IF
222          END IF
223       END IF
224 *
225       NORMA = SASUM( MN, S, 1 )
226       NORMB = CLANGE( 'One-norm', M, NRHS, B, LDB, DUMMY )
227 *
228       RETURN
229 *
230 *     End of CQRT15
231 *
232       END