ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pcpocon.f
Go to the documentation of this file.
1  SUBROUTINE pcpocon( UPLO, N, A, IA, JA, DESCA, ANORM, RCOND, WORK,
2  $ LWORK, RWORK, LRWORK, INFO )
3 *
4 * -- ScaLAPACK routine (version 1.7) --
5 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
6 * and University of California, Berkeley.
7 * May 25, 2001
8 *
9 * .. Scalar Arguments ..
10  CHARACTER UPLO
11  INTEGER IA, INFO, JA, LRWORK, LWORK, N
12  REAL ANORM, RCOND
13 * ..
14 * .. Array Arguments ..
15  INTEGER DESCA( * )
16  REAL RWORK( * )
17  COMPLEX A( * ), WORK( * )
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 * PCPOCON estimates the reciprocal of the condition number (in the
24 * 1-norm) of a complex Hermitian positive definite distributed matrix
25 * using the Cholesky factorization A = U**H*U or A = L*L**H computed by
26 * PCPOTRF.
27 *
28 * An estimate is obtained for norm(inv(A(IA:IA+N-1,JA:JA+N-1))), and
29 * the reciprocal of the condition number is computed as
30 * RCOND = 1 / ( norm( A(IA:IA+N-1,JA:JA+N-1) ) *
31 * norm( inv(A(IA:IA+N-1,JA:JA+N-1)) ) ).
32 *
33 * Notes
34 * =====
35 *
36 * Each global data object is described by an associated description
37 * vector. This vector stores the information required to establish
38 * the mapping between an object element and its corresponding process
39 * and memory location.
40 *
41 * Let A be a generic term for any 2D block cyclicly distributed array.
42 * Such a global array has an associated description vector DESCA.
43 * In the following comments, the character _ should be read as
44 * "of the global array".
45 *
46 * NOTATION STORED IN EXPLANATION
47 * --------------- -------------- --------------------------------------
48 * DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
49 * DTYPE_A = 1.
50 * CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
51 * the BLACS process grid A is distribu-
52 * ted over. The context itself is glo-
53 * bal, but the handle (the integer
54 * value) may vary.
55 * M_A (global) DESCA( M_ ) The number of rows in the global
56 * array A.
57 * N_A (global) DESCA( N_ ) The number of columns in the global
58 * array A.
59 * MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
60 * the rows of the array.
61 * NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
62 * the columns of the array.
63 * RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
64 * row of the array A is distributed.
65 * CSRC_A (global) DESCA( CSRC_ ) The process column over which the
66 * first column of the array A is
67 * distributed.
68 * LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
69 * array. LLD_A >= MAX(1,LOCr(M_A)).
70 *
71 * Let K be the number of rows or columns of a distributed matrix,
72 * and assume that its process grid has dimension p x q.
73 * LOCr( K ) denotes the number of elements of K that a process
74 * would receive if K were distributed over the p processes of its
75 * process column.
76 * Similarly, LOCc( K ) denotes the number of elements of K that a
77 * process would receive if K were distributed over the q processes of
78 * its process row.
79 * The values of LOCr() and LOCc() may be determined via a call to the
80 * ScaLAPACK tool function, NUMROC:
81 * LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
82 * LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
83 * An upper bound for these quantities may be computed by:
84 * LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
85 * LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
86 *
87 * Arguments
88 * =========
89 *
90 * UPLO (global input) CHARACTER
91 * Specifies whether the factor stored in
92 * A(IA:IA+N-1,JA:JA+N-1) is upper or lower triangular.
93 * = 'U': Upper triangular
94 * = 'L': Lower triangular
95 *
96 * N (global input) INTEGER
97 * The order of the distributed matrix A(IA:IA+N-1,JA:JA+N-1).
98 * N >= 0.
99 *
100 * A (local input) COMPLEX pointer into the local memory to
101 * an array of dimension ( LLD_A, LOCc(JA+N-1) ). On entry, this
102 * array contains the local pieces of the factors L or U from
103 * the Cholesky factorization A(IA:IA+N-1,JA:JA+N-1) = U'*U or
104 * L*L', as computed by PCPOTRF.
105 *
106 * IA (global input) INTEGER
107 * The row index in the global array A indicating the first
108 * row of sub( A ).
109 *
110 * JA (global input) INTEGER
111 * The column index in the global array A indicating the
112 * first column of sub( A ).
113 *
114 * DESCA (global and local input) INTEGER array of dimension DLEN_.
115 * The array descriptor for the distributed matrix A.
116 *
117 * ANORM (global input) REAL
118 * The 1-norm (or infinity-norm) of the hermitian distributed
119 * matrix A(IA:IA+N-1,JA:JA+N-1).
120 *
121 * RCOND (global output) REAL
122 * The reciprocal of the condition number of the distributed
123 * matrix A(IA:IA+N-1,JA:JA+N-1), computed as
124 * RCOND = 1 / ( norm( A(IA:IA+N-1,JA:JA+N-1) ) *
125 * norm( inv(A(IA:IA+N-1,JA:JA+N-1)) ) ).
126 *
127 * WORK (local workspace/local output) COMPLEX array,
128 * dimension (LWORK)
129 * On exit, WORK(1) returns the minimal and optimal LWORK.
130 *
131 * LWORK (local or global input) INTEGER
132 * The dimension of the array WORK.
133 * LWORK is local input and must be at least
134 * LWORK >= 2*LOCr(N+MOD(IA-1,MB_A)) +
135 * MAX( 2, MAX(NB_A*MAX(1,CEIL(P-1,Q)),LOCc(N+MOD(JA-1,NB_A)) +
136 * NB_A*MAX(1,CEIL(Q-1,P))) ).
137 *
138 * If LWORK = -1, then LWORK is global input and a workspace
139 * query is assumed; the routine only calculates the minimum
140 * and optimal size for all work arrays. Each of these
141 * values is returned in the first entry of the corresponding
142 * work array, and no error message is issued by PXERBLA.
143 *
144 * RWORK (local workspace/local output) REAL array,
145 * dimension (LRWORK)
146 * On exit, RWORK(1) returns the minimal and optimal LRWORK.
147 *
148 * LRWORK (local or global input) INTEGER
149 * The dimension of the array RWORK.
150 * LRWORK is local input and must be at least
151 * LRWORK >= 2*LOCc(N+MOD(JA-1,NB_A)).
152 *
153 * If LRWORK = -1, then LRWORK is global input and a workspace
154 * query is assumed; the routine only calculates the minimum
155 * and optimal size for all work arrays. Each of these
156 * values is returned in the first entry of the corresponding
157 * work array, and no error message is issued by PXERBLA.
158 *
159 *
160 * INFO (global output) INTEGER
161 * = 0: successful exit
162 * < 0: If the i-th argument is an array and the j-entry had
163 * an illegal value, then INFO = -(i*100+j), if the i-th
164 * argument is a scalar and had an illegal value, then
165 * INFO = -i.
166 *
167 * =====================================================================
168 *
169 * .. Parameters ..
170  INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
171  $ lld_, mb_, m_, nb_, n_, rsrc_
172  parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
173  $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
174  $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
175  REAL ONE, ZERO
176  parameter( one = 1.0e+0, zero = 0.0e+0 )
177 * ..
178 * .. Local Scalars ..
179  LOGICAL LQUERY, UPPER
180  CHARACTER CBTOP, COLCTOP, NORMIN, ROWCTOP
181  INTEGER IACOL, IAROW, ICOFF, ICTXT, IIA, IPNL, IPNU,
182  $ ipv, ipw, ipx, iroff, iv, ix, ixx, jja, jv,
183  $ jx, kase, lrwmin, lwmin, mycol, myrow, np,
184  $ npcol, nprow, npmod, nq, nqmod
185  REAL AINVNM, SCALE, SL, SU, SMLNUM
186  COMPLEX WMAX, ZDUM
187 * ..
188 * .. Local Arrays ..
189  INTEGER DESCV( DLEN_ ), DESCX( DLEN_ ), IDUM1( 3 ),
190  $ idum2( 3 )
191 * ..
192 * .. External Subroutines ..
193  EXTERNAL blacs_gridinfo, cgebr2d, cgebs2d, chk1mat,
194  $ descset, infog2l, pcamax, pchk1mat, pclatrs,
195  $ pclacon, pcsrscl, pb_topget, pb_topset,
196  $ pxerbla
197 * ..
198 * .. External Functions ..
199  LOGICAL LSAME
200  INTEGER ICEIL, INDXG2P, NUMROC
201  REAL PSLAMCH
202  EXTERNAL iceil, indxg2p, lsame, numroc, pslamch
203 * ..
204 * .. Intrinsic Functions ..
205  INTRINSIC abs, aimag, ichar, max, mod, real
206 * ..
207 * .. Statement Functions ..
208  REAL CABS1
209 * ..
210 * .. Statement Function definitions ..
211  cabs1( zdum ) = abs( real( zdum ) ) + abs( aimag( zdum ) )
212 * ..
213 * .. Executable Statements ..
214 *
215 * Get grid parameters
216 *
217  ictxt = desca( ctxt_ )
218  CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
219 *
220 * Test the input parameters
221 *
222  info = 0
223  IF( nprow.EQ.-1 ) THEN
224  info = -(600+ctxt_)
225  ELSE
226  CALL chk1mat( n, 2, n, 2, ia, ja, desca, 6, info )
227  IF( info.EQ.0 ) THEN
228  upper = lsame( uplo, 'U' )
229  iarow = indxg2p( ia, desca( mb_ ), myrow, desca( rsrc_ ),
230  $ nprow )
231  iacol = indxg2p( ja, desca( nb_ ), mycol, desca( csrc_ ),
232  $ npcol )
233  npmod = numroc( n + mod( ia-1, desca( mb_ ) ), desca( mb_ ),
234  $ myrow, iarow, nprow )
235  nqmod = numroc( n + mod( ja-1, desca( nb_ ) ), desca( nb_ ),
236  $ mycol, iacol, npcol )
237  lwmin = 2*npmod +
238  $ max( 2, max( desca( nb_ )*
239  $ max( 1, iceil( nprow-1, npcol ) ), nqmod +
240  $ desca( nb_ )*
241  $ max( 1, iceil( npcol-1, nprow ) ) ) )
242  work( 1 ) = real( lwmin )
243  lrwmin = 2*nqmod
244  rwork( 1 ) = real( lrwmin )
245  lquery = ( lwork.EQ.-1 .OR. lrwork.EQ.-1 )
246 *
247  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
248  info = -1
249  ELSE IF( anorm.LT.zero ) THEN
250  info = -7
251  ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
252  info = -10
253  ELSE IF( lrwork.LT.lrwmin .AND. .NOT.lquery ) THEN
254  info = -12
255  END IF
256  END IF
257 *
258  IF( upper ) THEN
259  idum1( 1 ) = ichar( 'U' )
260  ELSE
261  idum1( 1 ) = ichar( 'L' )
262  END IF
263  idum2( 1 ) = 1
264  IF( lwork.EQ.-1 ) THEN
265  idum1( 2 ) = -1
266  ELSE
267  idum1( 2 ) = 1
268  END IF
269  idum2( 2 ) = 10
270  IF( lrwork.EQ.-1 ) THEN
271  idum1( 3 ) = -1
272  ELSE
273  idum1( 3 ) = 1
274  END IF
275  idum2( 3 ) = 12
276  CALL pchk1mat( n, 2, n, 2, ia, ja, desca, 6, 3, idum1, idum2,
277  $ info )
278  END IF
279 *
280  IF( info.NE.0 ) THEN
281  CALL pxerbla( ictxt, 'PCPOCON', -info )
282  RETURN
283  ELSE IF( lquery ) THEN
284  RETURN
285  END IF
286 *
287 * Quick return if possible
288 *
289  rcond = zero
290  IF( n.EQ.0 ) THEN
291  rcond = one
292  RETURN
293  ELSE IF( anorm.EQ.zero ) THEN
294  RETURN
295  ELSE IF( n.EQ.1 ) THEN
296  rcond = one
297  RETURN
298  END IF
299 *
300  CALL pb_topget( ictxt, 'Combine', 'Columnwise', colctop )
301  CALL pb_topget( ictxt, 'Combine', 'Rowwise', rowctop )
302  CALL pb_topset( ictxt, 'Combine', 'Columnwise', '1-tree' )
303  CALL pb_topset( ictxt, 'Combine', 'Rowwise', '1-tree' )
304 *
305  smlnum = pslamch( ictxt, 'Safe minimum' )
306  iroff = mod( ia-1, desca( mb_ ) )
307  icoff = mod( ja-1, desca( nb_ ) )
308  CALL infog2l( ia, ja, desca, nprow, npcol, myrow, mycol, iia, jja,
309  $ iarow, iacol )
310  np = numroc( n+iroff, desca( mb_ ), myrow, iarow, nprow )
311  nq = numroc( n+icoff, desca( nb_ ), mycol, iacol, npcol )
312  iv = iroff + 1
313  ix = iv
314  jv = icoff + 1
315  jx = jv
316 *
317  ipx = 1
318  ipv = ipx + np
319  ipw = ipv + np
320  ipnl = 1
321  ipnu = ipnl + nq
322 *
323  CALL descset( descv, n+iroff, 1, desca( mb_ ), 1, iarow, mycol,
324  $ ictxt, max( 1, np ) )
325  CALL descset( descx, n+iroff, 1, desca( mb_ ), 1, iarow, mycol,
326  $ ictxt, max( 1, np ) )
327 *
328 * Estimate the 1-norm (or I-norm) of inv(A).
329 *
330  ainvnm = zero
331  kase = 0
332  normin = 'N'
333 *
334  10 CONTINUE
335  CALL pclacon( n, work( ipv ), iv, jv, descv, work( ipx ), ix, jx,
336  $ descx, ainvnm, kase )
337  IF( kase.NE.0 ) THEN
338  IF( upper ) THEN
339 *
340 * Multiply by inv(U').
341 *
342  descx( csrc_ ) = iacol
343  CALL pclatrs( 'Upper', 'Conjugate transpose', 'Non-unit',
344  $ normin, n, a, ia, ja, desca, work( ipx ),
345  $ ix, jx, descx, sl, rwork( ipnl ),
346  $ work( ipw ) )
347  descx( csrc_ ) = mycol
348  normin = 'Y'
349 *
350 * Multiply by inv(U).
351 *
352  descx( csrc_ ) = iacol
353  CALL pclatrs( 'Upper', 'No transpose', 'Non-unit', normin,
354  $ n, a, ia, ja, desca, work( ipx ), ix, jx,
355  $ descx, su, rwork( ipnu ), work( ipw ) )
356  descx( csrc_ ) = mycol
357  ELSE
358 *
359 * Multiply by inv(L).
360 *
361  descx( csrc_ ) = iacol
362  CALL pclatrs( 'Lower', 'No transpose', 'Non-unit', normin,
363  $ n, a, ia, ja, desca, work( ipx ), ix, jx,
364  $ descx, sl, rwork( ipnl ), work( ipw ) )
365  descx( csrc_ ) = mycol
366  normin = 'Y'
367 *
368 * Multiply by inv(L').
369 *
370  descx( csrc_ ) = iacol
371  CALL pclatrs( 'Lower', 'Conjugate transpose', 'Non-unit',
372  $ normin, n, a, ia, ja, desca, work( ipx ),
373  $ ix, jx, descx, su, rwork( ipnu ),
374  $ work( ipw ) )
375  descx( csrc_ ) = mycol
376  END IF
377 *
378 * Multiply by 1/SCALE if doing so will not cause overflow.
379 *
380  scale = sl*su
381  IF( scale.NE.one ) THEN
382  CALL pcamax( n, wmax, ixx, work( ipx ), ix, jx, descx, 1 )
383  IF( descx( m_ ).EQ.1 .AND. n.EQ.1 ) THEN
384  CALL pb_topget( ictxt, 'Broadcast', 'Columnwise', cbtop )
385  IF( myrow.EQ.iarow ) THEN
386  CALL cgebs2d( ictxt, 'Column', cbtop, 1, 1, wmax, 1 )
387  ELSE
388  CALL cgebr2d( ictxt, 'Column', cbtop, 1, 1, wmax, 1,
389  $ iarow, mycol )
390  END IF
391  END IF
392  IF( scale.LT.cabs1( wmax )*smlnum .OR. scale.EQ.zero )
393  $ GO TO 20
394  CALL pcsrscl( n, scale, work( ipx ), ix, jx, descx, 1 )
395  END IF
396  GO TO 10
397  END IF
398 *
399 * Compute the estimate of the reciprocal condition number.
400 *
401  IF( ainvnm.NE.zero )
402  $ rcond = ( one / ainvnm ) / anorm
403 *
404  20 CONTINUE
405 *
406  CALL pb_topset( ictxt, 'Combine', 'Columnwise', colctop )
407  CALL pb_topset( ictxt, 'Combine', 'Rowwise', rowctop )
408 *
409  RETURN
410 *
411 * End of PCPOCON
412 *
413  END
max
#define max(A, B)
Definition: pcgemr.c:180
pclacon
subroutine pclacon(N, V, IV, JV, DESCV, X, IX, JX, DESCX, EST, KASE)
Definition: pclacon.f:3
pclatrs
subroutine pclatrs(UPLO, TRANS, DIAG, NORMIN, N, A, IA, JA, DESCA, X, IX, JX, DESCX, SCALE, CNORM, WORK)
Definition: pclatrs.f:4
infog2l
subroutine infog2l(GRINDX, GCINDX, DESC, NPROW, NPCOL, MYROW, MYCOL, LRINDX, LCINDX, RSRC, CSRC)
Definition: infog2l.f:3
pcpocon
subroutine pcpocon(UPLO, N, A, IA, JA, DESCA, ANORM, RCOND, WORK, LWORK, RWORK, LRWORK, INFO)
Definition: pcpocon.f:3
pchk1mat
subroutine pchk1mat(MA, MAPOS0, NA, NAPOS0, IA, JA, DESCA, DESCAPOS0, NEXTRA, EX, EXPOS, INFO)
Definition: pchkxmat.f:3
pcsrscl
subroutine pcsrscl(N, SA, SX, IX, JX, DESCX, INCX)
Definition: pcsrscl.f:2
descset
subroutine descset(DESC, M, N, MB, NB, IRSRC, ICSRC, ICTXT, LLD)
Definition: descset.f:3
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