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