LAPACK: TESTING/LIN/cqrt12.f Source File

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