SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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 IF( itrnsw.EQ.0 ) THEN
240 wtemp = w( jcol )
241 ELSE
242 wtemp = dconjg( w( jcol ) )
243 END IF
244*
245 IF( itrnse.EQ.0 ) THEN
246 CALL pzaxpy( n, wtemp, e, 1, jcol, desce, 1, work, 1, jcol,
247 $ descw, 1 )
248 ELSE IF( itrnse.EQ.1 ) 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 IF( ( myrow.EQ.ii ) .AND. ( mycol.EQ.jj ) ) 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 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 $ ( dble( n )*ulp )
286*
287 RETURN
288*
289* End of PZGET22
290*
infog2l
subroutine infog2l(grindx, gcindx, desc, nprow, npcol, myrow, mycol, lrindx, lcindx, rsrc, csrc)
Definition infog2l.f:3
max
#define max(A, B)
Definition pcgemr.c:180
min
#define min(A, B)
Definition pcgemr.c:181
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
pzlange
double precision function pzlange(norm, m, n, a, ia, ja, desca, work)
Definition pzlange.f:3
lsame
logical function lsame(ca, cb)
Definition tools.f:1724
Here is the call graph for this function:
Here is the caller graph for this function: