ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pcget22.f
Go to the documentation of this file.
1  SUBROUTINE pcget22( TRANSA, TRANSE, TRANSW, N, A, DESCA, E, DESCE,
2  $ W, WORK, DESCW, RWORK, RESULT )
3 *
4 * -- ScaLAPACK testing routine (version 1.7) --
5 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
6 * and University of California, Berkeley.
7 * August 14, 2001
8 *
9 * .. Scalar Arguments ..
10  CHARACTER TRANSA, TRANSE, TRANSW
11  INTEGER N
12 * ..
13 * .. Array Arguments ..
14  INTEGER DESCA( * ), DESCE( * ), DESCW( * )
15  REAL RESULT( 2 ), RWORK( * )
16  COMPLEX A( * ), E( * ), W( * ), WORK( * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * PCGET22 does an eigenvector check.
23 *
24 * The basic test is:
25 *
26 * RESULT(1) = | A E - E W | / ( |A| |E| ulp )
27 *
28 * using the 1-norm. It also tests the normalization of E:
29 *
30 * RESULT(2) = max | m-norm(E(j)) - 1 | / ( n ulp )
31 * j
32 *
33 * where E(j) is the j-th eigenvector, and m-norm is the max-norm of a
34 * vector. The max-norm of a complex n-vector x in this case is the
35 * maximum of |re(x(i)| + |im(x(i)| over i = 1, ..., n.
36 *
37 * Arguments
38 * ==========
39 *
40 * TRANSA (input) CHARACTER*1
41 * Specifies whether or not A is transposed.
42 * = 'N': No transpose
43 * = 'T': Transpose
44 * = 'C': Conjugate transpose
45 *
46 * TRANSE (input) CHARACTER*1
47 * Specifies whether or not E is transposed.
48 * = 'N': No transpose, eigenvectors are in columns of E
49 * = 'T': Transpose, eigenvectors are in rows of E
50 * = 'C': Conjugate transpose, eigenvectors are in rows of E
51 *
52 * TRANSW (input) CHARACTER*1
53 * Specifies whether or not W is transposed.
54 * = 'N': No transpose
55 * = 'T': Transpose, same as TRANSW = 'N'
56 * = 'C': Conjugate transpose, use -WI(j) instead of WI(j)
57 *
58 * N (input) INTEGER
59 * The order of the matrix A. N >= 0.
60 *
61 * A (input) COMPLEX array, dimension (*)
62 * The matrix whose eigenvectors are in E.
63 *
64 * DESCA (input) INTEGER array, dimension(*)
65 *
66 * E (input) COMPLEX array, dimension (*)
67 * The matrix of eigenvectors. If TRANSE = 'N', the eigenvectors
68 * are stored in the columns of E, if TRANSE = 'T' or 'C', the
69 * eigenvectors are stored in the rows of E.
70 *
71 * DESCE (input) INTEGER array, dimension(*)
72 *
73 * W (input) COMPLEX array, dimension (N)
74 * The eigenvalues of A.
75 *
76 * WORK (workspace) COMPLEX array, dimension (*)
77 * DESCW (input) INTEGER array, dimension(*)
78 *
79 * RWORK (workspace) REAL array, dimension (N)
80 *
81 * RESULT (output) REAL array, dimension (2)
82 * RESULT(1) = | A E - E W | / ( |A| |E| ulp )
83 * RESULT(2) = max | m-norm(E(j)) - 1 | / ( n ulp )
84 * j
85 * Further Details
86 * ===============
87 *
88 * Contributed by Mark Fahey, June, 2000
89 *
90 * =====================================================================
91 *
92 * .. Parameters ..
93  INTEGER BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_,
94  $ mb_, nb_, rsrc_, csrc_, lld_
95  parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
96  $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
97  $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
98  REAL ZERO, ONE
99  parameter( zero = 0.0e+0, one = 1.0e+0 )
100  COMPLEX CZERO, CONE
101  parameter( czero = ( 0.0e+0, 0.0e+0 ),
102  $ cone = ( 1.0e+0, 0.0e+0 ) )
103 * ..
104 * .. Local Scalars ..
105  CHARACTER NORMA, NORME
106  INTEGER ICOL, II, IROW, ITRNSE, ITRNSW, J, JCOL, JJ,
107  $ jrow, jvec, lda, lde, ldw, mb, mycol, myrow,
108  $ nb, npcol, nprow, contxt, ra, ca, rsrc, csrc
109  REAL ANORM, ENORM, ENRMAX, ENRMIN, ERRNRM, TEMP1,
110  $ ulp, unfl
111  COMPLEX CDUM, WTEMP
112 * ..
113 * .. External Functions ..
114  LOGICAL LSAME
115  REAL PSLAMCH, PCLANGE
116  EXTERNAL lsame, pslamch, pclange
117 * ..
118 * .. External Subroutines ..
119  EXTERNAL blacs_gridinfo, sgamn2d, sgamx2d, infog2l,
120  $ pcaxpy, pcgemm, pclaset
121 * ..
122 * .. Intrinsic Functions ..
123  INTRINSIC abs, real, conjg, aimag, max, min
124 * ..
125 * .. Statement Functions ..
126  REAL CABS1
127 * ..
128 * .. Statement Function definitions ..
129  cabs1( cdum ) = abs( real( cdum ) ) + abs( aimag( cdum ) )
130 * ..
131 * .. Executable Statements ..
132 *
133 * Initialize RESULT (in case N=0)
134 *
135  result( 1 ) = zero
136  result( 2 ) = zero
137  IF( n.LE.0 )
138  $ RETURN
139 *
140  contxt = desca( ctxt_ )
141  rsrc = desca( rsrc_ )
142  csrc = desca( csrc_ )
143  nb = desca( nb_ )
144  mb = desca( mb_ )
145  lda = desca( lld_ )
146  lde = desce( lld_ )
147  ldw = descw( lld_ )
148  CALL blacs_gridinfo( contxt, nprow, npcol, myrow, mycol )
149 *
150  unfl = pslamch( contxt, 'Safe minimum' )
151  ulp = pslamch( contxt, 'Precision' )
152 *
153  itrnse = 0
154  itrnsw = 0
155  norma = 'O'
156  norme = 'O'
157 *
158  IF( lsame( transa, 'T' ) .OR. lsame( transa, 'C' ) ) THEN
159  norma = 'I'
160  END IF
161 *
162  IF( lsame( transe, 'T' ) ) THEN
163  itrnse = 1
164  norme = 'I'
165  ELSE IF( lsame( transe, 'C' ) ) THEN
166  itrnse = 2
167  norme = 'I'
168  END IF
169 *
170  IF( lsame( transw, 'C' ) ) THEN
171  itrnsw = 1
172  END IF
173 *
174 * Normalization of E:
175 *
176  enrmin = one / ulp
177  enrmax = zero
178  IF( itrnse.EQ.0 ) THEN
179  DO 20 jvec = 1, n
180  temp1 = zero
181  DO 10 j = 1, n
182  CALL infog2l( j, jvec, desce, nprow, npcol, myrow, mycol,
183  $ irow, icol, ii, jj )
184  IF( ( myrow.EQ.ii ) .AND. ( mycol.EQ.jj ) ) THEN
185  temp1 = max( temp1, cabs1( e( ( icol-1 )*lde+
186  $ irow ) ) )
187  END IF
188  10 CONTINUE
189  IF( mycol.EQ.jj ) THEN
190  CALL sgamx2d( contxt, 'Col', ' ', 1, 1, temp1, 1, ra, ca,
191  $ -1, -1, -1 )
192  enrmin = min( enrmin, temp1 )
193  enrmax = max( enrmax, temp1 )
194  END IF
195  20 CONTINUE
196  CALL sgamx2d( contxt, 'Row', ' ', 1, 1, enrmax, 1, ra, ca, -1,
197  $ -1, -1 )
198  CALL sgamn2d( contxt, 'Row', ' ', 1, 1, enrmin, 1, ra, ca, -1,
199  $ -1, -1 )
200  ELSE
201  DO 40 j = 1, n
202  temp1 = zero
203  DO 30 jvec = 1, n
204  CALL infog2l( j, jvec, desce, nprow, npcol, myrow, mycol,
205  $ irow, icol, ii, jj )
206  IF( ( myrow.EQ.ii ) .AND. ( mycol.EQ.jj ) ) THEN
207  temp1 = max( temp1, cabs1( e( ( icol-1 )*lde+
208  $ irow ) ) )
209  END IF
210  30 CONTINUE
211  IF( myrow.EQ.ii ) THEN
212  CALL sgamx2d( contxt, 'Row', ' ', 1, 1, temp1, 1, ra, ca,
213  $ -1, -1, -1 )
214  enrmin = min( enrmin, temp1 )
215  enrmax = max( enrmax, temp1 )
216  END IF
217  40 CONTINUE
218  CALL sgamx2d( contxt, 'Row', ' ', 1, 1, enrmax, 1, ra, ca, -1,
219  $ -1, -1 )
220  CALL sgamn2d( contxt, 'Row', ' ', 1, 1, enrmin, 1, ra, ca, -1,
221  $ -1, -1 )
222  END IF
223 *
224 * Norm of A:
225 *
226  anorm = max( pclange( norma, n, n, a, 1, 1, desca, rwork ), unfl )
227 *
228 * Norm of E:
229 *
230  enorm = max( pclange( norme, n, n, e, 1, 1, desce, rwork ), ulp )
231 *
232 * Norm of error:
233 *
234 * Error = AE - EW
235 *
236  CALL pclaset( 'Full', n, n, czero, czero, work, 1, 1, descw )
237 *
238  DO 60 jcol = 1, n
239  IF( itrnsw.EQ.0 ) THEN
240  wtemp = w( jcol )
241  ELSE
242  wtemp = conjg( w( jcol ) )
243  END IF
244 *
245  IF( itrnse.EQ.0 ) THEN
246  CALL pcaxpy( n, wtemp, e, 1, jcol, desce, 1, work, 1, jcol,
247  $ descw, 1 )
248  ELSE IF( itrnse.EQ.1 ) THEN
249  CALL pcaxpy( n, wtemp, e, jcol, 1, desce, n, work, 1, jcol,
250  $ descw, 1 )
251  ELSE
252  CALL pcaxpy( n, conjg( wtemp ), e, jcol, 1, desce, n, work,
253  $ 1, jcol, descw, 1 )
254  DO 50 jrow = 1, n
255  CALL infog2l( jrow, jcol, descw, nprow, npcol, myrow,
256  $ mycol, irow, icol, ii, jj )
257  IF( ( myrow.EQ.ii ) .AND. ( mycol.EQ.jj ) ) THEN
258  work( ( jcol-1 )*ldw+jrow )
259  $ = conjg( work( ( jcol-1 )*ldw+jrow ) )
260  END IF
261  50 CONTINUE
262  END IF
263  60 CONTINUE
264 *
265  CALL pcgemm( transa, transe, n, n, n, cone, a, 1, 1, desca, e, 1,
266  $ 1, desce, -cone, work, 1, 1, descw )
267 *
268  errnrm = pclange( 'One', n, n, work, 1, 1, descw, rwork ) / enorm
269 *
270 * Compute RESULT(1) (avoiding under/overflow)
271 *
272  IF( anorm.GT.errnrm ) THEN
273  result( 1 ) = ( errnrm / anorm ) / ulp
274  ELSE
275  IF( anorm.LT.one ) THEN
276  result( 1 ) = ( min( errnrm, anorm ) / anorm ) / ulp
277  ELSE
278  result( 1 ) = min( errnrm / anorm, one ) / ulp
279  END IF
280  END IF
281 *
282 * Compute RESULT(2) : the normalization error in E.
283 *
284  result( 2 ) = max( abs( enrmax-one ), abs( enrmin-one ) ) /
285  $ ( real( n )*ulp )
286 *
287  RETURN
288 *
289 * End of PCGET22
290 *
291  END
max
#define max(A, B)
Definition: pcgemr.c:180
infog2l
subroutine infog2l(GRINDX, GCINDX, DESC, NPROW, NPCOL, MYROW, MYCOL, LRINDX, LCINDX, RSRC, CSRC)
Definition: infog2l.f:3
pcget22
subroutine pcget22(TRANSA, TRANSE, TRANSW, N, A, DESCA, E, DESCE, W, WORK, DESCW, RWORK, RESULT)
Definition: pcget22.f:3
pclaset
subroutine pclaset(UPLO, M, N, ALPHA, BETA, A, IA, JA, DESCA)
Definition: pcblastst.f:7508
min
#define min(A, B)
Definition: pcgemr.c:181