OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
pzget22.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine pzget22 (transa, transe, transw, n, a, desca, e, desce, w, work, descw, rwork, result)

Function/Subroutine Documentation

◆ pzget22()

subroutine pzget22 ( character transa,
character transe,
character transw,
integer n,
complex*16, dimension( * ) a,
integer, dimension( * ) desca,
complex*16, dimension( * ) e,
integer, dimension( * ) desce,
complex*16, dimension( * ) w,
complex*16, dimension( * ) work,
integer, dimension( * ) descw,
double precision, dimension( * ) rwork,
double precision, dimension( 2 ) result )

Definition at line 1 of file pzget22.f.

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 DOUBLE PRECISION RESULT( 2 ), RWORK( * )
16 COMPLEX*16 A( * ), E( * ), W( * ), WORK( * )
17* ..
18*
19* Purpose
20* =======
21*
22* PZGET22 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*16 array, dimension (*)
62* The matrix whose eigenvectors are in E.
63*
64* DESCA (input) INTEGER array, dimension(*)
65*
66* E (input) COMPLEX*16 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*16 array, dimension (N)
74* The eigenvalues of A.
75*
76* WORK (workspace) COMPLEX*16 array, dimension (*)
77* DESCW (input) INTEGER array, dimension(*)
78*
79* RWORK (workspace) DOUBLE PRECISION array, dimension (N)
80*
81* RESULT (output) DOUBLE PRECISION 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 DOUBLE PRECISION ZERO, ONE
99 parameter( zero = 0.0d+0, one = 1.0d+0 )
100 COMPLEX*16 CZERO, CONE
101 parameter( czero = ( 0.0d+0, 0.0d+0 ),
102 $ cone = ( 1.0d+0, 0.0d+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, CA, CSRC, RA, RSRC
109 DOUBLE PRECISION ANORM, ENORM, ENRMAX, ENRMIN, ERRNRM, TEMP1,
110 $ ULP, UNFL
111 COMPLEX*16 CDUM, WTEMP
112* ..
113* .. External Functions ..
114 LOGICAL LSAME
115 DOUBLE PRECISION PDLAMCH, PZLANGE
116 EXTERNAL lsame, pdlamch, pzlange
117* ..
118* .. External Subroutines ..
119 EXTERNAL blacs_gridinfo, dgamn2d, dgamx2d, infog2l,
120 $ pzaxpy, pzgemm, pzlaset
121* ..
122* .. Intrinsic Functions ..
123 INTRINSIC abs, dble, dconjg, dimag, max, min
124* ..
125* .. Statement Functions ..
126 DOUBLE PRECISION CABS1
127* ..
128* .. Statement Function definitions ..
129 cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( 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 = pdlamch( contxt, 'Safe minimum' )
151 ulp = pdlamch( 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 dgamx2d( 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 dgamx2d( contxt, 'Row', ' ', 1, 1, enrmax, 1, ra, ca, -1,
197 $ -1, -1 )
198 CALL dgamn2d( 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 dgamx2d( 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 dgamx2d( contxt, 'Row', ' ', 1, 1, enrmax, 1, ra, ca, -1,
219 $ -1, -1 )
220 CALL dgamn2d( contxt, 'Row', ' ', 1, 1, enrmin, 1, ra, ca, -1,
221 $ -1, -1 )
222 END IF
223*
224* Norm of A:
225*
226 anorm = max( pzlange( norma, n, n, a, 1, 1, desca, rwork ), unfl )
227*
228* Norm of E:
229*
230 enorm = max( pzlange( norme, n, n, e, 1, 1, desce, rwork ), ulp )
231*
232* Norm of error:
233*
234* Error = AE - EW
235*
236 CALL pzlaset( '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 = DCONJG( W( JCOL ) )
243 END IF
244*
245.EQ. IF( ITRNSE0 ) THEN
246 CALL PZAXPY( N, WTEMP, E, 1, JCOL, DESCE, 1, WORK, 1, JCOL,
247 $ DESCW, 1 )
248.EQ. ELSE IF( ITRNSE1 ) THEN
249 CALL PZAXPY( N, WTEMP, E, JCOL, 1, DESCE, N, WORK, 1, JCOL,
250 $ DESCW, 1 )
251 ELSE
252 CALL PZAXPY( N, DCONJG( 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 $ = DCONJG( WORK( ( JCOL-1 )*LDW+JROW ) )
260 END IF
261 50 CONTINUE
262 END IF
263 60 CONTINUE
264*
265 CALL PZGEMM( TRANSA, TRANSE, N, N, N, CONE, A, 1, 1, DESCA, E, 1,
266 $ 1, DESCE, -CONE, WORK, 1, 1, DESCW )
267*
268 ERRNRM = PZLANGE( '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 $ ( DBLE( N )*ULP )
286*
287 RETURN
288*
289* End of PZGET22
290*
lde
logical function lde(ri, rj, lr)
Definition dblat2.f:2942
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:53
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
pzaxpy
subroutine pzaxpy(n, a, x, ix, jx, descx, incx, y, iy, jy, descy, incy)
Definition mpi.f:1436
infog2l
subroutine infog2l(grindx, gcindx, desc, nprow, npcol, myrow, mycol, lrindx, lcindx, rsrc, csrc)
Definition mpi.f:937
pzlange
double precision function pzlange(norm, m, n, a, ia, ja, desca, work)
Definition mpi.f:1287
blacs_gridinfo
subroutine blacs_gridinfo(cntxt, nprow, npcol, myrow, mycol)
Definition mpi.f:754
pdlamch
double precision function pdlamch(ictxt, cmach)
Definition pdblastst.f:6769
pzlaset
subroutine pzlaset(uplo, m, n, alpha, beta, a, ia, ja, desca)
Definition pzblastst.f:7509