LAPACK: TESTING/LIN/dqrt12.f Source File

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