LAPACK: TESTING/LIN/clatsp.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CLATSP( UPLO, N, X, 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            N
00010 *     ..
00011 *     .. Array Arguments ..
00012       INTEGER            ISEED( * )
00013       COMPLEX            X( * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  CLATSP generates a special test matrix for the complex symmetric
00020 *  (indefinite) factorization for packed matrices.  The pivot blocks of
00021 *  the generated matrix 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 (N*(N+1)/2)
00041 *          The generated matrix in packed storage format.  The matrix
00042 *          consists of 3x3 and 2x2 diagonal blocks which result in the
00043 *          pivot sequence given above.  The matrix outside these
00044 *          diagonal blocks is zero.
00045 *
00046 *  ISEED   (input/output) INTEGER array, dimension (4)
00047 *          On entry, the seed for the random number generator.  The last
00048 *          of the four integers must be odd.  (modified on exit)
00049 *
00050 *  =====================================================================
00051 *
00052 *     .. Parameters ..
00053       COMPLEX            EYE
00054       PARAMETER          ( EYE = ( 0.0, 1.0 ) )
00055 *     ..
00056 *     .. Local Scalars ..
00057       INTEGER            J, JJ, N5
00058       REAL               ALPHA, ALPHA3, BETA
00059       COMPLEX            A, B, C, R
00060 *     ..
00061 *     .. External Functions ..
00062       COMPLEX            CLARND
00063       EXTERNAL           CLARND
00064 *     ..
00065 *     .. Intrinsic Functions ..
00066       INTRINSIC          ABS, SQRT
00067 *     ..
00068 *     .. Executable Statements ..
00069 *
00070 *     Initialize constants
00071 *
00072       ALPHA = ( 1.+SQRT( 17. ) ) / 8.
00073       BETA = ALPHA - 1. / 1000.
00074       ALPHA3 = ALPHA*ALPHA*ALPHA
00075 *
00076 *     Fill the matrix with zeros.
00077 *
00078       DO 10 J = 1, N*( N+1 ) / 2
00079          X( J ) = 0.0
00080    10 CONTINUE
00081 *
00082 *     UPLO = 'U':  Upper triangular storage
00083 *
00084       IF( UPLO.EQ.'U' ) THEN
00085          N5 = N / 5
00086          N5 = N - 5*N5 + 1
00087 *
00088          JJ = N*( N+1 ) / 2
00089          DO 20 J = N, N5, -5
00090             A = ALPHA3*CLARND( 5, ISEED )
00091             B = CLARND( 5, ISEED ) / ALPHA
00092             C = A - 2.*B*EYE
00093             R = C / BETA
00094             X( JJ ) = A
00095             X( JJ-2 ) = B
00096             JJ = JJ - J
00097             X( JJ ) = CLARND( 2, ISEED )
00098             X( JJ-1 ) = R
00099             JJ = JJ - ( J-1 )
00100             X( JJ ) = C
00101             JJ = JJ - ( J-2 )
00102             X( JJ ) = CLARND( 2, ISEED )
00103             JJ = JJ - ( J-3 )
00104             X( JJ ) = CLARND( 2, ISEED )
00105             IF( ABS( X( JJ+( J-3 ) ) ).GT.ABS( X( JJ ) ) ) THEN
00106                X( JJ+( J-4 ) ) = 2.0*X( JJ+( J-3 ) )
00107             ELSE
00108                X( JJ+( J-4 ) ) = 2.0*X( JJ )
00109             END IF
00110             JJ = JJ - ( J-4 )
00111    20    CONTINUE
00112 *
00113 *        Clean-up for N not a multiple of 5.
00114 *
00115          J = N5 - 1
00116          IF( J.GT.2 ) THEN
00117             A = ALPHA3*CLARND( 5, ISEED )
00118             B = CLARND( 5, ISEED ) / ALPHA
00119             C = A - 2.*B*EYE
00120             R = C / BETA
00121             X( JJ ) = A
00122             X( JJ-2 ) = B
00123             JJ = JJ - J
00124             X( JJ ) = CLARND( 2, ISEED )
00125             X( JJ-1 ) = R
00126             JJ = JJ - ( J-1 )
00127             X( JJ ) = C
00128             JJ = JJ - ( J-2 )
00129             J = J - 3
00130          END IF
00131          IF( J.GT.1 ) THEN
00132             X( JJ ) = CLARND( 2, ISEED )
00133             X( JJ-J ) = CLARND( 2, ISEED )
00134             IF( ABS( X( JJ ) ).GT.ABS( X( JJ-J ) ) ) THEN
00135                X( JJ-1 ) = 2.0*X( JJ )
00136             ELSE
00137                X( JJ-1 ) = 2.0*X( JJ-J )
00138             END IF
00139             JJ = JJ - J - ( J-1 )
00140             J = J - 2
00141          ELSE IF( J.EQ.1 ) THEN
00142             X( JJ ) = CLARND( 2, ISEED )
00143             J = J - 1
00144          END IF
00145 *
00146 *     UPLO = 'L':  Lower triangular storage
00147 *
00148       ELSE
00149          N5 = N / 5
00150          N5 = N5*5
00151 *
00152          JJ = 1
00153          DO 30 J = 1, N5, 5
00154             A = ALPHA3*CLARND( 5, ISEED )
00155             B = CLARND( 5, ISEED ) / ALPHA
00156             C = A - 2.*B*EYE
00157             R = C / BETA
00158             X( JJ ) = A
00159             X( JJ+2 ) = B
00160             JJ = JJ + ( N-J+1 )
00161             X( JJ ) = CLARND( 2, ISEED )
00162             X( JJ+1 ) = R
00163             JJ = JJ + ( N-J )
00164             X( JJ ) = C
00165             JJ = JJ + ( N-J-1 )
00166             X( JJ ) = CLARND( 2, ISEED )
00167             JJ = JJ + ( N-J-2 )
00168             X( JJ ) = CLARND( 2, ISEED )
00169             IF( ABS( X( JJ-( N-J-2 ) ) ).GT.ABS( X( JJ ) ) ) THEN
00170                X( JJ-( N-J-2 )+1 ) = 2.0*X( JJ-( N-J-2 ) )
00171             ELSE
00172                X( JJ-( N-J-2 )+1 ) = 2.0*X( JJ )
00173             END IF
00174             JJ = JJ + ( N-J-3 )
00175    30    CONTINUE
00176 *
00177 *        Clean-up for N not a multiple of 5.
00178 *
00179          J = N5 + 1
00180          IF( J.LT.N-1 ) THEN
00181             A = ALPHA3*CLARND( 5, ISEED )
00182             B = CLARND( 5, ISEED ) / ALPHA
00183             C = A - 2.*B*EYE
00184             R = C / BETA
00185             X( JJ ) = A
00186             X( JJ+2 ) = B
00187             JJ = JJ + ( N-J+1 )
00188             X( JJ ) = CLARND( 2, ISEED )
00189             X( JJ+1 ) = R
00190             JJ = JJ + ( N-J )
00191             X( JJ ) = C
00192             JJ = JJ + ( N-J-1 )
00193             J = J + 3
00194          END IF
00195          IF( J.LT.N ) THEN
00196             X( JJ ) = CLARND( 2, ISEED )
00197             X( JJ+( N-J+1 ) ) = CLARND( 2, ISEED )
00198             IF( ABS( X( JJ ) ).GT.ABS( X( JJ+( N-J+1 ) ) ) ) THEN
00199                X( JJ+1 ) = 2.0*X( JJ )
00200             ELSE
00201                X( JJ+1 ) = 2.0*X( JJ+( N-J+1 ) )
00202             END IF
00203             JJ = JJ + ( N-J+1 ) + ( N-J )
00204             J = J + 2
00205          ELSE IF( J.EQ.N ) THEN
00206             X( JJ ) = CLARND( 2, ISEED )
00207             JJ = JJ + ( N-J+1 )
00208             J = J + 1
00209          END IF
00210       END IF
00211 *
00212       RETURN
00213 *
00214 *     End of CLATSP
00215 *
00216       END