ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pdgesvx.f
Go to the documentation of this file.
1  SUBROUTINE pdgesvx( 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, IWORK, LIWORK, 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, LIWORK,
14  $ LWORK, N, NRHS
15  DOUBLE PRECISION RCOND
16 * ..
17 * .. Array Arguments ..
18  INTEGER DESCA( * ), DESCAF( * ), DESCB( * ),
19  $ DESCX( * ), IPIV( * ), IWORK( * )
20  DOUBLE PRECISION A( * ), AF( * ), B( * ), BERR( * ), C( * ),
21  $ ferr( * ), r( * ), work( * ), x( * )
22 * ..
23 *
24 * Purpose
25 * =======
26 *
27 * PDGESVX uses the LU factorization to compute the solution to a real
28 * system of linear equations
29 *
30 * A(IA:IA+N-1,JA:JA+N-1) * X = B(IB:IB+N-1,JB:JB+NRHS-1),
31 *
32 * where A(IA:IA+N-1,JA:JA+N-1) is an N-by-N matrix and X and
33 * B(IB:IB+N-1,JB:JB+NRHS-1) are N-by-NRHS matrices.
34 *
35 * Error bounds on the solution and a condition estimate are also
36 * provided.
37 *
38 * Notes
39 * =====
40 *
41 * Each global data object is described by an associated description
42 * vector. This vector stores the information required to establish
43 * the mapping between an object element and its corresponding process
44 * and memory location.
45 *
46 * Let A be a generic term for any 2D block cyclicly distributed array.
47 * Such a global array has an associated description vector DESCA.
48 * In the following comments, the character _ should be read as
49 * "of the global array".
50 *
51 * NOTATION STORED IN EXPLANATION
52 * --------------- -------------- --------------------------------------
53 * DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
54 * DTYPE_A = 1.
55 * CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
56 * the BLACS process grid A is distribu-
57 * ted over. The context itself is glo-
58 * bal, but the handle (the integer
59 * value) may vary.
60 * M_A (global) DESCA( M_ ) The number of rows in the global
61 * array A.
62 * N_A (global) DESCA( N_ ) The number of columns in the global
63 * array A.
64 * MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
65 * the rows of the array.
66 * NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
67 * the columns of the array.
68 * RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
69 * row of the array A is distributed.
70 * CSRC_A (global) DESCA( CSRC_ ) The process column over which the
71 * first column of the array A is
72 * distributed.
73 * LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
74 * array. LLD_A >= MAX(1,LOCr(M_A)).
75 *
76 * Let K be the number of rows or columns of a distributed matrix,
77 * and assume that its process grid has dimension p x q.
78 * LOCr( K ) denotes the number of elements of K that a process
79 * would receive if K were distributed over the p processes of its
80 * process column.
81 * Similarly, LOCc( K ) denotes the number of elements of K that a
82 * process would receive if K were distributed over the q processes of
83 * its process row.
84 * The values of LOCr() and LOCc() may be determined via a call to the
85 * ScaLAPACK tool function, NUMROC:
86 * LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
87 * LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
88 * An upper bound for these quantities may be computed by:
89 * LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
90 * LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
91 *
92 * Description
93 * ===========
94 *
95 * In the following description, A denotes A(IA:IA+N-1,JA:JA+N-1),
96 * B denotes B(IB:IB+N-1,JB:JB+NRHS-1) and X denotes
97 * X(IX:IX+N-1,JX:JX+NRHS-1).
98 *
99 * The following steps are performed:
100 *
101 * 1. If FACT = 'E', real scaling factors are computed to equilibrate
102 * the system:
103 * TRANS = 'N': diag(R)*A*diag(C) *inv(diag(C))*X = diag(R)*B
104 * TRANS = 'T': (diag(R)*A*diag(C))**T *inv(diag(R))*X = diag(C)*B
105 * TRANS = 'C': (diag(R)*A*diag(C))**H *inv(diag(R))*X = diag(C)*B
106 * Whether or not the system will be equilibrated depends on the
107 * scaling of the matrix A, but if equilibration is used, A is
108 * overwritten by diag(R)*A*diag(C) and B by diag(R)*B (if TRANS='N')
109 * or diag(C)*B (if TRANS = 'T' or 'C').
110 *
111 * 2. If FACT = 'N' or 'E', the LU decomposition is used to factor the
112 * matrix A (after equilibration if FACT = 'E') as
113 * A = P * L * U,
114 * where P is a permutation matrix, L is a unit lower triangular
115 * matrix, and U is upper triangular.
116 *
117 * 3. The factored form of A is used to estimate the condition number
118 * of the matrix A. If the reciprocal of the condition number is
119 * less than machine precision, steps 4-6 are skipped.
120 *
121 * 4. The system of equations is solved for X using the factored form
122 * of A.
123 *
124 * 5. Iterative refinement is applied to improve the computed solution
125 * matrix and calculate error bounds and backward error estimates
126 * for it.
127 *
128 * 6. If FACT = 'E' and equilibration was used, the matrix X is
129 * premultiplied by diag(C) (if TRANS = 'N') or diag(R) (if
130 * TRANS = 'T' or 'C') so that it solves the original system
131 * before equilibration.
132 *
133 * Arguments
134 * =========
135 *
136 * FACT (global input) CHARACTER
137 * Specifies whether or not the factored form of the matrix
138 * A(IA:IA+N-1,JA:JA+N-1) is supplied on entry, and if not,
139 * whether the matrix A(IA:IA+N-1,JA:JA+N-1) should be
140 * equilibrated before it is factored.
141 * = 'F': On entry, AF(IAF:IAF+N-1,JAF:JAF+N-1) and IPIV con-
142 * tain the factored form of A(IA:IA+N-1,JA:JA+N-1).
143 * If EQUED is not 'N', the matrix
144 * A(IA:IA+N-1,JA:JA+N-1) has been equilibrated with
145 * scaling factors given by R and C.
146 * A(IA:IA+N-1,JA:JA+N-1), AF(IAF:IAF+N-1,JAF:JAF+N-1),
147 * and IPIV are not modified.
148 * = 'N': The matrix A(IA:IA+N-1,JA:JA+N-1) will be copied to
149 * AF(IAF:IAF+N-1,JAF:JAF+N-1) and factored.
150 * = 'E': The matrix A(IA:IA+N-1,JA:JA+N-1) will be equili-
151 * brated if necessary, then copied to
152 * AF(IAF:IAF+N-1,JAF:JAF+N-1) and factored.
153 *
154 * TRANS (global input) CHARACTER
155 * Specifies the form of the system of equations:
156 * = 'N': A(IA:IA+N-1,JA:JA+N-1) * X(IX:IX+N-1,JX:JX+NRHS-1)
157 * = B(IB:IB+N-1,JB:JB+NRHS-1) (No transpose)
158 * = 'T': A(IA:IA+N-1,JA:JA+N-1)**T * X(IX:IX+N-1,JX:JX+NRHS-1)
159 * = B(IB:IB+N-1,JB:JB+NRHS-1) (Transpose)
160 * = 'C': A(IA:IA+N-1,JA:JA+N-1)**H * X(IX:IX+N-1,JX:JX+NRHS-1)
161 * = B(IB:IB+N-1,JB:JB+NRHS-1) (Transpose)
162 *
163 * N (global input) INTEGER
164 * The number of rows and columns to be operated on, i.e. the
165 * order of the distributed submatrix A(IA:IA+N-1,JA:JA+N-1).
166 * N >= 0.
167 *
168 * NRHS (global input) INTEGER
169 * The number of right-hand sides, i.e., the number of columns
170 * of the distributed submatrices B(IB:IB+N-1,JB:JB+NRHS-1) and
171 * X(IX:IX+N-1,JX:JX+NRHS-1). NRHS >= 0.
172 *
173 * A (local input/local output) DOUBLE PRECISION pointer into
174 * the local memory to an array of local dimension
175 * (LLD_A,LOCc(JA+N-1)). On entry, the N-by-N matrix
176 * A(IA:IA+N-1,JA:JA+N-1). If FACT = 'F' and EQUED is not 'N',
177 * then A(IA:IA+N-1,JA:JA+N-1) must have been equilibrated by
178 * the scaling factors in R and/or C. A(IA:IA+N-1,JA:JA+N-1) is
179 * not modified if FACT = 'F' or 'N', or if FACT = 'E' and
180 * EQUED = 'N' on exit.
181 *
182 * On exit, if EQUED .ne. 'N', A(IA:IA+N-1,JA:JA+N-1) is scaled
183 * as follows:
184 * EQUED = 'R': A(IA:IA+N-1,JA:JA+N-1) :=
185 * diag(R) * A(IA:IA+N-1,JA:JA+N-1)
186 * EQUED = 'C': A(IA:IA+N-1,JA:JA+N-1) :=
187 * A(IA:IA+N-1,JA:JA+N-1) * diag(C)
188 * EQUED = 'B': A(IA:IA+N-1,JA:JA+N-1) :=
189 * diag(R) * A(IA:IA+N-1,JA:JA+N-1) * diag(C).
190 *
191 * IA (global input) INTEGER
192 * The row index in the global array A indicating the first
193 * row of sub( A ).
194 *
195 * JA (global input) INTEGER
196 * The column index in the global array A indicating the
197 * first column of sub( A ).
198 *
199 * DESCA (global and local input) INTEGER array of dimension DLEN_.
200 * The array descriptor for the distributed matrix A.
201 *
202 * AF (local input or local output) DOUBLE PRECISION pointer
203 * into the local memory to an array of local dimension
204 * (LLD_AF,LOCc(JA+N-1)). If FACT = 'F', then
205 * AF(IAF:IAF+N-1,JAF:JAF+N-1) is an input argument and on
206 * entry contains the factors L and U from the factorization
207 * A(IA:IA+N-1,JA:JA+N-1) = P*L*U as computed by PDGETRF.
208 * If EQUED .ne. 'N', then AF is the factored form of the
209 * equilibrated matrix A(IA:IA+N-1,JA:JA+N-1).
210 *
211 * If FACT = 'N', then AF(IAF:IAF+N-1,JAF:JAF+N-1) is an output
212 * argument and on exit returns the factors L and U from the
213 * factorization A(IA:IA+N-1,JA:JA+N-1) = P*L*U of the original
214 * matrix A(IA:IA+N-1,JA:JA+N-1).
215 *
216 * If FACT = 'E', then AF(IAF:IAF+N-1,JAF:JAF+N-1) is an output
217 * argument and on exit returns the factors L and U from the
218 * factorization A(IA:IA+N-1,JA:JA+N-1) = P*L*U of the equili-
219 * brated matrix A(IA:IA+N-1,JA:JA+N-1) (see the description of
220 * A(IA:IA+N-1,JA:JA+N-1) for the form of the equilibrated
221 * matrix).
222 *
223 * IAF (global input) INTEGER
224 * The row index in the global array AF indicating the first
225 * row of sub( AF ).
226 *
227 * JAF (global input) INTEGER
228 * The column index in the global array AF indicating the
229 * first column of sub( AF ).
230 *
231 * DESCAF (global and local input) INTEGER array of dimension DLEN_.
232 * The array descriptor for the distributed matrix AF.
233 *
234 * IPIV (local input or local output) INTEGER array, dimension
235 * LOCr(M_A)+MB_A. If FACT = 'F', then IPIV is an input argu-
236 * ment and on entry contains the pivot indices from the fac-
237 * torization A(IA:IA+N-1,JA:JA+N-1) = P*L*U as computed by
238 * PDGETRF; IPIV(i) -> The global row local row i was
239 * swapped with. This array must be aligned with
240 * A( IA:IA+N-1, * ).
241 *
242 * If FACT = 'N', then IPIV is an output argument and on exit
243 * contains the pivot indices from the factorization
244 * A(IA:IA+N-1,JA:JA+N-1) = P*L*U of the original matrix
245 * A(IA:IA+N-1,JA:JA+N-1).
246 *
247 * If FACT = 'E', then IPIV is an output argument and on exit
248 * contains the pivot indices from the factorization
249 * A(IA:IA+N-1,JA:JA+N-1) = P*L*U of the equilibrated matrix
250 * A(IA:IA+N-1,JA:JA+N-1).
251 *
252 * EQUED (global input or global output) CHARACTER
253 * Specifies the form of equilibration that was done.
254 * = 'N': No equilibration (always true if FACT = 'N').
255 * = 'R': Row equilibration, i.e., A(IA:IA+N-1,JA:JA+N-1) has
256 * been premultiplied by diag(R).
257 * = 'C': Column equilibration, i.e., A(IA:IA+N-1,JA:JA+N-1)
258 * has been postmultiplied by diag(C).
259 * = 'B': Both row and column equilibration, i.e.,
260 * A(IA:IA+N-1,JA:JA+N-1) has been replaced by
261 * diag(R) * A(IA:IA+N-1,JA:JA+N-1) * diag(C).
262 * EQUED is an input variable if FACT = 'F'; otherwise, it is an
263 * output variable.
264 *
265 * R (local input or local output) DOUBLE PRECISION array,
266 * dimension LOCr(M_A).
267 * The row scale factors for A(IA:IA+N-1,JA:JA+N-1).
268 * If EQUED = 'R' or 'B', A(IA:IA+N-1,JA:JA+N-1) is multiplied
269 * on the left by diag(R); if EQUED='N' or 'C', R is not acces-
270 * sed. R is an input variable if FACT = 'F'; otherwise, R is
271 * an output variable.
272 * If FACT = 'F' and EQUED = 'R' or 'B', each element of R must
273 * be positive.
274 * R is replicated in every process column, and is aligned
275 * with the distributed matrix A.
276 *
277 * C (local input or local output) DOUBLE PRECISION array,
278 * dimension LOCc(N_A).
279 * The column scale factors for A(IA:IA+N-1,JA:JA+N-1).
280 * If EQUED = 'C' or 'B', A(IA:IA+N-1,JA:JA+N-1) is multiplied
281 * on the right by diag(C); if EQUED = 'N' or 'R', C is not
282 * accessed. C is an input variable if FACT = 'F'; otherwise,
283 * C is an output variable. If FACT = 'F' and EQUED = 'C' or
284 * 'B', each element of C must be positive.
285 * C is replicated in every process row, and is aligned with
286 * the distributed matrix A.
287 *
288 * B (local input/local output) DOUBLE PRECISION pointer
289 * into the local memory to an array of local dimension
290 * (LLD_B,LOCc(JB+NRHS-1) ). On entry, the N-by-NRHS right-hand
291 * side matrix B(IB:IB+N-1,JB:JB+NRHS-1). On exit, if
292 * EQUED = 'N', B(IB:IB+N-1,JB:JB+NRHS-1) is not modified; if
293 * TRANS = 'N' and EQUED = 'R' or 'B', B is overwritten by
294 * diag(R)*B(IB:IB+N-1,JB:JB+NRHS-1); if TRANS = 'T' or 'C'
295 * and EQUED = 'C' or 'B', B(IB:IB+N-1,JB:JB+NRHS-1) is over-
296 * written by diag(C)*B(IB:IB+N-1,JB:JB+NRHS-1).
297 *
298 * IB (global input) INTEGER
299 * The row index in the global array B indicating the first
300 * row of sub( B ).
301 *
302 * JB (global input) INTEGER
303 * The column index in the global array B indicating the
304 * first column of sub( B ).
305 *
306 * DESCB (global and local input) INTEGER array of dimension DLEN_.
307 * The array descriptor for the distributed matrix B.
308 *
309 * X (local input/local output) DOUBLE PRECISION pointer
310 * into the local memory to an array of local dimension
311 * (LLD_X, LOCc(JX+NRHS-1)). If INFO = 0, the N-by-NRHS
312 * solution matrix X(IX:IX+N-1,JX:JX+NRHS-1) to the original
313 * system of equations. Note that A(IA:IA+N-1,JA:JA+N-1) and
314 * B(IB:IB+N-1,JB:JB+NRHS-1) are modified on exit if
315 * EQUED .ne. 'N', and the solution to the equilibrated system
316 * is inv(diag(C))*X(IX:IX+N-1,JX:JX+NRHS-1) if TRANS = 'N'
317 * and EQUED = 'C' or 'B', or
318 * inv(diag(R))*X(IX:IX+N-1,JX:JX+NRHS-1) if TRANS = 'T' or 'C'
319 * and EQUED = 'R' or 'B'.
320 *
321 * IX (global input) INTEGER
322 * The row index in the global array X indicating the first
323 * row of sub( X ).
324 *
325 * JX (global input) INTEGER
326 * The column index in the global array X indicating the
327 * first column of sub( X ).
328 *
329 * DESCX (global and local input) INTEGER array of dimension DLEN_.
330 * The array descriptor for the distributed matrix X.
331 *
332 * RCOND (global output) DOUBLE PRECISION
333 * The estimate of the reciprocal condition number of the matrix
334 * A(IA:IA+N-1,JA:JA+N-1) after equilibration (if done). If
335 * RCOND is less than the machine precision (in particular, if
336 * RCOND = 0), the matrix is singular to working precision.
337 * This condition is indicated by a return code of INFO > 0.
338 *
339 * FERR (local output) DOUBLE PRECISION array, dimension LOCc(N_B)
340 * The estimated forward error bounds for each solution vector
341 * X(j) (the j-th column of the solution matrix
342 * X(IX:IX+N-1,JX:JX+NRHS-1). If XTRUE is the true solution,
343 * FERR(j) bounds the magnitude of the largest entry in
344 * (X(j) - XTRUE) divided by the magnitude of the largest entry
345 * in X(j). The estimate is as reliable as the estimate for
346 * RCOND, and is almost always a slight overestimate of the
347 * true error. FERR is replicated in every process row, and is
348 * aligned with the matrices B and X.
349 *
350 * BERR (local output) DOUBLE PRECISION array, dimension LOCc(N_B).
351 * The componentwise relative backward error of each solution
352 * vector X(j) (i.e., the smallest relative change in
353 * any entry of A(IA:IA+N-1,JA:JA+N-1) or
354 * B(IB:IB+N-1,JB:JB+NRHS-1) that makes X(j) an exact solution).
355 * BERR is replicated in every process row, and is aligned
356 * with the matrices B and X.
357 *
358 * WORK (local workspace/local output) DOUBLE PRECISION array,
359 * dimension (LWORK)
360 * On exit, WORK(1) returns the minimal and optimal LWORK.
361 *
362 * LWORK (local or global input) INTEGER
363 * The dimension of the array WORK.
364 * LWORK is local input and must be at least
365 * LWORK = MAX( PDGECON( LWORK ), PDGERFS( LWORK ) )
366 * + LOCr( N_A ).
367 *
368 * If LWORK = -1, then LWORK is global input and a workspace
369 * query is assumed; the routine only calculates the minimum
370 * and optimal size for all work arrays. Each of these
371 * values is returned in the first entry of the corresponding
372 * work array, and no error message is issued by PXERBLA.
373 *
374 * IWORK (local workspace/local output) INTEGER array,
375 * dimension (LIWORK)
376 * On exit, IWORK(1) returns the minimal and optimal LIWORK.
377 *
378 * LIWORK (local or global input) INTEGER
379 * The dimension of the array IWORK.
380 * LIWORK is local input and must be at least
381 * LIWORK = LOCr(N_A).
382 *
383 * If LIWORK = -1, then LIWORK is global input and a workspace
384 * query is assumed; the routine only calculates the minimum
385 * and optimal size for all work arrays. Each of these
386 * values is returned in the first entry of the corresponding
387 * work array, and no error message is issued by PXERBLA.
388 *
389 *
390 * INFO (global output) INTEGER
391 * = 0: successful exit
392 * < 0: if INFO = -i, the i-th argument had an illegal value
393 * > 0: if INFO = i, and i is
394 * <= N: U(IA+I-1,IA+I-1) is exactly zero. The
395 * factorization has been completed, but the
396 * factor U is exactly singular, so the solution
397 * and error bounds could not be computed.
398 * = N+1: RCOND is less than machine precision. The
399 * factorization has been completed, but the
400 * matrix is singular to working precision, and
401 * the solution and error bounds have not been
402 * computed.
403 *
404 * =====================================================================
405 *
406 * .. Parameters ..
407  INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
408  $ LLD_, MB_, M_, NB_, N_, RSRC_
409  PARAMETER ( BLOCK_CYCLIC_2D = 1, dlen_ = 9, dtype_ = 1,
410  $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
411  $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
412  DOUBLE PRECISION ONE, ZERO
413  PARAMETER ( ONE = 1.0d+0, zero = 0.0d+0 )
414 * ..
415 * .. Local Scalars ..
416  LOGICAL COLEQU, EQUIL, LQUERY, NOFACT, NOTRAN, ROWEQU
417  CHARACTER NORM
418  INTEGER CONWRK, I, IACOL, IAROW, IAFROW, IBROW, IBCOL,
419  $ icoffa, icoffb, icoffx, ictxt, idumm,
420  $ iia, iib, iix,
421  $ infequ, iroffa, iroffaf, iroffb,
422  $ iroffx, ixcol, ixrow, j, jja, jjb, jjx,
423  $ lcm, lcmq,
424  $ liwmin, lwmin, mycol, myrow, np, npcol, nprow,
425  $ nq, nqb, nrhsq, rfswrk
426  DOUBLE PRECISION AMAX, ANORM, BIGNUM, COLCND, RCMAX, RCMIN,
427  $ ROWCND, SMLNUM
428 * ..
429 * .. Local Arrays ..
430  INTEGER CDESC( DLEN_ ), IDUM1( 5 ), IDUM2( 5 )
431 * ..
432 * .. External Subroutines ..
433  EXTERNAL blacs_gridinfo, chk1mat, descset, pchk2mat,
434  $ dgebr2d, dgebs2d, dgamn2d,
435  $ dgamx2d, infog2l, pdcopy, pdgecon,
436  $ pdgeequ, pdgerfs, pdgetrf, pdgetrs,
437  $ pdlacpy, pdlaqge, pxerbla
438 * ..
439 * .. External Functions ..
440  LOGICAL LSAME
441  INTEGER ICEIL, ILCM, INDXG2P, NUMROC
442  DOUBLE PRECISION PDLAMCH, PDLANGE
443  EXTERNAL iceil, ilcm, indxg2p, lsame, numroc, pdlange,
444  $ pdlamch
445 * ..
446 * .. Intrinsic Functions ..
447  INTRINSIC dble, ichar, max, min, mod
448 * ..
449 * .. Executable Statements ..
450 *
451 * Get grid parameters
452 *
453  ictxt = desca( ctxt_ )
454  CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
455 *
456 * Test the input parameters
457 *
458  info = 0
459  IF( nprow.EQ.-1 ) THEN
460  info = -(800+ctxt_)
461  ELSE
462  CALL chk1mat( n, 3, n, 3, ia, ja, desca, 8, info )
463  IF( lsame( fact, 'F' ) )
464  $ CALL chk1mat( n, 3, n, 3, iaf, jaf, descaf, 12, info )
465  CALL chk1mat( n, 3, nrhs, 4, ib, jb, descb, 20, info )
466  CALL chk1mat( n, 3, nrhs, 4, ix, jx, descx, 24, info )
467  nofact = lsame( fact, 'N' )
468  equil = lsame( fact, 'E' )
469  notran = lsame( trans, 'N' )
470  IF( nofact .OR. equil ) THEN
471  equed = 'N'
472  rowequ = .false.
473  colequ = .false.
474  ELSE
475  rowequ = lsame( equed, 'R' ) .OR. lsame( equed, 'B' )
476  colequ = lsame( equed, 'C' ) .OR. lsame( equed, 'B' )
477  smlnum = pdlamch( ictxt, 'Safe minimum' )
478  bignum = one / smlnum
479  END IF
480  IF( info.EQ.0 ) THEN
481  iarow = indxg2p( ia, desca( mb_ ), myrow, desca( rsrc_ ),
482  $ nprow )
483  iafrow = indxg2p( iaf, descaf( mb_ ), myrow,
484  $ descaf( rsrc_ ), nprow )
485  ibrow = indxg2p( ib, descb( mb_ ), myrow, descb( rsrc_ ),
486  $ nprow )
487  ixrow = indxg2p( ix, descx( mb_ ), myrow, descx( rsrc_ ),
488  $ nprow )
489  iroffa = mod( ia-1, desca( mb_ ) )
490  iroffaf = mod( iaf-1, descaf( mb_ ) )
491  icoffa = mod( ja-1, desca( nb_ ) )
492  iroffb = mod( ib-1, descb( mb_ ) )
493  icoffb = mod( jb-1, descb( nb_ ) )
494  iroffx = mod( ix-1, descx( mb_ ) )
495  icoffx = mod( jx-1, descx( nb_ ) )
496  CALL infog2l( ia, ja, desca, nprow, npcol, myrow,
497  $ mycol, iia, jja, iarow, iacol )
498  np = numroc( n+iroffa, desca( mb_ ), myrow, iarow,
499  $ nprow )
500  IF( myrow.EQ.iarow )
501  $ np = np-iroffa
502  nq = numroc( n+icoffa, desca( nb_ ), mycol, iacol,
503  $ npcol )
504  IF( mycol.EQ.iacol )
505  $ nq = nq-icoffa
506  nqb = iceil( n+iroffa, desca( nb_ )*npcol )
507  lcm = ilcm( nprow, npcol )
508  lcmq = lcm / npcol
509  conwrk = 2*np + 2*nq + max( 2, max( desca( nb_ )*
510  $ max( 1, iceil( nprow-1, npcol ) ), nq +
511  $ desca( nb_ )*
512  $ max( 1, iceil( npcol-1, nprow ) ) ) )
513  rfswrk = 3*np
514  IF( lsame( trans, 'N' ) ) THEN
515  rfswrk = rfswrk + np + nq +
516  $ iceil( nqb, lcmq )*desca( nb_ )
517  ELSE IF( lsame( trans, 'T' ).OR.lsame( trans, 'C' ) ) THEN
518  rfswrk = rfswrk + np + nq
519  END IF
520  lwmin = max( conwrk, rfswrk )
521  work( 1 ) = dble( lwmin )
522  liwmin = np
523  iwork( 1 ) = liwmin
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 dgamn2d( ictxt, 'Columnwise', ' ', 1, 1, rcmin,
554  $ 1, idumm, idumm, -1, -1, mycol )
555  CALL dgamx2d( 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 dgamn2d( ictxt, 'Rowwise', ' ', 1, 1, rcmin,
574  $ 1, idumm, idumm, -1, -1, mycol )
575  CALL dgamx2d( 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 ) = dble( lwmin )
590  iwork( 1 ) = liwmin
591  lquery = ( lwork.EQ.-1 .OR. liwork.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( liwork.LT.liwmin .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( liwork.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( liwork.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, 'PDGESVX', -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 pdgeequ( 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 pdlaqge( 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 pdcopy( n, c, 1, ja, cdesc, cdesc( lld_ ), work, ib, jb,
703  $ descb, 1 )
704  IF( mycol.EQ.ibcol ) THEN
705  CALL dgebs2d( ictxt, 'Rowwise', ' ', np, 1, work( iib ),
706  $ descb( lld_ ) )
707  ELSE
708  CALL dgebr2d( ictxt, 'Rowwise', ' ', np, 1, work( 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_ ) ) = work( 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 pdlacpy( 'Full', n, n, a, ia, ja, desca, af, iaf, jaf,
724  $ descaf )
725  CALL pdgetrf( 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 = pdlange( norm, n, n, a, ia, ja, desca, work )
744 *
745 * Compute the reciprocal of the condition number of A.
746 *
747  CALL pdgecon( norm, n, af, iaf, jaf, descaf, anorm, rcond, work,
748  $ lwork, iwork, liwork, info )
749 *
750 * Return if the matrix is singular to working precision.
751 *
752  IF( rcond.LT.pdlamch( ictxt, 'Epsilon' ) ) THEN
753  info = ia + n
754  RETURN
755  END IF
756 *
757 * Compute the solution matrix X.
758 *
759  CALL pdlacpy( 'Full', n, nrhs, b, ib, jb, descb, x, ix, jx,
760  $ descx )
761  CALL pdgetrs( 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 pdgerfs( 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, iwork, liwork, 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 pdcopy( n, c, 1, ja, cdesc, cdesc( lld_ ), work, ix,
791  $ jx, descx, 1 )
792  IF( mycol.EQ.ibcol ) THEN
793  CALL dgebs2d( ictxt, 'Rowwise', ' ', np, 1,
794  $ work( iix ), descx( lld_ ) )
795  ELSE
796  CALL dgebr2d( ictxt, 'Rowwise', ' ', np, 1,
797  $ work( iix ), descx( lld_ ), myrow, ibcol )
798  END IF
799 *
800  DO 80 j = jjx, jjx+nrhsq-1
801  DO 70 i = iix, iix+np-1
802  x( i+( j-1 )*descx( lld_ ) ) = work( i )*
803  $ x( i+( j-1 )*descx( lld_ ) )
804  70 CONTINUE
805  80 CONTINUE
806  DO 90 j = jjx, jjx+nrhsq-1
807  ferr( j ) = ferr( j ) / colcnd
808  90 CONTINUE
809  END IF
810  ELSE IF( rowequ ) THEN
811  DO 110 j = jjx, jjx+nrhsq-1
812  DO 100 i = iix, iix+np-1
813  x( i+( j-1 )*descx( lld_ ) ) = r( i )*
814  $ x( i+( j-1 )*descx( lld_ ) )
815  100 CONTINUE
816  110 CONTINUE
817  DO 120 j = jjx, jjx+nrhsq-1
818  ferr( j ) = ferr( j ) / rowcnd
819  120 CONTINUE
820  END IF
821 *
822  work( 1 ) = dble( lwmin )
823  iwork( 1 ) = liwmin
824 *
825  RETURN
826 *
827 * End of PDGESVX
828 *
829  END
max
#define max(A, B)
Definition: pcgemr.c:180
pdgesvx
subroutine pdgesvx(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, IWORK, LIWORK, INFO)
Definition: pdgesvx.f:5
infog2l
subroutine infog2l(GRINDX, GCINDX, DESC, NPROW, NPCOL, MYROW, MYCOL, LRINDX, LCINDX, RSRC, CSRC)
Definition: infog2l.f:3
pdgerfs
subroutine pdgerfs(TRANS, N, NRHS, A, IA, JA, DESCA, AF, IAF, JAF, DESCAF, IPIV, B, IB, JB, DESCB, X, IX, JX, DESCX, FERR, BERR, WORK, LWORK, IWORK, LIWORK, INFO)
Definition: pdgerfs.f:5
pdlaqge
subroutine pdlaqge(M, N, A, IA, JA, DESCA, R, C, ROWCND, COLCND, AMAX, EQUED)
Definition: pdlaqge.f:3
pdgecon
subroutine pdgecon(NORM, N, A, IA, JA, DESCA, ANORM, RCOND, WORK, LWORK, IWORK, LIWORK, INFO)
Definition: pdgecon.f:3
pdgetrs
subroutine pdgetrs(TRANS, N, NRHS, A, IA, JA, DESCA, IPIV, B, IB, JB, DESCB, INFO)
Definition: pdgetrs.f:3
pdgeequ
subroutine pdgeequ(M, N, A, IA, JA, DESCA, R, C, ROWCND, COLCND, AMAX, INFO)
Definition: pdgeequ.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
pdgetrf
subroutine pdgetrf(M, N, A, IA, JA, DESCA, IPIV, INFO)
Definition: pdgetrf.f:2
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
pdlacpy
subroutine pdlacpy(UPLO, M, N, A, IA, JA, DESCA, B, IB, JB, DESCB)
Definition: pdlacpy.f:3
min
#define min(A, B)
Definition: pcgemr.c:181