ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pspoequ.f
Go to the documentation of this file.
1  SUBROUTINE pspoequ( N, A, IA, JA, DESCA, SR, SC, SCOND, AMAX,
2  $ 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 1, 1997
8 *
9 * .. Scalar Arguments ..
10  INTEGER IA, INFO, JA, N
11  REAL AMAX, SCOND
12 * ..
13 * .. Array Arguments ..
14  INTEGER DESCA( * )
15  REAL A( * ), SC( * ), SR( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * PSPOEQU computes row and column scalings intended to
22 * equilibrate a distributed symmetric positive definite matrix
23 * sub( A ) = A(IA:IA+N-1,JA:JA+N-1) and reduce its condition number
24 * (with respect to the two-norm). SR and SC contain the scale
25 * factors, S(i) = 1/sqrt(A(i,i)), chosen so that the scaled distri-
26 * buted matrix B with elements B(i,j) = S(i)*A(i,j)*S(j) has ones on
27 * the diagonal. This choice of SR and SC puts the condition number
28 * of B within a factor N of the smallest possible condition number
29 * over all possible diagonal scalings.
30 *
31 * The scaling factor are stored along process rows in SR and along
32 * process columns in SC. The duplication of information simplifies
33 * greatly the application of the factors.
34 *
35 * Notes
36 * =====
37 *
38 * Each global data object is described by an associated description
39 * vector. This vector stores the information required to establish
40 * the mapping between an object element and its corresponding process
41 * and memory location.
42 *
43 * Let A be a generic term for any 2D block cyclicly distributed array.
44 * Such a global array has an associated description vector DESCA.
45 * In the following comments, the character _ should be read as
46 * "of the global array".
47 *
48 * NOTATION STORED IN EXPLANATION
49 * --------------- -------------- --------------------------------------
50 * DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
51 * DTYPE_A = 1.
52 * CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
53 * the BLACS process grid A is distribu-
54 * ted over. The context itself is glo-
55 * bal, but the handle (the integer
56 * value) may vary.
57 * M_A (global) DESCA( M_ ) The number of rows in the global
58 * array A.
59 * N_A (global) DESCA( N_ ) The number of columns in the global
60 * array A.
61 * MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
62 * the rows of the array.
63 * NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
64 * the columns of the array.
65 * RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
66 * row of the array A is distributed.
67 * CSRC_A (global) DESCA( CSRC_ ) The process column over which the
68 * first column of the array A is
69 * distributed.
70 * LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
71 * array. LLD_A >= MAX(1,LOCr(M_A)).
72 *
73 * Let K be the number of rows or columns of a distributed matrix,
74 * and assume that its process grid has dimension p x q.
75 * LOCr( K ) denotes the number of elements of K that a process
76 * would receive if K were distributed over the p processes of its
77 * process column.
78 * Similarly, LOCc( K ) denotes the number of elements of K that a
79 * process would receive if K were distributed over the q processes of
80 * its process row.
81 * The values of LOCr() and LOCc() may be determined via a call to the
82 * ScaLAPACK tool function, NUMROC:
83 * LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
84 * LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
85 * An upper bound for these quantities may be computed by:
86 * LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
87 * LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
88 *
89 * Arguments
90 * =========
91 *
92 * N (global input) INTEGER
93 * The number of rows and columns to be operated on i.e the
94 * order of the distributed submatrix sub( A ). N >= 0.
95 *
96 * A (local input) REAL pointer into the local memory to an
97 * array of local dimension ( LLD_A, LOCc(JA+N-1) ), the
98 * N-by-N symmetric positive definite distributed matrix
99 * sub( A ) whose scaling factors are to be computed. Only the
100 * diagonal elements of sub( A ) are referenced.
101 *
102 * IA (global input) INTEGER
103 * The row index in the global array A indicating the first
104 * row of sub( A ).
105 *
106 * JA (global input) INTEGER
107 * The column index in the global array A indicating the
108 * first column of sub( A ).
109 *
110 * DESCA (global and local input) INTEGER array of dimension DLEN_.
111 * The array descriptor for the distributed matrix A.
112 *
113 * SR (local output) REAL array, dimension LOCr(M_A)
114 * If INFO = 0, SR(IA:IA+N-1) contains the row scale factors
115 * for sub( A ). SR is aligned with the distributed matrix A,
116 * and replicated across every process column. SR is tied to the
117 * distributed matrix A.
118 *
119 * SC (local output) REAL array, dimension LOCc(N_A)
120 * If INFO = 0, SC(JA:JA+N-1) contains the column scale factors
121 * for A(IA:IA+M-1,JA:JA+N-1). SC is aligned with the distribu-
122 * ted matrix A, and replicated down every process row. SC is
123 * tied to the distributed matrix A.
124 *
125 * SCOND (global output) REAL
126 * If INFO = 0, SCOND contains the ratio of the smallest SR(i)
127 * (or SC(j)) to the largest SR(i) (or SC(j)), with
128 * IA <= i <= IA+N-1 and JA <= j <= JA+N-1. If SCOND >= 0.1
129 * and AMAX is neither too large nor too small, it is not worth
130 * scaling by SR (or SC).
131 *
132 * AMAX (global output) REAL
133 * Absolute value of largest matrix element. If AMAX is very
134 * close to overflow or very close to underflow, the matrix
135 * should be scaled.
136 *
137 * INFO (global output) INTEGER
138 * = 0: successful exit
139 * < 0: If the i-th argument is an array and the j-entry had
140 * an illegal value, then INFO = -(i*100+j), if the i-th
141 * argument is a scalar and had an illegal value, then
142 * INFO = -i.
143 * > 0: If INFO = K, the K-th diagonal entry of sub( A ) is
144 * nonpositive.
145 *
146 * =====================================================================
147 *
148 * .. Parameters ..
149  INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
150  $ lld_, mb_, m_, nb_, n_, rsrc_
151  parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
152  $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
153  $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
154  REAL ZERO, ONE
155  parameter( zero = 0.0e+0, one = 1.0e+0 )
156 * ..
157 * .. Local Scalars ..
158  CHARACTER ALLCTOP, COLCTOP, ROWCTOP
159  INTEGER IACOL, IAROW, ICOFF, ICTXT, ICURCOL, ICURROW,
160  $ idumm, ii, iia, ioffa, ioffd, iroff, j, jb, jj,
161  $ jja, jn, lda, ll, mycol, myrow, np, npcol,
162  $ nprow, nq
163  REAL AII, SMIN
164 * ..
165 * .. Local Arrays ..
166  INTEGER DESCSC( DLEN_ ), DESCSR( DLEN_ )
167 * ..
168 * .. External Subroutines ..
169  EXTERNAL blacs_gridinfo, chk1mat, descset, igamn2d,
170  $ infog2l, pchk1mat, pb_topget, pxerbla,
171  $ sgamn2d, sgamx2d, sgsum2d
172 * ..
173 * .. External Functions ..
174  INTEGER ICEIL, NUMROC
175  REAL PSLAMCH
176  EXTERNAL iceil, numroc, pslamch
177 * ..
178 * .. Intrinsic Functions ..
179  INTRINSIC max, min, mod, sqrt
180 * ..
181 * .. Executable Statements ..
182 *
183 * Get grid parameters
184 *
185  ictxt = desca( ctxt_ )
186  CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
187 *
188 * Test the input parameters.
189 *
190  info = 0
191  IF( nprow.EQ.-1 ) THEN
192  info = -(500+ctxt_)
193  ELSE
194  CALL chk1mat( n, 1, n, 1, ia, ja, desca, 5, info )
195  CALL pchk1mat( n, 1, n, 1, ia, ja, desca, 5, 0, idumm, idumm,
196  $ info )
197  END IF
198 *
199  IF( info.NE.0 ) THEN
200  CALL pxerbla( ictxt, 'PSPOEQU', -info )
201  RETURN
202  END IF
203 *
204 * Quick return if possible
205 *
206  IF( n.EQ.0 ) THEN
207  scond = one
208  amax = zero
209  RETURN
210  END IF
211 *
212  CALL pb_topget( ictxt, 'Combine', 'All', allctop )
213  CALL pb_topget( ictxt, 'Combine', 'Rowwise', rowctop )
214  CALL pb_topget( ictxt, 'Combine', 'Columnwise', colctop )
215 *
216 * Compute some local indexes
217 *
218  CALL infog2l( ia, ja, desca, nprow, npcol, myrow, mycol, iia, jja,
219  $ iarow, iacol )
220  iroff = mod( ia-1, desca( mb_ ) )
221  icoff = mod( ja-1, desca( nb_ ) )
222  np = numroc( n+iroff, desca( mb_ ), myrow, iarow, nprow )
223  nq = numroc( n+icoff, desca( nb_ ), mycol, iacol, npcol )
224  IF( myrow.EQ.iarow )
225  $ np = np - iroff
226  IF( mycol.EQ.iacol )
227  $ nq = nq - icoff
228  jn = min( iceil( ja, desca( nb_ ) ) * desca( nb_ ), ja+n-1 )
229  lda = desca( lld_ )
230 *
231 * Assign descriptors for SR and SC arrays
232 *
233  CALL descset( descsr, n, 1, desca( mb_ ), 1, 0, 0, ictxt,
234  $ max( 1, np ) )
235  CALL descset( descsc, 1, n, 1, desca( nb_ ), 0, 0, ictxt, 1 )
236 *
237 * Initialize the scaling factors to zero.
238 *
239  DO 10 ii = iia, iia+np-1
240  sr( ii ) = zero
241  10 CONTINUE
242 *
243  DO 20 jj = jja, jja+nq-1
244  sc( jj ) = zero
245  20 CONTINUE
246 *
247 * Find the minimum and maximum diagonal elements.
248 * Handle first block separately.
249 *
250  ii = iia
251  jj = jja
252  jb = jn-ja+1
253  smin = one / pslamch( ictxt, 'S' )
254  amax = zero
255 *
256  ioffa = ii+(jj-1)*lda
257  IF( myrow.EQ.iarow .AND. mycol.EQ.iacol ) THEN
258  ioffd = ioffa
259  DO 30 ll = 0, jb-1
260  aii = a( ioffd )
261  sr( ii+ll ) = aii
262  sc( jj+ll ) = aii
263  smin = min( smin, aii )
264  amax = max( amax, aii )
265  IF( aii.LE.zero .AND. info.EQ.0 )
266  $ info = ll + 1
267  ioffd = ioffd + lda + 1
268  30 CONTINUE
269  END IF
270 *
271  IF( myrow.EQ.iarow ) THEN
272  ii = ii + jb
273  ioffa = ioffa + jb
274  END IF
275  IF( mycol.EQ.iacol ) THEN
276  jj = jj + jb
277  ioffa = ioffa + jb*lda
278  END IF
279  icurrow = mod( iarow+1, nprow )
280  icurcol = mod( iacol+1, npcol )
281 *
282 * Loop over remaining blocks of columns
283 *
284  DO 50 j = jn+1, ja+n-1, desca( nb_ )
285  jb = min( n-j+ja, desca( nb_ ) )
286 *
287  IF( myrow.EQ.icurrow .AND. mycol.EQ.icurcol ) THEN
288  ioffd = ioffa
289  DO 40 ll = 0, jb-1
290  aii = a( ioffd )
291  sr( ii+ll ) = aii
292  sc( jj+ll ) = aii
293  smin = min( smin, aii )
294  amax = max( amax, aii )
295  IF( aii.LE.zero .AND. info.EQ.0 )
296  $ info = j + ll - ja + 1
297  ioffd = ioffd + lda + 1
298  40 CONTINUE
299  END IF
300 *
301  IF( myrow.EQ.icurrow ) THEN
302  ii = ii + jb
303  ioffa = ioffa + jb
304  END IF
305  IF( mycol.EQ.icurcol ) THEN
306  jj = jj + jb
307  ioffa = ioffa + jb*lda
308  END IF
309  icurrow = mod( icurrow+1, nprow )
310  icurcol = mod( icurcol+1, npcol )
311 *
312  50 CONTINUE
313 *
314 * Compute scaling factors
315 *
316  CALL sgsum2d( ictxt, 'Columnwise', colctop, 1, nq, sc( jja ),
317  $ 1, -1, mycol )
318  CALL sgsum2d( ictxt, 'Rowwise', rowctop, np, 1, sr( iia ),
319  $ max( 1, np ), -1, mycol )
320 *
321  CALL sgamx2d( ictxt, 'All', allctop, 1, 1, amax, 1, idumm, idumm,
322  $ -1, -1, mycol )
323  CALL sgamn2d( ictxt, 'All', allctop, 1, 1, smin, 1, idumm, idumm,
324  $ -1, -1, mycol )
325 *
326  IF( smin.LE.zero ) THEN
327 *
328 * Find the first non-positive diagonal element and return.
329 *
330  CALL igamn2d( ictxt, 'All', allctop, 1, 1, info, 1, ii, jj, -1,
331  $ -1, mycol )
332  RETURN
333 *
334  ELSE
335 *
336 * Set the scale factors to the reciprocals
337 * of the diagonal elements.
338 *
339  DO 60 ii = iia, iia+np-1
340  sr( ii ) = one / sqrt( sr( ii ) )
341  60 CONTINUE
342 *
343  DO 70 jj = jja, jja+nq-1
344  sc( jj ) = one / sqrt( sc( jj ) )
345  70 CONTINUE
346 *
347 * Compute SCOND = min(S(I)) / max(S(I))
348 *
349  scond = sqrt( smin ) / sqrt( amax )
350 *
351  END IF
352 *
353  RETURN
354 *
355 * End of PSPOEQU
356 *
357  END
pspoequ
subroutine pspoequ(N, A, IA, JA, DESCA, SR, SC, SCOND, AMAX, INFO)
Definition: pspoequ.f:3
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
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
min
#define min(A, B)
Definition: pcgemr.c:181