LAPACK: TESTING/LIN/zqrt12.f Source File

Go to the documentation of this file.
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