1 SUBROUTINE ZQRT15( 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 DOUBLE PRECISION NORMA, NORMB
11 * ..
12 * .. Array Arguments ..
13 INTEGER ISEED( 4 )
14 DOUBLE PRECISION S( * )
15 COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( LWORK )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * ZQRT15 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*16 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*16 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) DOUBLE PRECISION 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) DOUBLE PRECISION
64 * one-norm norm of A.
65 *
66 * NORMB (output) DOUBLE PRECISION
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*16 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 DOUBLE PRECISION ZERO, ONE, TWO, SVMIN
82 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0,
83 $ SVMIN = 0.1D+0 )
84 COMPLEX*16 CZERO, CONE
85 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
86 $ CONE = ( 1.0D+0, 0.0D+0 ) )
87 * ..
88 * .. Local Scalars ..
89 INTEGER INFO, J, MN
90 DOUBLE PRECISION BIGNUM, EPS, SMLNUM, TEMP
91 * ..
92 * .. Local Arrays ..
93 DOUBLE PRECISION DUMMY( 1 )
94 * ..
95 * .. External Functions ..
96 DOUBLE PRECISION DASUM, DLAMCH, DLARND, DZNRM2, ZLANGE
97 EXTERNAL DASUM, DLAMCH, DLARND, DZNRM2, ZLANGE
98 * ..
99 * .. External Subroutines ..
100 EXTERNAL DLABAD, DLAORD, DLASCL, XERBLA, ZDSCAL, ZGEMM,
101 $ ZLARF, ZLARNV, ZLAROR, ZLASCL, ZLASET
102 * ..
103 * .. Intrinsic Functions ..
104 INTRINSIC ABS, DCMPLX, MAX, MIN
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( 'ZQRT15', 16 )
111 RETURN
112 END IF
113 *
114 SMLNUM = DLAMCH( 'Safe minimum' )
115 BIGNUM = ONE / SMLNUM
116 CALL DLABAD( SMLNUM, BIGNUM )
117 EPS = DLAMCH( '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( 'ZQRT15', 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 = DLARND( 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 DLAORD( 'Decreasing', RANK, S, 1 )
149 *
150 * Generate 'rank' columns of a random orthogonal matrix in A
151 *
152 CALL ZLARNV( 2, ISEED, M, WORK )
153 CALL ZDSCAL( M, ONE / DZNRM2( M, WORK, 1 ), WORK, 1 )
154 CALL ZLASET( 'Full', M, RANK, CZERO, CONE, A, LDA )
155 CALL ZLARF( 'Left', M, RANK, WORK, 1, DCMPLX( 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 ZLARNV( 2, ISEED, RANK*NRHS, WORK )
163 CALL ZGEMM( '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 ZDSCAL( M, S( J ), A( 1, J ), 1 )
172 40 CONTINUE
173 IF( RANK.LT.N )
174 $ CALL ZLASET( 'Full', M, N-RANK, CZERO, CZERO,
175 $ A( 1, RANK+1 ), LDA )
176 CALL ZLAROR( '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 ZLASET( 'Full', M, N, CZERO, CZERO, A, LDA )
189 CALL ZLASET( 'Full', M, NRHS, CZERO, CZERO, B, LDB )
190 *
191 END IF
192 *
193 * Scale the matrix
194 *
195 IF( SCALE.NE.1 ) THEN
196 NORMA = ZLANGE( 'Max', M, N, A, LDA, DUMMY )
197 IF( NORMA.NE.ZERO ) THEN
198 IF( SCALE.EQ.2 ) THEN
199 *
200 * matrix scaled up
201 *
202 CALL ZLASCL( 'General', 0, 0, NORMA, BIGNUM, M, N, A,
203 $ LDA, INFO )
204 CALL DLASCL( 'General', 0, 0, NORMA, BIGNUM, MN, 1, S,
205 $ MN, INFO )
206 CALL ZLASCL( 'General', 0, 0, NORMA, BIGNUM, M, NRHS, B,
207 $ LDB, INFO )
208 ELSE IF( SCALE.EQ.3 ) THEN
209 *
210 * matrix scaled down
211 *
212 CALL ZLASCL( 'General', 0, 0, NORMA, SMLNUM, M, N, A,
213 $ LDA, INFO )
214 CALL DLASCL( 'General', 0, 0, NORMA, SMLNUM, MN, 1, S,
215 $ MN, INFO )
216 CALL ZLASCL( 'General', 0, 0, NORMA, SMLNUM, M, NRHS, B,
217 $ LDB, INFO )
218 ELSE
219 CALL XERBLA( 'ZQRT15', 1 )
220 RETURN
221 END IF
222 END IF
223 END IF
224 *
225 NORMA = DASUM( MN, S, 1 )
226 NORMB = ZLANGE( 'One-norm', M, NRHS, B, LDB, DUMMY )
227 *
228 RETURN
229 *
230 * End of ZQRT15
231 *
232 END
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 DOUBLE PRECISION NORMA, NORMB
11 * ..
12 * .. Array Arguments ..
13 INTEGER ISEED( 4 )
14 DOUBLE PRECISION S( * )
15 COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( LWORK )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * ZQRT15 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*16 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*16 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) DOUBLE PRECISION 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) DOUBLE PRECISION
64 * one-norm norm of A.
65 *
66 * NORMB (output) DOUBLE PRECISION
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*16 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 DOUBLE PRECISION ZERO, ONE, TWO, SVMIN
82 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0,
83 $ SVMIN = 0.1D+0 )
84 COMPLEX*16 CZERO, CONE
85 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
86 $ CONE = ( 1.0D+0, 0.0D+0 ) )
87 * ..
88 * .. Local Scalars ..
89 INTEGER INFO, J, MN
90 DOUBLE PRECISION BIGNUM, EPS, SMLNUM, TEMP
91 * ..
92 * .. Local Arrays ..
93 DOUBLE PRECISION DUMMY( 1 )
94 * ..
95 * .. External Functions ..
96 DOUBLE PRECISION DASUM, DLAMCH, DLARND, DZNRM2, ZLANGE
97 EXTERNAL DASUM, DLAMCH, DLARND, DZNRM2, ZLANGE
98 * ..
99 * .. External Subroutines ..
100 EXTERNAL DLABAD, DLAORD, DLASCL, XERBLA, ZDSCAL, ZGEMM,
101 $ ZLARF, ZLARNV, ZLAROR, ZLASCL, ZLASET
102 * ..
103 * .. Intrinsic Functions ..
104 INTRINSIC ABS, DCMPLX, MAX, MIN
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( 'ZQRT15', 16 )
111 RETURN
112 END IF
113 *
114 SMLNUM = DLAMCH( 'Safe minimum' )
115 BIGNUM = ONE / SMLNUM
116 CALL DLABAD( SMLNUM, BIGNUM )
117 EPS = DLAMCH( '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( 'ZQRT15', 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 = DLARND( 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 DLAORD( 'Decreasing', RANK, S, 1 )
149 *
150 * Generate 'rank' columns of a random orthogonal matrix in A
151 *
152 CALL ZLARNV( 2, ISEED, M, WORK )
153 CALL ZDSCAL( M, ONE / DZNRM2( M, WORK, 1 ), WORK, 1 )
154 CALL ZLASET( 'Full', M, RANK, CZERO, CONE, A, LDA )
155 CALL ZLARF( 'Left', M, RANK, WORK, 1, DCMPLX( 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 ZLARNV( 2, ISEED, RANK*NRHS, WORK )
163 CALL ZGEMM( '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 ZDSCAL( M, S( J ), A( 1, J ), 1 )
172 40 CONTINUE
173 IF( RANK.LT.N )
174 $ CALL ZLASET( 'Full', M, N-RANK, CZERO, CZERO,
175 $ A( 1, RANK+1 ), LDA )
176 CALL ZLAROR( '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 ZLASET( 'Full', M, N, CZERO, CZERO, A, LDA )
189 CALL ZLASET( 'Full', M, NRHS, CZERO, CZERO, B, LDB )
190 *
191 END IF
192 *
193 * Scale the matrix
194 *
195 IF( SCALE.NE.1 ) THEN
196 NORMA = ZLANGE( 'Max', M, N, A, LDA, DUMMY )
197 IF( NORMA.NE.ZERO ) THEN
198 IF( SCALE.EQ.2 ) THEN
199 *
200 * matrix scaled up
201 *
202 CALL ZLASCL( 'General', 0, 0, NORMA, BIGNUM, M, N, A,
203 $ LDA, INFO )
204 CALL DLASCL( 'General', 0, 0, NORMA, BIGNUM, MN, 1, S,
205 $ MN, INFO )
206 CALL ZLASCL( 'General', 0, 0, NORMA, BIGNUM, M, NRHS, B,
207 $ LDB, INFO )
208 ELSE IF( SCALE.EQ.3 ) THEN
209 *
210 * matrix scaled down
211 *
212 CALL ZLASCL( 'General', 0, 0, NORMA, SMLNUM, M, N, A,
213 $ LDA, INFO )
214 CALL DLASCL( 'General', 0, 0, NORMA, SMLNUM, MN, 1, S,
215 $ MN, INFO )
216 CALL ZLASCL( 'General', 0, 0, NORMA, SMLNUM, M, NRHS, B,
217 $ LDB, INFO )
218 ELSE
219 CALL XERBLA( 'ZQRT15', 1 )
220 RETURN
221 END IF
222 END IF
223 END IF
224 *
225 NORMA = DASUM( MN, S, 1 )
226 NORMB = ZLANGE( 'One-norm', M, NRHS, B, LDB, DUMMY )
227 *
228 RETURN
229 *
230 * End of ZQRT15
231 *
232 END