LAPACK: TESTING/EIG/dlatm4.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DLATM4( ITYPE, N, NZ1, NZ2, ISIGN, 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       INTEGER            IDIST, ISIGN, ITYPE, LDA, N, NZ1, NZ2
00010       DOUBLE PRECISION   AMAGN, RCOND, TRIANG
00011 *     ..
00012 *     .. Array Arguments ..
00013       INTEGER            ISEED( 4 )
00014       DOUBLE PRECISION   A( LDA, * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  DLATM4 generates basic square matrices, which may later be
00021 *  multiplied by others in order to produce test matrices.  It is
00022 *  intended mainly to be used to test the generalized eigenvalue
00023 *  routines.
00024 *
00025 *  It first generates the diagonal and (possibly) subdiagonal,
00026 *  according to the value of ITYPE, NZ1, NZ2, ISIGN, AMAGN, and RCOND.
00027 *  It then fills in the upper triangle with random numbers, if TRIANG is
00028 *  non-zero.
00029 *
00030 *  Arguments
00031 *  =========
00032 *
00033 *  ITYPE   (input) INTEGER
00034 *          The "type" of matrix on the diagonal and sub-diagonal.
00035 *          If ITYPE < 0, then type abs(ITYPE) is generated and then
00036 *             swapped end for end (A(I,J) := A'(N-J,N-I).)  See also
00037 *             the description of AMAGN and ISIGN.
00038 *
00039 *          Special types:
00040 *          = 0:  the zero matrix.
00041 *          = 1:  the identity.
00042 *          = 2:  a transposed Jordan block.
00043 *          = 3:  If N is odd, then a k+1 x k+1 transposed Jordan block
00044 *                followed by a k x k identity block, where k=(N-1)/2.
00045 *                If N is even, then k=(N-2)/2, and a zero diagonal entry
00046 *                is tacked onto the end.
00047 *
00048 *          Diagonal types.  The diagonal consists of NZ1 zeros, then
00049 *             k=N-NZ1-NZ2 nonzeros.  The subdiagonal is zero.  ITYPE
00050 *             specifies the nonzero diagonal entries as follows:
00051 *          = 4:  1, ..., k
00052 *          = 5:  1, RCOND, ..., RCOND
00053 *          = 6:  1, ..., 1, RCOND
00054 *          = 7:  1, a, a^2, ..., a^(k-1)=RCOND
00055 *          = 8:  1, 1-d, 1-2*d, ..., 1-(k-1)*d=RCOND
00056 *          = 9:  random numbers chosen from (RCOND,1)
00057 *          = 10: random numbers with distribution IDIST (see DLARND.)
00058 *
00059 *  N       (input) INTEGER
00060 *          The order of the matrix.
00061 *
00062 *  NZ1     (input) INTEGER
00063 *          If abs(ITYPE) > 3, then the first NZ1 diagonal entries will
00064 *          be zero.
00065 *
00066 *  NZ2     (input) INTEGER
00067 *          If abs(ITYPE) > 3, then the last NZ2 diagonal entries will
00068 *          be zero.
00069 *
00070 *  ISIGN   (input) INTEGER
00071 *          = 0: The sign of the diagonal and subdiagonal entries will
00072 *               be left unchanged.
00073 *          = 1: The diagonal and subdiagonal entries will have their
00074 *               sign changed at random.
00075 *          = 2: If ITYPE is 2 or 3, then the same as ISIGN=1.
00076 *               Otherwise, with probability 0.5, odd-even pairs of
00077 *               diagonal entries A(2*j-1,2*j-1), A(2*j,2*j) will be
00078 *               converted to a 2x2 block by pre- and post-multiplying
00079 *               by distinct random orthogonal rotations.  The remaining
00080 *               diagonal entries will have their sign changed at random.
00081 *
00082 *  AMAGN   (input) DOUBLE PRECISION
00083 *          The diagonal and subdiagonal entries will be multiplied by
00084 *          AMAGN.
00085 *
00086 *  RCOND   (input) DOUBLE PRECISION
00087 *          If abs(ITYPE) > 4, then the smallest diagonal entry will be
00088 *          entry will be RCOND.  RCOND must be between 0 and 1.
00089 *
00090 *  TRIANG  (input) DOUBLE PRECISION
00091 *          The entries above the diagonal will be random numbers with
00092 *          magnitude bounded by TRIANG (i.e., random numbers multiplied
00093 *          by TRIANG.)
00094 *
00095 *  IDIST   (input) INTEGER
00096 *          Specifies the type of distribution to be used to generate a
00097 *          random matrix.
00098 *          = 1:  UNIFORM( 0, 1 )
00099 *          = 2:  UNIFORM( -1, 1 )
00100 *          = 3:  NORMAL ( 0, 1 )
00101 *
00102 *  ISEED   (input/output) INTEGER array, dimension (4)
00103 *          On entry ISEED specifies the seed of the random number
00104 *          generator.  The values of ISEED are changed on exit, and can
00105 *          be used in the next call to DLATM4 to continue the same
00106 *          random number sequence.
00107 *          Note: ISEED(4) should be odd, for the random number generator
00108 *          used at present.
00109 *
00110 *  A       (output) DOUBLE PRECISION array, dimension (LDA, N)
00111 *          Array to be computed.
00112 *
00113 *  LDA     (input) INTEGER
00114 *          Leading dimension of A.  Must be at least 1 and at least N.
00115 *
00116 *  =====================================================================
00117 *
00118 *     .. Parameters ..
00119       DOUBLE PRECISION   ZERO, ONE, TWO
00120       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
00121       DOUBLE PRECISION   HALF
00122       PARAMETER          ( HALF = ONE / TWO )
00123 *     ..
00124 *     .. Local Scalars ..
00125       INTEGER            I, IOFF, ISDB, ISDE, JC, JD, JR, K, KBEG, KEND,
00126      $                   KLEN
00127       DOUBLE PRECISION   ALPHA, CL, CR, SAFMIN, SL, SR, SV1, SV2, TEMP
00128 *     ..
00129 *     .. External Functions ..
00130       DOUBLE PRECISION   DLAMCH, DLARAN, DLARND
00131       EXTERNAL           DLAMCH, DLARAN, DLARND
00132 *     ..
00133 *     .. External Subroutines ..
00134       EXTERNAL           DLASET
00135 *     ..
00136 *     .. Intrinsic Functions ..
00137       INTRINSIC          ABS, DBLE, EXP, LOG, MAX, MIN, MOD, SQRT
00138 *     ..
00139 *     .. Executable Statements ..
00140 *
00141       IF( N.LE.0 )
00142      $   RETURN
00143       CALL DLASET( 'Full', N, N, ZERO, ZERO, A, LDA )
00144 *
00145 *     Insure a correct ISEED
00146 *
00147       IF( MOD( ISEED( 4 ), 2 ).NE.1 )
00148      $   ISEED( 4 ) = ISEED( 4 ) + 1
00149 *
00150 *     Compute diagonal and subdiagonal according to ITYPE, NZ1, NZ2,
00151 *     and RCOND
00152 *
00153       IF( ITYPE.NE.0 ) THEN
00154          IF( ABS( ITYPE ).GE.4 ) THEN
00155             KBEG = MAX( 1, MIN( N, NZ1+1 ) )
00156             KEND = MAX( KBEG, MIN( N, N-NZ2 ) )
00157             KLEN = KEND + 1 - KBEG
00158          ELSE
00159             KBEG = 1
00160             KEND = N
00161             KLEN = N
00162          END IF
00163          ISDB = 1
00164          ISDE = 0
00165          GO TO ( 10, 30, 50, 80, 100, 120, 140, 160,
00166      $           180, 200 )ABS( ITYPE )
00167 *
00168 *        abs(ITYPE) = 1: Identity
00169 *
00170    10    CONTINUE
00171          DO 20 JD = 1, N
00172             A( JD, JD ) = ONE
00173    20    CONTINUE
00174          GO TO 220
00175 *
00176 *        abs(ITYPE) = 2: Transposed Jordan block
00177 *
00178    30    CONTINUE
00179          DO 40 JD = 1, N - 1
00180             A( JD+1, JD ) = ONE
00181    40    CONTINUE
00182          ISDB = 1
00183          ISDE = N - 1
00184          GO TO 220
00185 *
00186 *        abs(ITYPE) = 3: Transposed Jordan block, followed by the
00187 *                        identity.
00188 *
00189    50    CONTINUE
00190          K = ( N-1 ) / 2
00191          DO 60 JD = 1, K
00192             A( JD+1, JD ) = ONE
00193    60    CONTINUE
00194          ISDB = 1
00195          ISDE = K
00196          DO 70 JD = K + 2, 2*K + 1
00197             A( JD, JD ) = ONE
00198    70    CONTINUE
00199          GO TO 220
00200 *
00201 *        abs(ITYPE) = 4: 1,...,k
00202 *
00203    80    CONTINUE
00204          DO 90 JD = KBEG, KEND
00205             A( JD, JD ) = DBLE( JD-NZ1 )
00206    90    CONTINUE
00207          GO TO 220
00208 *
00209 *        abs(ITYPE) = 5: One large D value:
00210 *
00211   100    CONTINUE
00212          DO 110 JD = KBEG + 1, KEND
00213             A( JD, JD ) = RCOND
00214   110    CONTINUE
00215          A( KBEG, KBEG ) = ONE
00216          GO TO 220
00217 *
00218 *        abs(ITYPE) = 6: One small D value:
00219 *
00220   120    CONTINUE
00221          DO 130 JD = KBEG, KEND - 1
00222             A( JD, JD ) = ONE
00223   130    CONTINUE
00224          A( KEND, KEND ) = RCOND
00225          GO TO 220
00226 *
00227 *        abs(ITYPE) = 7: Exponentially distributed D values:
00228 *
00229   140    CONTINUE
00230          A( KBEG, KBEG ) = ONE
00231          IF( KLEN.GT.1 ) THEN
00232             ALPHA = RCOND**( ONE / DBLE( KLEN-1 ) )
00233             DO 150 I = 2, KLEN
00234                A( NZ1+I, NZ1+I ) = ALPHA**DBLE( I-1 )
00235   150       CONTINUE
00236          END IF
00237          GO TO 220
00238 *
00239 *        abs(ITYPE) = 8: Arithmetically distributed D values:
00240 *
00241   160    CONTINUE
00242          A( KBEG, KBEG ) = ONE
00243          IF( KLEN.GT.1 ) THEN
00244             ALPHA = ( ONE-RCOND ) / DBLE( KLEN-1 )
00245             DO 170 I = 2, KLEN
00246                A( NZ1+I, NZ1+I ) = DBLE( KLEN-I )*ALPHA + RCOND
00247   170       CONTINUE
00248          END IF
00249          GO TO 220
00250 *
00251 *        abs(ITYPE) = 9: Randomly distributed D values on ( RCOND, 1):
00252 *
00253   180    CONTINUE
00254          ALPHA = LOG( RCOND )
00255          DO 190 JD = KBEG, KEND
00256             A( JD, JD ) = EXP( ALPHA*DLARAN( ISEED ) )
00257   190    CONTINUE
00258          GO TO 220
00259 *
00260 *        abs(ITYPE) = 10: Randomly distributed D values from DIST
00261 *
00262   200    CONTINUE
00263          DO 210 JD = KBEG, KEND
00264             A( JD, JD ) = DLARND( IDIST, ISEED )
00265   210    CONTINUE
00266 *
00267   220    CONTINUE
00268 *
00269 *        Scale by AMAGN
00270 *
00271          DO 230 JD = KBEG, KEND
00272             A( JD, JD ) = AMAGN*DBLE( A( JD, JD ) )
00273   230    CONTINUE
00274          DO 240 JD = ISDB, ISDE
00275             A( JD+1, JD ) = AMAGN*DBLE( A( JD+1, JD ) )
00276   240    CONTINUE
00277 *
00278 *        If ISIGN = 1 or 2, assign random signs to diagonal and
00279 *        subdiagonal
00280 *
00281          IF( ISIGN.GT.0 ) THEN
00282             DO 250 JD = KBEG, KEND
00283                IF( DBLE( A( JD, JD ) ).NE.ZERO ) THEN
00284                   IF( DLARAN( ISEED ).GT.HALF )
00285      $               A( JD, JD ) = -A( JD, JD )
00286                END IF
00287   250       CONTINUE
00288             DO 260 JD = ISDB, ISDE
00289                IF( DBLE( A( JD+1, JD ) ).NE.ZERO ) THEN
00290                   IF( DLARAN( ISEED ).GT.HALF )
00291      $               A( JD+1, JD ) = -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                TEMP = A( JD, JD )
00301                A( JD, JD ) = A( KBEG+KEND-JD, KBEG+KEND-JD )
00302                A( KBEG+KEND-JD, KBEG+KEND-JD ) = TEMP
00303   270       CONTINUE
00304             DO 280 JD = 1, ( N-1 ) / 2
00305                TEMP = A( JD+1, JD )
00306                A( JD+1, JD ) = A( N+1-JD, N-JD )
00307                A( N+1-JD, N-JD ) = TEMP
00308   280       CONTINUE
00309          END IF
00310 *
00311 *        If ISIGN = 2, and no subdiagonals already, then apply
00312 *        random rotations to make 2x2 blocks.
00313 *
00314          IF( ISIGN.EQ.2 .AND. ITYPE.NE.2 .AND. ITYPE.NE.3 ) THEN
00315             SAFMIN = DLAMCH( 'S' )
00316             DO 290 JD = KBEG, KEND - 1, 2
00317                IF( DLARAN( ISEED ).GT.HALF ) THEN
00318 *
00319 *                 Rotation on left.
00320 *
00321                   CL = TWO*DLARAN( ISEED ) - ONE
00322                   SL = TWO*DLARAN( ISEED ) - ONE
00323                   TEMP = ONE / MAX( SAFMIN, SQRT( CL**2+SL**2 ) )
00324                   CL = CL*TEMP
00325                   SL = SL*TEMP
00326 *
00327 *                 Rotation on right.
00328 *
00329                   CR = TWO*DLARAN( ISEED ) - ONE
00330                   SR = TWO*DLARAN( ISEED ) - ONE
00331                   TEMP = ONE / MAX( SAFMIN, SQRT( CR**2+SR**2 ) )
00332                   CR = CR*TEMP
00333                   SR = SR*TEMP
00334 *
00335 *                 Apply
00336 *
00337                   SV1 = A( JD, JD )
00338                   SV2 = A( JD+1, JD+1 )
00339                   A( JD, JD ) = CL*CR*SV1 + SL*SR*SV2
00340                   A( JD+1, JD ) = -SL*CR*SV1 + CL*SR*SV2
00341                   A( JD, JD+1 ) = -CL*SR*SV1 + SL*CR*SV2
00342                   A( JD+1, JD+1 ) = SL*SR*SV1 + CL*CR*SV2
00343                END IF
00344   290       CONTINUE
00345          END IF
00346 *
00347       END IF
00348 *
00349 *     Fill in upper triangle (except for 2x2 blocks)
00350 *
00351       IF( TRIANG.NE.ZERO ) THEN
00352          IF( ISIGN.NE.2 .OR. ITYPE.EQ.2 .OR. ITYPE.EQ.3 ) THEN
00353             IOFF = 1
00354          ELSE
00355             IOFF = 2
00356             DO 300 JR = 1, N - 1
00357                IF( A( JR+1, JR ).EQ.ZERO )
00358      $            A( JR, JR+1 ) = TRIANG*DLARND( IDIST, ISEED )
00359   300       CONTINUE
00360          END IF
00361 *
00362          DO 320 JC = 2, N
00363             DO 310 JR = 1, JC - IOFF
00364                A( JR, JC ) = TRIANG*DLARND( IDIST, ISEED )
00365   310       CONTINUE
00366   320    CONTINUE
00367       END IF
00368 *
00369       RETURN
00370 *
00371 *     End of DLATM4
00372 *
00373       END