ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pssvdchk.f
Go to the documentation of this file.
1  SUBROUTINE pssvdchk( M, N, A, IA, JA, DESCA, U, IU, JU, DESCU, VT,
2  $ IVT, JVT, DESCVT, S, THRESH, WORK, LWORK,
3  $ RESULT, CHK, MTM )
4 *
5 * -- ScaLAPACK routine (version 1.7) --
6 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
7 * and University of California, Berkeley.
8 * May 1, 1997
9 *
10 * .. Scalar Arguments ..
11  INTEGER IA, IU, IVT, JA, JU, JVT, LWORK, M, N
12  REAL CHK, MTM, THRESH
13 * ..
14 * .. Array Arguments ..
15  INTEGER DESCA( * ), DESCU( * ), DESCVT( * ),
16  $ RESULT( * )
17  REAL A( * ), S( * ), U( * ), VT( * ), WORK( * )
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 * For given two-dimensional matrices A, U, VT, and one-dimensional
24 * array D compute the following four tests:
25 *
26 * (1) | A - U*diag(S) VT | / ( |A| max(M,N) ulp )
27 *
28 * (2) | I - U'*U | / ( M ulp )
29 *
30 * (3) | I - VT*VT' | / ( N ulp ),
31 *
32 * (4) S contains SIZE = MIN( M, N ) nonnegative values in
33 * decreasing order.
34 * It then compares result of computations (1)-(3)
35 * with TRESH and returns results of comparisons and test (4) in
36 * RESULT(I). When the i-th test fails, value of RESULT( I ) is set
37 * to 1.
38 *
39 * Notes
40 * =====
41 *
42 * Each global data object is described by an associated description
43 * vector. This vector stores the information required to establish
44 * the mapping between an object element and its corresponding process
45 * and memory location.
46 *
47 * Let A be a generic term for any 2D block cyclicly distributed array.
48 * Such a global array has an associated description vector DESCA.
49 * In the following comments, the character _ should be read as
50 * "of the global array".
51 *
52 * NOTATION STORED IN EXPLANATION
53 * --------------- -------------- --------------------------------------
54 * DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
55 * DTYPE_A = 1.
56 * CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
57 * the BLACS process grid A is distribu-
58 * ted over. The context itself is glo-
59 * bal, but the handle (the integer
60 * value) may vary.
61 * M_A (global) DESCA( M_ ) The number of rows in the global
62 * array A.
63 * N_A (global) DESCA( N_ ) The number of columns in the global
64 * array A.
65 * MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
66 * the rows of the array.
67 * NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
68 * the columns of the array.
69 * RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
70 * row of the array A is distributed.
71 * CSRC_A (global) DESCA( CSRC_ ) The process column over which the
72 * first column of the array A is
73 * distributed.
74 * LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
75 * array. LLD_A >= MAX(1,LOCr(M_A)).
76 *
77 * Let K be the number of rows or columns of a distributed matrix,
78 * and assume that its process grid has dimension p x q.
79 * LOCr( K ) denotes the number of elements of K that a process
80 * would receive if K were distributed over the p processes of its
81 * process column.
82 * Similarly, LOCc( K ) denotes the number of elements of K that a
83 * process would receive if K were distributed over the q processes of
84 * its process row.
85 * The values of LOCr() and LOCc() may be determined via a call to the
86 * ScaLAPACK tool function, NUMROC:
87 * LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
88 * LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
89 * An upper bound for these quantities may be computed by:
90 * LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
91 * LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
92 *
93 * Arguments
94 * =========
95 *
96 * MP = number of local rows in A and U
97 * NQ = number of local columns in A and VT
98 * SIZEP = number of local rows in VT
99 * SIZEQ = number of local columns in U
100 *
101 * M (global input) INTEGER
102 * Matrix size.
103 * The number of global rows in A and U and
104 *
105 * N (global input) INTEGER
106 * The number of global columns in A and VT.
107 *
108 * A (input) block cyclic distributed REAL array,
109 * global dimension (M, N), local dimension (DESCA( DLEN_ ), NQ)
110 * Contains the original test matrix.
111 *
112 * IA (global input) INTEGER
113 * The global row index of the submatrix of the distributed
114 * matrix A to operate on.
115 *
116 * JA (global input) INTEGER
117 * The global column index of the submatrix of the distributed
118 * matrix A to operate on.
119 *
120 * DESCA (global and local input) INTEGER array of dimension DLEN_
121 * The array descriptor for the distributed matrix A.
122 *
123 * U (local input) REAL array
124 * global dimension (M, SIZE), local dimension
125 * (DESCU( DLEN_ ), SIZEQ)
126 * Contains left singular vectors of matrix A.
127 *
128 * IU (global input) INTEGER
129 * The global row index of the submatrix of the distributed
130 * matrix U to operate on.
131 *
132 * JU (global input) INTEGER
133 * The global column index of the submatrix of the distributed
134 * matrix U to operate on.
135 *
136 * DESCU (global and local input) INTEGER array of dimension DLEN_
137 * The array descriptor for the distributed matrix U.
138 *
139 * VT (local input) REAL array
140 * global dimension (SIZE, N), local dimension
141 * (DESCVT( DLEN_ ), NQ)
142 * Contains right singular vectors of matrix A.
143 *
144 * IVT (global input) INTEGER
145 * The global row index of the submatrix of the distributed
146 * matrix VT to operate on.
147 *
148 * JVT (global input) INTEGER
149 * The global column index of the submatrix of the distributed
150 * matrix VT to operate on.
151 *
152 * DESCVT (global and local input) INTEGER array of dimension DLEN_
153 * The array descriptor for the distributed matrix VT.
154 *
155 * S (global input) REAL array, dimension (SIZE)
156 * Contains the computed singular values
157 *
158 * THRESH (input) REAL
159 * A test will count as "failed" if the "error", computed as
160 * described below, exceeds THRESH. Note that the error
161 * is scaled to be O(1), so THRESH should be a reasonably
162 * small multiple of 1, e.g., 10 or 100. In particular,
163 * it should not depend on the precision (single vs. double)
164 * or the size of the matrix. It must be at least zero.
165 *
166 * WORK (local workspace) REAL array, dimension (LWORK)
167 *
168 * LWORK (local input) INTEGER
169 * The length of the array WORK.
170 * LWORK >= 1 + SIZEQ*SIZEP + MAX[WORK(pdlange(size,size)),
171 * WORK(pdlange(m,n))],
172 * where
173 * SIZEQ = NUMROC( SIZE, DESCU( NB_ ), MYCOL, 0, NPCOL ),
174 * SIZEP = NUMROC( SIZE, DESCVT( MB_ ), MYROW, 0, NPROW ),
175 * and worekspaces required to call pdlange are
176 * WORK(pdlange(size,size)) < MAX(SIZEQ0,2) < SIZEB +2,
177 * WORK(pdlange(m,n)) < MAX(NQ0,2) < SIZEB +2,
178 * SIZEB = MAX(M, N)
179 * Finally, upper limit on required workspace is
180 * LWORK > 1 + SIZEQ*SIZEP + SIZEB + 2
181 *
182 * RESULT (global input/output) INTEGER array. Four first elements of
183 * the array are set to 0 or 1 depending on passing four
184 * respective tests ( see above in Purpose ). The elements of
185 * RESULT are set to
186 * 0 if the test passes i.e.
187 * | A - U*diag(S)*VT | / ( |A| max(M,N) ulp ) <= THRESH
188 * 1 if the test fails i.e.
189 * | A - U*diag(S)*VT | / ( |A| max(M,N) ulp ) > THRESH
190 *
191 * CHK (global output) REAL
192 * value of the | A - U*diag(S) VT | / ( |A| max(M,N) ulp )
193 *
194 * MTM (global output) REAL
195 * maximum of the two values:
196 * | I - U'*U | / ( M ulp ) and | I - VT*VT' | / ( N ulp )
197 *
198 * ======================================================================
199 *
200 * .. Parameters ..
201  INTEGER BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_,
202  $ MB_, NB_, RSRC_, CSRC_, LLD_
203  parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
204  $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
205  $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
206  REAL ZERO, ONE, MONE
207  PARAMETER ( ZERO = 0.0e+0, one = 1.0e+0, mone = -1.0e0 )
208 * ..
209 * .. Local Scalars ..
210  INTEGER I, INFO, LDR, LOCALCOL, LWMIN, MP, MX, MYCOL,
211  $ MYROW, NPCOL, NPROW, NQ, PCOL, PTRR, PTRWORK,
212  $ SIZE, sizep, sizepos, sizeq
213  REAL FIRST, NORMA, NORMAI, NORMU, NORMVT, SECOND,
214  $ THRESHA, ULP
215 * ..
216 * .. Local Arrays ..
217  INTEGER DESCR( DLEN_ )
218 * ..
219 * .. External Functions ..
220  INTEGER INDXG2L, INDXG2P, NUMROC
221  REAL PSLAMCH, PSLANGE
222  EXTERNAL indxg2l, indxg2p, numroc, pslamch, pslange
223 * ..
224 * .. External Subroutines ..
225  EXTERNAL blacs_gridinfo, chk1mat, descinit, sscal,
226  $ pselset, psgemm, pslaset, pxerbla
227 * ..
228 * .. Intrinsic Functions ..
229  INTRINSIC max, min
230 * ..
231 * .. Executable Statements ..
232 *
233 * Test the input parameters.
234 *
235  CALL blacs_gridinfo( desca( ctxt_ ), nprow, npcol, myrow, mycol )
236  info = 0
237  SIZE = min( m, n )
238 *
239 * Sizepos is a number of parameters to pdsvdchk plus one. It's used
240 * for the error reporting.
241 *
242  sizepos = 22
243  IF( nprow.EQ.-1 ) THEN
244  info = -607
245  ELSE
246  CALL chk1mat( m, 1, n, 2, ia, ja, desca, 6, info )
247  CALL chk1mat( m, 1, SIZE, sizepos, iu, ju, descu, 10, info )
248  CALL chk1mat( SIZE, sizepos, n, 2, ivt, jvt, descvt, 14, info )
249  END IF
250 *
251  IF( info.EQ.0 ) THEN
252 *
253 * Calculate workspace
254 *
255  mp = numroc( m, desca( mb_ ), myrow, 0, nprow )
256  nq = numroc( n, desca( nb_ ), mycol, 0, npcol )
257  sizep = numroc( SIZE, descvt( mb_ ), myrow, 0, nprow )
258  sizeq = numroc( SIZE, descu( nb_ ), mycol, 0, npcol )
259  mx = max( sizeq, nq )
260  lwmin = 2 + sizeq*sizep + max( 2, mx )
261  work( 1 ) = lwmin
262  IF( lwork.EQ.-1 )
263  $ GO TO 40
264  IF( lwork.LT.lwmin ) THEN
265  info = -18
266  ELSE IF( thresh.LE.0 ) THEN
267  info = -16
268  END IF
269  END IF
270  IF( info.NE.0 ) THEN
271  CALL pxerbla( desca( ctxt_ ), 'PSSVDCHK', -info )
272  RETURN
273  END IF
274 *
275  ldr = max( 1, sizep )
276  ulp = pslamch( desca( ctxt_ ), 'P' )
277  normai = pslange( '1', m, n, a, 1, 1, desca, work )
278 *
279 * Allocate array R of global dimension SIZE x SIZE for testing
280 *
281  ptrr = 2
282  ptrwork = ptrr + sizeq*sizep
283 *
284  CALL descinit( descr, SIZE, SIZE, descvt( mb_ ), descu( nb_ ), 0,
285  $ 0, desca( ctxt_ ), ldr, info )
286 *
287 * Test 2. Form identity matrix R and make check norm(U'*U - I )
288 *
289  CALL pslaset( 'Full', SIZE, SIZE, zero, one, work( ptrr ), 1, 1,
290  $ descr )
291  CALL psgemm( 'T', 'N', SIZE, SIZE, m, one, u, 1, 1, descu, u, 1,
292  $ 1, descu, mone, work( ptrr ), 1, 1, descr )
293 *
294  normu = pslange( '1', SIZE, SIZE, work( ptrr ), 1, 1, descr,
295  $ work( ptrwork ) )
296 *
297  normu = normu / ulp / SIZE / thresh
298  IF( normu.GT.1. )
299  $ result( 2 ) = 1
300 *
301 * Test3. Form identity matrix R and check norm(VT*VT' - I )
302 *
303  CALL pslaset( 'Full', SIZE, SIZE, zero, one, work( ptrr ), 1, 1,
304  $ descr )
305  CALL psgemm( 'N', 'T', SIZE, SIZE, n, one, vt, 1, 1, descvt, vt,
306  $ 1, 1, descvt, mone, work( ptrr ), 1, 1, descr )
307  normvt = pslange( '1', SIZE, SIZE, work( ptrr ), 1, 1, descr,
308  $ work( ptrwork ) )
309 *
310  normvt = normvt / ulp / SIZE / thresh
311  IF( normvt.GT.1. )
312  $ result( 3 ) = 1
313 *
314  mtm = max( normvt, normu )*thresh
315 *
316 * Test 1.
317 * Initialize R = diag( S )
318 *
319  CALL pslaset( 'Full', SIZE, SIZE, zero, zero, work( ptrr ), 1, 1,
320  $ descr )
321 *
322  DO 10 i = 1, SIZE
323  CALL pselset( work( ptrr ), i, i, descr, s( i ) )
324  10 CONTINUE
325 *
326 * Calculate U = U*R
327 *
328  DO 20 i = 1, SIZE
329  pcol = indxg2p( i, descu( nb_ ), 0, 0, npcol )
330  localcol = indxg2l( i, descu( nb_ ), 0, 0, npcol )
331  IF( mycol.EQ.pcol ) THEN
332  CALL sscal( mp, s( i ), u( ( localcol-1 )*descu( lld_ )+1 ),
333  $ 1 )
334  END IF
335  20 CONTINUE
336 *
337 * Calculate A = U*VT - A
338 *
339  CALL psgemm( 'N', 'N', m, n, SIZE, one, u, 1, 1, descu, vt, 1, 1,
340  $ descvt, mone, a, 1, 1, desca )
341 *
342  norma = pslange( '1', m, n, a, 1, 1, desca, work( ptrwork ) )
343  thresha = normai*max( m, n )*ulp*thresh
344 *
345  IF( norma.GT.thresha )
346  $ result( 1 ) = 1
347 *
348  IF( thresha.EQ.0 ) THEN
349  chk = 0.0e0
350  ELSE
351  chk = norma / thresha*thresh
352  END IF
353 *
354 * Test 4.
355 *
356  DO 30 i = 1, SIZE - 1
357  first = s( i )
358  second = s( i+1 )
359  IF( first.LT.second )
360  $ result( 4 ) = 1
361  30 CONTINUE
362  40 CONTINUE
363  RETURN
364  END
max
#define max(A, B)
Definition: pcgemr.c:180
pselset
subroutine pselset(A, IA, JA, DESCA, ALPHA)
Definition: pselset.f:2
descinit
subroutine descinit(DESC, M, N, MB, NB, IRSRC, ICSRC, ICTXT, LLD, INFO)
Definition: descinit.f:3
pssvdchk
subroutine pssvdchk(M, N, A, IA, JA, DESCA, U, IU, JU, DESCU, VT, IVT, JVT, DESCVT, S, THRESH, WORK, LWORK, RESULT, CHK, MTM)
Definition: pssvdchk.f:4
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
pslaset
subroutine pslaset(UPLO, M, N, ALPHA, BETA, A, IA, JA, DESCA)
Definition: psblastst.f:6863
min
#define min(A, B)
Definition: pcgemr.c:181