LAPACK: TESTING/LIN/cqrt15.f Source File

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