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