ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pdsvdcmp.f
Go to the documentation of this file.
1  SUBROUTINE pdsvdcmp( M, N, JOBTYPE, S, SC, U, UC, IU, JU, DESCU,
2  $ VT, VTC, IVT, JVT, DESCVT, THRESH, RESULT,
3  $ DELTA, WORK, LWORK )
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 IU, IVT, JOBTYPE, JU, JVT, LWORK, M, N
12  DOUBLE PRECISION DELTA, THRESH
13 * ..
14 * .. Array Arguments ..
15  INTEGER DESCU( * ), DESCVT( * ), RESULT( * )
16  DOUBLE PRECISION S( * ), SC( * ), U( * ), UC( * ), VT( * ),
17  $ vtc( * ), work( * )
18 * ..
19 *
20 * Purpose
21 * ========
22 * Testing how accurately "full" and "partial" decomposition options
23 * provided by PDGESVD correspond to each other.
24 *
25 * Notes
26 * =====
27 *
28 * Each global data object is described by an associated description
29 * vector. This vector stores the information required to establish
30 * the mapping between an object element and its corresponding process
31 * and memory location.
32 *
33 * Let A be a generic term for any 2D block cyclicly distributed array.
34 * Such a global array has an associated description vector DESCA.
35 * In the following comments, the character _ should be read as
36 * "of the global array".
37 *
38 * NOTATION STORED IN EXPLANATION
39 * --------------- -------------- --------------------------------------
40 * DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
41 * DTYPE_A = 1.
42 * CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
43 * the BLACS process grid A is distribu-
44 * ted over. The context itself is glo-
45 * bal, but the handle (the integer
46 * value) may vary.
47 * M_A (global) DESCA( M_ ) The number of rows in the global
48 * array A.
49 * N_A (global) DESCA( N_ ) The number of columns in the global
50 * array A.
51 * MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
52 * the rows of the array.
53 * NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
54 * the columns of the array.
55 * RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
56 * row of the array A is distributed.
57 * CSRC_A (global) DESCA( CSRC_ ) The process column over which the
58 * first column of the array A is
59 * distributed.
60 * LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
61 * array. LLD_A >= MAX(1,LOCr(M_A)).
62 *
63 * Let K be the number of rows or columns of a distributed matrix,
64 * and assume that its process grid has dimension p x q.
65 * LOCr( K ) denotes the number of elements of K that a process
66 * would receive if K were distributed over the p processes of its
67 * process column.
68 * Similarly, LOCc( K ) denotes the number of elements of K that a
69 * process would receive if K were distributed over the q processes of
70 * its process row.
71 * The values of LOCr() and LOCc() may be determined via a call to the
72 * ScaLAPACK tool function, NUMROC:
73 * LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
74 * LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
75 * An upper bound for these quantities may be computed by:
76 * LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
77 * LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
78 *
79 * Arguments
80 * ==========
81 *
82 * M (global input) INTEGER
83 * Number of rows of the distributed matrix, for which
84 * SVD was calculated
85 *
86 * N (global input) INTEGER
87 * Number of columns of the distributed matrix, for which
88 * SVD was calculated
89 *
90 * JOBTYPE (global input) INTEGER
91 * Depending on the value of this parameter,
92 * the following comparisons are performed:
93 *
94 * JOBTYPE | COMPARISON
95 * -------------------------------------------
96 * 2 | | U - UC | / ( M ulp ) > THRESH,
97 * 3 | | VT - VTC | / ( N ulp ) > THRESH
98 *
99 * In addition, for JOBTYPE = 2:4 comparison
100 * | S1 - S2 | / ( SIZE ulp |S| ) > THRESH
101 * is performed. Positive result of any of the comparisons
102 * typically indicates erroneous computations and sets
103 * to one corresponding element of array RESULT
104 *
105 * S (global input) DOUBLE PRECISION array of singular values
106 * calculated for JOBTYPE equal to 1
107 *
108 * SC (global input) DOUBLE PRECISION array of singular values
109 * calculated for JOBTYPE nonequal to 1
110 *
111 * U (local input) DOUBLE PRECISION array of left singular
112 * vectors calculated for JOBTYPE equal to 1, local
113 * dimension (MP, SIZEQ), global dimension (M, SIZE)
114 *
115 * UC (local input) DOUBLE PRECISION array of left singular
116 * vectors calculated for JOBTYPE non equal to 1, local
117 * dimension (MP, SIZEQ), global dimension (M, SIZE)
118 *
119 * IU (global input) INTEGER
120 * The row index in the global array U indicating the first
121 * row of sub( U ).
122 *
123 * JU (global input) INTEGER
124 * The column index in the global array U indicating the
125 * first column of sub( U ).
126 *
127 * DESCU (global input) INTEGER array of dimension DLEN_
128 * The array descriptor for the distributed matrix U and UC
129 *
130 * V (local input) DOUBLE PRECISION array of right singular
131 * vectors calculated for JOBTYPE equal to 1, local
132 * dimension (SIZEP, NQ), global dimension (SIZE, N)
133 *
134 * VC (local input) DOUBLE PRECISION array of right singular
135 * vectors calculated for JOBTYPE non equal to 1, local
136 * dimension (SIZEP, NQ), global dimension (SIZE, N)
137 *
138 * IVT (global input) INTEGER
139 * The row index in the global array VT indicating the first
140 * row of sub( VT ).
141 *
142 * JVT (global input) INTEGER
143 * The column index in the global array VT indicating the
144 * first column of sub( VT ).
145 *
146 * DESCVT (global input) INTEGER array of dimension DLEN_
147 * The array descriptor for the distributed matrix VT and
148 * VTC
149 *
150 * THRESH (global input) DOUBLE PRECISION
151 * The threshold value for the test ratios. A result is
152 * included in the output file if RESULT >= THRESH. The test
153 * ratios are scaled to be O(1), so THRESH should be a small
154 * multiple of 1, e.g., 10 or 100. To have every test ratio
155 * printed, use THRESH = 0.
156 *
157 * RESULT (global input/output) INTEGER array.
158 * Every nonzero entry corresponds to erroneous computation.
159 *
160 * DELTA (global output) DOUBLE PRECISION
161 * maximum of the available of the following three values
162 * | U - UC | / ( M ulp THRESH ),
163 * | VT - VT | / ( N ulp THRESH ),
164 * | S1 - S2 | / ( SIZE ulp |S| THRESH )
165 *
166 * WORK (local workspace/output) DOUBLE PRECISION array,
167 * dimension (LWORK)
168 * On exit, WORK(1) returns the optimal LWORK.
169 *
170 * LWORK (local input) INTEGER
171 * The dimension of the array WORK.
172 *
173 * ======================================================================
174 *
175 * .. Parameters ..
176  INTEGER BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_,
177  $ MB_, NB_, RSRC_, CSRC_, LLD_
178  parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
179  $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
180  $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
181 * ..
182 * .. Local Scalars ..
183  INTEGER COLPTR, I, INFO, J, LWMIN, MYCOL, MYROW, NPCOL,
184  $ NPROW, NQ, RESULTS, SIZE, SIZEPOS, SIZEQ
185  DOUBLE PRECISION ACCUR, CMP, NORMDIFS, NORMDIFU, NORMDIFV,
186  $ NORMS, ULP
187 * ..
188 * .. External Functions ..
189  INTEGER NUMROC
190  DOUBLE PRECISION DLANGE, PDLAMCH, PDLANGE
191  EXTERNAL numroc, dlange, pdlamch, pdlange
192 * ..
193 * .. External Subroutines ..
194  EXTERNAL blacs_gridinfo, chk1mat, pxerbla
195 * ..
196 * .. Intrinsic Functions ..
197  INTRINSIC max, min
198 * ..
199 * .. Executable Statements ..
200 * This is just to keep ftnchek happy
201  IF( block_cyclic_2d*csrc_*dlen_*dtype_*mb_*m_*n_*rsrc_.LT.0 )
202  $ RETURN
203 *
204  results = 0
205  normdifs = 0
206  normdifu = 0
207  normdifv = 0
208  SIZE = min( m, n )
209 *
210 * Sizepos is a number of parameters to pdsvdcmp plus one. It's used
211 * for the error reporting.
212 *
213  sizepos = 17
214  info = 0
215  CALL blacs_gridinfo( descu( ctxt_ ), nprow, npcol, myrow, mycol )
216  IF( nprow.EQ.-1 ) THEN
217  info = -607
218  ELSE
219  CALL chk1mat( m, 1, SIZE, sizepos, 1, 1, descu, 8, info )
220  CALL chk1mat( SIZE, sizepos, n, 2, 1, 1, descvt, 11, info )
221  END IF
222 *
223  IF( info.EQ.0 ) THEN
224 *
225 * Calculate workspace.
226 *
227  sizeq = numroc( SIZE, descu( nb_ ), mycol, 0, npcol )
228  nq = numroc( n, descvt( nb_ ), mycol, 0, npcol )
229  lwmin = max( sizeq, nq ) + 4
230  work( 1 ) = lwmin
231  IF( lwork.EQ.-1 )
232  $ GO TO 60
233  IF( lwork.LT.lwmin ) THEN
234  info = -16
235  ELSE IF( thresh.LE.0 ) THEN
236  info = -12
237  END IF
238  END IF
239 *
240  IF( info.NE.0 ) THEN
241  CALL pxerbla( descu( ctxt_ ), 'PDSVDCMP', -info )
242  RETURN
243  END IF
244 *
245  ulp = pdlamch( descu( ctxt_ ), 'P' )
246 *
247 * Make comparison of singular values.
248 *
249  norms = dlange( '1', SIZE, 1, s, SIZE, work )
250  DO 10 i = 1, SIZE
251  sc( i ) = s( i ) - sc( i )
252  10 CONTINUE
253 *
254  normdifs = dlange( '1', SIZE, 1, sc, SIZE, work )
255  accur = ulp*size*norms*thresh
256 *
257  IF( normdifs.GT.accur )
258  $ results = 1
259  IF( normdifs.EQ.0 .AND. accur.EQ.0 ) THEN
260  normdifs = 0
261  ELSE
262  normdifs = normdifs / accur
263  END IF
264 *
265  IF( jobtype.EQ.2 ) THEN
266 *
267  result( 5 ) = results
268  accur = ulp*m*thresh
269  DO 30 j = 1, sizeq
270  colptr = descu( lld_ )*( j-1 )
271  DO 20 i = 1, descu( lld_ )
272  uc( i+colptr ) = u( i+colptr ) - uc( i+colptr )
273  20 CONTINUE
274  30 CONTINUE
275 *
276  normdifu = pdlange( '1', m, SIZE, uc, iu, ju, descu, work )
277 *
278  IF( normdifu.GE.accur )
279  $ result( 6 ) = 1
280  IF( normdifu.EQ.0 .AND. accur.EQ.0 ) THEN
281  normdifu = 0
282  ELSE
283  normdifu = normdifu / accur
284  END IF
285 *
286  ELSE IF( jobtype.EQ.3 ) THEN
287 *
288  result( 7 ) = results
289  accur = ulp*n*thresh
290  DO 50 j = 1, nq
291  colptr = descvt( lld_ )*( j-1 )
292  DO 40 i = 1, descvt( lld_ )
293  vtc( i+colptr ) = vt( i+colptr ) - vtc( i+colptr )
294  40 CONTINUE
295  50 CONTINUE
296 *
297  normdifv = pdlange( '1', SIZE, n, vtc, ivt, jvt, descvt, work )
298 *
299  IF( normdifv.GE.accur )
300  $ result( 8 ) = 1
301 *
302  IF( normdifv.EQ.0 .AND. accur.EQ.0 ) THEN
303  normdifv = 0
304  ELSE
305  normdifv = normdifv / accur
306  END IF
307 *
308  ELSE IF( jobtype.EQ.4 ) THEN
309 *
310  result( 9 ) = results
311 *
312  END IF
313 *
314  cmp = max( normdifv, normdifu )
315  delta = max( cmp, normdifs )
316 *
317  60 CONTINUE
318 *
319 * End of PDSVDCMP
320 *
321  RETURN
322  END
max
#define max(A, B)
Definition: pcgemr.c:180
pdsvdcmp
subroutine pdsvdcmp(M, N, JOBTYPE, S, SC, U, UC, IU, JU, DESCU, VT, VTC, IVT, JVT, DESCVT, THRESH, RESULT, DELTA, WORK, LWORK)
Definition: pdsvdcmp.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
min
#define min(A, B)
Definition: pcgemr.c:181