![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DQRT12 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 * Definition: 00009 * =========== 00010 * 00011 * DOUBLE PRECISION FUNCTION DQRT12( M, N, A, LDA, S, WORK, LWORK ) 00012 * 00013 * .. Scalar Arguments .. 00014 * INTEGER LDA, LWORK, M, N 00015 * .. 00016 * .. Array Arguments .. 00017 * DOUBLE PRECISION A( LDA, * ), S( * ), WORK( LWORK ) 00018 * .. 00019 * 00020 * 00021 *> \par Purpose: 00022 * ============= 00023 *> 00024 *> \verbatim 00025 *> 00026 *> DQRT12 computes the singular values `svlues' of the upper trapezoid 00027 *> of A(1:M,1:N) and returns the ratio 00028 *> 00029 *> || s - svlues||/(||svlues||*eps*max(M,N)) 00030 *> \endverbatim 00031 * 00032 * Arguments: 00033 * ========== 00034 * 00035 *> \param[in] M 00036 *> \verbatim 00037 *> M is INTEGER 00038 *> The number of rows of the matrix A. 00039 *> \endverbatim 00040 *> 00041 *> \param[in] N 00042 *> \verbatim 00043 *> N is INTEGER 00044 *> The number of columns of the matrix A. 00045 *> \endverbatim 00046 *> 00047 *> \param[in] A 00048 *> \verbatim 00049 *> A is DOUBLE PRECISION array, dimension (LDA,N) 00050 *> The M-by-N matrix A. Only the upper trapezoid is referenced. 00051 *> \endverbatim 00052 *> 00053 *> \param[in] LDA 00054 *> \verbatim 00055 *> LDA is INTEGER 00056 *> The leading dimension of the array A. 00057 *> \endverbatim 00058 *> 00059 *> \param[in] S 00060 *> \verbatim 00061 *> S is DOUBLE PRECISION array, dimension (min(M,N)) 00062 *> The singular values of the matrix A. 00063 *> \endverbatim 00064 *> 00065 *> \param[out] WORK 00066 *> \verbatim 00067 *> WORK is DOUBLE PRECISION array, dimension (LWORK) 00068 *> \endverbatim 00069 *> 00070 *> \param[in] LWORK 00071 *> \verbatim 00072 *> LWORK is INTEGER 00073 *> The length of the array WORK. LWORK >= max(M*N + 4*min(M,N) + 00074 *> max(M,N), M*N+2*MIN( M, N )+4*N). 00075 *> \endverbatim 00076 * 00077 * Authors: 00078 * ======== 00079 * 00080 *> \author Univ. of Tennessee 00081 *> \author Univ. of California Berkeley 00082 *> \author Univ. of Colorado Denver 00083 *> \author NAG Ltd. 00084 * 00085 *> \date November 2011 00086 * 00087 *> \ingroup double_lin 00088 * 00089 * ===================================================================== 00090 DOUBLE PRECISION FUNCTION DQRT12( M, N, A, LDA, S, WORK, LWORK ) 00091 * 00092 * -- LAPACK test routine (version 3.4.0) -- 00093 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00094 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00095 * November 2011 00096 * 00097 * .. Scalar Arguments .. 00098 INTEGER LDA, LWORK, M, N 00099 * .. 00100 * .. Array Arguments .. 00101 DOUBLE PRECISION A( LDA, * ), S( * ), WORK( LWORK ) 00102 * .. 00103 * 00104 * ===================================================================== 00105 * 00106 * .. Parameters .. 00107 DOUBLE PRECISION ZERO, ONE 00108 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 00109 * .. 00110 * .. Local Scalars .. 00111 INTEGER I, INFO, ISCL, J, MN 00112 DOUBLE PRECISION ANRM, BIGNUM, NRMSVL, SMLNUM 00113 * .. 00114 * .. External Functions .. 00115 DOUBLE PRECISION DASUM, DLAMCH, DLANGE, DNRM2 00116 EXTERNAL DASUM, DLAMCH, DLANGE, DNRM2 00117 * .. 00118 * .. External Subroutines .. 00119 EXTERNAL DAXPY, DBDSQR, DGEBD2, DLABAD, DLASCL, DLASET, 00120 $ XERBLA 00121 * .. 00122 * .. Intrinsic Functions .. 00123 INTRINSIC DBLE, MAX, MIN 00124 * .. 00125 * .. Local Arrays .. 00126 DOUBLE PRECISION DUMMY( 1 ) 00127 * .. 00128 * .. Executable Statements .. 00129 * 00130 DQRT12 = ZERO 00131 * 00132 * Test that enough workspace is supplied 00133 * 00134 IF( LWORK.LT.MAX( M*N+4*MIN( M, N )+MAX( M, N ), 00135 $ M*N+2*MIN( M, N )+4*N) ) THEN 00136 CALL XERBLA( 'DQRT12', 7 ) 00137 RETURN 00138 END IF 00139 * 00140 * Quick return if possible 00141 * 00142 MN = MIN( M, N ) 00143 IF( MN.LE.ZERO ) 00144 $ RETURN 00145 * 00146 NRMSVL = DNRM2( MN, S, 1 ) 00147 * 00148 * Copy upper triangle of A into work 00149 * 00150 CALL DLASET( 'Full', M, N, ZERO, ZERO, WORK, M ) 00151 DO 20 J = 1, N 00152 DO 10 I = 1, MIN( J, M ) 00153 WORK( ( J-1 )*M+I ) = A( I, J ) 00154 10 CONTINUE 00155 20 CONTINUE 00156 * 00157 * Get machine parameters 00158 * 00159 SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' ) 00160 BIGNUM = ONE / SMLNUM 00161 CALL DLABAD( SMLNUM, BIGNUM ) 00162 * 00163 * Scale work if max entry outside range [SMLNUM,BIGNUM] 00164 * 00165 ANRM = DLANGE( 'M', M, N, WORK, M, DUMMY ) 00166 ISCL = 0 00167 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 00168 * 00169 * Scale matrix norm up to SMLNUM 00170 * 00171 CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, WORK, M, INFO ) 00172 ISCL = 1 00173 ELSE IF( ANRM.GT.BIGNUM ) THEN 00174 * 00175 * Scale matrix norm down to BIGNUM 00176 * 00177 CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, WORK, M, INFO ) 00178 ISCL = 1 00179 END IF 00180 * 00181 IF( ANRM.NE.ZERO ) THEN 00182 * 00183 * Compute SVD of work 00184 * 00185 CALL DGEBD2( M, N, WORK, M, WORK( M*N+1 ), WORK( M*N+MN+1 ), 00186 $ WORK( M*N+2*MN+1 ), WORK( M*N+3*MN+1 ), 00187 $ WORK( M*N+4*MN+1 ), INFO ) 00188 CALL DBDSQR( 'Upper', MN, 0, 0, 0, WORK( M*N+1 ), 00189 $ WORK( M*N+MN+1 ), DUMMY, MN, DUMMY, 1, DUMMY, MN, 00190 $ WORK( M*N+2*MN+1 ), INFO ) 00191 * 00192 IF( ISCL.EQ.1 ) THEN 00193 IF( ANRM.GT.BIGNUM ) THEN 00194 CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MN, 1, 00195 $ WORK( M*N+1 ), MN, INFO ) 00196 END IF 00197 IF( ANRM.LT.SMLNUM ) THEN 00198 CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MN, 1, 00199 $ WORK( M*N+1 ), MN, INFO ) 00200 END IF 00201 END IF 00202 * 00203 ELSE 00204 * 00205 DO 30 I = 1, MN 00206 WORK( M*N+I ) = ZERO 00207 30 CONTINUE 00208 END IF 00209 * 00210 * Compare s and singular values of work 00211 * 00212 CALL DAXPY( MN, -ONE, S, 1, WORK( M*N+1 ), 1 ) 00213 DQRT12 = DASUM( MN, WORK( M*N+1 ), 1 ) / 00214 $ ( DLAMCH( 'Epsilon' )*DBLE( MAX( M, N ) ) ) 00215 IF( NRMSVL.NE.ZERO ) 00216 $ DQRT12 = DQRT12 / NRMSVL 00217 * 00218 RETURN 00219 * 00220 * End of DQRT12 00221 * 00222 END