OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
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.EQ..AND..EQ. IF( ( MYROWII ) ( MYCOLJJ ) ) THEN
207 TEMP1 = MAX( TEMP1, CABS1( E( ( ICOL-1 )*LDE+
208 $ IROW ) ) )
209 END IF
210 30 CONTINUE
211.EQ. IF( MYROWII ) 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.EQ. IF( ITRNSW0 ) THEN
240 WTEMP = W( JCOL )
241 ELSE
242 WTEMP = CONJG( W( JCOL ) )
243 END IF
244*
245.EQ. IF( ITRNSE0 ) THEN
246 CALL PCAXPY( N, WTEMP, E, 1, JCOL, DESCE, 1, WORK, 1, JCOL,
247 $ DESCW, 1 )
248.EQ. ELSE IF( ITRNSE1 ) 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.EQ..AND..EQ. IF( ( MYROWII ) ( MYCOLJJ ) ) 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.GT. IF( ANORMERRNRM ) THEN
273 RESULT( 1 ) = ( ERRNRM / ANORM ) / ULP
274 ELSE
275.LT. IF( ANORMONE ) 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
lde
logical function lde(ri, rj, lr)
Definition dblat2.f:2942
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
infog2l
subroutine infog2l(grindx, gcindx, desc, nprow, npcol, myrow, mycol, lrindx, lcindx, rsrc, csrc)
Definition mpi.f:937
pcaxpy
subroutine pcaxpy(n, a, x, ix, jx, descx, incx, y, iy, jy, descy, incy)
Definition mpi.f:1425
blacs_gridinfo
subroutine blacs_gridinfo(cntxt, nprow, npcol, myrow, mycol)
Definition mpi.f:754
pclaset
subroutine pclaset(uplo, m, n, alpha, beta, a, ia, ja, desca)
Definition pcblastst.f:7508
pcget22
subroutine pcget22(transa, transe, transw, n, a, desca, e, desce, w, work, descw, rwork, result)
Definition pcget22.f:3