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