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