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