LAPACK: TESTING/LIN/clatsy.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CLATSY( UPLO, N, X, LDX, ISEED )
00002 *
00003 *  -- LAPACK auxiliary test routine (version 3.1) --
00004 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00005 *     November 2006
00006 *
00007 *     .. Scalar Arguments ..
00008       CHARACTER          UPLO
00009       INTEGER            LDX, N
00010 *     ..
00011 *     .. Array Arguments ..
00012       INTEGER            ISEED( * )
00013       COMPLEX            X( LDX, * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  CLATSY generates a special test matrix for the complex symmetric
00020 *  (indefinite) factorization.  The pivot blocks of the generated matrix
00021 *  will be in the following order:
00022 *     2x2 pivot block, non diagonalizable
00023 *     1x1 pivot block
00024 *     2x2 pivot block, diagonalizable
00025 *     (cycle repeats)
00026 *  A row interchange is required for each non-diagonalizable 2x2 block.
00027 *
00028 *  Arguments
00029 *  =========
00030 *
00031 *  UPLO    (input) CHARACTER
00032 *          Specifies whether the generated matrix is to be upper or
00033 *          lower triangular.
00034 *          = 'U':  Upper triangular
00035 *          = 'L':  Lower triangular
00036 *
00037 *  N       (input) INTEGER
00038 *          The dimension of the matrix to be generated.
00039 *
00040 *  X       (output) COMPLEX array, dimension (LDX,N)
00041 *          The generated matrix, consisting of 3x3 and 2x2 diagonal
00042 *          blocks which result in the pivot sequence given above.
00043 *          The matrix outside of these diagonal blocks is zero.
00044 *
00045 *  LDX     (input) INTEGER
00046 *          The leading dimension of the array X.
00047 *
00048 *  ISEED   (input/output) INTEGER array, dimension (4)
00049 *          On entry, the seed for the random number generator.  The last
00050 *          of the four integers must be odd.  (modified on exit)
00051 *
00052 *  =====================================================================
00053 *
00054 *     .. Parameters ..
00055       COMPLEX            EYE
00056       PARAMETER          ( EYE = ( 0.0, 1.0 ) )
00057 *     ..
00058 *     .. Local Scalars ..
00059       INTEGER            I, J, N5
00060       REAL               ALPHA, ALPHA3, BETA
00061       COMPLEX            A, B, C, R
00062 *     ..
00063 *     .. External Functions ..
00064       COMPLEX            CLARND
00065       EXTERNAL           CLARND
00066 *     ..
00067 *     .. Intrinsic Functions ..
00068       INTRINSIC          ABS, SQRT
00069 *     ..
00070 *     .. Executable Statements ..
00071 *
00072 *     Initialize constants
00073 *
00074       ALPHA = ( 1.+SQRT( 17. ) ) / 8.
00075       BETA = ALPHA - 1. / 1000.
00076       ALPHA3 = ALPHA*ALPHA*ALPHA
00077 *
00078 *     UPLO = 'U':  Upper triangular storage
00079 *
00080       IF( UPLO.EQ.'U' ) THEN
00081 *
00082 *        Fill the upper triangle of the matrix with zeros.
00083 *
00084          DO 20 J = 1, N
00085             DO 10 I = 1, J
00086                X( I, J ) = 0.0
00087    10       CONTINUE
00088    20    CONTINUE
00089          N5 = N / 5
00090          N5 = N - 5*N5 + 1
00091 *
00092          DO 30 I = N, N5, -5
00093             A = ALPHA3*CLARND( 5, ISEED )
00094             B = CLARND( 5, ISEED ) / ALPHA
00095             C = A - 2.*B*EYE
00096             R = C / BETA
00097             X( I, I ) = A
00098             X( I-2, I ) = B
00099             X( I-2, I-1 ) = R
00100             X( I-2, I-2 ) = C
00101             X( I-1, I-1 ) = CLARND( 2, ISEED )
00102             X( I-3, I-3 ) = CLARND( 2, ISEED )
00103             X( I-4, I-4 ) = CLARND( 2, ISEED )
00104             IF( ABS( X( I-3, I-3 ) ).GT.ABS( X( I-4, I-4 ) ) ) THEN
00105                X( I-4, I-3 ) = 2.0*X( I-3, I-3 )
00106             ELSE
00107                X( I-4, I-3 ) = 2.0*X( I-4, I-4 )
00108             END IF
00109    30    CONTINUE
00110 *
00111 *        Clean-up for N not a multiple of 5.
00112 *
00113          I = N5 - 1
00114          IF( I.GT.2 ) THEN
00115             A = ALPHA3*CLARND( 5, ISEED )
00116             B = CLARND( 5, ISEED ) / ALPHA
00117             C = A - 2.*B*EYE
00118             R = C / BETA
00119             X( I, I ) = A
00120             X( I-2, I ) = B
00121             X( I-2, I-1 ) = R
00122             X( I-2, I-2 ) = C
00123             X( I-1, I-1 ) = CLARND( 2, ISEED )
00124             I = I - 3
00125          END IF
00126          IF( I.GT.1 ) THEN
00127             X( I, I ) = CLARND( 2, ISEED )
00128             X( I-1, I-1 ) = CLARND( 2, ISEED )
00129             IF( ABS( X( I, I ) ).GT.ABS( X( I-1, I-1 ) ) ) THEN
00130                X( I-1, I ) = 2.0*X( I, I )
00131             ELSE
00132                X( I-1, I ) = 2.0*X( I-1, I-1 )
00133             END IF
00134             I = I - 2
00135          ELSE IF( I.EQ.1 ) THEN
00136             X( I, I ) = CLARND( 2, ISEED )
00137             I = I - 1
00138          END IF
00139 *
00140 *     UPLO = 'L':  Lower triangular storage
00141 *
00142       ELSE
00143 *
00144 *        Fill the lower triangle of the matrix with zeros.
00145 *
00146          DO 50 J = 1, N
00147             DO 40 I = J, N
00148                X( I, J ) = 0.0
00149    40       CONTINUE
00150    50    CONTINUE
00151          N5 = N / 5
00152          N5 = N5*5
00153 *
00154          DO 60 I = 1, N5, 5
00155             A = ALPHA3*CLARND( 5, ISEED )
00156             B = CLARND( 5, ISEED ) / ALPHA
00157             C = A - 2.*B*EYE
00158             R = C / BETA
00159             X( I, I ) = A
00160             X( I+2, I ) = B
00161             X( I+2, I+1 ) = R
00162             X( I+2, I+2 ) = C
00163             X( I+1, I+1 ) = CLARND( 2, ISEED )
00164             X( I+3, I+3 ) = CLARND( 2, ISEED )
00165             X( I+4, I+4 ) = CLARND( 2, ISEED )
00166             IF( ABS( X( I+3, I+3 ) ).GT.ABS( X( I+4, I+4 ) ) ) THEN
00167                X( I+4, I+3 ) = 2.0*X( I+3, I+3 )
00168             ELSE
00169                X( I+4, I+3 ) = 2.0*X( I+4, I+4 )
00170             END IF
00171    60    CONTINUE
00172 *
00173 *        Clean-up for N not a multiple of 5.
00174 *
00175          I = N5 + 1
00176          IF( I.LT.N-1 ) THEN
00177             A = ALPHA3*CLARND( 5, ISEED )
00178             B = CLARND( 5, ISEED ) / ALPHA
00179             C = A - 2.*B*EYE
00180             R = C / BETA
00181             X( I, I ) = A
00182             X( I+2, I ) = B
00183             X( I+2, I+1 ) = R
00184             X( I+2, I+2 ) = C
00185             X( I+1, I+1 ) = CLARND( 2, ISEED )
00186             I = I + 3
00187          END IF
00188          IF( I.LT.N ) THEN
00189             X( I, I ) = CLARND( 2, ISEED )
00190             X( I+1, I+1 ) = CLARND( 2, ISEED )
00191             IF( ABS( X( I, I ) ).GT.ABS( X( I+1, I+1 ) ) ) THEN
00192                X( I+1, I ) = 2.0*X( I, I )
00193             ELSE
00194                X( I+1, I ) = 2.0*X( I+1, I+1 )
00195             END IF
00196             I = I + 2
00197          ELSE IF( I.EQ.N ) THEN
00198             X( I, I ) = CLARND( 2, ISEED )
00199             I = I + 1
00200          END IF
00201       END IF
00202 *
00203       RETURN
00204 *
00205 *     End of CLATSY
00206 *
00207       END