|
| 1 | +* |
| 2 | +* |
| 3 | + PROGRAM SAMPLE_PDSYEV_CALL |
| 4 | +* |
| 5 | +* |
| 6 | +* -- ScaLAPACK routine (version 1.2) -- |
| 7 | +* University of Tennessee, Knoxville, Oak Ridge National Laboratory, |
| 8 | +* and University of California, Berkeley. |
| 9 | +* May 10, 1996 |
| 10 | +* |
| 11 | +* This routine contains a sample call to PDSYEV. |
| 12 | +* When compiled and run, it produces output which can be |
| 13 | +* pasted directly into matlab. |
| 14 | +* |
| 15 | +* .. Parameters .. |
| 16 | + INTEGER LWORK, MAXN |
| 17 | + DOUBLE PRECISION ZERO |
| 18 | + PARAMETER ( LWORK = 264, MAXN = 100, ZERO = 0.0D+0 ) |
| 19 | + INTEGER LDA |
| 20 | + DOUBLE PRECISION MONE |
| 21 | + INTEGER MAXPROCS |
| 22 | + PARAMETER ( LDA = MAXN, MONE = -1.0D+0, MAXPROCS = 512 ) |
| 23 | +* .. |
| 24 | +* .. Local Scalars .. |
| 25 | + INTEGER CONTEXT, I, IAM, INFO, MYCOL, MYROW, N, NB, |
| 26 | + $ NPCOL, NPROCS, NPROW |
| 27 | +* .. |
| 28 | +* .. Local Arrays .. |
| 29 | + INTEGER DESCA( 50 ), DESCZ( 50 ) |
| 30 | + DOUBLE PRECISION A( LDA, LDA ), W( MAXN ), |
| 31 | + $ WORK( LWORK ), Z( LDA, LDA ) |
| 32 | +* .. |
| 33 | +* .. External Subroutines .. |
| 34 | + EXTERNAL BLACS_EXIT, BLACS_GET, BLACS_GRIDEXIT, |
| 35 | + $ BLACS_GRIDINFO, BLACS_GRIDINIT, BLACS_PINFO, |
| 36 | + $ BLACS_SETUP, DESCINIT, PDLAMODHILB, PDLAPRNT, |
| 37 | + $ PDSYEV |
| 38 | +* .. |
| 39 | +* .. Executable Statements .. |
| 40 | +* |
| 41 | +* |
| 42 | +* Set up the problem |
| 43 | +* |
| 44 | + N = 4 |
| 45 | + NB = 1 |
| 46 | + NPROW = 2 |
| 47 | + NPCOL = 2 |
| 48 | +* |
| 49 | +* |
| 50 | +* Initialize the BLACS |
| 51 | +* |
| 52 | + CALL BLACS_PINFO( IAM, NPROCS ) |
| 53 | + IF( ( NPROCS.LT.1 ) ) THEN |
| 54 | + CALL BLACS_SETUP( IAM, NPROW*NPCOL ) |
| 55 | + END IF |
| 56 | +* |
| 57 | +* |
| 58 | +* Initialize a single BLACS context |
| 59 | +* |
| 60 | + CALL BLACS_GET( -1, 0, CONTEXT ) |
| 61 | + CALL BLACS_GRIDINIT( CONTEXT, 'R', NPROW, NPCOL ) |
| 62 | + CALL BLACS_GRIDINFO( CONTEXT, NPROW, NPCOL, MYROW, MYCOL ) |
| 63 | +* |
| 64 | +* Bail out if this process is not a part of this context. |
| 65 | +* |
| 66 | + IF( MYROW.EQ.-1 ) |
| 67 | + $ GO TO 20 |
| 68 | +* |
| 69 | +* |
| 70 | +* These are basic array descriptors |
| 71 | +* |
| 72 | + CALL DESCINIT( DESCA, N, N, NB, NB, 0, 0, CONTEXT, LDA, INFO ) |
| 73 | + CALL DESCINIT( DESCZ, N, N, NB, NB, 0, 0, CONTEXT, LDA, INFO ) |
| 74 | +* |
| 75 | +* Build a matrix that you can create with |
| 76 | +* a one line matlab command: hilb(n) + diag([1:-1/n:1/n]) |
| 77 | +* |
| 78 | + CALL PDLAMODHILB( N, A, 1, 1, DESCA, INFO ) |
| 79 | +* |
| 80 | +* Ask PDSYEV to compute the entire eigendecomposition |
| 81 | +* |
| 82 | + CALL PDSYEV( 'V', 'U', N, A, 1, 1, DESCA, W, Z, 1, 1, |
| 83 | + $ DESCZ, WORK, LWORK, INFO ) |
| 84 | +* |
| 85 | +* Print out the eigenvectors |
| 86 | +* |
| 87 | + CALL PDLAPRNT( N, N, Z, 1, 1, DESCZ, 0, 0, 'Z', 6, WORK ) |
| 88 | +* |
| 89 | +* Print out matlab code which will check the residual |
| 90 | +* |
| 91 | + IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN |
| 92 | + PRINT *, ' N =', N |
| 93 | + PRINT *, ' A = hilb(N) + diag([1:-1/N:1/N])' |
| 94 | + DO 10 I = 1, N |
| 95 | + PRINT *, ' W(', I, ')=', W( I ), ';' |
| 96 | + 10 CONTINUE |
| 97 | + PRINT *, ' backerror = A - Z * diag(W) * Z'' ' |
| 98 | + PRINT *, ' resid = A * Z - Z * diag(W)' |
| 99 | + PRINT *, ' ortho = Z'' * Z - eye(N)' |
| 100 | + PRINT *, ' norm(backerror)' |
| 101 | + PRINT *, ' norm(resid)' |
| 102 | + PRINT *, ' norm(ortho)' |
| 103 | + END IF |
| 104 | +* |
| 105 | + CALL BLACS_GRIDEXIT( CONTEXT ) |
| 106 | +* |
| 107 | + 20 CONTINUE |
| 108 | +* |
| 109 | + CALL BLACS_EXIT( 0 ) |
| 110 | +* |
| 111 | +* |
| 112 | +* Uncomment this line on SUN systems to avoid the useless print out |
| 113 | +* |
| 114 | +* CALL IEEE_FLAGS( 'clear', 'exception', 'underflow', '') |
| 115 | +* |
| 116 | +* |
| 117 | + 9999 FORMAT( 'W=diag([', 4D16.12, ']);' ) |
| 118 | +* |
| 119 | + STOP |
| 120 | + END |
| 121 | +* |
| 122 | + SUBROUTINE PDLAMODHILB( N, A, IA, JA, DESCA, INFO ) |
| 123 | +* |
| 124 | +* -- ScaLAPACK routine (version 1.2) -- |
| 125 | +* University of Tennessee, Knoxville, Oak Ridge National Laboratory, |
| 126 | +* and University of California, Berkeley. |
| 127 | +* May 10, 1996 |
| 128 | +* |
| 129 | +* |
| 130 | +* |
| 131 | +* |
| 132 | +* .. Parameters .. |
| 133 | + INTEGER BLOCK_CYCLIC_2D, DLEN_, DT_, CTXT_, M_, N_, |
| 134 | + $ MB_, NB_, RSRC_, CSRC_, LLD_ |
| 135 | + PARAMETER ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DT_ = 1, |
| 136 | + $ CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6, |
| 137 | + $ RSRC_ = 7, CSRC_ = 8, LLD_ = 9 ) |
| 138 | + DOUBLE PRECISION ONE |
| 139 | + PARAMETER ( ONE = 1.0D+0 ) |
| 140 | +* .. |
| 141 | +* .. Scalar Arguments .. |
| 142 | + INTEGER IA, INFO, JA, N |
| 143 | +* .. |
| 144 | +* .. Array Arguments .. |
| 145 | + INTEGER DESCA( * ) |
| 146 | + DOUBLE PRECISION A( * ) |
| 147 | +* .. |
| 148 | +* .. Local Scalars .. |
| 149 | + INTEGER I, J, MYCOL, MYROW, NPCOL, NPROW |
| 150 | +* .. |
| 151 | +* .. External Subroutines .. |
| 152 | + EXTERNAL BLACS_GRIDINFO, PDELSET |
| 153 | +* .. |
| 154 | +* .. Intrinsic Functions .. |
| 155 | + INTRINSIC DBLE |
| 156 | +* .. |
| 157 | +* .. Executable Statements .. |
| 158 | +* |
| 159 | +* |
| 160 | +* The matlab code for a real matrix is: |
| 161 | +* hilb(n) + diag( [ 1:-1/n:1/n ] ) |
| 162 | +* The matlab code for a complex matrix is: |
| 163 | +* hilb(N) + toeplitz( [ 1 (1:(N-1))*i ] ) |
| 164 | +* |
| 165 | +* This is just to keep ftnchek happy |
| 166 | + IF( BLOCK_CYCLIC_2D*CSRC_*CTXT_*DLEN_*DT_*LLD_*MB_*M_*NB_*N_* |
| 167 | + $ RSRC_.LT.0 )RETURN |
| 168 | +* |
| 169 | + INFO = 0 |
| 170 | +* |
| 171 | + CALL BLACS_GRIDINFO( DESCA( CTXT_ ), NPROW, NPCOL, MYROW, MYCOL ) |
| 172 | +* |
| 173 | +* |
| 174 | + IF( IA.NE.1 ) THEN |
| 175 | + INFO = -3 |
| 176 | + ELSE IF( JA.NE.1 ) THEN |
| 177 | + INFO = -4 |
| 178 | + END IF |
| 179 | +* |
| 180 | + DO 20 J = 1, N |
| 181 | + DO 10 I = 1, N |
| 182 | + IF( I.EQ.J ) THEN |
| 183 | + CALL PDELSET( A, I, J, DESCA, |
| 184 | + $ ( DBLE( N-I+1 ) ) / DBLE( N )+ONE / |
| 185 | + $ ( DBLE( I+J )-ONE ) ) |
| 186 | + ELSE |
| 187 | + CALL PDELSET( A, I, J, DESCA, ONE / ( DBLE( I+J )-ONE ) ) |
| 188 | + END IF |
| 189 | + 10 CONTINUE |
| 190 | + 20 CONTINUE |
| 191 | +* |
| 192 | +* |
| 193 | + RETURN |
| 194 | +* |
| 195 | +* End of PDLAMODHLIB |
| 196 | +* |
| 197 | + END |
0 commit comments