TESTING/LIN/sqrt12.f

   1 *> \brief \b SQRT12
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *  Definition:
   9 *  ===========
  10 *
  11 *       REAL             FUNCTION SQRT12( M, N, A, LDA, S, WORK, LWORK )
  12 *
  13 *       .. Scalar Arguments ..
  14 *       INTEGER            LDA, LWORK, M, N
  15 *       ..
  16 *       .. Array Arguments ..
  17 *       REAL               A( LDA, * ), S( * ), WORK( LWORK )
  18 *       ..
  19 *
  20 *
  21 *> \par Purpose:
  22 *  =============
  23 *>
  24 *> \verbatim
  25 *>
  26 *> SQRT12 computes the singular values `svlues' of the upper trapezoid
  27 *> of A(1:M,1:N) and returns the ratio
  28 *>
  29 *>      || s - svlues||/(||svlues||*eps*max(M,N))
  30 *> \endverbatim
  31 *
  32 *  Arguments:
  33 *  ==========
  34 *
  35 *> \param[in] M
  36 *> \verbatim
  37 *>          M is INTEGER
  38 *>          The number of rows of the matrix A.
  39 *> \endverbatim
  40 *>
  41 *> \param[in] N
  42 *> \verbatim
  43 *>          N is INTEGER
  44 *>          The number of columns of the matrix A.
  45 *> \endverbatim
  46 *>
  47 *> \param[in] A
  48 *> \verbatim
  49 *>          A is REAL array, dimension (LDA,N)
  50 *>          The M-by-N matrix A. Only the upper trapezoid is referenced.
  51 *> \endverbatim
  52 *>
  53 *> \param[in] LDA
  54 *> \verbatim
  55 *>          LDA is INTEGER
  56 *>          The leading dimension of the array A.
  57 *> \endverbatim
  58 *>
  59 *> \param[in] S
  60 *> \verbatim
  61 *>          S is REAL array, dimension (min(M,N))
  62 *>          The singular values of the matrix A.
  63 *> \endverbatim
  64 *>
  65 *> \param[out] WORK
  66 *> \verbatim
  67 *>          WORK is REAL array, dimension (LWORK)
  68 *> \endverbatim
  69 *>
  70 *> \param[in] LWORK
  71 *> \verbatim
  72 *>          LWORK is INTEGER
  73 *>          The length of the array WORK. LWORK >= max(M*N + 4*min(M,N) +
  74 *>          max(M,N), M*N+2*MIN( M, N )+4*N).
  75 *> \endverbatim
  76 *
  77 *  Authors:
  78 *  ========
  79 *
  80 *> \author Univ. of Tennessee
  81 *> \author Univ. of California Berkeley
  82 *> \author Univ. of Colorado Denver
  83 *> \author NAG Ltd.
  84 *
  85 *> \date November 2011
  86 *
  87 *> \ingroup single_lin
  88 *
  89 *  =====================================================================
  90       REAL             FUNCTION SQRT12( M, N, A, LDA, S, WORK, LWORK )
  91 *
  92 *  -- LAPACK test routine (version 3.4.0) --
  93 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  94 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  95 *     November 2011
  96 *
  97 *     .. Scalar Arguments ..
  98       INTEGER            LDA, LWORK, M, N
  99 *     ..
 100 *     .. Array Arguments ..
 101       REAL               A( LDA, * ), S( * ), WORK( LWORK )
 102 *     ..
 103 *
 104 *  =====================================================================
 105 *
 106 *     .. Parameters ..
 107       REAL               ZERO, ONE
 108       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
 109 *     ..
 110 *     .. Local Scalars ..
 111       INTEGER            I, INFO, ISCL, J, MN
 112       REAL               ANRM, BIGNUM, NRMSVL, SMLNUM
 113 *     ..
 114 *     .. External Functions ..
 115       REAL               SASUM, SLAMCH, SLANGE, SNRM2
 116       EXTERNAL           SASUM, SLAMCH, SLANGE, SNRM2
 117 *     ..
 118 *     .. External Subroutines ..
 119       EXTERNAL           SAXPY, SBDSQR, SGEBD2, SLABAD, SLASCL, SLASET,
 120      $                   XERBLA
 121 *     ..
 122 *     .. Intrinsic Functions ..
 123       INTRINSIC          MAX, MIN, REAL
 124 *     ..
 125 *     .. Local Arrays ..
 126       REAL               DUMMY( 1 )
 127 *     ..
 128 *     .. Executable Statements ..
 129 *
 130       SQRT12 = ZERO
 131 *
 132 *     Test that enough workspace is supplied
 133 *
 134       IF( LWORK.LT.MAX( M*N+4*MIN( M, N )+MAX( M, N ),
 135      $                  M*N+2*MIN( M, N )+4*N) ) THEN
 136          CALL XERBLA( 'SQRT12', 7 )
 137          RETURN
 138       END IF
 139 *
 140 *     Quick return if possible
 141 *
 142       MN = MIN( M, N )
 143       IF( MN.LE.ZERO )
 144      $   RETURN
 145 *
 146       NRMSVL = SNRM2( MN, S, 1 )
 147 *
 148 *     Copy upper triangle of A into work
 149 *
 150       CALL SLASET( 'Full', M, N, ZERO, ZERO, WORK, M )
 151       DO 20 J = 1, N
 152          DO 10 I = 1, MIN( J, M )
 153             WORK( ( J-1 )*M+I ) = A( I, J )
 154    10    CONTINUE
 155    20 CONTINUE
 156 *
 157 *     Get machine parameters
 158 *
 159       SMLNUM = SLAMCH( 'S' ) / SLAMCH( 'P' )
 160       BIGNUM = ONE / SMLNUM
 161       CALL SLABAD( SMLNUM, BIGNUM )
 162 *
 163 *     Scale work if max entry outside range [SMLNUM,BIGNUM]
 164 *
 165       ANRM = SLANGE( 'M', M, N, WORK, M, DUMMY )
 166       ISCL = 0
 167       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
 168 *
 169 *        Scale matrix norm up to SMLNUM
 170 *
 171          CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, WORK, M, INFO )
 172          ISCL = 1
 173       ELSE IF( ANRM.GT.BIGNUM ) THEN
 174 *
 175 *        Scale matrix norm down to BIGNUM
 176 *
 177          CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, WORK, M, INFO )
 178          ISCL = 1
 179       END IF
 180 *
 181       IF( ANRM.NE.ZERO ) THEN
 182 *
 183 *        Compute SVD of work
 184 *
 185          CALL SGEBD2( M, N, WORK, M, WORK( M*N+1 ), WORK( M*N+MN+1 ),
 186      $                WORK( M*N+2*MN+1 ), WORK( M*N+3*MN+1 ),
 187      $                WORK( M*N+4*MN+1 ), INFO )
 188          CALL SBDSQR( 'Upper', MN, 0, 0, 0, WORK( M*N+1 ),
 189      $                WORK( M*N+MN+1 ), DUMMY, MN, DUMMY, 1, DUMMY, MN,
 190      $                WORK( M*N+2*MN+1 ), INFO )
 191 *
 192          IF( ISCL.EQ.1 ) THEN
 193             IF( ANRM.GT.BIGNUM ) THEN
 194                CALL SLASCL( 'G', 0, 0, BIGNUM, ANRM, MN, 1,
 195      $                      WORK( M*N+1 ), MN, INFO )
 196             END IF
 197             IF( ANRM.LT.SMLNUM ) THEN
 198                CALL SLASCL( 'G', 0, 0, SMLNUM, ANRM, MN, 1,
 199      $                      WORK( M*N+1 ), MN, INFO )
 200             END IF
 201          END IF
 202 *
 203       ELSE
 204 *
 205          DO 30 I = 1, MN
 206             WORK( M*N+I ) = ZERO
 207    30    CONTINUE
 208       END IF
 209 *
 210 *     Compare s and singular values of work
 211 *
 212       CALL SAXPY( MN, -ONE, S, 1, WORK( M*N+1 ), 1 )
 213       SQRT12 = SASUM( MN, WORK( M*N+1 ), 1 ) /
 214      $         ( SLAMCH( 'Epsilon' )*REAL( MAX( M, N ) ) )
 215       IF( NRMSVL.NE.ZERO )
 216      $   SQRT12 = SQRT12 / NRMSVL
 217 *
 218       RETURN
 219 *
 220 *     End of SQRT12
 221 *
 222       END