ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pzgsepchk.f
Go to the documentation of this file.
1 *
2 *
3  SUBROUTINE pzgsepchk( IBTYPE, MS, NV, A, IA, JA, DESCA, B, IB, JB,
4  $ DESCB, THRESH, Q, IQ, JQ, DESCQ, C, IC, JC,
5  $ DESCC, W, WORK, LWORK, TSTNRM, RESULT )
6 *
7 * -- ScaLAPACK routine (version 1.7) --
8 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
9 * and University of California, Berkeley.
10 * April 15, 1997
11 *
12 * .. Scalar Arguments ..
13  INTEGER IA, IB, IBTYPE, IC, IQ, JA, JB, JC, JQ, LWORK,
14  $ MS, NV, RESULT
15  DOUBLE PRECISION THRESH, TSTNRM
16 * ..
17 * .. Array Arguments ..
18 *
19  INTEGER DESCA( * ), DESCB( * ), DESCC( * ), DESCQ( * )
20  DOUBLE PRECISION W( * ), WORK( * )
21  COMPLEX*16 A( * ), B( * ), C( * ), Q( * )
22 * ..
23 *
24 *
25 * Purpose
26 * =======
27 *
28 * PZGSEPCHK checks a decomposition of the form
29 *
30 * A Q = B Q D or
31 * A B Q = Q D or
32 * B A Q = Q D
33 *
34 * where A is a symmetric matrix, B is
35 * symmetric positive definite, Q is orthogonal, and D is diagonal.
36 *
37 * One of the following test ratios is computed:
38 *
39 * IBTYPE = 1: TSTNRM = | A Q - B Q D | / ( |A| |Q| n ulp )
40 *
41 * IBTYPE = 2: TSTNRM = | A B Q - Q D | / ( |A| |Q| n ulp )
42 *
43 * IBTYPE = 3: TSTNRM = | B A Q - Q D | / ( |A| |Q| n ulp )
44 *
45 *
46 * Notes
47 * =====
48 *
49 *
50 * Each global data object is described by an associated description
51 * vector. This vector stores the information required to establish
52 * the mapping between an object element and its corresponding process
53 * and memory location.
54 *
55 * Let A be a generic term for any 2D block cyclicly distributed array.
56 * Such a global array has an associated description vector DESCA.
57 * In the following comments, the character _ should be read as
58 * "of the global array".
59 *
60 * NOTATION STORED IN EXPLANATION
61 * --------------- -------------- --------------------------------------
62 * DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
63 * DTYPE_A = 1.
64 * CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
65 * the BLACS process grid A is distribu-
66 * ted over. The context itself is glo-
67 * bal, but the handle (the integer
68 * value) may vary.
69 * M_A (global) DESCA( M_ ) The number of rows in the global
70 * array A.
71 * N_A (global) DESCA( N_ ) The number of columns in the global
72 * array A.
73 * MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
74 * the rows of the array.
75 * NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
76 * the columns of the array.
77 * RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
78 * row of the array A is distributed.
79 * CSRC_A (global) DESCA( CSRC_ ) The process column over which the
80 * first column of the array A is
81 * distributed.
82 * LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
83 * array. LLD_A >= MAX(1,LOCr(M_A)).
84 *
85 * Let K be the number of rows or columns of a distributed matrix,
86 * and assume that its process grid has dimension p x q.
87 * LOCr( K ) denotes the number of elements of K that a process
88 * would receive if K were distributed over the p processes of its
89 * process column.
90 * Similarly, LOCc( K ) denotes the number of elements of K that a
91 * process would receive if K were distributed over the q processes of
92 * its process row.
93 * The values of LOCr() and LOCc() may be determined via a call to the
94 * ScaLAPACK tool function, NUMROC:
95 * LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
96 * LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
97 * An upper bound for these quantities may be computed by:
98 * LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
99 * LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
100 *
101 *
102 * Arguments
103 * =========
104 *
105 * MP = number of local rows in A, B and Q
106 * MQ = number of local columns in A
107 * NQ = number of local columns in B and Q
108 *
109 * IBTYPE (input) INTEGER
110 * The form of the symmetric generalized eigenproblem.
111 * = 1: A*Q = (lambda)*B*Q
112 * = 2: A*B*Q = (lambda)*Q
113 * = 3: B*A*Q = (lambda)*Q
114 *
115 * MS (global input) INTEGER
116 * Matrix size.
117 * The number of global rows in A, B, C and Q
118 * Also, the number of columns in A
119 *
120 * NV (global input) INTEGER
121 * Number of eigenvectors
122 * The number of global columns in C and Q.
123 *
124 * A (local input) REAL pointer to an
125 * array in local memory of dimension (LLD_A, LOCc(JA+N-1)).
126 * This array contains the local pieces of the M-by-N
127 * distributed test matrix A
128 *
129 * IA (global input) INTEGER
130 * A's global row index, which points to the beginning of the
131 * submatrix which is to be operated on.
132 *
133 * JA (global input) INTEGER
134 * A's global column index, which points to the beginning of
135 * the submatrix which is to be operated on.
136 *
137 * DESCA (global and local input) INTEGER array of dimension 8
138 * The array descriptor for the distributed matrix A.
139 *
140 * B (local input) REAL pointer to an
141 * array in local memory of dimension (LLD_B, LOCc(JB+N-1)).
142 * This array contains the local pieces of the M-by-N
143 * distributed test matrix B
144 *
145 * IB (global input) INTEGER
146 * B's global row index, which points to the beginning of the
147 * submatrix which is to be operated on.
148 *
149 * JB (global input) INTEGER
150 * B's global column index, which points to the beginning of
151 * the submatrix which is to be operated on.
152 *
153 * DESCB (global and local input) INTEGER array of dimension 8
154 * The array descriptor for the distributed matrix B.
155 *
156 * THRESH (input) REAL
157 * A test will count as "failed" if the "error", computed as
158 * described below, exceeds THRESH. Note that the error
159 * is scaled to be O(1), so THRESH should be a reasonably
160 * small multiple of 1, e.g., 10 or 100. In particular,
161 * it should not depend on the precision (single vs. double)
162 * or the size of the matrix. It must be at least zero.
163 *
164 * Q (local input) REAL array
165 * global dimension (MS, NV),
166 * local dimension (DESCA( DLEN_ ), NQ)
167 *
168 * Contains the eigenvectors as computed by PSSYEVX
169 *
170 * IQ (global input) INTEGER
171 * Q's global row index, which points to the beginning of the
172 * submatrix which is to be operated on.
173 *
174 * JQ (global input) INTEGER
175 * Q's global column index, which points to the beginning of
176 * the submatrix which is to be operated on.
177 *
178 * DESCQ (global and local input) INTEGER array of dimension 8
179 * The array descriptor for the distributed matrix Q.
180 *
181 * C (local workspace) REAL array,
182 * global dimension (MS, NV),
183 * local dimension (DESCA( DLEN_ ), MQ)
184 *
185 * Accumulator for computing AQ -QL
186 *
187 * IC (global input) INTEGER
188 * C's global row index, which points to the beginning of the
189 * submatrix which is to be operated on.
190 *
191 * JC (global input) INTEGER
192 * C's global column index, which points to the beginning of
193 * the submatrix which is to be operated on.
194 *
195 * DESCC (global and local input) INTEGER array of dimension 8
196 * The array descriptor for the distributed matrix C.
197 *
198 * W (global input) REAL array, dimension (NV)
199 *
200 * Contains the computed eigenvalues
201 *
202 * WORK (local workspace) REAL array,
203 * dimension (LWORK)
204 *
205 * LWORK (local input) INTEGER
206 * The length of the array WORK.
207 * LWORK >= NUMROC( NV, DESCA( NB_ ), MYCOL, 0, NPCOL )
208 *
209 * TSTNRM (global output) REAL
210 *
211 * RESULT (global output) INTEGER
212 * 0 if the test passes
213 * 1 if the test fails
214 *
215 * .. Local Scalars ..
216 *
217  INTEGER I, INFO, MYCOL, MYROW, NPCOL, NPROW, NQ
218  DOUBLE PRECISION ANORM, ULP
219 * ..
220 * .. Parameters ..
221  INTEGER BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_,
222  $ MB_, NB_, RSRC_, CSRC_, LLD_
223  parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
224  $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
225  $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
226  DOUBLE PRECISION ONE, ZERO
227  PARAMETER ( ONE = 1.0d+0, zero = 0.0d+0 )
228  COMPLEX*16 CONE, CNEGONE, CZERO
229  parameter( cone = 1.0d+0, cnegone = -1.0d+0,
230  $ czero = 0.0d+0 )
231 * ..
232 * .. External Functions ..
233  INTEGER NUMROC
234  DOUBLE PRECISION DLAMCH, PZLANGE
235  EXTERNAL numroc, dlamch, pzlange
236 * ..
237 * .. External Subroutines ..
238  EXTERNAL blacs_gridinfo, chk1mat, pxerbla, pzdscal,
239  $ pzgemm
240 * ..
241 * .. Intrinsic Functions ..
242  INTRINSIC max
243 * ..
244 * .. Executable Statements ..
245 * This is just to keep ftnchek happy
246  IF( block_cyclic_2d*csrc_*ctxt_*dlen_*dtype_*lld_*mb_*m_*nb_*n_*
247  $ rsrc_.LT.0 )RETURN
248 *
249  result = 0
250 *
251  CALL blacs_gridinfo( desca( ctxt_ ), nprow, npcol, myrow, mycol )
252 *
253  info = 0
254  CALL chk1mat( ms, 1, ms, 2, ia, ja, desca, 7, info )
255  CALL chk1mat( ms, 1, ms, 2, ib, jb, descb, 11, info )
256  CALL chk1mat( ms, 1, nv, 2, iq, jq, descq, 16, info )
257  CALL chk1mat( ms, 1, nv, 2, ib, jb, descb, 20, info )
258 *
259  IF( info.EQ.0 ) THEN
260 *
261  nq = numroc( nv, desca( nb_ ), mycol, 0, npcol )
262 *
263  IF( iq.NE.1 ) THEN
264  info = -14
265  ELSE IF( jq.NE.1 ) THEN
266  info = -15
267  ELSE IF( ia.NE.1 ) THEN
268  info = -5
269  ELSE IF( ja.NE.1 ) THEN
270  info = -6
271  ELSE IF( ib.NE.1 ) THEN
272  info = -9
273  ELSE IF( jb.NE.1 ) THEN
274  info = -10
275  ELSE IF( lwork.LT.nq ) THEN
276  info = -23
277  END IF
278  END IF
279 *
280  IF( info.NE.0 ) THEN
281  CALL pxerbla( desca( ctxt_ ), 'PZGSEPCHK', -info )
282  RETURN
283  END IF
284 *
285  result = 0
286  ulp = dlamch( 'Epsilon' )
287 *
288 * Compute product of Max-norms of A and Q.
289 *
290  anorm = pzlange( 'M', ms, ms, a, ia, ja, desca, work )*
291  $ pzlange( 'M', ms, nv, q, iq, jq, descq, work )
292  IF( anorm.EQ.zero )
293  $ anorm = one
294 *
295  IF( ibtype.EQ.1 ) THEN
296 *
297 * Norm of AQ - BQD
298 *
299 * C = AQ
300 *
301  CALL pzgemm( 'N', 'N', ms, nv, ms, cone, a, ia, ja, desca, q,
302  $ iq, jq, descq, czero, c, ic, jc, descc )
303 *
304 * Q = QD
305 *
306  DO 10 i = 1, nv
307  CALL pzdscal( ms, w( i ), q, iq, jq+i-1, descq, 1 )
308  10 CONTINUE
309 *
310 * C = C - BQ (i.e. AQ-BQD)
311 *
312  CALL pzgemm( 'N', 'N', ms, nv, ms, cone, b, ib, jb, descb, q,
313  $ iq, jq, descq, cnegone, c, ic, jc, descc )
314 *
315  tstnrm = ( pzlange( 'M', ms, nv, c, ic, jc, descc, work ) /
316  $ anorm ) / ( max( ms, 1 )*ulp )
317 *
318 *
319  ELSE IF( ibtype.EQ.2 ) THEN
320 *
321 * Norm of ABQ - QD
322 *
323 *
324 * C = BQ
325 *
326  CALL pzgemm( 'N', 'N', ms, nv, ms, cone, b, ib, jb, descb, q,
327  $ iq, jq, descq, czero, c, ic, jc, descc )
328 *
329 * Q = QD
330 *
331  DO 20 i = 1, nv
332  CALL pzdscal( ms, w( i ), q, iq, jq+i-1, descq, 1 )
333  20 CONTINUE
334 *
335 * Q = AC - Q
336 *
337  CALL pzgemm( 'N', 'N', ms, nv, ms, cone, a, ia, ja, desca, c,
338  $ ic, jc, descc, cnegone, q, iq, jq, descq )
339 *
340  tstnrm = ( pzlange( 'M', ms, nv, q, iq, jq, descq, work ) /
341  $ anorm ) / ( max( ms, 1 )*ulp )
342 *
343  ELSE IF( ibtype.EQ.3 ) THEN
344 *
345 * Norm of BAQ - QD
346 *
347 *
348 * C = AQ
349 *
350  CALL pzgemm( 'N', 'N', ms, nv, ms, cone, a, ia, ja, desca, q,
351  $ iq, jq, descq, czero, c, ic, jc, descc )
352 *
353 * Q = QD
354 *
355  DO 30 i = 1, nv
356  CALL pzdscal( ms, w( i ), q, iq, jq+i-1, descq, 1 )
357  30 CONTINUE
358 *
359 * Q = BC - Q
360 *
361  CALL pzgemm( 'N', 'N', ms, nv, ms, cone, b, ib, jb, descb, c,
362  $ ic, jc, descc, cnegone, q, iq, jq, descq )
363 *
364  tstnrm = ( pzlange( 'M', ms, nv, q, iq, jq, descq, work ) /
365  $ anorm ) / ( max( ms, 1 )*ulp )
366 *
367  END IF
368 *
369  IF( tstnrm.GT.thresh .OR. ( tstnrm-tstnrm.NE.0.0d0 ) ) THEN
370  result = 1
371  END IF
372  RETURN
373 *
374 * End of PZGSEPCHK
375 *
376  END
pzgsepchk
subroutine pzgsepchk(IBTYPE, MS, NV, A, IA, JA, DESCA, B, IB, JB, DESCB, THRESH, Q, IQ, JQ, DESCQ, C, IC, JC, DESCC, W, WORK, LWORK, TSTNRM, RESULT)
Definition: pzgsepchk.f:6
max
#define max(A, B)
Definition: pcgemr.c:180
chk1mat
subroutine chk1mat(MA, MAPOS0, NA, NAPOS0, IA, JA, DESCA, DESCAPOS0, INFO)
Definition: chk1mat.f:3
pxerbla
subroutine pxerbla(ICTXT, SRNAME, INFO)
Definition: pxerbla.f:2