SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
pcgecon.f
Go to the documentation of this file.
1 SUBROUTINE pcgecon( NORM, 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 NORM
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* PCGECON estimates the reciprocal of the condition number of a general
24* distributed complex matrix A(IA:IA+N-1,JA:JA+N-1), in either the
25* 1-norm or the infinity-norm, using the LU factorization computed by
26* PCGETRF.
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* NORM (global input) CHARACTER
91* Specifies whether the 1-norm condition number or the
92* infinity-norm condition number is required:
93* = '1' or 'O': 1-norm
94* = 'I': Infinity-norm
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
101* to an array of dimension ( LLD_A, LOCc(JA+N-1) ). On entry,
102* this array contains the local pieces of the factors L and U
103* from the factorization A(IA:IA+N-1,JA:JA+N-1) = P*L*U; the
104* unit diagonal elements of L are not stored.
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* If NORM = '1' or 'O', the 1-norm of the original distributed
119* matrix A(IA:IA+N-1,JA:JA+N-1).
120* If NORM = 'I', the infinity-norm of the original distributed
121* matrix A(IA:IA+N-1,JA:JA+N-1).
122*
123* RCOND (global output) REAL
124* The reciprocal of the condition number of the distributed
125* matrix A(IA:IA+N-1,JA:JA+N-1), computed as
126* RCOND = 1 / ( norm( A(IA:IA+N-1,JA:JA+N-1) ) *
127* norm( inv(A(IA:IA+N-1,JA:JA+N-1)) ) ).
128*
129* WORK (local workspace/local output) COMPLEX array,
130* dimension (LWORK)
131* On exit, WORK(1) returns the minimal and optimal LWORK.
132*
133* LWORK (local or global input) INTEGER
134* The dimension of the array WORK.
135* LWORK is local input and must be at least
136* LWORK >= 2*LOCr(N+MOD(IA-1,MB_A)) +
137* MAX( 2, MAX(NB_A*CEIL(NPROW-1,NPCOL),LOCc(N+MOD(JA-1,NB_A)) +
138* NB_A*CEIL(NPCOL-1,NPROW)) ).
139*
140* LOCr and LOCc values can be computed using the ScaLAPACK
141* tool function NUMROC; NPROW and NPCOL can be determined by
142* calling the subroutine BLACS_GRIDINFO.
143*
144* If LWORK = -1, then LWORK is global input and a workspace
145* query is assumed; the routine only calculates the minimum
146* and optimal size for all work arrays. Each of these
147* values is returned in the first entry of the corresponding
148* work array, and no error message is issued by PXERBLA.
149*
150* RWORK (local workspace/local output) REAL array,
151* dimension (LRWORK)
152* On exit, RWORK(1) returns the minimal and optimal LRWORK.
153*
154* LRWORK (local or global input) INTEGER
155* The dimension of the array RWORK.
156* LRWORK is local input and must be at least
157* LRWORK >= MAX( 1, 2*LOCc(N+MOD(JA-1,NB_A)) ).
158*
159* If LRWORK = -1, then LRWORK is global input and a workspace
160* query is assumed; the routine only calculates the minimum
161* and optimal size for all work arrays. Each of these
162* values is returned in the first entry of the corresponding
163* work array, and no error message is issued by PXERBLA.
164*
165*
166* INFO (global output) INTEGER
167* = 0: successful exit
168* < 0: If the i-th argument is an array and the j-entry had
169* an illegal value, then INFO = -(i*100+j), if the i-th
170* argument is a scalar and had an illegal value, then
171* INFO = -i.
172*
173* =====================================================================
174*
175* .. Parameters ..
176 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
177 $ lld_, mb_, m_, nb_, n_, rsrc_
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 REAL ONE, ZERO
182 parameter( one = 1.0e+0, zero = 0.0e+0 )
183* ..
184* .. Local Scalars ..
185 LOGICAL LQUERY, ONENRM
186 CHARACTER CBTOP, COLCTOP, NORMIN, ROWCTOP
187 INTEGER IACOL, IAROW, ICOFF, ICTXT, IIA, IPNL, IPNU,
188 $ ipv, ipw, ipx, iroff, iv, ix, ixx, jja, jv, jx,
189 $ kase, kase1, lrwmin, lwmin, mycol, myrow, np,
190 $ npcol, npmod, nprow, nq, nqmod
191 REAL AINVNM, SCALE, SL, SMLNUM, SU
192 COMPLEX WMAX, ZDUM
193* ..
194* .. Local Arrays ..
195 INTEGER DESCV( DLEN_ ), DESCX( DLEN_ ), IDUM1( 3 ),
196 $ idum2( 3 )
197* ..
198* .. External Subroutines ..
199 EXTERNAL blacs_gridinfo, chk1mat, descset, cgebr2d,
200 $ cgebs2d, infog2l, pcamax, pchk1mat,
201 $ pclatrs, pclacon, pcsrscl, pb_topget,
202 $ pb_topset, pxerbla
203* ..
204* .. External Functions ..
205 LOGICAL LSAME
206 INTEGER ICEIL, INDXG2P, NUMROC
207 REAL PSLAMCH
208 EXTERNAL iceil, indxg2p, lsame, numroc, pslamch
209* ..
210* .. Intrinsic Functions ..
211 INTRINSIC abs, aimag, ichar, max, mod, real
212* ..
213* .. Statement Functions ..
214 REAL CABS1
215* ..
216* .. Statement Function definitions ..
217 cabs1( zdum ) = abs( real( zdum ) ) + abs( aimag( zdum ) )
218* ..
219* .. Executable Statements ..
220*
221* Get grid parameters
222*
223 ictxt = desca( ctxt_ )
224 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
225*
226* Test the input parameters
227*
228 info = 0
229 IF( nprow.EQ.-1 ) THEN
230 info = -( 600 + ctxt_ )
231 ELSE
232 CALL chk1mat( n, 2, n, 2, ia, ja, desca, 6, info )
233 IF( info.EQ.0 ) THEN
234 onenrm = norm.EQ.'1' .OR. lsame( norm, 'O' )
235 iarow = indxg2p( ia, desca( mb_ ), myrow, desca( rsrc_ ),
236 $ nprow )
237 iacol = indxg2p( ja, desca( nb_ ), mycol, desca( csrc_ ),
238 $ npcol )
239 npmod = numroc( n + mod( ia-1, desca( mb_ ) ), desca( mb_ ),
240 $ myrow, iarow, nprow )
241 nqmod = numroc( n + mod( ja-1, desca( nb_ ) ), desca( nb_ ),
242 $ mycol, iacol, npcol )
243 lwmin = 2*npmod +
244 $ max( 2, max( desca( nb_ )*
245 $ max( 1, iceil( nprow-1, npcol ) ), nqmod +
246 $ desca( nb_ )*
247 $ max( 1, iceil( npcol-1, nprow ) ) ) )
248 work( 1 ) = real( lwmin )
249 lrwmin = max( 1, 2*nqmod )
250 rwork( 1 ) = real( lrwmin )
251 lquery = ( lwork.EQ.-1 .OR. lrwork.EQ.-1 )
252*
253 IF( .NOT.onenrm .AND. .NOT.lsame( norm, 'I' ) ) THEN
254 info = -1
255 ELSE IF( anorm.LT.zero ) THEN
256 info = -7
257 ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
258 info = -10
259 ELSE IF( lrwork.LT.lrwmin .AND. .NOT.lquery ) THEN
260 info = -12
261 END IF
262 END IF
263*
264 IF( onenrm ) THEN
265 idum1( 1 ) = ichar( '1' )
266 ELSE
267 idum1( 1 ) = ichar( 'I' )
268 END IF
269 idum2( 1 ) = 1
270 IF( lwork.EQ.-1 ) THEN
271 idum1( 2 ) = -1
272 ELSE
273 idum1( 2 ) = 1
274 END IF
275 idum2( 2 ) = 10
276 IF( lrwork.EQ.-1 ) THEN
277 idum1( 3 ) = -1
278 ELSE
279 idum1( 3 ) = 1
280 END IF
281 idum2( 3 ) = 12
282 CALL pchk1mat( n, 2, n, 2, ia, ja, desca, 6, 3, idum1, idum2,
283 $ info )
284 END IF
285*
286 IF( info.NE.0 ) THEN
287 CALL pxerbla( ictxt, 'PCGECON', -info )
288 RETURN
289 ELSE IF( lquery ) THEN
290 RETURN
291 END IF
292*
293* Quick return if possible
294*
295 rcond = zero
296 IF( n.EQ.0 ) THEN
297 rcond = one
298 RETURN
299 ELSE IF( anorm.EQ.zero ) THEN
300 RETURN
301 ELSE IF( n.EQ.1 ) THEN
302 rcond = one
303 RETURN
304 END IF
305*
306 CALL pb_topget( ictxt, 'Combine', 'Columnwise', colctop )
307 CALL pb_topget( ictxt, 'Combine', 'Rowwise', rowctop )
308 CALL pb_topset( ictxt, 'Combine', 'Columnwise', '1-tree' )
309 CALL pb_topset( ictxt, 'Combine', 'Rowwise', '1-tree' )
310*
311 smlnum = pslamch( ictxt, 'Safe minimum' )
312 iroff = mod( ia-1, desca( mb_ ) )
313 icoff = mod( ja-1, desca( nb_ ) )
314 CALL infog2l( ia, ja, desca, nprow, npcol, myrow, mycol, iia, jja,
315 $ iarow, iacol )
316 np = numroc( n+iroff, desca( mb_ ), myrow, iarow, nprow )
317 nq = numroc( n+icoff, desca( nb_ ), mycol, iacol, npcol )
318 iv = iroff + 1
319 ix = iv
320 jv = icoff + 1
321 jx = jv
322*
323 ipx = 1
324 ipv = ipx + np
325 ipw = ipv + np
326 ipnl = 1
327 ipnu = ipnl + nq
328*
329 CALL descset( descv, n+iroff, 1, desca( mb_ ), 1, iarow, mycol,
330 $ ictxt, max( 1, np ) )
331 CALL descset( descx, n+iroff, 1, desca( mb_ ), 1, iarow, mycol,
332 $ ictxt, max( 1, np ) )
333*
334* Estimate the norm of inv(A).
335*
336 ainvnm = zero
337 normin = 'N'
338 IF( onenrm ) THEN
339 kase1 = 1
340 ELSE
341 kase1 = 2
342 END IF
343 kase = 0
344*
345 10 CONTINUE
346 CALL pclacon( n, work( ipv ), iv, jv, descv, work( ipx ), ix, jx,
347 $ descx, ainvnm, kase )
348 IF( kase.NE.0 ) THEN
349 IF( kase.EQ.kase1 ) THEN
350*
351* Multiply by inv(L).
352*
353 descx( csrc_ ) = iacol
354 CALL pclatrs( 'Lower', 'No transpose', 'Unit', normin,
355 $ n, a, ia, ja, desca, work( ipx ), ix, jx,
356 $ descx, sl, rwork( ipnl ), work( ipw ) )
357 descx( csrc_ ) = mycol
358*
359* Multiply by inv(U).
360*
361 descx( csrc_ ) = iacol
362 CALL pclatrs( 'Upper', 'No transpose', 'Non-unit', normin,
363 $ n, a, ia, ja, desca, work( ipx ), ix, jx,
364 $ descx, su, rwork( ipnu ), work( ipw ) )
365 descx( csrc_ ) = mycol
366 ELSE
367*
368* Multiply by inv(U').
369*
370 descx( csrc_ ) = iacol
371 CALL pclatrs( 'Upper', 'Conjugate transpose', 'Non-unit',
372 $ normin, n, a, ia, ja, desca, work( ipx ), ix,
373 $ jx, descx, su, rwork( ipnu ), work( ipw ) )
374 descx( csrc_ ) = mycol
375*
376* Multiply by inv(L').
377*
378 descx( csrc_ ) = iacol
379 CALL pclatrs( 'Lower', 'Conjugate transpose', 'Unit',
380 $ normin, n, a, ia, ja, desca, work( ipx ),
381 $ ix, jx, descx, sl, rwork( ipnl ),
382 $ work( ipw ) )
383 descx( csrc_ ) = mycol
384 END IF
385*
386* Divide X by 1/(SL*SU) if doing so will not cause overflow.
387*
388 scale = sl*su
389 normin = 'Y'
390 IF( scale.NE.one ) THEN
391 CALL pcamax( n, wmax, ixx, work( ipx ), ix, jx, descx, 1 )
392 IF( descx( m_ ).EQ.1 .AND. n.EQ.1 ) THEN
393 CALL pb_topget( ictxt, 'Broadcast', 'Columnwise', cbtop )
394 IF( myrow.EQ.iarow ) THEN
395 CALL cgebs2d( ictxt, 'Column', cbtop, 1, 1, wmax, 1 )
396 ELSE
397 CALL cgebr2d( ictxt, 'Column', cbtop, 1, 1, wmax, 1,
398 $ iarow, mycol )
399 END IF
400 END IF
401 IF( scale.LT.cabs1( wmax )*smlnum .OR. scale.EQ.zero )
402 $ GO TO 20
403 CALL pcsrscl( n, scale, work( ipx ), ix, jx, descx, 1 )
404 END IF
405 GO TO 10
406 END IF
407*
408* Compute the estimate of the reciprocal condition number.
409*
410 IF( ainvnm.NE.zero )
411 $ rcond = ( one / ainvnm ) / anorm
412*
413 20 CONTINUE
414*
415 CALL pb_topset( ictxt, 'Combine', 'Columnwise', colctop )
416 CALL pb_topset( ictxt, 'Combine', 'Rowwise', rowctop )
417*
418 RETURN
419*
420* End of PCGECON
421*
422 END
subroutine chk1mat(ma, mapos0, na, napos0, ia, ja, desca, descapos0, info)
Definition chk1mat.f:3
subroutine descset(desc, m, n, mb, nb, irsrc, icsrc, ictxt, lld)
Definition descset.f:3
subroutine infog2l(grindx, gcindx, desc, nprow, npcol, myrow, mycol, lrindx, lcindx, rsrc, csrc)
Definition infog2l.f:3
subroutine pcgecon(norm, n, a, ia, ja, desca, anorm, rcond, work, lwork, rwork, lrwork, info)
Definition pcgecon.f:3
#define max(A, B)
Definition pcgemr.c:180
subroutine pchk1mat(ma, mapos0, na, napos0, ia, ja, desca, descapos0, nextra, ex, expos, info)
Definition pchkxmat.f:3
subroutine pclacon(n, v, iv, jv, descv, x, ix, jx, descx, est, kase)
Definition pclacon.f:3
subroutine pclatrs(uplo, trans, diag, normin, n, a, ia, ja, desca, x, ix, jx, descx, scale, cnorm, work)
Definition pclatrs.f:4
subroutine pcsrscl(n, sa, sx, ix, jx, descx, incx)
Definition pcsrscl.f:2
subroutine pxerbla(ictxt, srname, info)
Definition pxerbla.f:2