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*> \date November 2011
86*
87*> \ingroup double_lin
88*
89*  =====================================================================
90      DOUBLE PRECISION FUNCTION DQRT12( M, N, A, LDA, S, WORK, LWORK )
91*
92*  -- LAPACK test routine (version 3.4.0) --
93*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
94*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
95*     November 2011
96*
97*     .. Scalar Arguments ..
98      INTEGER            LDA, LWORK, M, N
99*     ..
100*     .. Array Arguments ..
101      DOUBLE PRECISION   A( LDA, * ), S( * ), WORK( LWORK )
102*     ..
103*
104*  =====================================================================
105*
106*     .. Parameters ..
107      DOUBLE PRECISION   ZERO, ONE
108      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
109*     ..
110*     .. Local Scalars ..
111      INTEGER            I, INFO, ISCL, J, MN
112      DOUBLE PRECISION   ANRM, BIGNUM, NRMSVL, SMLNUM
113*     ..
114*     .. External Functions ..
115      DOUBLE PRECISION   DASUM, DLAMCH, DLANGE, DNRM2
116      EXTERNAL           DASUM, DLAMCH, DLANGE, DNRM2
117*     ..
118*     .. External Subroutines ..
119      EXTERNAL           DAXPY, DBDSQR, DGEBD2, DLABAD, DLASCL, DLASET,
120     $                   XERBLA
121*     ..
122*     .. Intrinsic Functions ..
123      INTRINSIC          DBLE, MAX, MIN
124*     ..
125*     .. Local Arrays ..
126      DOUBLE PRECISION   DUMMY( 1 )
127*     ..
128*     .. Executable Statements ..
129*
130      DQRT12 = ZERO
131*
132*     Test that enough workspace is supplied
133*
134      IF( LWORK.LT.MAX( M*N+4*MIN( M, N )+MAX( M, N ),
135     $                  M*N+2*MIN( M, N )+4*N) ) THEN
136         CALL XERBLA( 'DQRT12', 7 )
137         RETURN
138      END IF
139*
140*     Quick return if possible
141*
142      MN = MIN( M, N )
143      IF( MN.LE.ZERO )
144     $   RETURN
145*
146      NRMSVL = DNRM2( MN, S, 1 )
147*
148*     Copy upper triangle of A into work
149*
150      CALL DLASET( 'Full', M, N, ZERO, ZERO, WORK, M )
151      DO 20 J = 1, N
152         DO 10 I = 1, MIN( J, M )
153            WORK( ( J-1 )*M+I ) = A( I, J )
154   10    CONTINUE
155   20 CONTINUE
156*
157*     Get machine parameters
158*
159      SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' )
160      BIGNUM = ONE / SMLNUM
161      CALL DLABAD( SMLNUM, BIGNUM )
162*
163*     Scale work if max entry outside range [SMLNUM,BIGNUM]
164*
165      ANRM = DLANGE( 'M', M, N, WORK, M, DUMMY )
166      ISCL = 0
167      IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
168*
169*        Scale matrix norm up to SMLNUM
170*
171         CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, WORK, M, INFO )
172         ISCL = 1
173      ELSE IF( ANRM.GT.BIGNUM ) THEN
174*
175*        Scale matrix norm down to BIGNUM
176*
177         CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, WORK, M, INFO )
178         ISCL = 1
179      END IF
180*
181      IF( ANRM.NE.ZERO ) THEN
182*
183*        Compute SVD of work
184*
185         CALL DGEBD2( M, N, WORK, M, WORK( M*N+1 ), WORK( M*N+MN+1 ),
186     $                WORK( M*N+2*MN+1 ), WORK( M*N+3*MN+1 ),
187     $                WORK( M*N+4*MN+1 ), INFO )
188         CALL DBDSQR( 'Upper', MN, 0, 0, 0, WORK( M*N+1 ),
189     $                WORK( M*N+MN+1 ), DUMMY, MN, DUMMY, 1, DUMMY, MN,
190     $                WORK( M*N+2*MN+1 ), INFO )
191*
192         IF( ISCL.EQ.1 ) THEN
193            IF( ANRM.GT.BIGNUM ) THEN
194               CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MN, 1,
195     $                      WORK( M*N+1 ), MN, INFO )
196            END IF
197            IF( ANRM.LT.SMLNUM ) THEN
198               CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MN, 1,
199     $                      WORK( M*N+1 ), MN, INFO )
200            END IF
201         END IF
202*
203      ELSE
204*
205         DO 30 I = 1, MN
206            WORK( M*N+I ) = ZERO
207   30    CONTINUE
208      END IF
209*
210*     Compare s and singular values of work
211*
212      CALL DAXPY( MN, -ONE, S, 1, WORK( M*N+1 ), 1 )
213      DQRT12 = DASUM( MN, WORK( M*N+1 ), 1 ) /
214     $         ( DLAMCH( 'Epsilon' )*DBLE( MAX( M, N ) ) )
215      IF( NRMSVL.NE.ZERO )
216     $   DQRT12 = DQRT12 / NRMSVL
217*
218      RETURN
219*
220*     End of DQRT12
221*
222      END
223