ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pdsvdchk.f
Go to the documentation of this file.
1  SUBROUTINE pdsvdchk( 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  DOUBLE PRECISION CHK, MTM, THRESH
13 * ..
14 * .. Array Arguments ..
15  INTEGER DESCA( * ), DESCU( * ), DESCVT( * ),
16  $ RESULT( * )
17  DOUBLE PRECISION 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 DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension (SIZE)
156 * Contains the computed singular values
157 *
158 * THRESH (input) DOUBLE PRECISION
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) DOUBLE PRECISION 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) DOUBLE PRECISION
192 * value of the | A - U*diag(S) VT | / ( |A| max(M,N) ulp )
193 *
194 * MTM (global output) DOUBLE PRECISION
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  DOUBLE PRECISION ZERO, ONE, MONE
207  PARAMETER ( ZERO = 0.0d+0, one = 1.0d+0, mone = -1.0d0 )
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  DOUBLE PRECISION 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  DOUBLE PRECISION PDLAMCH, PDLANGE
222  EXTERNAL indxg2l, indxg2p, numroc, pdlamch, pdlange
223 * ..
224 * .. External Subroutines ..
225  EXTERNAL blacs_gridinfo, chk1mat, descinit, dscal,
226  $ pdelset, pdgemm, pdlaset, pxerbla
227 * ..
228 * .. Intrinsic Functions ..
229  INTRINSIC max, min
230 * ..
231 * .. Executable Statements ..
232 * This is just to keep ftnchek happy
233  IF( block_cyclic_2d*csrc_*dtype_*m_*n_*rsrc_.LT.0 ) RETURN
234 *
235 * Test the input parameters.
236 *
237  CALL blacs_gridinfo( desca( ctxt_ ), nprow, npcol, myrow, mycol )
238  info = 0
239  SIZE = min( m, n )
240 *
241 * Sizepos is a number of parameters to pdsvdchk plus one. It's used
242 * for the error reporting.
243 *
244  sizepos = 22
245  IF( nprow.EQ.-1 ) THEN
246  info = -607
247  ELSE
248  CALL chk1mat( m, 1, n, 2, ia, ja, desca, 6, info )
249  CALL chk1mat( m, 1, SIZE, sizepos, iu, ju, descu, 10, info )
250  CALL chk1mat( SIZE, sizepos, n, 2, ivt, jvt, descvt, 14, info )
251  END IF
252 *
253  IF( info.EQ.0 ) THEN
254 *
255 * Calculate workspace
256 *
257  mp = numroc( m, desca( mb_ ), myrow, 0, nprow )
258  nq = numroc( n, desca( nb_ ), mycol, 0, npcol )
259  sizep = numroc( SIZE, descvt( mb_ ), myrow, 0, nprow )
260  sizeq = numroc( SIZE, descu( nb_ ), mycol, 0, npcol )
261  mx = max( sizeq, nq )
262  lwmin = 2 + sizeq*sizep + max( 2, mx )
263  work( 1 ) = lwmin
264  IF( lwork.EQ.-1 )
265  $ GO TO 40
266  IF( lwork.LT.lwmin ) THEN
267  info = -18
268  ELSE IF( thresh.LE.0 ) THEN
269  info = -16
270  END IF
271  END IF
272  IF( info.NE.0 ) THEN
273  CALL pxerbla( desca( ctxt_ ), 'PDSVDCHK', -info )
274  RETURN
275  END IF
276 *
277  ldr = max( 1, sizep )
278  ulp = pdlamch( desca( ctxt_ ), 'P' )
279  normai = pdlange( '1', m, n, a, ia, ja, desca, work )
280 *
281 * Allocate array R of global dimension SIZE x SIZE for testing
282 *
283  ptrr = 2
284  ptrwork = ptrr + sizeq*sizep
285 *
286  CALL descinit( descr, SIZE, SIZE, descvt( mb_ ), descu( nb_ ), 0,
287  $ 0, desca( ctxt_ ), ldr, info )
288 *
289 * Test 2. Form identity matrix R and make check norm(U'*U - I )
290 *
291  CALL pdlaset( 'Full', SIZE, SIZE, zero, one, work( ptrr ), 1, 1,
292  $ descr )
293  CALL pdgemm( 'T', 'N', SIZE, SIZE, m, one, u, iu, ju, descu, u,
294  $ iu, ju, descu, mone, work( ptrr ), 1, 1, descr )
295 *
296  normu = pdlange( '1', SIZE, SIZE, work( ptrr ), 1, 1, descr,
297  $ work( ptrwork ) )
298 *
299  normu = normu / ulp / SIZE / thresh
300  IF( normu.GT.1. )
301  $ result( 2 ) = 1
302 *
303 * Test3. Form identity matrix R and check norm(VT*VT' - I )
304 *
305  CALL pdlaset( 'Full', SIZE, SIZE, zero, one, work( ptrr ), 1, 1,
306  $ descr )
307  CALL pdgemm( 'N', 'T', SIZE, SIZE, n, one, vt, ivt, jvt, descvt,
308  $ vt, ivt, jvt, descvt, mone, work( ptrr ),
309  $ 1, 1, descr )
310  normvt = pdlange( '1', SIZE, SIZE, work( ptrr ), 1, 1, descr,
311  $ work( ptrwork ) )
312 *
313  normvt = normvt / ulp / SIZE / thresh
314  IF( normvt.GT.1. )
315  $ result( 3 ) = 1
316 *
317  mtm = max( normvt, normu )*thresh
318 *
319 * Test 1.
320 * Initialize R = diag( S )
321 *
322  CALL pdlaset( 'Full', SIZE, SIZE, zero, zero, work( ptrr ), 1, 1,
323  $ descr )
324 *
325  DO 10 i = 1, SIZE
326  CALL pdelset( work( ptrr ), i, i, descr, s( i ) )
327  10 CONTINUE
328 *
329 * Calculate U = U*R
330 *
331  DO 20 i = 1, SIZE
332  pcol = indxg2p( i, descu( nb_ ), 0, 0, npcol )
333  localcol = indxg2l( i, descu( nb_ ), 0, 0, npcol )
334  IF( mycol.EQ.pcol ) THEN
335  CALL dscal( mp, s( i ), u( ( localcol-1 )*descu( lld_ )+1 ),
336  $ 1 )
337  END IF
338  20 CONTINUE
339 *
340 * Calculate A = U*VT - A
341 *
342  CALL pdgemm( 'N', 'N', m, n, SIZE, one, u, iu, ju, descu, vt,
343  $ ivt, jvt, descvt, mone, a, ia, ja, desca )
344 *
345  norma = pdlange( '1', m, n, a, ia, ja, desca, work( ptrwork ) )
346  thresha = normai*max( m, n )*ulp*thresh
347 *
348  IF( norma.GT.thresha )
349  $ result( 1 ) = 1
350 *
351  IF( thresha.EQ.0 ) THEN
352  chk = 0.0d0
353  ELSE
354  chk = norma / thresha*thresh
355  END IF
356 *
357 * Test 4.
358 *
359  DO 30 i = 1, SIZE - 1
360  first = s( i )
361  second = s( i+1 )
362  IF( first.LT.second )
363  $ result( 4 ) = 1
364  30 CONTINUE
365  40 CONTINUE
366  RETURN
367  END
pdsvdchk
subroutine pdsvdchk(M, N, A, IA, JA, DESCA, U, IU, JU, DESCU, VT, IVT, JVT, DESCVT, S, THRESH, WORK, LWORK, RESULT, CHK, MTM)
Definition: pdsvdchk.f:4
max
#define max(A, B)
Definition: pcgemr.c:180
descinit
subroutine descinit(DESC, M, N, MB, NB, IRSRC, ICSRC, ICTXT, LLD, INFO)
Definition: descinit.f:3
pdlaset
subroutine pdlaset(UPLO, M, N, ALPHA, BETA, A, IA, JA, DESCA)
Definition: pdblastst.f:6862
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
pdelset
subroutine pdelset(A, IA, JA, DESCA, ALPHA)
Definition: pdelset.f:2
min
#define min(A, B)
Definition: pcgemr.c:181