 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 DOUBLE PRECISION FUNCTION DQRT12( M, N, A, LDA, S, WORK, LWORK ) 00002 * 00003 * -- LAPACK test routine (version 3.1.1) -- 00004 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00005 * January 2007 00006 * 00007 * .. Scalar Arguments .. 00008 INTEGER LDA, LWORK, M, N 00009 * .. 00010 * .. Array Arguments .. 00011 DOUBLE PRECISION A( LDA, * ), S( * ), WORK( LWORK ) 00012 * .. 00013 * 00014 * Purpose 00015 * ======= 00016 * 00017 * DQRT12 computes the singular values `svlues' of the upper trapezoid 00018 * of A(1:M,1:N) and returns the ratio 00019 * 00020 * || s - svlues||/(||svlues||*eps*max(M,N)) 00021 * 00022 * Arguments 00023 * ========= 00024 * 00025 * M (input) INTEGER 00026 * The number of rows of the matrix A. 00027 * 00028 * N (input) INTEGER 00029 * The number of columns of the matrix A. 00030 * 00031 * A (input) DOUBLE PRECISION array, dimension (LDA,N) 00032 * The M-by-N matrix A. Only the upper trapezoid is referenced. 00033 * 00034 * LDA (input) INTEGER 00035 * The leading dimension of the array A. 00036 * 00037 * S (input) DOUBLE PRECISION array, dimension (min(M,N)) 00038 * The singular values of the matrix A. 00039 * 00040 * WORK (workspace) DOUBLE PRECISION array, dimension (LWORK) 00041 * 00042 * LWORK (input) INTEGER 00043 * The length of the array WORK. LWORK >= max(M*N + 4*min(M,N) + 00044 * max(M,N), M*N+2*MIN( M, N )+4*N). 00045 * 00046 * ===================================================================== 00047 * 00048 * .. Parameters .. 00049 DOUBLE PRECISION ZERO, ONE 00050 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 00051 * .. 00052 * .. Local Scalars .. 00053 INTEGER I, INFO, ISCL, J, MN 00054 DOUBLE PRECISION ANRM, BIGNUM, NRMSVL, SMLNUM 00055 * .. 00056 * .. External Functions .. 00057 DOUBLE PRECISION DASUM, DLAMCH, DLANGE, DNRM2 00058 EXTERNAL DASUM, DLAMCH, DLANGE, DNRM2 00059 * .. 00060 * .. External Subroutines .. 00061 EXTERNAL DAXPY, DBDSQR, DGEBD2, DLABAD, DLASCL, DLASET, 00062 $ XERBLA 00063 * .. 00064 * .. Intrinsic Functions .. 00065 INTRINSIC DBLE, MAX, MIN 00066 * .. 00067 * .. Local Arrays .. 00068 DOUBLE PRECISION DUMMY( 1 ) 00069 * .. 00070 * .. Executable Statements .. 00071 * 00072 DQRT12 = ZERO 00073 * 00074 * Test that enough workspace is supplied 00075 * 00076 IF( LWORK.LT.MAX( M*N+4*MIN( M, N )+MAX( M, N ), 00077 $ M*N+2*MIN( M, N )+4*N) ) THEN 00078 CALL XERBLA( 'DQRT12', 7 ) 00079 RETURN 00080 END IF 00081 * 00082 * Quick return if possible 00083 * 00084 MN = MIN( M, N ) 00085 IF( MN.LE.ZERO ) 00086 $ RETURN 00087 * 00088 NRMSVL = DNRM2( MN, S, 1 ) 00089 * 00090 * Copy upper triangle of A into work 00091 * 00092 CALL DLASET( 'Full', M, N, ZERO, ZERO, WORK, M ) 00093 DO 20 J = 1, N 00094 DO 10 I = 1, MIN( J, M ) 00095 WORK( ( J-1 )*M+I ) = A( I, J ) 00096 10 CONTINUE 00097 20 CONTINUE 00098 * 00099 * Get machine parameters 00100 * 00101 SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' ) 00102 BIGNUM = ONE / SMLNUM 00103 CALL DLABAD( SMLNUM, BIGNUM ) 00104 * 00105 * Scale work if max entry outside range [SMLNUM,BIGNUM] 00106 * 00107 ANRM = DLANGE( 'M', M, N, WORK, M, DUMMY ) 00108 ISCL = 0 00109 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 00110 * 00111 * Scale matrix norm up to SMLNUM 00112 * 00113 CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, WORK, M, INFO ) 00114 ISCL = 1 00115 ELSE IF( ANRM.GT.BIGNUM ) THEN 00116 * 00117 * Scale matrix norm down to BIGNUM 00118 * 00119 CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, WORK, M, INFO ) 00120 ISCL = 1 00121 END IF 00122 * 00123 IF( ANRM.NE.ZERO ) THEN 00124 * 00125 * Compute SVD of work 00126 * 00127 CALL DGEBD2( M, N, WORK, M, WORK( M*N+1 ), WORK( M*N+MN+1 ), 00128 $ WORK( M*N+2*MN+1 ), WORK( M*N+3*MN+1 ), 00129 $ WORK( M*N+4*MN+1 ), INFO ) 00130 CALL DBDSQR( 'Upper', MN, 0, 0, 0, WORK( M*N+1 ), 00131 $ WORK( M*N+MN+1 ), DUMMY, MN, DUMMY, 1, DUMMY, MN, 00132 $ WORK( M*N+2*MN+1 ), INFO ) 00133 * 00134 IF( ISCL.EQ.1 ) THEN 00135 IF( ANRM.GT.BIGNUM ) THEN 00136 CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MN, 1, 00137 $ WORK( M*N+1 ), MN, INFO ) 00138 END IF 00139 IF( ANRM.LT.SMLNUM ) THEN 00140 CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MN, 1, 00141 $ WORK( M*N+1 ), MN, INFO ) 00142 END IF 00143 END IF 00144 * 00145 ELSE 00146 * 00147 DO 30 I = 1, MN 00148 WORK( M*N+I ) = ZERO 00149 30 CONTINUE 00150 END IF 00151 * 00152 * Compare s and singular values of work 00153 * 00154 CALL DAXPY( MN, -ONE, S, 1, WORK( M*N+1 ), 1 ) 00155 DQRT12 = DASUM( MN, WORK( M*N+1 ), 1 ) / 00156 $ ( DLAMCH( 'Epsilon' )*DBLE( MAX( M, N ) ) ) 00157 IF( NRMSVL.NE.ZERO ) 00158 $ DQRT12 = DQRT12 / NRMSVL 00159 * 00160 RETURN 00161 * 00162 * End of DQRT12 00163 * 00164 END
1.7.3 
