LAPACK: TESTING/EIG/zlatm4.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZLATM4( ITYPE, N, NZ1, NZ2, RSIGN, AMAGN, RCOND,
00002      $                   TRIANG, IDIST, ISEED, A, LDA )
00003 *
00004 *  -- LAPACK auxiliary test routine (version 3.1) --
00005 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       LOGICAL            RSIGN
00010       INTEGER            IDIST, ITYPE, LDA, N, NZ1, NZ2
00011       DOUBLE PRECISION   AMAGN, RCOND, TRIANG
00012 *     ..
00013 *     .. Array Arguments ..
00014       INTEGER            ISEED( 4 )
00015       COMPLEX*16         A( LDA, * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  ZLATM4 generates basic square matrices, which may later be
00022 *  multiplied by others in order to produce test matrices.  It is
00023 *  intended mainly to be used to test the generalized eigenvalue
00024 *  routines.
00025 *
00026 *  It first generates the diagonal and (possibly) subdiagonal,
00027 *  according to the value of ITYPE, NZ1, NZ2, RSIGN, AMAGN, and RCOND.
00028 *  It then fills in the upper triangle with random numbers, if TRIANG is
00029 *  non-zero.
00030 *
00031 *  Arguments
00032 *  =========
00033 *
00034 *  ITYPE   (input) INTEGER
00035 *          The "type" of matrix on the diagonal and sub-diagonal.
00036 *          If ITYPE < 0, then type abs(ITYPE) is generated and then
00037 *             swapped end for end (A(I,J) := A'(N-J,N-I).)  See also
00038 *             the description of AMAGN and RSIGN.
00039 *
00040 *          Special types:
00041 *          = 0:  the zero matrix.
00042 *          = 1:  the identity.
00043 *          = 2:  a transposed Jordan block.
00044 *          = 3:  If N is odd, then a k+1 x k+1 transposed Jordan block
00045 *                followed by a k x k identity block, where k=(N-1)/2.
00046 *                If N is even, then k=(N-2)/2, and a zero diagonal entry
00047 *                is tacked onto the end.
00048 *
00049 *          Diagonal types.  The diagonal consists of NZ1 zeros, then
00050 *             k=N-NZ1-NZ2 nonzeros.  The subdiagonal is zero.  ITYPE
00051 *             specifies the nonzero diagonal entries as follows:
00052 *          = 4:  1, ..., k
00053 *          = 5:  1, RCOND, ..., RCOND
00054 *          = 6:  1, ..., 1, RCOND
00055 *          = 7:  1, a, a^2, ..., a^(k-1)=RCOND
00056 *          = 8:  1, 1-d, 1-2*d, ..., 1-(k-1)*d=RCOND
00057 *          = 9:  random numbers chosen from (RCOND,1)
00058 *          = 10: random numbers with distribution IDIST (see ZLARND.)
00059 *
00060 *  N       (input) INTEGER
00061 *          The order of the matrix.
00062 *
00063 *  NZ1     (input) INTEGER
00064 *          If abs(ITYPE) > 3, then the first NZ1 diagonal entries will
00065 *          be zero.
00066 *
00067 *  NZ2     (input) INTEGER
00068 *          If abs(ITYPE) > 3, then the last NZ2 diagonal entries will
00069 *          be zero.
00070 *
00071 *  RSIGN   (input) LOGICAL
00072 *          = .TRUE.:  The diagonal and subdiagonal entries will be
00073 *                     multiplied by random numbers of magnitude 1.
00074 *          = .FALSE.: The diagonal and subdiagonal entries will be
00075 *                     left as they are (usually non-negative real.)
00076 *
00077 *  AMAGN   (input) DOUBLE PRECISION
00078 *          The diagonal and subdiagonal entries will be multiplied by
00079 *          AMAGN.
00080 *
00081 *  RCOND   (input) DOUBLE PRECISION
00082 *          If abs(ITYPE) > 4, then the smallest diagonal entry will be
00083 *          RCOND.  RCOND must be between 0 and 1.
00084 *
00085 *  TRIANG  (input) DOUBLE PRECISION
00086 *          The entries above the diagonal will be random numbers with
00087 *          magnitude bounded by TRIANG (i.e., random numbers multiplied
00088 *          by TRIANG.)
00089 *
00090 *  IDIST   (input) INTEGER
00091 *          On entry, DIST specifies the type of distribution to be used
00092 *          to generate a random matrix .
00093 *          = 1: real and imaginary parts each UNIFORM( 0, 1 )
00094 *          = 2: real and imaginary parts each UNIFORM( -1, 1 )
00095 *          = 3: real and imaginary parts each NORMAL( 0, 1 )
00096 *          = 4: complex number uniform in DISK( 0, 1 )
00097 *
00098 *  ISEED   (input/output) INTEGER array, dimension (4)
00099 *          On entry ISEED specifies the seed of the random number
00100 *          generator.  The values of ISEED are changed on exit, and can
00101 *          be used in the next call to ZLATM4 to continue the same
00102 *          random number sequence.
00103 *          Note: ISEED(4) should be odd, for the random number generator
00104 *          used at present.
00105 *
00106 *  A       (output) COMPLEX*16 array, dimension (LDA, N)
00107 *          Array to be computed.
00108 *
00109 *  LDA     (input) INTEGER
00110 *          Leading dimension of A.  Must be at least 1 and at least N.
00111 *
00112 *  =====================================================================
00113 *
00114 *     .. Parameters ..
00115       DOUBLE PRECISION   ZERO, ONE
00116       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00117       COMPLEX*16         CZERO, CONE
00118       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
00119      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
00120 *     ..
00121 *     .. Local Scalars ..
00122       INTEGER            I, ISDB, ISDE, JC, JD, JR, K, KBEG, KEND, KLEN
00123       DOUBLE PRECISION   ALPHA
00124       COMPLEX*16         CTEMP
00125 *     ..
00126 *     .. External Functions ..
00127       DOUBLE PRECISION   DLARAN
00128       COMPLEX*16         ZLARND
00129       EXTERNAL           DLARAN, ZLARND
00130 *     ..
00131 *     .. External Subroutines ..
00132       EXTERNAL           ZLASET
00133 *     ..
00134 *     .. Intrinsic Functions ..
00135       INTRINSIC          ABS, DBLE, DCMPLX, EXP, LOG, MAX, MIN, MOD
00136 *     ..
00137 *     .. Executable Statements ..
00138 *
00139       IF( N.LE.0 )
00140      $   RETURN
00141       CALL ZLASET( 'Full', N, N, CZERO, CZERO, A, LDA )
00142 *
00143 *     Insure a correct ISEED
00144 *
00145       IF( MOD( ISEED( 4 ), 2 ).NE.1 )
00146      $   ISEED( 4 ) = ISEED( 4 ) + 1
00147 *
00148 *     Compute diagonal and subdiagonal according to ITYPE, NZ1, NZ2,
00149 *     and RCOND
00150 *
00151       IF( ITYPE.NE.0 ) THEN
00152          IF( ABS( ITYPE ).GE.4 ) THEN
00153             KBEG = MAX( 1, MIN( N, NZ1+1 ) )
00154             KEND = MAX( KBEG, MIN( N, N-NZ2 ) )
00155             KLEN = KEND + 1 - KBEG
00156          ELSE
00157             KBEG = 1
00158             KEND = N
00159             KLEN = N
00160          END IF
00161          ISDB = 1
00162          ISDE = 0
00163          GO TO ( 10, 30, 50, 80, 100, 120, 140, 160,
00164      $           180, 200 )ABS( ITYPE )
00165 *
00166 *        abs(ITYPE) = 1: Identity
00167 *
00168    10    CONTINUE
00169          DO 20 JD = 1, N
00170             A( JD, JD ) = CONE
00171    20    CONTINUE
00172          GO TO 220
00173 *
00174 *        abs(ITYPE) = 2: Transposed Jordan block
00175 *
00176    30    CONTINUE
00177          DO 40 JD = 1, N - 1
00178             A( JD+1, JD ) = CONE
00179    40    CONTINUE
00180          ISDB = 1
00181          ISDE = N - 1
00182          GO TO 220
00183 *
00184 *        abs(ITYPE) = 3: Transposed Jordan block, followed by the
00185 *                        identity.
00186 *
00187    50    CONTINUE
00188          K = ( N-1 ) / 2
00189          DO 60 JD = 1, K
00190             A( JD+1, JD ) = CONE
00191    60    CONTINUE
00192          ISDB = 1
00193          ISDE = K
00194          DO 70 JD = K + 2, 2*K + 1
00195             A( JD, JD ) = CONE
00196    70    CONTINUE
00197          GO TO 220
00198 *
00199 *        abs(ITYPE) = 4: 1,...,k
00200 *
00201    80    CONTINUE
00202          DO 90 JD = KBEG, KEND
00203             A( JD, JD ) = DCMPLX( JD-NZ1 )
00204    90    CONTINUE
00205          GO TO 220
00206 *
00207 *        abs(ITYPE) = 5: One large D value:
00208 *
00209   100    CONTINUE
00210          DO 110 JD = KBEG + 1, KEND
00211             A( JD, JD ) = DCMPLX( RCOND )
00212   110    CONTINUE
00213          A( KBEG, KBEG ) = CONE
00214          GO TO 220
00215 *
00216 *        abs(ITYPE) = 6: One small D value:
00217 *
00218   120    CONTINUE
00219          DO 130 JD = KBEG, KEND - 1
00220             A( JD, JD ) = CONE
00221   130    CONTINUE
00222          A( KEND, KEND ) = DCMPLX( RCOND )
00223          GO TO 220
00224 *
00225 *        abs(ITYPE) = 7: Exponentially distributed D values:
00226 *
00227   140    CONTINUE
00228          A( KBEG, KBEG ) = CONE
00229          IF( KLEN.GT.1 ) THEN
00230             ALPHA = RCOND**( ONE / DBLE( KLEN-1 ) )
00231             DO 150 I = 2, KLEN
00232                A( NZ1+I, NZ1+I ) = DCMPLX( ALPHA**DBLE( I-1 ) )
00233   150       CONTINUE
00234          END IF
00235          GO TO 220
00236 *
00237 *        abs(ITYPE) = 8: Arithmetically distributed D values:
00238 *
00239   160    CONTINUE
00240          A( KBEG, KBEG ) = CONE
00241          IF( KLEN.GT.1 ) THEN
00242             ALPHA = ( ONE-RCOND ) / DBLE( KLEN-1 )
00243             DO 170 I = 2, KLEN
00244                A( NZ1+I, NZ1+I ) = DCMPLX( DBLE( KLEN-I )*ALPHA+RCOND )
00245   170       CONTINUE
00246          END IF
00247          GO TO 220
00248 *
00249 *        abs(ITYPE) = 9: Randomly distributed D values on ( RCOND, 1):
00250 *
00251   180    CONTINUE
00252          ALPHA = LOG( RCOND )
00253          DO 190 JD = KBEG, KEND
00254             A( JD, JD ) = EXP( ALPHA*DLARAN( ISEED ) )
00255   190    CONTINUE
00256          GO TO 220
00257 *
00258 *        abs(ITYPE) = 10: Randomly distributed D values from DIST
00259 *
00260   200    CONTINUE
00261          DO 210 JD = KBEG, KEND
00262             A( JD, JD ) = ZLARND( IDIST, ISEED )
00263   210    CONTINUE
00264 *
00265   220    CONTINUE
00266 *
00267 *        Scale by AMAGN
00268 *
00269          DO 230 JD = KBEG, KEND
00270             A( JD, JD ) = AMAGN*DBLE( A( JD, JD ) )
00271   230    CONTINUE
00272          DO 240 JD = ISDB, ISDE
00273             A( JD+1, JD ) = AMAGN*DBLE( A( JD+1, JD ) )
00274   240    CONTINUE
00275 *
00276 *        If RSIGN = .TRUE., assign random signs to diagonal and
00277 *        subdiagonal
00278 *
00279          IF( RSIGN ) THEN
00280             DO 250 JD = KBEG, KEND
00281                IF( DBLE( A( JD, JD ) ).NE.ZERO ) THEN
00282                   CTEMP = ZLARND( 3, ISEED )
00283                   CTEMP = CTEMP / ABS( CTEMP )
00284                   A( JD, JD ) = CTEMP*DBLE( A( JD, JD ) )
00285                END IF
00286   250       CONTINUE
00287             DO 260 JD = ISDB, ISDE
00288                IF( DBLE( A( JD+1, JD ) ).NE.ZERO ) THEN
00289                   CTEMP = ZLARND( 3, ISEED )
00290                   CTEMP = CTEMP / ABS( CTEMP )
00291                   A( JD+1, JD ) = CTEMP*DBLE( A( JD+1, JD ) )
00292                END IF
00293   260       CONTINUE
00294          END IF
00295 *
00296 *        Reverse if ITYPE < 0
00297 *
00298          IF( ITYPE.LT.0 ) THEN
00299             DO 270 JD = KBEG, ( KBEG+KEND-1 ) / 2
00300                CTEMP = A( JD, JD )
00301                A( JD, JD ) = A( KBEG+KEND-JD, KBEG+KEND-JD )
00302                A( KBEG+KEND-JD, KBEG+KEND-JD ) = CTEMP
00303   270       CONTINUE
00304             DO 280 JD = 1, ( N-1 ) / 2
00305                CTEMP = A( JD+1, JD )
00306                A( JD+1, JD ) = A( N+1-JD, N-JD )
00307                A( N+1-JD, N-JD ) = CTEMP
00308   280       CONTINUE
00309          END IF
00310 *
00311       END IF
00312 *
00313 *     Fill in upper triangle
00314 *
00315       IF( TRIANG.NE.ZERO ) THEN
00316          DO 300 JC = 2, N
00317             DO 290 JR = 1, JC - 1
00318                A( JR, JC ) = TRIANG*ZLARND( IDIST, ISEED )
00319   290       CONTINUE
00320   300    CONTINUE
00321       END IF
00322 *
00323       RETURN
00324 *
00325 *     End of ZLATM4
00326 *
00327       END