 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 DOUBLE PRECISION FUNCTION ZQRT12( 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 DOUBLE PRECISION RWORK( * ), S( * ) 00013 COMPLEX*16 A( LDA, * ), WORK( LWORK ) 00014 * .. 00015 * 00016 * Purpose 00017 * ======= 00018 * 00019 * ZQRT12 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*16 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) DOUBLE PRECISION array, dimension (min(M,N)) 00040 * The singular values of the matrix A. 00041 * 00042 * WORK (workspace) COMPLEX*16 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) DOUBLE PRECISION array, dimension (2*min(M,N)) 00049 * 00050 * ===================================================================== 00051 * 00052 * .. Parameters .. 00053 DOUBLE PRECISION ZERO, ONE 00054 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 00055 * .. 00056 * .. Local Scalars .. 00057 INTEGER I, INFO, ISCL, J, MN 00058 DOUBLE PRECISION ANRM, BIGNUM, NRMSVL, SMLNUM 00059 * .. 00060 * .. Local Arrays .. 00061 DOUBLE PRECISION DUMMY( 1 ) 00062 * .. 00063 * .. External Functions .. 00064 DOUBLE PRECISION DASUM, DLAMCH, DNRM2, ZLANGE 00065 EXTERNAL DASUM, DLAMCH, DNRM2, ZLANGE 00066 * .. 00067 * .. External Subroutines .. 00068 EXTERNAL DAXPY, DBDSQR, DLABAD, DLASCL, XERBLA, ZGEBD2, 00069 $ ZLASCL, ZLASET 00070 * .. 00071 * .. Intrinsic Functions .. 00072 INTRINSIC DBLE, DCMPLX, MAX, MIN 00073 * .. 00074 * .. Executable Statements .. 00075 * 00076 ZQRT12 = 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( 'ZQRT12', 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 = DNRM2( MN, S, 1 ) 00092 * 00093 * Copy upper triangle of A into work 00094 * 00095 CALL ZLASET( 'Full', M, N, DCMPLX( ZERO ), DCMPLX( ZERO ), WORK, 00096 $ M ) 00097 DO 20 J = 1, N 00098 DO 10 I = 1, MIN( J, M ) 00099 WORK( ( J-1 )*M+I ) = A( I, J ) 00100 10 CONTINUE 00101 20 CONTINUE 00102 * 00103 * Get machine parameters 00104 * 00105 SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' ) 00106 BIGNUM = ONE / SMLNUM 00107 CALL DLABAD( SMLNUM, BIGNUM ) 00108 * 00109 * Scale work if max entry outside range [SMLNUM,BIGNUM] 00110 * 00111 ANRM = ZLANGE( 'M', M, N, WORK, M, DUMMY ) 00112 ISCL = 0 00113 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 00114 * 00115 * Scale matrix norm up to SMLNUM 00116 * 00117 CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, WORK, M, INFO ) 00118 ISCL = 1 00119 ELSE IF( ANRM.GT.BIGNUM ) THEN 00120 * 00121 * Scale matrix norm down to BIGNUM 00122 * 00123 CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, WORK, M, INFO ) 00124 ISCL = 1 00125 END IF 00126 * 00127 IF( ANRM.NE.ZERO ) THEN 00128 * 00129 * Compute SVD of work 00130 * 00131 CALL ZGEBD2( M, N, WORK, M, RWORK( 1 ), RWORK( MN+1 ), 00132 $ WORK( M*N+1 ), WORK( M*N+MN+1 ), 00133 $ WORK( M*N+2*MN+1 ), INFO ) 00134 CALL DBDSQR( 'Upper', MN, 0, 0, 0, RWORK( 1 ), RWORK( MN+1 ), 00135 $ DUMMY, MN, DUMMY, 1, DUMMY, MN, RWORK( 2*MN+1 ), 00136 $ INFO ) 00137 * 00138 IF( ISCL.EQ.1 ) THEN 00139 IF( ANRM.GT.BIGNUM ) THEN 00140 CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MN, 1, RWORK( 1 ), 00141 $ MN, INFO ) 00142 END IF 00143 IF( ANRM.LT.SMLNUM ) THEN 00144 CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MN, 1, RWORK( 1 ), 00145 $ MN, INFO ) 00146 END IF 00147 END IF 00148 * 00149 ELSE 00150 * 00151 DO 30 I = 1, MN 00152 RWORK( I ) = ZERO 00153 30 CONTINUE 00154 END IF 00155 * 00156 * Compare s and singular values of work 00157 * 00158 CALL DAXPY( MN, -ONE, S, 1, RWORK( 1 ), 1 ) 00159 ZQRT12 = DASUM( MN, RWORK( 1 ), 1 ) / 00160 $ ( DLAMCH( 'Epsilon' )*DBLE( MAX( M, N ) ) ) 00161 IF( NRMSVL.NE.ZERO ) 00162 $ ZQRT12 = ZQRT12 / NRMSVL 00163 * 00164 RETURN 00165 * 00166 * End of ZQRT12 00167 * 00168 END
1.7.3 
