LAPACK: TESTING/LIN/sqrt15.f Source File

Go to the documentation of this file.
00001       SUBROUTINE SQRT15( SCALE, RKSEL, M, N, NRHS, A, LDA, B, LDB, S,
00002      $                   RANK, NORMA, NORMB, ISEED, WORK, LWORK )
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, LDB, LWORK, M, N, NRHS, RANK, RKSEL, SCALE
00010       REAL               NORMA, NORMB
00011 *     ..
00012 *     .. Array Arguments ..
00013       INTEGER            ISEED( 4 )
00014       REAL               A( LDA, * ), B( LDB, * ), S( * ), WORK( LWORK )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  SQRT15 generates a matrix with full or deficient rank and of various
00021 *  norms.
00022 *
00023 *  Arguments
00024 *  =========
00025 *
00026 *  SCALE   (input) INTEGER
00027 *          SCALE = 1: normally scaled matrix
00028 *          SCALE = 2: matrix scaled up
00029 *          SCALE = 3: matrix scaled down
00030 *
00031 *  RKSEL   (input) INTEGER
00032 *          RKSEL = 1: full rank matrix
00033 *          RKSEL = 2: rank-deficient matrix
00034 *
00035 *  M       (input) INTEGER
00036 *          The number of rows of the matrix A.
00037 *
00038 *  N       (input) INTEGER
00039 *          The number of columns of A.
00040 *
00041 *  NRHS    (input) INTEGER
00042 *          The number of columns of B.
00043 *
00044 *  A       (output) REAL array, dimension (LDA,N)
00045 *          The M-by-N matrix A.
00046 *
00047 *  LDA     (input) INTEGER
00048 *          The leading dimension of the array A.
00049 *
00050 *  B       (output) REAL array, dimension (LDB, NRHS)
00051 *          A matrix that is in the range space of matrix A.
00052 *
00053 *  LDB     (input) INTEGER
00054 *          The leading dimension of the array B.
00055 *
00056 *  S       (output) REAL array, dimension MIN(M,N)
00057 *          Singular values of A.
00058 *
00059 *  RANK    (output) INTEGER
00060 *          number of nonzero singular values of A.
00061 *
00062 *  NORMA   (output) REAL
00063 *          one-norm of A.
00064 *
00065 *  NORMB   (output) REAL
00066 *          one-norm of B.
00067 *
00068 *  ISEED   (input/output) integer array, dimension (4)
00069 *          seed for random number generator.
00070 *
00071 *  WORK    (workspace) REAL array, dimension (LWORK)
00072 *
00073 *  LWORK   (input) INTEGER
00074 *          length of work space required.
00075 *          LWORK >= MAX(M+MIN(M,N),NRHS*MIN(M,N),2*N+M)
00076 *
00077 *  =====================================================================
00078 *
00079 *     .. Parameters ..
00080       REAL               ZERO, ONE, TWO, SVMIN
00081       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0,
00082      $                   SVMIN = 0.1E0 )
00083 *     ..
00084 *     .. Local Scalars ..
00085       INTEGER            INFO, J, MN
00086       REAL               BIGNUM, EPS, SMLNUM, TEMP
00087 *     ..
00088 *     .. Local Arrays ..
00089       REAL               DUMMY( 1 )
00090 *     ..
00091 *     .. External Functions ..
00092       REAL               SASUM, SLAMCH, SLANGE, SLARND, SNRM2
00093       EXTERNAL           SASUM, SLAMCH, SLANGE, SLARND, SNRM2
00094 *     ..
00095 *     .. External Subroutines ..
00096       EXTERNAL           SGEMM, SLAORD, SLARF, SLARNV, SLAROR, SLASCL,
00097      $                   SLASET, SSCAL, XERBLA
00098 *     ..
00099 *     .. Intrinsic Functions ..
00100       INTRINSIC          ABS, MAX, MIN
00101 *     ..
00102 *     .. Executable Statements ..
00103 *
00104       MN = MIN( M, N )
00105       IF( LWORK.LT.MAX( M+MN, MN*NRHS, 2*N+M ) ) THEN
00106          CALL XERBLA( 'SQRT15', 16 )
00107          RETURN
00108       END IF
00109 *
00110       SMLNUM = SLAMCH( 'Safe minimum' )
00111       BIGNUM = ONE / SMLNUM
00112       EPS = SLAMCH( 'Epsilon' )
00113       SMLNUM = ( SMLNUM / EPS ) / EPS
00114       BIGNUM = ONE / SMLNUM
00115 *
00116 *     Determine rank and (unscaled) singular values
00117 *
00118       IF( RKSEL.EQ.1 ) THEN
00119          RANK = MN
00120       ELSE IF( RKSEL.EQ.2 ) THEN
00121          RANK = ( 3*MN ) / 4
00122          DO 10 J = RANK + 1, MN
00123             S( J ) = ZERO
00124    10    CONTINUE
00125       ELSE
00126          CALL XERBLA( 'SQRT15', 2 )
00127       END IF
00128 *
00129       IF( RANK.GT.0 ) THEN
00130 *
00131 *        Nontrivial case
00132 *
00133          S( 1 ) = ONE
00134          DO 30 J = 2, RANK
00135    20       CONTINUE
00136             TEMP = SLARND( 1, ISEED )
00137             IF( TEMP.GT.SVMIN ) THEN
00138                S( J ) = ABS( TEMP )
00139             ELSE
00140                GO TO 20
00141             END IF
00142    30    CONTINUE
00143          CALL SLAORD( 'Decreasing', RANK, S, 1 )
00144 *
00145 *        Generate 'rank' columns of a random orthogonal matrix in A
00146 *
00147          CALL SLARNV( 2, ISEED, M, WORK )
00148          CALL SSCAL( M, ONE / SNRM2( M, WORK, 1 ), WORK, 1 )
00149          CALL SLASET( 'Full', M, RANK, ZERO, ONE, A, LDA )
00150          CALL SLARF( 'Left', M, RANK, WORK, 1, TWO, A, LDA,
00151      $               WORK( M+1 ) )
00152 *
00153 *        workspace used: m+mn
00154 *
00155 *        Generate consistent rhs in the range space of A
00156 *
00157          CALL SLARNV( 2, ISEED, RANK*NRHS, WORK )
00158          CALL SGEMM( 'No transpose', 'No transpose', M, NRHS, RANK, ONE,
00159      $               A, LDA, WORK, RANK, ZERO, B, LDB )
00160 *
00161 *        work space used: <= mn *nrhs
00162 *
00163 *        generate (unscaled) matrix A
00164 *
00165          DO 40 J = 1, RANK
00166             CALL SSCAL( M, S( J ), A( 1, J ), 1 )
00167    40    CONTINUE
00168          IF( RANK.LT.N )
00169      $      CALL SLASET( 'Full', M, N-RANK, ZERO, ZERO, A( 1, RANK+1 ),
00170      $                   LDA )
00171          CALL SLAROR( 'Right', 'No initialization', M, N, A, LDA, ISEED,
00172      $                WORK, INFO )
00173 *
00174       ELSE
00175 *
00176 *        work space used 2*n+m
00177 *
00178 *        Generate null matrix and rhs
00179 *
00180          DO 50 J = 1, MN
00181             S( J ) = ZERO
00182    50    CONTINUE
00183          CALL SLASET( 'Full', M, N, ZERO, ZERO, A, LDA )
00184          CALL SLASET( 'Full', M, NRHS, ZERO, ZERO, B, LDB )
00185 *
00186       END IF
00187 *
00188 *     Scale the matrix
00189 *
00190       IF( SCALE.NE.1 ) THEN
00191          NORMA = SLANGE( 'Max', M, N, A, LDA, DUMMY )
00192          IF( NORMA.NE.ZERO ) THEN
00193             IF( SCALE.EQ.2 ) THEN
00194 *
00195 *              matrix scaled up
00196 *
00197                CALL SLASCL( 'General', 0, 0, NORMA, BIGNUM, M, N, A,
00198      $                      LDA, INFO )
00199                CALL SLASCL( 'General', 0, 0, NORMA, BIGNUM, MN, 1, S,
00200      $                      MN, INFO )
00201                CALL SLASCL( 'General', 0, 0, NORMA, BIGNUM, M, NRHS, B,
00202      $                      LDB, INFO )
00203             ELSE IF( SCALE.EQ.3 ) THEN
00204 *
00205 *              matrix scaled down
00206 *
00207                CALL SLASCL( 'General', 0, 0, NORMA, SMLNUM, M, N, A,
00208      $                      LDA, INFO )
00209                CALL SLASCL( 'General', 0, 0, NORMA, SMLNUM, MN, 1, S,
00210      $                      MN, INFO )
00211                CALL SLASCL( 'General', 0, 0, NORMA, SMLNUM, M, NRHS, B,
00212      $                      LDB, INFO )
00213             ELSE
00214                CALL XERBLA( 'SQRT15', 1 )
00215                RETURN
00216             END IF
00217          END IF
00218       END IF
00219 *
00220       NORMA = SASUM( MN, S, 1 )
00221       NORMB = SLANGE( 'One-norm', M, NRHS, B, LDB, DUMMY )
00222 *
00223       RETURN
00224 *
00225 *     End of SQRT15
00226 *
00227       END