ScaLAPACK  2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
pclatms.f
Go to the documentation of this file.
00001 *
00002 *
00003       SUBROUTINE PCLATMS( M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX,
00004      $                    KL, KU, PACK, A, IA, JA, DESCA, ORDER, WORK,
00005      $                    LWORK, INFO )
00006 *
00007 *  -- ScaLAPACK routine (version 1.7) --
00008 *     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
00009 *     and University of California, Berkeley.
00010 *     May 1, 1997
00011 *
00012 *     .. Scalar Arguments ..
00013       CHARACTER          DIST, PACK, SYM
00014       INTEGER            IA, INFO, JA, KL, KU, LWORK, M, MODE, N, ORDER
00015       REAL               COND, DMAX
00016 *     ..
00017 *     .. Array Arguments ..
00018       INTEGER            DESCA( * ), ISEED( 4 )
00019       REAL               D( * )
00020       COMPLEX            A( * ), WORK( * )
00021 *     ..
00022 *
00023 *  Purpose
00024 *  =======
00025 *
00026 *     PCLATMS generates random Hermitian matrices with specified
00027 *     eigenvalues for testing SCALAPACK programs.
00028 *
00029 *     PCLATMS operates by applying the following sequence of
00030 *     operations:
00031 *
00032 *       Set the diagonal to D, where D may be input or
00033 *          computed according to MODE, COND, DMAX, and SYM
00034 *          as described below.
00035 *
00036 *           Generate a dense M x N matrix by multiplying D on the left
00037 *               and the right by random unitary matrices, then:
00038 *
00039 *           Reduce the bandwidth according to KL and KU, using
00040 *           Householder transformations.
00041 *           ### bandwidth reduction NOT SUPPORTED ###
00042 *
00043 *  Arguments
00044 *  =========
00045 *
00046 *  M      - (global input) INTEGER
00047 *           The number of rows of A. Not modified.
00048 *
00049 *  N      - (global input) INTEGER
00050 *           The number of columns of A. Not modified.
00051 *           ### M .ne. N unsupported
00052 *
00053 *  DIST   - (global input) CHARACTER*1
00054 *           On entry, DIST specifies the type of distribution to be used
00055 *           to generate the random eigen-/singular values.
00056 *           'U' => UNIFORM( 0, 1 )  ( 'U' for uniform )
00057 *           'S' => UNIFORM( -1, 1 ) ( 'S' for symmetric )
00058 *           'N' => NORMAL( 0, 1 )   ( 'N' for normal )
00059 *           Not modified.
00060 *
00061 *  ISEED  - (global input) INTEGER array, dimension ( 4 )
00062 *           On entry ISEED specifies the seed of the random number
00063 *           generator. They should lie between 0 and 4095 inclusive,
00064 *           and ISEED(4) should be odd. The random number generator
00065 *           uses a linear congruential sequence limited to small
00066 *           integers, and so should produce machine independent
00067 *           random numbers. The values of ISEED are changed on
00068 *           exit, and can be used in the next call to CLATMS
00069 *           to continue the same random number sequence.
00070 *           Changed on exit.
00071 *
00072 *  SYM    - (global input) CHARACTER*1
00073 *           If SYM='S' or 'H', the generated matrix is Hermitian, with
00074 *             eigenvalues specified by D, COND, MODE, and DMAX; they
00075 *             may be positive, negative, or zero.
00076 *           If SYM='P', the generated matrix is Hermitian, with
00077 *             eigenvalues (= singular values) specified by D, COND,
00078 *             MODE, and DMAX; they will not be negative.
00079 *           If SYM='N', the generated matrix is nonsymmetric, with
00080 *             singular values specified by D, COND, MODE, and DMAX;
00081 *             they will not be negative.
00082 *           ### SYM = 'N' NOT SUPPORTED ###
00083 *           Not modified.
00084 *
00085 *  D      - (local input/output) REAL array,
00086 *           dimension ( MIN( M , N ) )
00087 *           This array is used to specify the singular values or
00088 *           eigenvalues of A (see SYM, above.)  If MODE=0, then D is
00089 *           assumed to contain the singular/eigenvalues, otherwise
00090 *           they will be computed according to MODE, COND, and DMAX,
00091 *           and placed in D.
00092 *           Modified if MODE is nonzero.
00093 *
00094 *  MODE   - (global input) INTEGER
00095 *           On entry this describes how the singular/eigenvalues are to
00096 *           be specified:
00097 *           MODE = 0 means use D as input
00098 *           MODE = 1 sets D(1)=1 and D(2:N)=1.0/COND
00099 *           MODE = 2 sets D(1:N-1)=1 and D(N)=1.0/COND
00100 *           MODE = 3 sets D(I)=COND**(-(I-1)/(N-1))
00101 *           MODE = 4 sets D(i)=1 - (i-1)/(N-1)*(1 - 1/COND)
00102 *           MODE = 5 sets D to random numbers in the range
00103 *                    ( 1/COND , 1 ) such that their logarithms
00104 *                    are uniformly distributed.
00105 *           MODE = 6 set D to random numbers from same distribution
00106 *                    as the rest of the matrix.
00107 *           MODE < 0 has the same meaning as ABS(MODE), except that
00108 *              the order of the elements of D is reversed.
00109 *           Thus if MODE is positive, D has entries ranging from
00110 *              1 to 1/COND, if negative, from 1/COND to 1,
00111 *           If SYM='S' or 'H', and MODE is neither 0, 6, nor -6, then
00112 *              the elements of D will also be multiplied by a random
00113 *              sign (i.e., +1 or -1.)
00114 *           Not modified.
00115 *
00116 *  COND   - (global input) REAL
00117 *           On entry, this is used as described under MODE above.
00118 *           If used, it must be >= 1. Not modified.
00119 *
00120 *  DMAX   - (global input) REAL
00121 *           If MODE is neither -6, 0 nor 6, the contents of D, as
00122 *           computed according to MODE and COND, will be scaled by
00123 *           DMAX / max(abs(D(i))); thus, the maximum absolute eigen- or
00124 *           singular value (which is to say the norm) will be abs(DMAX).
00125 *           Note that DMAX need not be positive: if DMAX is negative
00126 *           (or zero), D will be scaled by a negative number (or zero).
00127 *           Not modified.
00128 *
00129 *  KL     - (global input) INTEGER
00130 *           This specifies the lower bandwidth of the  matrix. For
00131 *           example, KL=0 implies upper triangular, KL=1 implies upper
00132 *           Hessenberg, and KL being at least M-1 means that the matrix
00133 *           has full lower bandwidth.  KL must equal KU if the matrix
00134 *           is Hermitian.
00135 *           Not modified.
00136 *           ### 1 <= KL < N-1 is NOT SUPPORTED ###
00137 *
00138 *  KU     - (global input) INTEGER
00139 *           This specifies the upper bandwidth of the  matrix. For
00140 *           example, KU=0 implies lower triangular, KU=1 implies lower
00141 *           Hessenberg, and KU being at least N-1 means that the matrix
00142 *           has full upper bandwidth.  KL must equal KU if the matrix
00143 *           is Hermitian.
00144 *           Not modified.
00145 *           ### 1 <= KU < N-1 is NOT SUPPORTED ###
00146 *
00147 *  PACK   - (global input) CHARACTER*1
00148 *           This specifies packing of matrix as follows:
00149 *           'N' => no packing
00150 *           ### PACK must be 'N'  all other options NOT SUPPORTED ###
00151 *
00152 *  A      - (local output) COMPLEX array
00153 *           Global dimension (M, N), local dimension (MP, NQ)
00154 *           On exit A is the desired test matrix.
00155 *
00156 *  IA      (global input) INTEGER
00157 *          A's global row index, which points to the beginning of the
00158 *          submatrix which is to be operated on.
00159 *
00160 *  JA      (global input) INTEGER
00161 *          A's global column index, which points to the beginning of
00162 *          the submatrix which is to be operated on.
00163 *
00164 *  DESCA   (global and local input) INTEGER array of dimension DLEN_.
00165 *          The array descriptor for the distributed matrix A.
00166 *
00167 *  ORDER  - (input) INTEGER
00168 *           The number of reflectors used to define the orthogonal
00169 *           matrix Q.  A = Q * D * Q'
00170 *           Higher ORDER requires more computation and communication.
00171 *
00172 *  WORK   - (local input/output) COMPLEX array,
00173 *           dimension (LWORK)
00174 *
00175 *  LWORK  - (local input) INTEGER dimension of WORK
00176 *           LWORK >= SIZETMS as returned by PCLASIZESEP
00177 *
00178 *  INFO   - (global output) INTEGER
00179 *           Error code.  On exit, INFO will be set to one of the
00180 *           following values:
00181 *             0 => normal return
00182 *            -1 => M negative or unequal to N and SYM='S', 'H', or 'P'
00183 *            -2 => N negative
00184 *            -3 => DIST illegal string
00185 *            -5 => SYM illegal string
00186 *            -7 => MODE not in range -6 to 6
00187 *            -8 => COND less than 1.0, and MODE neither -6, 0 nor 6
00188 *           -10 => KL negative
00189 *           -11 => KU negative, or SYM='S' or 'H' and KU not equal to KL
00190 *           -16 => DESCA is inconsistent
00191 *           -17 => ORDER not in the range 0 to N inclusive
00192 *            1  => Error return from SLATM1
00193 *            2  => Cannot scale to DMAX (max. sing. value is 0)
00194 *            3  => Error return from PCLAGHE
00195 *
00196 *-----------------------------------------------------------------------
00197 *
00198 *
00199 *     .. Parameters ..
00200       INTEGER            BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_,
00201      $                   MB_, NB_, RSRC_, CSRC_, LLD_
00202       PARAMETER          ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1,
00203      $                   CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6,
00204      $                   RSRC_ = 7, CSRC_ = 8, LLD_ = 9 )
00205       REAL               ZERO, ONE
00206       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00207       COMPLEX            CZERO
00208       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ) )
00209 *     ..
00210 *     .. Local Scalars ..
00211       INTEGER            I, IDIST, IINFO, IPACK, IRSIGN, ISYM, LLB,
00212      $                   MNMIN, MYCOL, MYROW, NP, NPCOL, NPROW, NQ
00213       REAL               ALPHA, TEMP
00214 *     ..
00215 *     .. Local Arrays ..
00216       INTEGER            IDUM1( 1 ), IDUM2( 1 )
00217 *     ..
00218 *     .. External Functions ..
00219       LOGICAL            LSAME
00220       INTEGER            NUMROC
00221       EXTERNAL           LSAME, NUMROC
00222 *     ..
00223 *     .. External Subroutines ..
00224       EXTERNAL           BLACS_GRIDINFO, CHK1MAT, CLASET, PCHK1MAT,
00225      $                   PCLAGHE, PXERBLA, SLATM1, SSCAL
00226 *     ..
00227 *     .. Intrinsic Functions ..
00228       INTRINSIC          ABS, MAX, MIN, MOD
00229 *     ..
00230 *     .. Executable Statements ..
00231 *       This is just to keep ftnchek happy
00232       IF( BLOCK_CYCLIC_2D*CSRC_*CTXT_*DLEN_*DTYPE_*LLD_*MB_*M_*NB_*N_*
00233      $    RSRC_.LT.0 )RETURN
00234 *
00235 *     1)      Decode and Test the input parameters.
00236 *             Initialize flags & seed.
00237 *
00238 *
00239       INFO = 0
00240 *
00241       CALL BLACS_GRIDINFO( DESCA( CTXT_ ), NPROW, NPCOL, MYROW, MYCOL )
00242       IF( ( MYROW.GE.NPROW .OR. MYROW.LT.0 ) .OR.
00243      $    ( MYCOL.GE.NPCOL .OR. MYCOL.LT.0 ) )RETURN
00244 *
00245       NP = NUMROC( N, DESCA( MB_ ), MYROW, 0, NPROW )
00246       NQ = NUMROC( N, DESCA( NB_ ), MYCOL, 0, NPCOL )
00247 *
00248 *     Quick return if possible
00249 *
00250       IF( M.EQ.0 .OR. N.EQ.0 )
00251      $   RETURN
00252 *
00253 *     Decode DIST
00254 *
00255       IF( LSAME( DIST, 'U' ) ) THEN
00256          IDIST = 1
00257       ELSE IF( LSAME( DIST, 'S' ) ) THEN
00258          IDIST = 2
00259       ELSE IF( LSAME( DIST, 'N' ) ) THEN
00260          IDIST = 3
00261       ELSE
00262          IDIST = -1
00263       END IF
00264 *
00265 *     Decode SYM
00266 *
00267       IF( LSAME( SYM, 'N' ) ) THEN
00268          ISYM = 1
00269          IRSIGN = 0
00270       ELSE IF( LSAME( SYM, 'P' ) ) THEN
00271          ISYM = 2
00272          IRSIGN = 0
00273       ELSE IF( LSAME( SYM, 'S' ) ) THEN
00274          ISYM = 2
00275          IRSIGN = 1
00276       ELSE IF( LSAME( SYM, 'H' ) ) THEN
00277          ISYM = 2
00278          IRSIGN = 1
00279       ELSE
00280          ISYM = -1
00281       END IF
00282 *
00283 *     Decode PACK
00284 *
00285       IF( LSAME( PACK, 'N' ) ) THEN
00286          IPACK = 0
00287       ELSE
00288          IPACK = 1
00289       END IF
00290 *
00291 *     Set certain internal parameters
00292 *
00293       MNMIN = MIN( M, N )
00294       LLB = MIN( KL, M-1 )
00295 *
00296       IF( ORDER.EQ.0 )
00297      $   ORDER = N
00298 *
00299 *     Set INFO if an error
00300 *
00301       IF( NPROW.EQ.-1 ) THEN
00302          INFO = -( 1600+CTXT_ )
00303       ELSE
00304          CALL CHK1MAT( M, 1, N, 2, IA, JA, DESCA, 16, INFO )
00305          IF( INFO.EQ.0 ) THEN
00306             IF( M.NE.N .AND. ISYM.NE.1 ) THEN
00307                INFO = -2
00308             ELSE IF( IDIST.EQ.-1 ) THEN
00309                INFO = -3
00310             ELSE IF( ISYM.EQ.-1 ) THEN
00311                INFO = -5
00312             ELSE IF( ABS( MODE ).GT.6 ) THEN
00313                INFO = -7
00314             ELSE IF( ( MODE.NE.0 .AND. ABS( MODE ).NE.6 ) .AND. COND.LT.
00315      $               ONE ) THEN
00316                INFO = -8
00317             ELSE IF( KL.LT.0 ) THEN
00318                INFO = -10
00319             ELSE IF( KU.LT.0 .OR. ( ISYM.NE.1 .AND. KL.NE.KU ) ) THEN
00320                INFO = -11
00321             ELSE IF( ( ORDER.LT.0 ) .OR. ( ORDER.GT.N ) ) THEN
00322                INFO = -17
00323             END IF
00324          END IF
00325          CALL PCHK1MAT( M, 1, N, 2, IA, JA, DESCA, 16, 0, IDUM1, IDUM2,
00326      $                  INFO )
00327       END IF
00328 *
00329 *     Check for unsupported features
00330 *
00331       IF( ISYM.NE.2 ) THEN
00332          INFO = -5
00333       ELSE IF( IPACK.NE.0 ) THEN
00334          INFO = -12
00335       ELSE IF( KL.GT.0 .AND. KL.LT.M-1 ) THEN
00336          INFO = -10
00337       ELSE IF( KU.GT.0 .AND. KU.LT.N-1 ) THEN
00338          INFO = -11
00339       ELSE IF( LLB.NE.0 .AND. LLB.NE.M-1 ) THEN
00340          INFO = -10
00341       END IF
00342       IF( INFO.NE.0 ) THEN
00343          CALL PXERBLA( DESCA( CTXT_ ), 'PCLATMS', -INFO )
00344          RETURN
00345       END IF
00346 *
00347 *     Initialize random number generator
00348 *
00349       DO 10 I = 1, 4
00350          ISEED( I ) = MOD( ABS( ISEED( I ) ), 4096 )
00351    10 CONTINUE
00352 *
00353       IF( MOD( ISEED( 4 ), 2 ).NE.1 )
00354      $   ISEED( 4 ) = ISEED( 4 ) + 1
00355 *
00356 *     2)      Set up D  if indicated.
00357 *
00358 *             Compute D according to COND and MODE
00359 *
00360       CALL SLATM1( MODE, COND, IRSIGN, IDIST, ISEED, D, MNMIN, IINFO )
00361 *
00362       IF( IINFO.NE.0 ) THEN
00363          INFO = 1
00364          RETURN
00365       END IF
00366 *
00367 *
00368       IF( MODE.NE.0 .AND. ABS( MODE ).NE.6 ) THEN
00369 *
00370 *        Scale by DMAX
00371 *
00372          TEMP = ABS( D( 1 ) )
00373          DO 20 I = 2, MNMIN
00374             TEMP = MAX( TEMP, ABS( D( I ) ) )
00375    20    CONTINUE
00376 *
00377          IF( TEMP.GT.ZERO ) THEN
00378             ALPHA = DMAX / TEMP
00379          ELSE
00380             INFO = 2
00381             RETURN
00382          END IF
00383 *
00384          CALL SSCAL( MNMIN, ALPHA, D, 1 )
00385 *
00386       END IF
00387 *
00388       CALL CLASET( 'A', NP, NQ, CZERO, CZERO, A, DESCA( LLD_ ) )
00389 *
00390 *     Hermitian -- A = U D U'
00391 *
00392       CALL PCLAGHE( M, LLB, D, A, IA, JA, DESCA, ISEED, ORDER, WORK,
00393      $              LWORK, IINFO )
00394 *
00395       RETURN
00396 *
00397 *     End of PCLATMS
00398 *
00399       END