ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pzgecon.f
Go to the documentation of this file.
1  SUBROUTINE pzgecon( 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  DOUBLE PRECISION ANORM, RCOND
13 * ..
14 * .. Array Arguments ..
15  INTEGER DESCA( * )
16  DOUBLE PRECISION RWORK( * )
17  COMPLEX*16 A( * ), WORK( * )
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 * PZGECON 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 * PZGETRF.
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*16 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) DOUBLE PRECISION
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) DOUBLE PRECISION
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*16 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) DOUBLE PRECISION 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  DOUBLE PRECISION ONE, ZERO
182  parameter( one = 1.0d+0, zero = 0.0d+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  DOUBLE PRECISION AINVNM, SCALE, SL, SMLNUM, SU
192  COMPLEX*16 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, infog2l,
200  $ pchk1mat, pb_topget, pb_topset, pxerbla,
201  $ pzamax, pzlatrs, pzlacon, pzdrscl,
202  $ zgebr2d, zgebs2d
203 * ..
204 * .. External Functions ..
205  LOGICAL LSAME
206  INTEGER ICEIL, INDXG2P, NUMROC
207  DOUBLE PRECISION PDLAMCH
208  EXTERNAL iceil, indxg2p, lsame, numroc, pdlamch
209 * ..
210 * .. Intrinsic Functions ..
211  INTRINSIC abs, dble, dimag, ichar, max, mod
212 * ..
213 * .. Statement Functions ..
214  DOUBLE PRECISION CABS1
215 * ..
216 * .. Statement Function definitions ..
217  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( 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 ) = dble( lwmin )
249  lrwmin = max( 1, 2*nqmod )
250  rwork( 1 ) = dble( 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, 'PZGECON', -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 = pdlamch( 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 pzlacon( 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 pzlatrs( '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 pzlatrs( '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 pzlatrs( '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 pzlatrs( '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 pzamax( 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 zgebs2d( ictxt, 'Column', cbtop, 1, 1, wmax, 1 )
396  ELSE
397  CALL zgebr2d( 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 pzdrscl( 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 PZGECON
421 *
422  END
max
#define max(A, B)
Definition: pcgemr.c:180
infog2l
subroutine infog2l(GRINDX, GCINDX, DESC, NPROW, NPCOL, MYROW, MYCOL, LRINDX, LCINDX, RSRC, CSRC)
Definition: infog2l.f:3
pzgecon
subroutine pzgecon(NORM, N, A, IA, JA, DESCA, ANORM, RCOND, WORK, LWORK, RWORK, LRWORK, INFO)
Definition: pzgecon.f:3
pzdrscl
subroutine pzdrscl(N, SA, SX, IX, JX, DESCX, INCX)
Definition: pzdrscl.f:2
pchk1mat
subroutine pchk1mat(MA, MAPOS0, NA, NAPOS0, IA, JA, DESCA, DESCAPOS0, NEXTRA, EX, EXPOS, INFO)
Definition: pchkxmat.f:3
descset
subroutine descset(DESC, M, N, MB, NB, IRSRC, ICSRC, ICTXT, LLD)
Definition: descset.f:3
pzlacon
subroutine pzlacon(N, V, IV, JV, DESCV, X, IX, JX, DESCX, EST, KASE)
Definition: pzlacon.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
pzlatrs
subroutine pzlatrs(UPLO, TRANS, DIAG, NORMIN, N, A, IA, JA, DESCA, X, IX, JX, DESCX, SCALE, CNORM, WORK)
Definition: pzlatrs.f:4