1 DOUBLE PRECISION FUNCTION ZQRT12( 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 DOUBLE PRECISION RWORK( * ), S( * )
13 COMPLEX*16 A( LDA, * ), WORK( LWORK )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * ZQRT12 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*16 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) DOUBLE PRECISION array, dimension (min(M,N))
40 * The singular values of the matrix A.
41 *
42 * WORK (workspace) COMPLEX*16 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) DOUBLE PRECISION array, dimension (2*min(M,N))
49 *
50 * =====================================================================
51 *
52 * .. Parameters ..
53 DOUBLE PRECISION ZERO, ONE
54 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
55 * ..
56 * .. Local Scalars ..
57 INTEGER I, INFO, ISCL, J, MN
58 DOUBLE PRECISION ANRM, BIGNUM, NRMSVL, SMLNUM
59 * ..
60 * .. Local Arrays ..
61 DOUBLE PRECISION DUMMY( 1 )
62 * ..
63 * .. External Functions ..
64 DOUBLE PRECISION DASUM, DLAMCH, DNRM2, ZLANGE
65 EXTERNAL DASUM, DLAMCH, DNRM2, ZLANGE
66 * ..
67 * .. External Subroutines ..
68 EXTERNAL DAXPY, DBDSQR, DLABAD, DLASCL, XERBLA, ZGEBD2,
69 $ ZLASCL, ZLASET
70 * ..
71 * .. Intrinsic Functions ..
72 INTRINSIC DBLE, DCMPLX, MAX, MIN
73 * ..
74 * .. Executable Statements ..
75 *
76 ZQRT12 = 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( 'ZQRT12', 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 = DNRM2( MN, S, 1 )
92 *
93 * Copy upper triangle of A into work
94 *
95 CALL ZLASET( 'Full', M, N, DCMPLX( ZERO ), DCMPLX( ZERO ), WORK,
96 $ M )
97 DO 20 J = 1, N
98 DO 10 I = 1, MIN( J, M )
99 WORK( ( J-1 )*M+I ) = A( I, J )
100 10 CONTINUE
101 20 CONTINUE
102 *
103 * Get machine parameters
104 *
105 SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' )
106 BIGNUM = ONE / SMLNUM
107 CALL DLABAD( SMLNUM, BIGNUM )
108 *
109 * Scale work if max entry outside range [SMLNUM,BIGNUM]
110 *
111 ANRM = ZLANGE( 'M', M, N, WORK, M, DUMMY )
112 ISCL = 0
113 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
114 *
115 * Scale matrix norm up to SMLNUM
116 *
117 CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, WORK, M, INFO )
118 ISCL = 1
119 ELSE IF( ANRM.GT.BIGNUM ) THEN
120 *
121 * Scale matrix norm down to BIGNUM
122 *
123 CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, WORK, M, INFO )
124 ISCL = 1
125 END IF
126 *
127 IF( ANRM.NE.ZERO ) THEN
128 *
129 * Compute SVD of work
130 *
131 CALL ZGEBD2( M, N, WORK, M, RWORK( 1 ), RWORK( MN+1 ),
132 $ WORK( M*N+1 ), WORK( M*N+MN+1 ),
133 $ WORK( M*N+2*MN+1 ), INFO )
134 CALL DBDSQR( 'Upper', MN, 0, 0, 0, RWORK( 1 ), RWORK( MN+1 ),
135 $ DUMMY, MN, DUMMY, 1, DUMMY, MN, RWORK( 2*MN+1 ),
136 $ INFO )
137 *
138 IF( ISCL.EQ.1 ) THEN
139 IF( ANRM.GT.BIGNUM ) THEN
140 CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MN, 1, RWORK( 1 ),
141 $ MN, INFO )
142 END IF
143 IF( ANRM.LT.SMLNUM ) THEN
144 CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MN, 1, RWORK( 1 ),
145 $ MN, INFO )
146 END IF
147 END IF
148 *
149 ELSE
150 *
151 DO 30 I = 1, MN
152 RWORK( I ) = ZERO
153 30 CONTINUE
154 END IF
155 *
156 * Compare s and singular values of work
157 *
158 CALL DAXPY( MN, -ONE, S, 1, RWORK( 1 ), 1 )
159 ZQRT12 = DASUM( MN, RWORK( 1 ), 1 ) /
160 $ ( DLAMCH( 'Epsilon' )*DBLE( MAX( M, N ) ) )
161 IF( NRMSVL.NE.ZERO )
162 $ ZQRT12 = ZQRT12 / NRMSVL
163 *
164 RETURN
165 *
166 * End of ZQRT12
167 *
168 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 DOUBLE PRECISION RWORK( * ), S( * )
13 COMPLEX*16 A( LDA, * ), WORK( LWORK )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * ZQRT12 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*16 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) DOUBLE PRECISION array, dimension (min(M,N))
40 * The singular values of the matrix A.
41 *
42 * WORK (workspace) COMPLEX*16 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) DOUBLE PRECISION array, dimension (2*min(M,N))
49 *
50 * =====================================================================
51 *
52 * .. Parameters ..
53 DOUBLE PRECISION ZERO, ONE
54 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
55 * ..
56 * .. Local Scalars ..
57 INTEGER I, INFO, ISCL, J, MN
58 DOUBLE PRECISION ANRM, BIGNUM, NRMSVL, SMLNUM
59 * ..
60 * .. Local Arrays ..
61 DOUBLE PRECISION DUMMY( 1 )
62 * ..
63 * .. External Functions ..
64 DOUBLE PRECISION DASUM, DLAMCH, DNRM2, ZLANGE
65 EXTERNAL DASUM, DLAMCH, DNRM2, ZLANGE
66 * ..
67 * .. External Subroutines ..
68 EXTERNAL DAXPY, DBDSQR, DLABAD, DLASCL, XERBLA, ZGEBD2,
69 $ ZLASCL, ZLASET
70 * ..
71 * .. Intrinsic Functions ..
72 INTRINSIC DBLE, DCMPLX, MAX, MIN
73 * ..
74 * .. Executable Statements ..
75 *
76 ZQRT12 = 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( 'ZQRT12', 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 = DNRM2( MN, S, 1 )
92 *
93 * Copy upper triangle of A into work
94 *
95 CALL ZLASET( 'Full', M, N, DCMPLX( ZERO ), DCMPLX( ZERO ), WORK,
96 $ M )
97 DO 20 J = 1, N
98 DO 10 I = 1, MIN( J, M )
99 WORK( ( J-1 )*M+I ) = A( I, J )
100 10 CONTINUE
101 20 CONTINUE
102 *
103 * Get machine parameters
104 *
105 SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' )
106 BIGNUM = ONE / SMLNUM
107 CALL DLABAD( SMLNUM, BIGNUM )
108 *
109 * Scale work if max entry outside range [SMLNUM,BIGNUM]
110 *
111 ANRM = ZLANGE( 'M', M, N, WORK, M, DUMMY )
112 ISCL = 0
113 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
114 *
115 * Scale matrix norm up to SMLNUM
116 *
117 CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, WORK, M, INFO )
118 ISCL = 1
119 ELSE IF( ANRM.GT.BIGNUM ) THEN
120 *
121 * Scale matrix norm down to BIGNUM
122 *
123 CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, WORK, M, INFO )
124 ISCL = 1
125 END IF
126 *
127 IF( ANRM.NE.ZERO ) THEN
128 *
129 * Compute SVD of work
130 *
131 CALL ZGEBD2( M, N, WORK, M, RWORK( 1 ), RWORK( MN+1 ),
132 $ WORK( M*N+1 ), WORK( M*N+MN+1 ),
133 $ WORK( M*N+2*MN+1 ), INFO )
134 CALL DBDSQR( 'Upper', MN, 0, 0, 0, RWORK( 1 ), RWORK( MN+1 ),
135 $ DUMMY, MN, DUMMY, 1, DUMMY, MN, RWORK( 2*MN+1 ),
136 $ INFO )
137 *
138 IF( ISCL.EQ.1 ) THEN
139 IF( ANRM.GT.BIGNUM ) THEN
140 CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MN, 1, RWORK( 1 ),
141 $ MN, INFO )
142 END IF
143 IF( ANRM.LT.SMLNUM ) THEN
144 CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MN, 1, RWORK( 1 ),
145 $ MN, INFO )
146 END IF
147 END IF
148 *
149 ELSE
150 *
151 DO 30 I = 1, MN
152 RWORK( I ) = ZERO
153 30 CONTINUE
154 END IF
155 *
156 * Compare s and singular values of work
157 *
158 CALL DAXPY( MN, -ONE, S, 1, RWORK( 1 ), 1 )
159 ZQRT12 = DASUM( MN, RWORK( 1 ), 1 ) /
160 $ ( DLAMCH( 'Epsilon' )*DBLE( MAX( M, N ) ) )
161 IF( NRMSVL.NE.ZERO )
162 $ ZQRT12 = ZQRT12 / NRMSVL
163 *
164 RETURN
165 *
166 * End of ZQRT12
167 *
168 END