ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pstrcon.f
Go to the documentation of this file.
1  SUBROUTINE pstrcon( NORM, UPLO, DIAG, N, A, IA, JA, DESCA, RCOND,
2  $ WORK, 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 *
10 * .. Scalar Arguments ..
11  CHARACTER DIAG, NORM, UPLO
12  INTEGER IA, JA, INFO, LIWORK, LWORK, N
13  REAL RCOND
14 * ..
15 * .. Array Arguments ..
16  INTEGER DESCA( * ), IWORK( * )
17  REAL A( * ), WORK( * )
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 * PSTRCON estimates the reciprocal of the condition number of a
24 * triangular distributed matrix A(IA:IA+N-1,JA:JA+N-1), in either the
25 * 1-norm or the infinity-norm.
26 *
27 * The norm of A(IA:IA+N-1,JA:JA+N-1) is computed and an estimate is
28 * obtained for norm(inv(A(IA:IA+N-1,JA:JA+N-1))), then the reciprocal
29 * 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 * UPLO (global input) CHARACTER
97 * = 'U': A(IA:IA+N-1,JA:JA+N-1) is upper triangular;
98 * = 'L': A(IA:IA+N-1,JA:JA+N-1) is lower triangular.
99 *
100 * DIAG (global input) CHARACTER
101 * = 'N': A(IA:IA+N-1,JA:JA+N-1) is non-unit triangular;
102 * = 'U': A(IA:IA+N-1,JA:JA+N-1) is unit triangular.
103 *
104 * N (global input) INTEGER
105 * The order of the distributed matrix A(IA:IA+N-1,JA:JA+N-1).
106 * N >= 0.
107 *
108 * A (local input) REAL pointer into the local memory
109 * to an array of dimension ( LLD_A, LOCc(JA+N-1) ). This array
110 * contains the local pieces of the triangular distributed
111 * matrix A(IA:IA+N-1,JA:JA+N-1). If UPLO = 'U', the leading
112 * N-by-N upper triangular part of this distributed matrix con-
113 * tains the upper triangular matrix, and its strictly lower
114 * triangular part is not referenced. If UPLO = 'L', the
115 * leading N-by-N lower triangular part of this ditributed
116 * matrix contains the lower triangular matrix, and the strictly
117 * upper triangular part is not referenced. If DIAG = 'U', the
118 * diagonal elements of A(IA:IA+N-1,JA:JA+N-1) are also not
119 * referenced and are assumed to be 1.
120 *
121 * IA (global input) INTEGER
122 * The row index in the global array A indicating the first
123 * row of sub( A ).
124 *
125 * JA (global input) INTEGER
126 * The column index in the global array A indicating the
127 * first column of sub( A ).
128 *
129 * DESCA (global and local input) INTEGER array of dimension DLEN_.
130 * The array descriptor for the distributed matrix A.
131 *
132 * RCOND (global output) REAL
133 * The reciprocal of the condition number of the distributed
134 * matrix A(IA:IA+N-1,JA:JA+N-1), computed as
135 * RCOND = 1 / ( norm( A(IA:IA+N-1,JA:JA+N-1) ) *
136 * norm( inv(A(IA:IA+N-1,JA:JA+N-1)) ) ).
137 *
138 * WORK (local workspace/local output) REAL array,
139 * dimension (LWORK)
140 * On exit, WORK(1) returns the minimal and optimal LWORK.
141 *
142 * LWORK (local or global input) INTEGER
143 * The dimension of the array WORK.
144 * LWORK is local input and must be at least
145 * LWORK >= 2*LOCr(N+MOD(IA-1,MB_A)) + LOCc(N+MOD(JA-1,NB_A))
146 * + MAX( 2, MAX( NB_A*MAX( 1, CEIL(NPROW-1,NPCOL) ),
147 * LOCc(N+MOD(JA-1,NB_A)) +
148 * NB_A*MAX( 1, CEIL(NPCOL-1,NPROW) ) ).
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 * IWORK (local workspace/local output) INTEGER array,
157 * dimension (LIWORK)
158 * On exit, IWORK(1) returns the minimal and optimal LIWORK.
159 *
160 * LIWORK (local or global input) INTEGER
161 * The dimension of the array IWORK.
162 * LIWORK is local input and must be at least
163 * LIWORK >= LOCr(N+MOD(IA-1,MB_A)).
164 *
165 * If LIWORK = -1, then LIWORK 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  REAL ONE, ZERO
188  parameter( one = 1.0e+0, zero = 0.0e+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, liwmin, lwmin, mycol, myrow, np, npcol,
196  $ npmod, nprow, nq, nqmod
197  REAL AINVNM, ANORM, SCALE, SMLNUM
198  REAL WMAX
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  $ pchk1mat, psamax, pslatrs, pslacon,
207  $ psrscl, pb_topget, pb_topset, pxerbla, sgebr2d,
208  $ sgebs2d
209 * ..
210 * .. External Functions ..
211  LOGICAL LSAME
212  INTEGER ICEIL, INDXG2P, NUMROC
213  REAL PSLAMCH, PSLANTR
214  EXTERNAL iceil, indxg2p, lsame, numroc, pslamch,
215  $ pslantr
216 * ..
217 * .. Intrinsic Functions ..
218  INTRINSIC abs, ichar, max, mod, real
219 * ..
220 * .. Executable Statements ..
221 *
222 * Get grid parameters
223 *
224  ictxt = desca( ctxt_ )
225  CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
226 *
227 * Test the input parameters
228 *
229  info = 0
230  IF( nprow.EQ.-1 ) THEN
231  info = -( 800 + ctxt_ )
232  ELSE
233  CALL chk1mat( n, 4, n, 4, ia, ja, desca, 8, info )
234  IF( info.EQ.0 ) THEN
235  upper = lsame( uplo, 'U' )
236  onenrm = norm.EQ.'1' .OR. lsame( norm, 'O' )
237  nounit = lsame( diag, 'N' )
238  iarow = indxg2p( ia, desca( mb_ ), myrow, desca( rsrc_ ),
239  $ nprow )
240  iacol = indxg2p( ja, desca( nb_ ), mycol, desca( csrc_ ),
241  $ npcol )
242  npmod = numroc( n + mod( ia-1, desca( mb_ ) ), desca( mb_ ),
243  $ myrow, iarow, nprow )
244  nqmod = numroc( n + mod( ja-1, desca( nb_ ) ), desca( nb_ ),
245  $ mycol, iacol, npcol )
246  lwmin = 2*npmod + nqmod +
247  $ max( 2, max( desca( nb_ )*
248  $ max( 1, iceil( nprow-1, npcol ) ), nqmod +
249  $ desca( nb_ )*
250  $ max( 1, iceil( npcol-1, nprow ) ) ) )
251  work( 1 ) = real( lwmin )
252  liwmin = npmod
253  iwork( 1 ) = liwmin
254  lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
255 *
256  IF( .NOT.onenrm .AND. .NOT.lsame( norm, 'I' ) ) THEN
257  info = -1
258  ELSE IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
259  info = -2
260  ELSE IF( .NOT.nounit .AND. .NOT.lsame( diag, 'U' ) ) THEN
261  info = -3
262  ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
263  info = -11
264  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
265  info = -13
266  END IF
267  END IF
268 *
269  IF( onenrm ) THEN
270  idum1( 1 ) = ichar( '1' )
271  ELSE
272  idum1( 1 ) = ichar( 'I' )
273  END IF
274  idum2( 1 ) = 1
275  IF( upper ) THEN
276  idum1( 2 ) = ichar( 'U' )
277  ELSE
278  idum1( 2 ) = ichar( 'L' )
279  END IF
280  idum2( 2 ) = 2
281  IF( nounit ) THEN
282  idum1( 3 ) = ichar( 'N' )
283  ELSE
284  idum1( 3 ) = ichar( 'U' )
285  END IF
286  idum2( 3 ) = 3
287  IF( lwork.EQ.-1 ) THEN
288  idum1( 4 ) = -1
289  ELSE
290  idum1( 4 ) = 1
291  END IF
292  idum2( 4 ) = 11
293  IF( liwork.EQ.-1 ) THEN
294  idum1( 5 ) = -1
295  ELSE
296  idum1( 5 ) = 1
297  END IF
298  idum2( 5 ) = 13
299  CALL pchk1mat( n, 4, n, 4, ia, ja, desca, 8, 5, idum1, idum2,
300  $ info )
301  END IF
302 *
303  IF( info.NE.0 ) THEN
304  CALL pxerbla( ictxt, 'PSTRCON', -info )
305  RETURN
306  ELSE IF( lquery ) THEN
307  RETURN
308  END IF
309 *
310 * Quick return if possible
311 *
312  IF( n.EQ.0 ) THEN
313  rcond = one
314  RETURN
315  END IF
316 *
317  CALL pb_topget( ictxt, 'Combine', 'Columnwise', colctop )
318  CALL pb_topget( ictxt, 'Combine', 'Rowwise', rowctop )
319  CALL pb_topset( ictxt, 'Combine', 'Columnwise', '1-tree' )
320  CALL pb_topset( ictxt, 'Combine', 'Rowwise', '1-tree' )
321 *
322  rcond = zero
323  smlnum = pslamch( ictxt, 'Safe minimum' )*real( max( 1, n ) )
324  CALL infog2l( ia, ja, desca, nprow, npcol, myrow, mycol, iia, jja,
325  $ iarow, iacol )
326  iroff = mod( ia-1, desca( mb_ ) )
327  icoff = mod( ja-1, desca( nb_ ) )
328  np = numroc( n+iroff, desca( mb_ ), myrow, iarow, nprow )
329  nq = numroc( n+icoff, desca( nb_ ), mycol, iacol, npcol )
330  iv = iroff + 1
331  ix = iv
332  jv = icoff + 1
333  jx = jv
334 *
335  ipx = 1
336  ipv = ipx + np
337  ipn = ipv + np
338  ipw = ipn + nq
339 *
340  CALL descset( descv, n+iroff, 1, desca( mb_ ), 1, iarow, mycol,
341  $ ictxt, max( 1, np ) )
342  CALL descset( descx, n+iroff, 1, desca( mb_ ), 1, iarow, mycol,
343  $ ictxt, max( 1, np ) )
344 *
345 * Compute the norm of the triangular matrix A.
346 *
347  anorm = pslantr( norm, uplo, diag, n, n, a, ia, ja, desca, work )
348 *
349 * Continue only if ANORM > 0.
350 *
351  IF( anorm.GT.zero ) THEN
352 *
353 * Estimate the norm of the inverse of A.
354 *
355  ainvnm = zero
356  normin = 'N'
357  IF( onenrm ) THEN
358  kase1 = 1
359  ELSE
360  kase1 = 2
361  END IF
362  kase = 0
363  10 CONTINUE
364  CALL pslacon( n, work( ipv ), iv, jv, descv, work( ipx ),
365  $ ix, jx, descx, iwork, ainvnm, kase )
366  IF( kase.NE.0 ) THEN
367  IF( kase.EQ.kase1 ) THEN
368 *
369 * Multiply by inv(A).
370 *
371  descx( csrc_ ) = iacol
372  CALL pslatrs( uplo, 'No transpose', diag, normin,
373  $ n, a, ia, ja, desca, work( ipx ), ix, jx,
374  $ descx, scale, work( ipn ), work( ipw ) )
375  descx( csrc_ ) = mycol
376  ELSE
377 *
378 * Multiply by inv(A').
379 *
380  descx( csrc_ ) = iacol
381  CALL pslatrs( uplo, 'Transpose', diag, normin,
382  $ n, a, ia, ja, desca, work( ipx ), ix, jx,
383  $ descx, scale, work( ipn ), work( ipw ) )
384  descx( csrc_ ) = mycol
385  END IF
386  normin = 'Y'
387 *
388 * Multiply by 1/SCALE if doing so will not cause overflow.
389 *
390  IF( scale.NE.one ) THEN
391  CALL psamax( n, wmax, ixx, work( ipx ), ix, jx,
392  $ descx, 1 )
393  IF( descx( m_ ).EQ.1 .AND. n.EQ.1 ) THEN
394  CALL pb_topget( ictxt, 'Broadcast', 'Columnwise',
395  $ cbtop )
396  IF( myrow.EQ.iarow ) THEN
397  CALL sgebs2d( ictxt, 'Column', cbtop, 1, 1, wmax,
398  $ 1 )
399  ELSE
400  CALL sgebr2d( ictxt, 'Column', cbtop, 1, 1, wmax,
401  $ 1, iarow, mycol )
402  END IF
403  END IF
404  IF( scale.LT.abs( wmax )*smlnum .OR. scale.EQ.zero )
405  $ GO TO 20
406  CALL psrscl( n, scale, work( ipx ), ix, jx, descx, 1 )
407  END IF
408  GO TO 10
409  END IF
410 *
411 * Compute the estimate of the reciprocal condition number.
412 *
413  IF( ainvnm.NE.zero )
414  $ rcond = ( one / anorm ) / ainvnm
415  END IF
416 *
417  20 CONTINUE
418 *
419  CALL pb_topset( ictxt, 'Combine', 'Columnwise', colctop )
420  CALL pb_topset( ictxt, 'Combine', 'Rowwise', rowctop )
421 *
422  RETURN
423 *
424 * End of PSTRCON
425 *
426  END
pslatrs
subroutine pslatrs(UPLO, TRANS, DIAG, NORMIN, N, A, IA, JA, DESCA, X, IX, JX, DESCX, SCALE, CNORM, WORK)
Definition: pslatrs.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
pstrcon
subroutine pstrcon(NORM, UPLO, DIAG, N, A, IA, JA, DESCA, RCOND, WORK, LWORK, IWORK, LIWORK, INFO)
Definition: pstrcon.f:3
psrscl
subroutine psrscl(N, SA, SX, IX, JX, DESCX, INCX)
Definition: psrscl.f:2
pslacon
subroutine pslacon(N, V, IV, JV, DESCV, X, IX, JX, DESCX, ISGN, EST, KASE)
Definition: pslacon.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
pxerbla
subroutine pxerbla(ICTXT, SRNAME, INFO)
Definition: pxerbla.f:2