1*> \brief \b DQRT12
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*  Definition:
9*  ===========
10*
11*       DOUBLE PRECISION FUNCTION DQRT12( M, N, A, LDA, S, WORK, LWORK )
12*
13*       .. Scalar Arguments ..
14*       INTEGER            LDA, LWORK, M, N
15*       ..
16*       .. Array Arguments ..
17*       DOUBLE PRECISION   A( LDA, * ), S( * ), WORK( LWORK )
18*       ..
19*
20*
21*> \par Purpose:
22*  =============
23*>
24*> \verbatim
25*>
26*> DQRT12 computes the singular values `svlues' of the upper trapezoid
27*> of A(1:M,1:N) and returns the ratio
28*>
29*>      || s - svlues||/(||svlues||*eps*max(M,N))
30*> \endverbatim
31*
32*  Arguments:
33*  ==========
34*
35*> \param[in] M
36*> \verbatim
37*>          M is INTEGER
38*>          The number of rows of the matrix A.
39*> \endverbatim
40*>
41*> \param[in] N
42*> \verbatim
43*>          N is INTEGER
44*>          The number of columns of the matrix A.
45*> \endverbatim
46*>
47*> \param[in] A
48*> \verbatim
49*>          A is DOUBLE PRECISION array, dimension (LDA,N)
50*>          The M-by-N matrix A. Only the upper trapezoid is referenced.
51*> \endverbatim
52*>
53*> \param[in] LDA
54*> \verbatim
55*>          LDA is INTEGER
56*>          The leading dimension of the array A.
57*> \endverbatim
58*>
59*> \param[in] S
60*> \verbatim
61*>          S is DOUBLE PRECISION array, dimension (min(M,N))
62*>          The singular values of the matrix A.
63*> \endverbatim
64*>
65*> \param[out] WORK
66*> \verbatim
67*>          WORK is DOUBLE PRECISION array, dimension (LWORK)
68*> \endverbatim
69*>
70*> \param[in] LWORK
71*> \verbatim
72*>          LWORK is INTEGER
73*>          The length of the array WORK. LWORK >= max(M*N + 4*min(M,N) +
74*>          max(M,N), M*N+2*MIN( M, N )+4*N).
75*> \endverbatim
76*
77*  Authors:
78*  ========
79*
80*> \author Univ. of Tennessee
81*> \author Univ. of California Berkeley
82*> \author Univ. of Colorado Denver
83*> \author NAG Ltd.
84*
85*> \ingroup double_lin
86*
87*  =====================================================================
88      DOUBLE PRECISION FUNCTION DQRT12( M, N, A, LDA, S, WORK, LWORK )
89*
90*  -- LAPACK test routine --
91*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
92*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
93*
94*     .. Scalar Arguments ..
95      INTEGER            LDA, LWORK, M, N
96*     ..
97*     .. Array Arguments ..
98      DOUBLE PRECISION   A( LDA, * ), S( * ), WORK( LWORK )
99*     ..
100*
101*  =====================================================================
102*
103*     .. Parameters ..
104      DOUBLE PRECISION   ZERO, ONE
105      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
106*     ..
107*     .. Local Scalars ..
108      INTEGER            I, INFO, ISCL, J, MN
109      DOUBLE PRECISION   ANRM, BIGNUM, NRMSVL, SMLNUM
110*     ..
111*     .. External Functions ..
112      DOUBLE PRECISION   DASUM, DLAMCH, DLANGE, DNRM2
113      EXTERNAL           DASUM, DLAMCH, DLANGE, DNRM2
114*     ..
115*     .. External Subroutines ..
116      EXTERNAL           DAXPY, DBDSQR, DGEBD2, DLABAD, DLASCL, DLASET,
117     $                   XERBLA
118*     ..
119*     .. Intrinsic Functions ..
120      INTRINSIC          DBLE, MAX, MIN
121*     ..
122*     .. Local Arrays ..
123      DOUBLE PRECISION   DUMMY( 1 )
124*     ..
125*     .. Executable Statements ..
126*
127      DQRT12 = ZERO
128*
129*     Test that enough workspace is supplied
130*
131      IF( LWORK.LT.MAX( M*N+4*MIN( M, N )+MAX( M, N ),
132     $                  M*N+2*MIN( M, N )+4*N) ) THEN
133         CALL XERBLA( 'DQRT12', 7 )
134         RETURN
135      END IF
136*
137*     Quick return if possible
138*
139      MN = MIN( M, N )
140      IF( MN.LE.ZERO )
141     $   RETURN
142*
143      NRMSVL = DNRM2( MN, S, 1 )
144*
145*     Copy upper triangle of A into work
146*
147      CALL DLASET( 'Full', M, N, ZERO, ZERO, WORK, M )
148      DO 20 J = 1, N
149         DO 10 I = 1, MIN( J, M )
150            WORK( ( J-1 )*M+I ) = A( I, J )
151   10    CONTINUE
152   20 CONTINUE
153*
154*     Get machine parameters
155*
156      SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' )
157      BIGNUM = ONE / SMLNUM
158      CALL DLABAD( SMLNUM, BIGNUM )
159*
160*     Scale work if max entry outside range [SMLNUM,BIGNUM]
161*
162      ANRM = DLANGE( 'M', M, N, WORK, M, DUMMY )
163      ISCL = 0
164      IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
165*
166*        Scale matrix norm up to SMLNUM
167*
168         CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, WORK, M, INFO )
169         ISCL = 1
170      ELSE IF( ANRM.GT.BIGNUM ) THEN
171*
172*        Scale matrix norm down to BIGNUM
173*
174         CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, WORK, M, INFO )
175         ISCL = 1
176      END IF
177*
178      IF( ANRM.NE.ZERO ) THEN
179*
180*        Compute SVD of work
181*
182         CALL DGEBD2( M, N, WORK, M, WORK( M*N+1 ), WORK( M*N+MN+1 ),
183     $                WORK( M*N+2*MN+1 ), WORK( M*N+3*MN+1 ),
184     $                WORK( M*N+4*MN+1 ), INFO )
185         CALL DBDSQR( 'Upper', MN, 0, 0, 0, WORK( M*N+1 ),
186     $                WORK( M*N+MN+1 ), DUMMY, MN, DUMMY, 1, DUMMY, MN,
187     $                WORK( M*N+2*MN+1 ), INFO )
188*
189         IF( ISCL.EQ.1 ) THEN
190            IF( ANRM.GT.BIGNUM ) THEN
191               CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MN, 1,
192     $                      WORK( M*N+1 ), MN, INFO )
193            END IF
194            IF( ANRM.LT.SMLNUM ) THEN
195               CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MN, 1,
196     $                      WORK( M*N+1 ), MN, INFO )
197            END IF
198         END IF
199*
200      ELSE
201*
202         DO 30 I = 1, MN
203            WORK( M*N+I ) = ZERO
204   30    CONTINUE
205      END IF
206*
207*     Compare s and singular values of work
208*
209      CALL DAXPY( MN, -ONE, S, 1, WORK( M*N+1 ), 1 )
210      DQRT12 = DASUM( MN, WORK( M*N+1 ), 1 ) /
211     $         ( DLAMCH( 'Epsilon' )*DBLE( MAX( M, N ) ) )
212      IF( NRMSVL.NE.ZERO )
213     $   DQRT12 = DQRT12 / NRMSVL
214*
215      RETURN
216*
217*     End of DQRT12
218*
219      END
220