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