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