ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pcgesvx.f
Go to the documentation of this file.
1  SUBROUTINE pcgesvx( FACT, TRANS, N, NRHS, A, IA, JA, DESCA, AF,
2  $ IAF, JAF, DESCAF, IPIV, EQUED, R, C, B, IB,
3  $ JB, DESCB, X, IX, JX, DESCX, RCOND, FERR,
4  $ BERR, WORK, LWORK, RWORK, LRWORK, INFO )
5 *
6 * -- ScaLAPACK routine (version 1.7) --
7 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
8 * and University of California, Berkeley.
9 * December 31, 1998
10 *
11 * .. Scalar Arguments ..
12  CHARACTER EQUED, FACT, TRANS
13  INTEGER IA, IAF, IB, INFO, IX, JA, JAF, JB, JX, LRWORK,
14  $ LWORK, N, NRHS
15  REAL RCOND
16 * ..
17 * .. Array Arguments ..
18  INTEGER DESCA( * ), DESCAF( * ), DESCB( * ),
19  $ DESCX( * ), IPIV( * )
20  REAL BERR( * ), C( * ), FERR( * ), R( * ),
21  $ rwork( * )
22  COMPLEX A( * ), AF( * ), B( * ), WORK( * ), X( * )
23 * ..
24 *
25 * Purpose
26 * =======
27 *
28 * PCGESVX uses the LU factorization to compute the solution to a
29 * complex system of linear equations
30 *
31 * A(IA:IA+N-1,JA:JA+N-1) * X = B(IB:IB+N-1,JB:JB+NRHS-1),
32 *
33 * where A(IA:IA+N-1,JA:JA+N-1) is an N-by-N matrix and X and
34 * B(IB:IB+N-1,JB:JB+NRHS-1) are N-by-NRHS matrices.
35 *
36 * Error bounds on the solution and a condition estimate are also
37 * provided.
38 *
39 * Notes
40 * =====
41 *
42 * Each global data object is described by an associated description
43 * vector. This vector stores the information required to establish
44 * the mapping between an object element and its corresponding process
45 * and memory location.
46 *
47 * Let A be a generic term for any 2D block cyclicly distributed array.
48 * Such a global array has an associated description vector DESCA.
49 * In the following comments, the character _ should be read as
50 * "of the global array".
51 *
52 * NOTATION STORED IN EXPLANATION
53 * --------------- -------------- --------------------------------------
54 * DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
55 * DTYPE_A = 1.
56 * CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
57 * the BLACS process grid A is distribu-
58 * ted over. The context itself is glo-
59 * bal, but the handle (the integer
60 * value) may vary.
61 * M_A (global) DESCA( M_ ) The number of rows in the global
62 * array A.
63 * N_A (global) DESCA( N_ ) The number of columns in the global
64 * array A.
65 * MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
66 * the rows of the array.
67 * NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
68 * the columns of the array.
69 * RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
70 * row of the array A is distributed.
71 * CSRC_A (global) DESCA( CSRC_ ) The process column over which the
72 * first column of the array A is
73 * distributed.
74 * LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
75 * array. LLD_A >= MAX(1,LOCr(M_A)).
76 *
77 * Let K be the number of rows or columns of a distributed matrix,
78 * and assume that its process grid has dimension p x q.
79 * LOCr( K ) denotes the number of elements of K that a process
80 * would receive if K were distributed over the p processes of its
81 * process column.
82 * Similarly, LOCc( K ) denotes the number of elements of K that a
83 * process would receive if K were distributed over the q processes of
84 * its process row.
85 * The values of LOCr() and LOCc() may be determined via a call to the
86 * ScaLAPACK tool function, NUMROC:
87 * LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
88 * LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
89 * An upper bound for these quantities may be computed by:
90 * LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
91 * LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
92 *
93 * Description
94 * ===========
95 *
96 * In the following description, A denotes A(IA:IA+N-1,JA:JA+N-1),
97 * B denotes B(IB:IB+N-1,JB:JB+NRHS-1) and X denotes
98 * X(IX:IX+N-1,JX:JX+NRHS-1).
99 *
100 * The following steps are performed:
101 *
102 * 1. If FACT = 'E', real scaling factors are computed to equilibrate
103 * the system:
104 * TRANS = 'N': diag(R)*A*diag(C) *inv(diag(C))*X = diag(R)*B
105 * TRANS = 'T': (diag(R)*A*diag(C))**T *inv(diag(R))*X = diag(C)*B
106 * TRANS = 'C': (diag(R)*A*diag(C))**H *inv(diag(R))*X = diag(C)*B
107 * Whether or not the system will be equilibrated depends on the
108 * scaling of the matrix A, but if equilibration is used, A is
109 * overwritten by diag(R)*A*diag(C) and B by diag(R)*B (if TRANS='N')
110 * or diag(C)*B (if TRANS = 'T' or 'C').
111 *
112 * 2. If FACT = 'N' or 'E', the LU decomposition is used to factor the
113 * matrix A (after equilibration if FACT = 'E') as
114 * A = P * L * U,
115 * where P is a permutation matrix, L is a unit lower triangular
116 * matrix, and U is upper triangular.
117 *
118 * 3. The factored form of A is used to estimate the condition number
119 * of the matrix A. If the reciprocal of the condition number is
120 * less than machine precision, steps 4-6 are skipped.
121 *
122 * 4. The system of equations is solved for X using the factored form
123 * of A.
124 *
125 * 5. Iterative refinement is applied to improve the computed solution
126 * matrix and calculate error bounds and backward error estimates
127 * for it.
128 *
129 * 6. If FACT = 'E' and equilibration was used, the matrix X is
130 * premultiplied by diag(C) (if TRANS = 'N') or diag(R) (if
131 * TRANS = 'T' or 'C') so that it solves the original system
132 * before equilibration.
133 *
134 * Arguments
135 * =========
136 *
137 * FACT (global input) CHARACTER
138 * Specifies whether or not the factored form of the matrix
139 * A(IA:IA+N-1,JA:JA+N-1) is supplied on entry, and if not,
140 * whether the matrix A(IA:IA+N-1,JA:JA+N-1) should be
141 * equilibrated before it is factored.
142 * = 'F': On entry, AF(IAF:IAF+N-1,JAF:JAF+N-1) and IPIV con-
143 * tain the factored form of A(IA:IA+N-1,JA:JA+N-1).
144 * If EQUED is not 'N', the matrix
145 * A(IA:IA+N-1,JA:JA+N-1) has been equilibrated with
146 * scaling factors given by R and C.
147 * A(IA:IA+N-1,JA:JA+N-1), AF(IAF:IAF+N-1,JAF:JAF+N-1),
148 * and IPIV are not modified.
149 * = 'N': The matrix A(IA:IA+N-1,JA:JA+N-1) will be copied to
150 * AF(IAF:IAF+N-1,JAF:JAF+N-1) and factored.
151 * = 'E': The matrix A(IA:IA+N-1,JA:JA+N-1) will be equili-
152 * brated if necessary, then copied to
153 * AF(IAF:IAF+N-1,JAF:JAF+N-1) and factored.
154 *
155 * TRANS (global input) CHARACTER
156 * Specifies the form of the system of equations:
157 * = 'N': A(IA:IA+N-1,JA:JA+N-1) * X(IX:IX+N-1,JX:JX+NRHS-1)
158 * = B(IB:IB+N-1,JB:JB+NRHS-1) (No transpose)
159 * = 'T': A(IA:IA+N-1,JA:JA+N-1)**T * X(IX:IX+N-1,JX:JX+NRHS-1)
160 * = B(IB:IB+N-1,JB:JB+NRHS-1) (Transpose)
161 * = 'C': A(IA:IA+N-1,JA:JA+N-1)**H * X(IX:IX+N-1,JX:JX+NRHS-1)
162 * = B(IB:IB+N-1,JB:JB+NRHS-1) (Conjugate transpose)
163 *
164 * N (global input) INTEGER
165 * The number of rows and columns to be operated on, i.e. the
166 * order of the distributed submatrix A(IA:IA+N-1,JA:JA+N-1).
167 * N >= 0.
168 *
169 * NRHS (global input) INTEGER
170 * The number of right-hand sides, i.e., the number of columns
171 * of the distributed submatrices B(IB:IB+N-1,JB:JB+NRHS-1) and
172 * X(IX:IX+N-1,JX:JX+NRHS-1). NRHS >= 0.
173 *
174 * A (local input/local output) COMPLEX pointer into
175 * the local memory to an array of local dimension
176 * (LLD_A,LOCc(JA+N-1)). On entry, the N-by-N matrix
177 * A(IA:IA+N-1,JA:JA+N-1). If FACT = 'F' and EQUED is not 'N',
178 * then A(IA:IA+N-1,JA:JA+N-1) must have been equilibrated by
179 * the scaling factors in R and/or C. A(IA:IA+N-1,JA:JA+N-1) is
180 * not modified if FACT = 'F' or 'N', or if FACT = 'E' and
181 * EQUED = 'N' on exit.
182 *
183 * On exit, if EQUED .ne. 'N', A(IA:IA+N-1,JA:JA+N-1) is scaled
184 * as follows:
185 * EQUED = 'R': A(IA:IA+N-1,JA:JA+N-1) :=
186 * diag(R) * A(IA:IA+N-1,JA:JA+N-1)
187 * EQUED = 'C': A(IA:IA+N-1,JA:JA+N-1) :=
188 * A(IA:IA+N-1,JA:JA+N-1) * diag(C)
189 * EQUED = 'B': A(IA:IA+N-1,JA:JA+N-1) :=
190 * diag(R) * A(IA:IA+N-1,JA:JA+N-1) * diag(C).
191 *
192 * IA (global input) INTEGER
193 * The row index in the global array A indicating the first
194 * row of sub( A ).
195 *
196 * JA (global input) INTEGER
197 * The column index in the global array A indicating the
198 * first column of sub( A ).
199 *
200 * DESCA (global and local input) INTEGER array of dimension DLEN_.
201 * The array descriptor for the distributed matrix A.
202 *
203 * AF (local input or local output) COMPLEX pointer
204 * into the local memory to an array of local dimension
205 * (LLD_AF,LOCc(JA+N-1)). If FACT = 'F', then
206 * AF(IAF:IAF+N-1,JAF:JAF+N-1) is an input argument and on
207 * entry contains the factors L and U from the factorization
208 * A(IA:IA+N-1,JA:JA+N-1) = P*L*U as computed by PCGETRF.
209 * If EQUED .ne. 'N', then AF is the factored form of the
210 * equilibrated matrix A(IA:IA+N-1,JA:JA+N-1).
211 *
212 * If FACT = 'N', then AF(IAF:IAF+N-1,JAF:JAF+N-1) is an output
213 * argument and on exit returns the factors L and U from the
214 * factorization A(IA:IA+N-1,JA:JA+N-1) = P*L*U of the original
215 * matrix A(IA:IA+N-1,JA:JA+N-1).
216 *
217 * If FACT = 'E', then AF(IAF:IAF+N-1,JAF:JAF+N-1) is an output
218 * argument and on exit returns the factors L and U from the
219 * factorization A(IA:IA+N-1,JA:JA+N-1) = P*L*U of the equili-
220 * brated matrix A(IA:IA+N-1,JA:JA+N-1) (see the description of
221 * A(IA:IA+N-1,JA:JA+N-1) for the form of the equilibrated
222 * matrix).
223 *
224 * IAF (global input) INTEGER
225 * The row index in the global array AF indicating the first
226 * row of sub( AF ).
227 *
228 * JAF (global input) INTEGER
229 * The column index in the global array AF indicating the
230 * first column of sub( AF ).
231 *
232 * DESCAF (global and local input) INTEGER array of dimension DLEN_.
233 * The array descriptor for the distributed matrix AF.
234 *
235 * IPIV (local input or local output) INTEGER array, dimension
236 * LOCr(M_A)+MB_A. If FACT = 'F', then IPIV is an input argu-
237 * ment and on entry contains the pivot indices from the fac-
238 * torization A(IA:IA+N-1,JA:JA+N-1) = P*L*U as computed by
239 * PCGETRF; IPIV(i) -> The global row local row i was
240 * swapped with. This array must be aligned with
241 * A( IA:IA+N-1, * ).
242 *
243 * If FACT = 'N', then IPIV is an output argument and on exit
244 * contains the pivot indices from the factorization
245 * A(IA:IA+N-1,JA:JA+N-1) = P*L*U of the original matrix
246 * A(IA:IA+N-1,JA:JA+N-1).
247 *
248 * If FACT = 'E', then IPIV is an output argument and on exit
249 * contains the pivot indices from the factorization
250 * A(IA:IA+N-1,JA:JA+N-1) = P*L*U of the equilibrated matrix
251 * A(IA:IA+N-1,JA:JA+N-1).
252 *
253 * EQUED (global input or global output) CHARACTER
254 * Specifies the form of equilibration that was done.
255 * = 'N': No equilibration (always true if FACT = 'N').
256 * = 'R': Row equilibration, i.e., A(IA:IA+N-1,JA:JA+N-1) has
257 * been premultiplied by diag(R).
258 * = 'C': Column equilibration, i.e., A(IA:IA+N-1,JA:JA+N-1)
259 * has been postmultiplied by diag(C).
260 * = 'B': Both row and column equilibration, i.e.,
261 * A(IA:IA+N-1,JA:JA+N-1) has been replaced by
262 * diag(R) * A(IA:IA+N-1,JA:JA+N-1) * diag(C).
263 * EQUED is an input variable if FACT = 'F'; otherwise, it is an
264 * output variable.
265 *
266 * R (local input or local output) REAL array,
267 * dimension LOCr(M_A).
268 * The row scale factors for A(IA:IA+N-1,JA:JA+N-1).
269 * If EQUED = 'R' or 'B', A(IA:IA+N-1,JA:JA+N-1) is multiplied
270 * on the left by diag(R); if EQUED='N' or 'C', R is not acces-
271 * sed. R is an input variable if FACT = 'F'; otherwise, R is
272 * an output variable.
273 * If FACT = 'F' and EQUED = 'R' or 'B', each element of R must
274 * be positive.
275 * R is replicated in every process column, and is aligned
276 * with the distributed matrix A.
277 *
278 * C (local input or local output) REAL array,
279 * dimension LOCc(N_A).
280 * The column scale factors for A(IA:IA+N-1,JA:JA+N-1).
281 * If EQUED = 'C' or 'B', A(IA:IA+N-1,JA:JA+N-1) is multiplied
282 * on the right by diag(C); if EQUED = 'N' or 'R', C is not
283 * accessed. C is an input variable if FACT = 'F'; otherwise,
284 * C is an output variable. If FACT = 'F' and EQUED = 'C' or
285 * 'B', each element of C must be positive.
286 * C is replicated in every process row, and is aligned with
287 * the distributed matrix A.
288 *
289 * B (local input/local output) COMPLEX pointer
290 * into the local memory to an array of local dimension
291 * (LLD_B,LOCc(JB+NRHS-1) ). On entry, the N-by-NRHS right-hand
292 * side matrix B(IB:IB+N-1,JB:JB+NRHS-1). On exit, if
293 * EQUED = 'N', B(IB:IB+N-1,JB:JB+NRHS-1) is not modified; if
294 * TRANS = 'N' and EQUED = 'R' or 'B', B is overwritten by
295 * diag(R)*B(IB:IB+N-1,JB:JB+NRHS-1); if TRANS = 'T' or 'C'
296 * and EQUED = 'C' or 'B', B(IB:IB+N-1,JB:JB+NRHS-1) is over-
297 * written by diag(C)*B(IB:IB+N-1,JB:JB+NRHS-1).
298 *
299 * IB (global input) INTEGER
300 * The row index in the global array B indicating the first
301 * row of sub( B ).
302 *
303 * JB (global input) INTEGER
304 * The column index in the global array B indicating the
305 * first column of sub( B ).
306 *
307 * DESCB (global and local input) INTEGER array of dimension DLEN_.
308 * The array descriptor for the distributed matrix B.
309 *
310 * X (local input/local output) COMPLEX pointer
311 * into the local memory to an array of local dimension
312 * (LLD_X, LOCc(JX+NRHS-1)). If INFO = 0, the N-by-NRHS
313 * solution matrix X(IX:IX+N-1,JX:JX+NRHS-1) to the original
314 * system of equations. Note that A(IA:IA+N-1,JA:JA+N-1) and
315 * B(IB:IB+N-1,JB:JB+NRHS-1) are modified on exit if
316 * EQUED .ne. 'N', and the solution to the equilibrated system
317 * is inv(diag(C))*X(IX:IX+N-1,JX:JX+NRHS-1) if TRANS = 'N'
318 * and EQUED = 'C' or 'B', or
319 * inv(diag(R))*X(IX:IX+N-1,JX:JX+NRHS-1) if TRANS = 'T' or 'C'
320 * and EQUED = 'R' or 'B'.
321 *
322 * IX (global input) INTEGER
323 * The row index in the global array X indicating the first
324 * row of sub( X ).
325 *
326 * JX (global input) INTEGER
327 * The column index in the global array X indicating the
328 * first column of sub( X ).
329 *
330 * DESCX (global and local input) INTEGER array of dimension DLEN_.
331 * The array descriptor for the distributed matrix X.
332 *
333 * RCOND (global output) REAL
334 * The estimate of the reciprocal condition number of the matrix
335 * A(IA:IA+N-1,JA:JA+N-1) after equilibration (if done). If
336 * RCOND is less than the machine precision (in particular, if
337 * RCOND = 0), the matrix is singular to working precision.
338 * This condition is indicated by a return code of INFO > 0.
339 *
340 * FERR (local output) REAL array, dimension LOCc(N_B)
341 * The estimated forward error bounds for each solution vector
342 * X(j) (the j-th column of the solution matrix
343 * X(IX:IX+N-1,JX:JX+NRHS-1). If XTRUE is the true solution,
344 * FERR(j) bounds the magnitude of the largest entry in
345 * (X(j) - XTRUE) divided by the magnitude of the largest entry
346 * in X(j). The estimate is as reliable as the estimate for
347 * RCOND, and is almost always a slight overestimate of the
348 * true error. FERR is replicated in every process row, and is
349 * aligned with the matrices B and X.
350 *
351 * BERR (local output) REAL array, dimension LOCc(N_B).
352 * The componentwise relative backward error of each solution
353 * vector X(j) (i.e., the smallest relative change in
354 * any entry of A(IA:IA+N-1,JA:JA+N-1) or
355 * B(IB:IB+N-1,JB:JB+NRHS-1) that makes X(j) an exact solution).
356 * BERR is replicated in every process row, and is aligned
357 * with the matrices B and X.
358 *
359 * WORK (local workspace/local output) COMPLEX array,
360 * dimension (LWORK)
361 * On exit, WORK(1) returns the minimal and optimal LWORK.
362 *
363 * LWORK (local or global input) INTEGER
364 * The dimension of the array WORK.
365 * LWORK is local input and must be at least
366 * LWORK = MAX( PCGECON( LWORK ), PCGERFS( LWORK ) )
367 * + LOCr( N_A ).
368 *
369 * If LWORK = -1, then LWORK is global input and a workspace
370 * query is assumed; the routine only calculates the minimum
371 * and optimal size for all work arrays. Each of these
372 * values is returned in the first entry of the corresponding
373 * work array, and no error message is issued by PXERBLA.
374 *
375 * RWORK (local workspace/local output) REAL array,
376 * dimension (LRWORK)
377 * On exit, RWORK(1) returns the minimal and optimal LRWORK.
378 *
379 * LRWORK (local or global input) INTEGER
380 * The dimension of the array RWORK.
381 * LRWORK is local input and must be at least
382 * LRWORK = 2*LOCc(N_A).
383 *
384 * If LRWORK = -1, then LRWORK is global input and a workspace
385 * query is assumed; the routine only calculates the minimum
386 * and optimal size for all work arrays. Each of these
387 * values is returned in the first entry of the corresponding
388 * work array, and no error message is issued by PXERBLA.
389 *
390 *
391 * INFO (global output) INTEGER
392 * = 0: successful exit
393 * < 0: if INFO = -i, the i-th argument had an illegal value
394 * > 0: if INFO = i, and i is
395 * <= N: U(IA+I-1,IA+I-1) is exactly zero. The
396 * factorization has been completed, but the
397 * factor U is exactly singular, so the solution
398 * and error bounds could not be computed.
399 * = N+1: RCOND is less than machine precision. The
400 * factorization has been completed, but the
401 * matrix is singular to working precision, and
402 * the solution and error bounds have not been
403 * computed.
404 *
405 * =====================================================================
406 *
407 * .. Parameters ..
408  INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
409  $ LLD_, MB_, M_, NB_, N_, RSRC_
410  PARAMETER ( BLOCK_CYCLIC_2D = 1, dlen_ = 9, dtype_ = 1,
411  $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
412  $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
413  REAL ONE, ZERO
414  PARAMETER ( ONE = 1.0e+0, zero = 0.0e+0 )
415 * ..
416 * .. Local Scalars ..
417  LOGICAL COLEQU, EQUIL, LQUERY, NOFACT, NOTRAN, ROWEQU
418  CHARACTER NORM
419  INTEGER CONWRK, I, IACOL, IAROW, IAFROW, IBROW, IBCOL,
420  $ icoffa, icoffb, icoffx, ictxt, idumm,
421  $ iia, iib, iix,
422  $ infequ, iroffa, iroffaf, iroffb,
423  $ iroffx, ixcol, ixrow, j, jja, jjb, jjx,
424  $ lcm, lcmq,
425  $ lrwmin, lwmin, mycol, myrow, np, npcol, nprow,
426  $ nq, nqb, nrhsq, rfswrk
427  REAL AMAX, ANORM, BIGNUM, COLCND, RCMAX, RCMIN,
428  $ ROWCND, SMLNUM
429 * ..
430 * .. Local Arrays ..
431  INTEGER CDESC( DLEN_ ), IDUM1( 5 ), IDUM2( 5 )
432 * ..
433 * .. External Subroutines ..
434  EXTERNAL blacs_gridinfo, chk1mat, descset, pchk2mat,
435  $ infog2l, pcgecon, pcgeequ, pcgerfs,
436  $ pcgetrf, pcgetrs, pclacpy,
437  $ pclaqge, pscopy, pxerbla, sgebr2d,
438  $ sgebs2d, sgamn2d, sgamx2d
439 * ..
440 * .. External Functions ..
441  LOGICAL LSAME
442  INTEGER ICEIL, ILCM, INDXG2P, NUMROC
443  REAL PSLAMCH, PCLANGE
444  EXTERNAL iceil, ilcm, indxg2p, lsame, numroc, pclange,
445  $ pslamch
446 * ..
447 * .. Intrinsic Functions ..
448  INTRINSIC ichar, max, min, mod, real
449 * ..
450 * .. Executable Statements ..
451 *
452 * Get grid parameters
453 *
454  ictxt = desca( ctxt_ )
455  CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
456 *
457 * Test the input parameters
458 *
459  info = 0
460  IF( nprow.EQ.-1 ) THEN
461  info = -(800+ctxt_)
462  ELSE
463  CALL chk1mat( n, 3, n, 3, ia, ja, desca, 8, info )
464  IF( lsame( fact, 'F' ) )
465  $ CALL chk1mat( n, 3, n, 3, iaf, jaf, descaf, 12, info )
466  CALL chk1mat( n, 3, nrhs, 4, ib, jb, descb, 20, info )
467  CALL chk1mat( n, 3, nrhs, 4, ix, jx, descx, 24, info )
468  nofact = lsame( fact, 'N' )
469  equil = lsame( fact, 'E' )
470  notran = lsame( trans, 'N' )
471  IF( nofact .OR. equil ) THEN
472  equed = 'N'
473  rowequ = .false.
474  colequ = .false.
475  ELSE
476  rowequ = lsame( equed, 'R' ) .OR. lsame( equed, 'B' )
477  colequ = lsame( equed, 'C' ) .OR. lsame( equed, 'B' )
478  smlnum = pslamch( ictxt, 'Safe minimum' )
479  bignum = one / smlnum
480  END IF
481  IF( info.EQ.0 ) THEN
482  iarow = indxg2p( ia, desca( mb_ ), myrow, desca( rsrc_ ),
483  $ nprow )
484  iafrow = indxg2p( iaf, descaf( mb_ ), myrow,
485  $ descaf( rsrc_ ), nprow )
486  ibrow = indxg2p( ib, descb( mb_ ), myrow, descb( rsrc_ ),
487  $ nprow )
488  ixrow = indxg2p( ix, descx( mb_ ), myrow, descx( rsrc_ ),
489  $ nprow )
490  iroffa = mod( ia-1, desca( mb_ ) )
491  iroffaf = mod( iaf-1, descaf( mb_ ) )
492  icoffa = mod( ja-1, desca( nb_ ) )
493  iroffb = mod( ib-1, descb( mb_ ) )
494  icoffb = mod( jb-1, descb( nb_ ) )
495  iroffx = mod( ix-1, descx( mb_ ) )
496  icoffx = mod( jx-1, descx( nb_ ) )
497  CALL infog2l( ia, ja, desca, nprow, npcol, myrow,
498  $ mycol, iia, jja, iarow, iacol )
499  np = numroc( n+iroffa, desca( mb_ ), myrow, iarow,
500  $ nprow )
501  IF( myrow.EQ.iarow )
502  $ np = np-iroffa
503  nq = numroc( n+icoffa, desca( nb_ ), mycol, iacol,
504  $ npcol )
505  IF( mycol.EQ.iacol )
506  $ nq = nq-icoffa
507  nqb = iceil( n+iroffa, desca( nb_ )*npcol )
508  lcm = ilcm( nprow, npcol )
509  lcmq = lcm / npcol
510  conwrk = 2*np + 2*nq + max( 2, max( desca( nb_ )*
511  $ max( 1, iceil( nprow-1, npcol ) ), nq +
512  $ desca( nb_ )*
513  $ max( 1, iceil( npcol-1, nprow ) ) ) )
514  rfswrk = 3*np
515  IF( lsame( trans, 'N' ) ) THEN
516  rfswrk = rfswrk + np + nq +
517  $ iceil( nqb, lcmq )*desca( nb_ )
518  ELSE IF( lsame( trans, 'T' ).OR.lsame( trans, 'C' ) ) THEN
519  rfswrk = rfswrk + np + nq
520  END IF
521  lwmin = max( conwrk, rfswrk )
522  lrwmin = max( 2*nq, np )
523  rwork( 1 ) = real( lrwmin )
524  IF( .NOT.nofact .AND. .NOT.equil .AND.
525  $ .NOT.lsame( fact, 'F' ) ) THEN
526  info = -1
527  ELSE IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND.
528  $ .NOT. lsame( trans, 'C' ) ) THEN
529  info = -2
530  ELSE IF( iroffa.NE.0 ) THEN
531  info = -6
532  ELSE IF( icoffa.NE.0 .OR. iroffa.NE.icoffa ) THEN
533  info = -7
534  ELSE IF( desca( mb_ ).NE.desca( nb_ ) ) THEN
535  info = -(800+nb_)
536  ELSE IF( iafrow.NE.iarow ) THEN
537  info = -10
538  ELSE IF( iroffaf.NE.0 ) THEN
539  info = -10
540  ELSE IF( ictxt.NE.descaf( ctxt_ ) ) THEN
541  info = -(1200+ctxt_)
542  ELSE IF( lsame( fact, 'F' ) .AND. .NOT.
543  $ ( rowequ .OR. colequ .OR. lsame( equed, 'N' ) ) ) THEN
544  info = -13
545  ELSE
546  IF( rowequ ) THEN
547  rcmin = bignum
548  rcmax = zero
549  DO 10 j = iia, iia + np - 1
550  rcmin = min( rcmin, r( j ) )
551  rcmax = max( rcmax, r( j ) )
552  10 CONTINUE
553  CALL sgamn2d( ictxt, 'Columnwise', ' ', 1, 1, rcmin,
554  $ 1, idumm, idumm, -1, -1, mycol )
555  CALL sgamx2d( ictxt, 'Columnwise', ' ', 1, 1, rcmax,
556  $ 1, idumm, idumm, -1, -1, mycol )
557  IF( rcmin.LE.zero ) THEN
558  info = -14
559  ELSE IF( n.GT.0 ) THEN
560  rowcnd = max( rcmin, smlnum ) /
561  $ min( rcmax, bignum )
562  ELSE
563  rowcnd = one
564  END IF
565  END IF
566  IF( colequ .AND. info.EQ.0 ) THEN
567  rcmin = bignum
568  rcmax = zero
569  DO 20 j = jja, jja+nq-1
570  rcmin = min( rcmin, c( j ) )
571  rcmax = max( rcmax, c( j ) )
572  20 CONTINUE
573  CALL sgamn2d( ictxt, 'Rowwise', ' ', 1, 1, rcmin,
574  $ 1, idumm, idumm, -1, -1, mycol )
575  CALL sgamx2d( ictxt, 'Rowwise', ' ', 1, 1, rcmax,
576  $ 1, idumm, idumm, -1, -1, mycol )
577  IF( rcmin.LE.zero ) THEN
578  info = -15
579  ELSE IF( n.GT.0 ) THEN
580  colcnd = max( rcmin, smlnum ) /
581  $ min( rcmax, bignum )
582  ELSE
583  colcnd = one
584  END IF
585  END IF
586  END IF
587  END IF
588 *
589  work( 1 ) = real( lwmin )
590  rwork( 1 ) = real( lrwmin )
591  lquery = ( lwork.EQ.-1 .OR. lrwork.EQ.-1 )
592  IF( info.EQ.0 ) THEN
593  IF( ibrow.NE.iarow ) THEN
594  info = -18
595  ELSE IF( ixrow.NE.ibrow ) THEN
596  info = -22
597  ELSE IF( descb( mb_ ).NE.desca( nb_ ) ) THEN
598  info = -(2000+nb_)
599  ELSE IF( ictxt.NE.descb( ctxt_ ) ) THEN
600  info = -(2000+ctxt_)
601  ELSE IF( descx( mb_ ).NE.desca( nb_ ) ) THEN
602  info = -(2400+nb_)
603  ELSE IF( ictxt.NE.descx( ctxt_ ) ) THEN
604  info = -(2400+ctxt_)
605  ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
606  info = -29
607  ELSE IF( lrwork.LT.lrwmin .AND. .NOT.lquery ) THEN
608  info = -31
609  END IF
610  idum1( 1 ) = ichar( fact )
611  idum2( 1 ) = 1
612  idum1( 2 ) = ichar( trans )
613  idum2( 2 ) = 2
614  IF( lsame( fact, 'F' ) ) THEN
615  idum1( 3 ) = ichar( equed )
616  idum2( 3 ) = 14
617  IF( lwork.EQ.-1 ) THEN
618  idum1( 4 ) = -1
619  ELSE
620  idum1( 4 ) = 1
621  END IF
622  idum2( 4 ) = 29
623  IF( lrwork.EQ.-1 ) THEN
624  idum1( 5 ) = -1
625  ELSE
626  idum1( 5 ) = 1
627  END IF
628  idum2( 5 ) = 31
629  CALL pchk2mat( n, 3, n, 3, ia, ja, desca, 8, n, 3,
630  $ nrhs, 4, ib, jb, descb, 20, 5, idum1,
631  $ idum2, info )
632  ELSE
633  IF( lwork.EQ.-1 ) THEN
634  idum1( 3 ) = -1
635  ELSE
636  idum1( 3 ) = 1
637  END IF
638  idum2( 3 ) = 29
639  IF( lrwork.EQ.-1 ) THEN
640  idum1( 4 ) = -1
641  ELSE
642  idum1( 4 ) = 1
643  END IF
644  idum2( 4 ) = 31
645  CALL pchk2mat( n, 3, n, 3, ia, ja, desca, 8, n, 3,
646  $ nrhs, 4, ib, jb, descb, 20, 4, idum1,
647  $ idum2, info )
648  END IF
649  END IF
650  END IF
651 *
652  IF( info.NE.0 ) THEN
653  CALL pxerbla( ictxt, 'PCGESVX', -info )
654  RETURN
655  ELSE IF( lquery ) THEN
656  RETURN
657  END IF
658 *
659  IF( equil ) THEN
660 *
661 * Compute row and column scalings to equilibrate the matrix A.
662 *
663  CALL pcgeequ( n, n, a, ia, ja, desca, r, c, rowcnd, colcnd,
664  $ amax, infequ )
665  IF( infequ.EQ.0 ) THEN
666 *
667 * Equilibrate the matrix.
668 *
669  CALL pclaqge( n, n, a, ia, ja, desca, r, c, rowcnd, colcnd,
670  $ amax, equed )
671  rowequ = lsame( equed, 'R' ) .OR. lsame( equed, 'B' )
672  colequ = lsame( equed, 'C' ) .OR. lsame( equed, 'B' )
673  END IF
674  END IF
675 *
676 * Scale the right-hand side.
677 *
678  CALL infog2l( ib, jb, descb, nprow, npcol, myrow, mycol, iib,
679  $ jjb, ibrow, ibcol )
680  np = numroc( n+iroffb, descb( mb_ ), myrow, ibrow, nprow )
681  nrhsq = numroc( nrhs+icoffb, descb( nb_ ), mycol, ibcol, npcol )
682  IF( myrow.EQ.ibrow )
683  $ np = np-iroffb
684  IF( mycol.EQ.ibcol )
685  $ nrhsq = nrhsq-icoffb
686 *
687  IF( notran ) THEN
688  IF( rowequ ) THEN
689  DO 40 j = jjb, jjb+nrhsq-1
690  DO 30 i = iib, iib+np-1
691  b( i+( j-1 )*descb( lld_ ) ) = r( i )*
692  $ b( i+( j-1 )*descb( lld_ ) )
693  30 CONTINUE
694  40 CONTINUE
695  END IF
696  ELSE IF( colequ ) THEN
697 *
698 * Transpose the Column scale factors
699 *
700  CALL descset( cdesc, 1, n+icoffa, 1, desca( nb_ ), myrow,
701  $ iacol, ictxt, 1 )
702  CALL pscopy( n, c, 1, ja, cdesc, cdesc( lld_ ), rwork, ib, jb,
703  $ descb, 1 )
704  IF( mycol.EQ.ibcol ) THEN
705  CALL sgebs2d( ictxt, 'Rowwise', ' ', np, 1, rwork( iib ),
706  $ descb( lld_ ) )
707  ELSE
708  CALL sgebr2d( ictxt, 'Rowwise', ' ', np, 1, rwork( iib ),
709  $ descb( lld_ ), myrow, ibcol )
710  END IF
711  DO 60 j = jjb, jjb+nrhsq-1
712  DO 50 i = iib, iib+np-1
713  b( i+( j-1 )*descb( lld_ ) ) = rwork( i )*
714  $ b( i+( j-1 )*descb( lld_ ) )
715  50 CONTINUE
716  60 CONTINUE
717  END IF
718 *
719  IF( nofact.OR.equil ) THEN
720 *
721 * Compute the LU factorization of A.
722 *
723  CALL pclacpy( 'Full', n, n, a, ia, ja, desca, af, iaf, jaf,
724  $ descaf )
725  CALL pcgetrf( n, n, af, iaf, jaf, descaf, ipiv, info )
726 *
727 * Return if INFO is non-zero.
728 *
729  IF( info.NE.0 ) THEN
730  IF( info.GT.0 )
731  $ rcond = zero
732  RETURN
733  END IF
734  END IF
735 *
736 * Compute the norm of the matrix A.
737 *
738  IF( notran ) THEN
739  norm = '1'
740  ELSE
741  norm = 'I'
742  END IF
743  anorm = pclange( norm, n, n, a, ia, ja, desca, rwork )
744 *
745 * Compute the reciprocal of the condition number of A.
746 *
747  CALL pcgecon( norm, n, af, iaf, jaf, descaf, anorm, rcond, work,
748  $ lwork, rwork, lrwork, info )
749 *
750 * Return if the matrix is singular to working precision.
751 *
752  IF( rcond.LT.pslamch( ictxt, 'Epsilon' ) ) THEN
753  info = ia + n
754  RETURN
755  END IF
756 *
757 * Compute the solution matrix X.
758 *
759  CALL pclacpy( 'Full', n, nrhs, b, ib, jb, descb, x, ix, jx,
760  $ descx )
761  CALL pcgetrs( trans, n, nrhs, af, iaf, jaf, descaf, ipiv, x, ix,
762  $ jx, descx, info )
763 *
764 * Use iterative refinement to improve the computed solution and
765 * compute error bounds and backward error estimates for it.
766 *
767  CALL pcgerfs( trans, n, nrhs, a, ia, ja, desca, af, iaf, jaf,
768  $ descaf, ipiv, b, ib, jb, descb, x, ix, jx, descx,
769  $ ferr, berr, work, lwork, rwork, lrwork, info )
770 *
771 * Transform the solution matrix X to a solution of the original
772 * system.
773 *
774  CALL infog2l( ix, jx, descx, nprow, npcol, myrow, mycol, iix,
775  $ jjx, ixrow, ixcol )
776  np = numroc( n+iroffx, descx( mb_ ), myrow, ixrow, nprow )
777  nrhsq = numroc( nrhs+icoffx, descx( nb_ ), mycol, ixcol, npcol )
778  IF( myrow.EQ.ibrow )
779  $ np = np-iroffx
780  IF( mycol.EQ.ibcol )
781  $ nrhsq = nrhsq-icoffx
782 *
783  IF( notran ) THEN
784  IF( colequ ) THEN
785 *
786 * Transpose the column scaling factors
787 *
788  CALL descset( cdesc, 1, n+icoffa, 1, desca( nb_ ), myrow,
789  $ iacol, ictxt, 1 )
790  CALL pscopy( n, c, 1, ja, cdesc, cdesc( lld_ ), rwork, ix,
791  $ jx, descx, 1 )
792  IF( mycol.EQ.ibcol ) THEN
793  CALL sgebs2d( ictxt, 'Rowwise', ' ', np, 1,
794  $ rwork( iix ), descx( lld_ ) )
795  ELSE
796  CALL sgebr2d( ictxt, 'Rowwise', ' ', np, 1,
797  $ rwork( iix ), descx( lld_ ), myrow,
798  $ ibcol )
799  END IF
800 *
801  DO 80 j = jjx, jjx+nrhsq-1
802  DO 70 i = iix, iix+np-1
803  x( i+( j-1 )*descx( lld_ ) ) = rwork( i )*
804  $ x( i+( j-1 )*descx( lld_ ) )
805  70 CONTINUE
806  80 CONTINUE
807  DO 90 j = jjx, jjx+nrhsq-1
808  ferr( j ) = ferr( j ) / colcnd
809  90 CONTINUE
810  END IF
811  ELSE IF( rowequ ) THEN
812  DO 110 j = jjx, jjx+nrhsq-1
813  DO 100 i = iix, iix+np-1
814  x( i+( j-1 )*descx( lld_ ) ) = r( i )*
815  $ x( i+( j-1 )*descx( lld_ ) )
816  100 CONTINUE
817  110 CONTINUE
818  DO 120 j = jjx, jjx+nrhsq-1
819  ferr( j ) = ferr( j ) / rowcnd
820  120 CONTINUE
821  END IF
822 *
823  work( 1 ) = real( lwmin )
824  rwork( 1 ) = real( lrwmin )
825 *
826  RETURN
827 *
828 * End of PCGESVX
829 *
830  END
max
#define max(A, B)
Definition: pcgemr.c:180
pcgerfs
subroutine pcgerfs(TRANS, N, NRHS, A, IA, JA, DESCA, AF, IAF, JAF, DESCAF, IPIV, B, IB, JB, DESCB, X, IX, JX, DESCX, FERR, BERR, WORK, LWORK, RWORK, LRWORK, INFO)
Definition: pcgerfs.f:5
pcgetrf
subroutine pcgetrf(M, N, A, IA, JA, DESCA, IPIV, INFO)
Definition: pcgetrf.f:2
pcgeequ
subroutine pcgeequ(M, N, A, IA, JA, DESCA, R, C, ROWCND, COLCND, AMAX, INFO)
Definition: pcgeequ.f:3
infog2l
subroutine infog2l(GRINDX, GCINDX, DESC, NPROW, NPCOL, MYROW, MYCOL, LRINDX, LCINDX, RSRC, CSRC)
Definition: infog2l.f:3
pcgetrs
subroutine pcgetrs(TRANS, N, NRHS, A, IA, JA, DESCA, IPIV, B, IB, JB, DESCB, INFO)
Definition: pcgetrs.f:3
pchk2mat
subroutine pchk2mat(MA, MAPOS0, NA, NAPOS0, IA, JA, DESCA, DESCAPOS0, MB, MBPOS0, NB, NBPOS0, IB, JB, DESCB, DESCBPOS0, NEXTRA, EX, EXPOS, INFO)
Definition: pchkxmat.f:175
descset
subroutine descset(DESC, M, N, MB, NB, IRSRC, ICSRC, ICTXT, LLD)
Definition: descset.f:3
pclacpy
subroutine pclacpy(UPLO, M, N, A, IA, JA, DESCA, B, IB, JB, DESCB)
Definition: pclacpy.f:3
pclaqge
subroutine pclaqge(M, N, A, IA, JA, DESCA, R, C, ROWCND, COLCND, AMAX, EQUED)
Definition: pclaqge.f:3
pcgesvx
subroutine pcgesvx(FACT, TRANS, N, NRHS, A, IA, JA, DESCA, AF, IAF, JAF, DESCAF, IPIV, EQUED, R, C, B, IB, JB, DESCB, X, IX, JX, DESCX, RCOND, FERR, BERR, WORK, LWORK, RWORK, LRWORK, INFO)
Definition: pcgesvx.f:5
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
pcgecon
subroutine pcgecon(NORM, N, A, IA, JA, DESCA, ANORM, RCOND, WORK, LWORK, RWORK, LRWORK, INFO)
Definition: pcgecon.f:3
min
#define min(A, B)
Definition: pcgemr.c:181