SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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,
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
subroutine chk1mat(ma, mapos0, na, napos0, ia, ja, desca, descapos0, info)
Definition chk1mat.f:3
subroutine descset(desc, m, n, mb, nb, irsrc, icsrc, ictxt, lld)
Definition descset.f:3
subroutine infog2l(grindx, gcindx, desc, nprow, npcol, myrow, mycol, lrindx, lcindx, rsrc, csrc)
Definition infog2l.f:3
subroutine pcgecon(norm, n, a, ia, ja, desca, anorm, rcond, work, lwork, rwork, lrwork, info)
Definition pcgecon.f:3
subroutine pcgeequ(m, n, a, ia, ja, desca, r, c, rowcnd, colcnd, amax, info)
Definition pcgeequ.f:3
#define max(A, B)
Definition pcgemr.c:180
#define min(A, B)
Definition pcgemr.c:181
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
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
subroutine pcgetrf(m, n, a, ia, ja, desca, ipiv, info)
Definition pcgetrf.f:2
subroutine pcgetrs(trans, n, nrhs, a, ia, ja, desca, ipiv, b, ib, jb, descb, info)
Definition pcgetrs.f:3
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
subroutine pclacpy(uplo, m, n, a, ia, ja, desca, b, ib, jb, descb)
Definition pclacpy.f:3
subroutine pclaqge(m, n, a, ia, ja, desca, r, c, rowcnd, colcnd, amax, equed)
Definition pclaqge.f:3
subroutine pxerbla(ictxt, srname, info)
Definition pxerbla.f:2