1 REAL FUNCTION CQRT12( M, N, A, LDA, S, WORK, LWORK,
2 $ RWORK )
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, LWORK, M, N
10 * ..
11 * .. Array Arguments ..
12 REAL RWORK( * ), S( * )
13 COMPLEX A( LDA, * ), WORK( LWORK )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * CQRT12 computes the singular values `svlues' of the upper trapezoid
20 * of A(1:M,1:N) and returns the ratio
21 *
22 * || s - svlues||/(||svlues||*eps*max(M,N))
23 *
24 * Arguments
25 * =========
26 *
27 * M (input) INTEGER
28 * The number of rows of the matrix A.
29 *
30 * N (input) INTEGER
31 * The number of columns of the matrix A.
32 *
33 * A (input) COMPLEX array, dimension (LDA,N)
34 * The M-by-N matrix A. Only the upper trapezoid is referenced.
35 *
36 * LDA (input) INTEGER
37 * The leading dimension of the array A.
38 *
39 * S (input) REAL array, dimension (min(M,N))
40 * The singular values of the matrix A.
41 *
42 * WORK (workspace) COMPLEX array, dimension (LWORK)
43 *
44 * LWORK (input) INTEGER
45 * The length of the array WORK. LWORK >= M*N + 2*min(M,N) +
46 * max(M,N).
47 *
48 * RWORK (workspace) REAL array, dimension (4*min(M,N))
49 *
50 * =====================================================================
51 *
52 * .. Parameters ..
53 REAL ZERO, ONE
54 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 )
55 * ..
56 * .. Local Scalars ..
57 INTEGER I, INFO, ISCL, J, MN
58 REAL ANRM, BIGNUM, NRMSVL, SMLNUM
59 * ..
60 * .. Local Arrays ..
61 REAL DUMMY( 1 )
62 * ..
63 * .. External Functions ..
64 REAL CLANGE, SASUM, SLAMCH, SNRM2
65 EXTERNAL CLANGE, SASUM, SLAMCH, SNRM2
66 * ..
67 * .. External Subroutines ..
68 EXTERNAL CGEBD2, CLASCL, CLASET, SAXPY, SBDSQR, SLABAD,
69 $ SLASCL, XERBLA
70 * ..
71 * .. Intrinsic Functions ..
72 INTRINSIC CMPLX, MAX, MIN, REAL
73 * ..
74 * .. Executable Statements ..
75 *
76 CQRT12 = ZERO
77 *
78 * Test that enough workspace is supplied
79 *
80 IF( LWORK.LT.M*N+2*MIN( M, N )+MAX( M, N ) ) THEN
81 CALL XERBLA( 'CQRT12', 7 )
82 RETURN
83 END IF
84 *
85 * Quick return if possible
86 *
87 MN = MIN( M, N )
88 IF( MN.LE.ZERO )
89 $ RETURN
90 *
91 NRMSVL = SNRM2( MN, S, 1 )
92 *
93 * Copy upper triangle of A into work
94 *
95 CALL CLASET( 'Full', M, N, CMPLX( ZERO ), CMPLX( ZERO ), WORK, M )
96 DO 20 J = 1, N
97 DO 10 I = 1, MIN( J, M )
98 WORK( ( J-1 )*M+I ) = A( I, J )
99 10 CONTINUE
100 20 CONTINUE
101 *
102 * Get machine parameters
103 *
104 SMLNUM = SLAMCH( 'S' ) / SLAMCH( 'P' )
105 BIGNUM = ONE / SMLNUM
106 CALL SLABAD( SMLNUM, BIGNUM )
107 *
108 * Scale work if max entry outside range [SMLNUM,BIGNUM]
109 *
110 ANRM = CLANGE( 'M', M, N, WORK, M, DUMMY )
111 ISCL = 0
112 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
113 *
114 * Scale matrix norm up to SMLNUM
115 *
116 CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, WORK, M, INFO )
117 ISCL = 1
118 ELSE IF( ANRM.GT.BIGNUM ) THEN
119 *
120 * Scale matrix norm down to BIGNUM
121 *
122 CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, WORK, M, INFO )
123 ISCL = 1
124 END IF
125 *
126 IF( ANRM.NE.ZERO ) THEN
127 *
128 * Compute SVD of work
129 *
130 CALL CGEBD2( M, N, WORK, M, RWORK( 1 ), RWORK( MN+1 ),
131 $ WORK( M*N+1 ), WORK( M*N+MN+1 ),
132 $ WORK( M*N+2*MN+1 ), INFO )
133 CALL SBDSQR( 'Upper', MN, 0, 0, 0, RWORK( 1 ), RWORK( MN+1 ),
134 $ DUMMY, MN, DUMMY, 1, DUMMY, MN, RWORK( 2*MN+1 ),
135 $ INFO )
136 *
137 IF( ISCL.EQ.1 ) THEN
138 IF( ANRM.GT.BIGNUM ) THEN
139 CALL SLASCL( 'G', 0, 0, BIGNUM, ANRM, MN, 1, RWORK( 1 ),
140 $ MN, INFO )
141 END IF
142 IF( ANRM.LT.SMLNUM ) THEN
143 CALL SLASCL( 'G', 0, 0, SMLNUM, ANRM, MN, 1, RWORK( 1 ),
144 $ MN, INFO )
145 END IF
146 END IF
147 *
148 ELSE
149 *
150 DO 30 I = 1, MN
151 RWORK( I ) = ZERO
152 30 CONTINUE
153 END IF
154 *
155 * Compare s and singular values of work
156 *
157 CALL SAXPY( MN, -ONE, S, 1, RWORK( 1 ), 1 )
158 CQRT12 = SASUM( MN, RWORK( 1 ), 1 ) /
159 $ ( SLAMCH( 'Epsilon' )*REAL( MAX( M, N ) ) )
160 IF( NRMSVL.NE.ZERO )
161 $ CQRT12 = CQRT12 / NRMSVL
162 *
163 RETURN
164 *
165 * End of CQRT12
166 *
167 END
2 $ RWORK )
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, LWORK, M, N
10 * ..
11 * .. Array Arguments ..
12 REAL RWORK( * ), S( * )
13 COMPLEX A( LDA, * ), WORK( LWORK )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * CQRT12 computes the singular values `svlues' of the upper trapezoid
20 * of A(1:M,1:N) and returns the ratio
21 *
22 * || s - svlues||/(||svlues||*eps*max(M,N))
23 *
24 * Arguments
25 * =========
26 *
27 * M (input) INTEGER
28 * The number of rows of the matrix A.
29 *
30 * N (input) INTEGER
31 * The number of columns of the matrix A.
32 *
33 * A (input) COMPLEX array, dimension (LDA,N)
34 * The M-by-N matrix A. Only the upper trapezoid is referenced.
35 *
36 * LDA (input) INTEGER
37 * The leading dimension of the array A.
38 *
39 * S (input) REAL array, dimension (min(M,N))
40 * The singular values of the matrix A.
41 *
42 * WORK (workspace) COMPLEX array, dimension (LWORK)
43 *
44 * LWORK (input) INTEGER
45 * The length of the array WORK. LWORK >= M*N + 2*min(M,N) +
46 * max(M,N).
47 *
48 * RWORK (workspace) REAL array, dimension (4*min(M,N))
49 *
50 * =====================================================================
51 *
52 * .. Parameters ..
53 REAL ZERO, ONE
54 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 )
55 * ..
56 * .. Local Scalars ..
57 INTEGER I, INFO, ISCL, J, MN
58 REAL ANRM, BIGNUM, NRMSVL, SMLNUM
59 * ..
60 * .. Local Arrays ..
61 REAL DUMMY( 1 )
62 * ..
63 * .. External Functions ..
64 REAL CLANGE, SASUM, SLAMCH, SNRM2
65 EXTERNAL CLANGE, SASUM, SLAMCH, SNRM2
66 * ..
67 * .. External Subroutines ..
68 EXTERNAL CGEBD2, CLASCL, CLASET, SAXPY, SBDSQR, SLABAD,
69 $ SLASCL, XERBLA
70 * ..
71 * .. Intrinsic Functions ..
72 INTRINSIC CMPLX, MAX, MIN, REAL
73 * ..
74 * .. Executable Statements ..
75 *
76 CQRT12 = ZERO
77 *
78 * Test that enough workspace is supplied
79 *
80 IF( LWORK.LT.M*N+2*MIN( M, N )+MAX( M, N ) ) THEN
81 CALL XERBLA( 'CQRT12', 7 )
82 RETURN
83 END IF
84 *
85 * Quick return if possible
86 *
87 MN = MIN( M, N )
88 IF( MN.LE.ZERO )
89 $ RETURN
90 *
91 NRMSVL = SNRM2( MN, S, 1 )
92 *
93 * Copy upper triangle of A into work
94 *
95 CALL CLASET( 'Full', M, N, CMPLX( ZERO ), CMPLX( ZERO ), WORK, M )
96 DO 20 J = 1, N
97 DO 10 I = 1, MIN( J, M )
98 WORK( ( J-1 )*M+I ) = A( I, J )
99 10 CONTINUE
100 20 CONTINUE
101 *
102 * Get machine parameters
103 *
104 SMLNUM = SLAMCH( 'S' ) / SLAMCH( 'P' )
105 BIGNUM = ONE / SMLNUM
106 CALL SLABAD( SMLNUM, BIGNUM )
107 *
108 * Scale work if max entry outside range [SMLNUM,BIGNUM]
109 *
110 ANRM = CLANGE( 'M', M, N, WORK, M, DUMMY )
111 ISCL = 0
112 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
113 *
114 * Scale matrix norm up to SMLNUM
115 *
116 CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, WORK, M, INFO )
117 ISCL = 1
118 ELSE IF( ANRM.GT.BIGNUM ) THEN
119 *
120 * Scale matrix norm down to BIGNUM
121 *
122 CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, WORK, M, INFO )
123 ISCL = 1
124 END IF
125 *
126 IF( ANRM.NE.ZERO ) THEN
127 *
128 * Compute SVD of work
129 *
130 CALL CGEBD2( M, N, WORK, M, RWORK( 1 ), RWORK( MN+1 ),
131 $ WORK( M*N+1 ), WORK( M*N+MN+1 ),
132 $ WORK( M*N+2*MN+1 ), INFO )
133 CALL SBDSQR( 'Upper', MN, 0, 0, 0, RWORK( 1 ), RWORK( MN+1 ),
134 $ DUMMY, MN, DUMMY, 1, DUMMY, MN, RWORK( 2*MN+1 ),
135 $ INFO )
136 *
137 IF( ISCL.EQ.1 ) THEN
138 IF( ANRM.GT.BIGNUM ) THEN
139 CALL SLASCL( 'G', 0, 0, BIGNUM, ANRM, MN, 1, RWORK( 1 ),
140 $ MN, INFO )
141 END IF
142 IF( ANRM.LT.SMLNUM ) THEN
143 CALL SLASCL( 'G', 0, 0, SMLNUM, ANRM, MN, 1, RWORK( 1 ),
144 $ MN, INFO )
145 END IF
146 END IF
147 *
148 ELSE
149 *
150 DO 30 I = 1, MN
151 RWORK( I ) = ZERO
152 30 CONTINUE
153 END IF
154 *
155 * Compare s and singular values of work
156 *
157 CALL SAXPY( MN, -ONE, S, 1, RWORK( 1 ), 1 )
158 CQRT12 = SASUM( MN, RWORK( 1 ), 1 ) /
159 $ ( SLAMCH( 'Epsilon' )*REAL( MAX( M, N ) ) )
160 IF( NRMSVL.NE.ZERO )
161 $ CQRT12 = CQRT12 / NRMSVL
162 *
163 RETURN
164 *
165 * End of CQRT12
166 *
167 END