ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pcposvx.f
Go to the documentation of this file.
1  SUBROUTINE pcposvx( FACT, UPLO, N, NRHS, A, IA, JA, DESCA, AF,
2  $ IAF, JAF, DESCAF, EQUED, SR, SC, B, IB, JB,
3  $ DESCB, X, IX, JX, DESCX, RCOND, FERR, BERR,
4  $ 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, UPLO
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( * ), DESCX( * )
19  REAL BERR( * ), FERR( * ), SC( * ),
20  $ SR( * ), RWORK( * )
21  COMPLEX A( * ), AF( * ),
22  $ b( * ), work( * ), x( * )
23 * ..
24 *
25 * Purpose
26 * =======
27 *
28 * PCPOSVX uses the Cholesky factorization A = U**H*U or A = L*L**H to
29 * compute the solution to a 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. In the following comments Y denotes Y(IY:IY+M-1,JY:JY+K-1)
38 * a M-by-K matrix where Y can be A, AF, B and X.
39 *
40 * Notes
41 * =====
42 *
43 * Each global data object is described by an associated description
44 * vector. This vector stores the information required to establish
45 * the mapping between an object element and its corresponding process
46 * and memory location.
47 *
48 * Let A be a generic term for any 2D block cyclicly distributed array.
49 * Such a global array has an associated description vector DESCA.
50 * In the following comments, the character _ should be read as
51 * "of the global array".
52 *
53 * NOTATION STORED IN EXPLANATION
54 * --------------- -------------- --------------------------------------
55 * DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
56 * DTYPE_A = 1.
57 * CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
58 * the BLACS process grid A is distribu-
59 * ted over. The context itself is glo-
60 * bal, but the handle (the integer
61 * value) may vary.
62 * M_A (global) DESCA( M_ ) The number of rows in the global
63 * array A.
64 * N_A (global) DESCA( N_ ) The number of columns in the global
65 * array A.
66 * MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
67 * the rows of the array.
68 * NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
69 * the columns of the array.
70 * RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
71 * row of the array A is distributed.
72 * CSRC_A (global) DESCA( CSRC_ ) The process column over which the
73 * first column of the array A is
74 * distributed.
75 * LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
76 * array. LLD_A >= MAX(1,LOCr(M_A)).
77 *
78 * Let K be the number of rows or columns of a distributed matrix,
79 * and assume that its process grid has dimension p x q.
80 * LOCr( K ) denotes the number of elements of K that a process
81 * would receive if K were distributed over the p processes of its
82 * process column.
83 * Similarly, LOCc( K ) denotes the number of elements of K that a
84 * process would receive if K were distributed over the q processes of
85 * its process row.
86 * The values of LOCr() and LOCc() may be determined via a call to the
87 * ScaLAPACK tool function, NUMROC:
88 * LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
89 * LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
90 * An upper bound for these quantities may be computed by:
91 * LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
92 * LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
93 *
94 * Description
95 * ===========
96 *
97 * The following steps are performed:
98 *
99 * 1. If FACT = 'E', real scaling factors are computed to equilibrate
100 * the system:
101 * diag(SR) * A * diag(SC) * inv(diag(SC)) * X = diag(SR) * B
102 * Whether or not the system will be equilibrated depends on the
103 * scaling of the matrix A, but if equilibration is used, A is
104 * overwritten by diag(SR)*A*diag(SC) and B by diag(SR)*B.
105 *
106 * 2. If FACT = 'N' or 'E', the Cholesky decomposition is used to
107 * factor the matrix A (after equilibration if FACT = 'E') as
108 * A = U**T* U, if UPLO = 'U', or
109 * A = L * L**T, if UPLO = 'L',
110 * where U is an upper triangular matrix and L is a lower triangular
111 * matrix.
112 *
113 * 3. The factored form of A is used to estimate the condition number
114 * of the matrix A. If the reciprocal of the condition number is
115 * less than machine precision, steps 4-6 are skipped.
116 *
117 * 4. The system of equations is solved for X using the factored form
118 * of A.
119 *
120 * 5. Iterative refinement is applied to improve the computed solution
121 * matrix and calculate error bounds and backward error estimates
122 * for it.
123 *
124 * 6. If equilibration was used, the matrix X is premultiplied by
125 * diag(SR) so that it solves the original system before
126 * equilibration.
127 *
128 * Arguments
129 * =========
130 *
131 * FACT (global input) CHARACTER
132 * Specifies whether or not the factored form of the matrix A is
133 * supplied on entry, and if not, whether the matrix A should be
134 * equilibrated before it is factored.
135 * = 'F': On entry, AF contains the factored form of A.
136 * If EQUED = 'Y', the matrix A has been equilibrated
137 * with scaling factors given by S. A and AF will not
138 * be modified.
139 * = 'N': The matrix A will be copied to AF and factored.
140 * = 'E': The matrix A will be equilibrated if necessary, then
141 * copied to AF and factored.
142 *
143 * UPLO (global input) CHARACTER
144 * = 'U': Upper triangle of A is stored;
145 * = 'L': Lower triangle of A is stored.
146 *
147 * N (global input) INTEGER
148 * The number of rows and columns to be operated on, i.e. the
149 * order of the distributed submatrix A(IA:IA+N-1,JA:JA+N-1).
150 * N >= 0.
151 *
152 * NRHS (global input) INTEGER
153 * The number of right hand sides, i.e., the number of columns
154 * of the distributed submatrices B and X. NRHS >= 0.
155 *
156 * A (local input/local output) COMPLEX pointer into
157 * the local memory to an array of local dimension
158 * ( LLD_A, LOCc(JA+N-1) ).
159 * On entry, the Hermitian matrix A, except if FACT = 'F' and
160 * EQUED = 'Y', then A must contain the equilibrated matrix
161 * diag(SR)*A*diag(SC). If UPLO = 'U', the leading
162 * N-by-N upper triangular part of A contains the upper
163 * triangular part of the matrix A, and the strictly lower
164 * triangular part of A is not referenced. If UPLO = 'L', the
165 * leading N-by-N lower triangular part of A contains the lower
166 * triangular part of the matrix A, and the strictly upper
167 * triangular part of A is not referenced. A is not modified if
168 * FACT = 'F' or 'N', or if FACT = 'E' and EQUED = 'N' on exit.
169 *
170 * On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by
171 * diag(SR)*A*diag(SC).
172 *
173 * IA (global input) INTEGER
174 * The row index in the global array A indicating the first
175 * row of sub( A ).
176 *
177 * JA (global input) INTEGER
178 * The column index in the global array A indicating the
179 * first column of sub( A ).
180 *
181 * DESCA (global and local input) INTEGER array of dimension DLEN_.
182 * The array descriptor for the distributed matrix A.
183 *
184 * AF (local input or local output) COMPLEX pointer
185 * into the local memory to an array of local dimension
186 * ( LLD_AF, LOCc(JA+N-1)).
187 * If FACT = 'F', then AF is an input argument and on entry
188 * contains the triangular factor U or L from the Cholesky
189 * factorization A = U**T*U or A = L*L**T, in the same storage
190 * format as A. If EQUED .ne. 'N', then AF is the factored form
191 * of the equilibrated matrix diag(SR)*A*diag(SC).
192 *
193 * If FACT = 'N', then AF is an output argument and on exit
194 * returns the triangular factor U or L from the Cholesky
195 * factorization A = U**T*U or A = L*L**T of the original
196 * matrix A.
197 *
198 * If FACT = 'E', then AF is an output argument and on exit
199 * returns the triangular factor U or L from the Cholesky
200 * factorization A = U**T*U or A = L*L**T of the equilibrated
201 * matrix A (see the description of A for the form of the
202 * equilibrated matrix).
203 *
204 * IAF (global input) INTEGER
205 * The row index in the global array AF indicating the first
206 * row of sub( AF ).
207 *
208 * JAF (global input) INTEGER
209 * The column index in the global array AF indicating the
210 * first column of sub( AF ).
211 *
212 * DESCAF (global and local input) INTEGER array of dimension DLEN_.
213 * The array descriptor for the distributed matrix AF.
214 *
215 * EQUED (global input/global output) CHARACTER
216 * Specifies the form of equilibration that was done.
217 * = 'N': No equilibration (always true if FACT = 'N').
218 * = 'Y': Equilibration was done, i.e., A has been replaced by
219 * diag(SR) * A * diag(SC).
220 * EQUED is an input variable if FACT = 'F'; otherwise, it is an
221 * output variable.
222 *
223 * SR (local input/local output) COMPLEX array,
224 * dimension (LLD_A)
225 * The scale factors for A distributed across process rows;
226 * not accessed if EQUED = 'N'. SR is an input variable if
227 * FACT = 'F'; otherwise, SR is an output variable.
228 * If FACT = 'F' and EQUED = 'Y', each element of SR must be
229 * positive.
230 *
231 * SC (local input/local output) COMPLEX array,
232 * dimension (LOC(N_A))
233 * The scale factors for A distributed across
234 * process columns; not accessed if EQUED = 'N'. SC is an input
235 * variable if FACT = 'F'; otherwise, SC is an output variable.
236 * If FACT = 'F' and EQUED = 'Y', each element of SC must be
237 * positive.
238 *
239 * B (local input/local output) COMPLEX pointer into
240 * the local memory to an array of local dimension
241 * ( LLD_B, LOCc(JB+NRHS-1) ).
242 * On entry, the N-by-NRHS right-hand side matrix B.
243 * On exit, if EQUED = 'N', B is not modified; if TRANS = 'N'
244 * and EQUED = 'R' or 'B', B is overwritten by diag(R)*B; if
245 * TRANS = 'T' or 'C' and EQUED = 'C' or 'B', B is overwritten
246 * by diag(C)*B.
247 *
248 * IB (global input) INTEGER
249 * The row index in the global array B indicating the first
250 * row of sub( B ).
251 *
252 * JB (global input) INTEGER
253 * The column index in the global array B indicating the
254 * first column of sub( B ).
255 *
256 * DESCB (global and local input) INTEGER array of dimension DLEN_.
257 * The array descriptor for the distributed matrix B.
258 *
259 * X (local input/local output) COMPLEX pointer into
260 * the local memory to an array of local dimension
261 * ( LLD_X, LOCc(JX+NRHS-1) ).
262 * If INFO = 0, the N-by-NRHS solution matrix X to the original
263 * system of equations. Note that A and B are modified on exit
264 * if EQUED .ne. 'N', and the solution to the equilibrated
265 * system is inv(diag(SC))*X if TRANS = 'N' and EQUED = 'C' or
266 * 'B', or inv(diag(SR))*X if TRANS = 'T' or 'C' and EQUED = 'R'
267 * or 'B'.
268 *
269 * IX (global input) INTEGER
270 * The row index in the global array X indicating the first
271 * row of sub( X ).
272 *
273 * JX (global input) INTEGER
274 * The column index in the global array X indicating the
275 * first column of sub( X ).
276 *
277 * DESCX (global and local input) INTEGER array of dimension DLEN_.
278 * The array descriptor for the distributed matrix X.
279 *
280 * RCOND (global output) REAL
281 * The estimate of the reciprocal condition number of the matrix
282 * A after equilibration (if done). If RCOND is less than the
283 * machine precision (in particular, if RCOND = 0), the matrix
284 * is singular to working precision. This condition is
285 * indicated by a return code of INFO > 0, and the solution and
286 * error bounds are not computed.
287 *
288 * FERR (local output) REAL array, dimension (LOC(N_B))
289 * The estimated forward error bounds for each solution vector
290 * X(j) (the j-th column of the solution matrix X).
291 * If XTRUE is the true solution, FERR(j) bounds the magnitude
292 * of the largest entry in (X(j) - XTRUE) divided by
293 * the magnitude of the largest entry in X(j). The quality of
294 * the error bound depends on the quality of the estimate of
295 * norm(inv(A)) computed in the code; if the estimate of
296 * norm(inv(A)) is accurate, the error bound is guaranteed.
297 *
298 * BERR (local output) REAL array, dimension (LOC(N_B))
299 * The componentwise relative backward error of each solution
300 * vector X(j) (i.e., the smallest relative change in
301 * any entry of A or B that makes X(j) an exact solution).
302 *
303 * WORK (local workspace/local output) COMPLEX array,
304 * dimension (LWORK)
305 * On exit, WORK(1) returns the minimal and optimal LWORK.
306 *
307 * LWORK (local or global input) INTEGER
308 * The dimension of the array WORK.
309 * LWORK is local input and must be at least
310 * LWORK = MAX( PCPOCON( LWORK ), PCPORFS( LWORK ) )
311 * + LOCr( N_A ).
312 * LWORK = 3*DESCA( LLD_ )
313 *
314 * If LWORK = -1, then LWORK is global input and a workspace
315 * query is assumed; the routine only calculates the minimum
316 * and optimal size for all work arrays. Each of these
317 * values is returned in the first entry of the corresponding
318 * work array, and no error message is issued by PXERBLA.
319 *
320 * RWORK (local workspace/local output) REAL array,
321 * dimension (LRWORK)
322 * On exit, RWORK(1) returns the minimal and optimal LRWORK.
323 *
324 * LRWORK (local or global input) INTEGER
325 * The dimension of the array RWORK.
326 * LRWORK is local input and must be at least
327 * LRWORK = 2*LOCc(N_A).
328 *
329 * If LRWORK = -1, then LRWORK is global input and a workspace
330 * query is assumed; the routine only calculates the minimum
331 * and optimal size for all work arrays. Each of these
332 * values is returned in the first entry of the corresponding
333 * work array, and no error message is issued by PXERBLA.
334 *
335 *
336 * INFO (global output) INTEGER
337 * = 0: successful exit
338 * < 0: if INFO = -i, the i-th argument had an illegal value
339 * > 0: if INFO = i, and i is
340 * <= N: if INFO = i, the leading minor of order i of A
341 * is not positive definite, so the factorization
342 * could not be completed, and the solution and error
343 * bounds could not be computed.
344 * = N+1: RCOND is less than machine precision. The
345 * factorization has been completed, but the matrix
346 * is singular to working precision, and the solution
347 * and error bounds have not been computed.
348 *
349 * =====================================================================
350 *
351 * .. Parameters ..
352  INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
353  $ LLD_, MB_, M_, NB_, N_, RSRC_
354  PARAMETER ( BLOCK_CYCLIC_2D = 1, dlen_ = 9, dtype_ = 1,
355  $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
356  $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
357  REAL ONE, ZERO
358  PARAMETER ( ONE = 1.0e+0, zero = 0.0e+0 )
359 * ..
360 * .. Local Scalars ..
361  LOGICAL EQUIL, LQUERY, NOFACT, RCEQU
362  INTEGER I, IACOL, IAROW, IAFROW, IBROW, IBCOL, ICOFF,
363  $ ICOFFA, ICTXT, IDUMM, IIA, IIB, IIX, INFEQU,
364  $ iroff, iroffa, iroffaf, iroffb, iroffx, ixcol,
365  $ ixrow, j, jja, jjb, jjx, ldb, ldx, lrwmin,
366  $ lwmin, mycol, myrow, np, npcol, nprow, nrhsq,
367  $ nq
368  REAL AMAX, ANORM, BIGNUM, SCOND, SMAX, SMIN, SMLNUM
369 * ..
370 * .. Local Arrays ..
371  INTEGER IDUM1( 5 ), IDUM2( 5 )
372 * ..
373 * .. External Subroutines ..
374  EXTERNAL blacs_gridinfo, chk1mat, pchk2mat, infog2l,
375  $ pcpocon, pcpoequ,
376  $ pcporfs, pcpotrf, pcpotrs,
377  $ pclacpy, pclaqsy, pxerbla,
378  $ sgamn2d, sgamx2d
379 * ..
380 * .. External Functions ..
381  LOGICAL LSAME
382  INTEGER INDXG2P, NUMROC
383  REAL PCLANHE, PSLAMCH
384  EXTERNAL indxg2p, lsame, numroc, pclanhe, pslamch
385 * ..
386 * .. Intrinsic Functions ..
387  INTRINSIC ichar, max, min, mod
388 * ..
389 * .. Executable Statements ..
390 *
391 * Get grid parameters
392 *
393  ictxt = desca( ctxt_ )
394  CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
395 *
396 * Test the input parameters
397 *
398  info = 0
399  IF( nprow.EQ.-1 ) THEN
400  info = -(800+ctxt_)
401  ELSE
402  CALL chk1mat( n, 3, n, 3, ia, ja, desca, 8, info )
403  IF( lsame( fact, 'F' ) )
404  $ CALL chk1mat( n, 3, n, 3, iaf, jaf, descaf, 12, info )
405  CALL chk1mat( n, 3, nrhs, 4, ib, jb, descb, 20, info )
406  IF( info.EQ.0 ) THEN
407  iarow = indxg2p( ia, desca( mb_ ), myrow, desca( rsrc_ ),
408  $ nprow )
409  iafrow = indxg2p( iaf, descaf( mb_ ), myrow,
410  $ descaf( rsrc_ ), nprow )
411  ibrow = indxg2p( ib, descb( mb_ ), myrow, descb( rsrc_ ),
412  $ nprow )
413  ixrow = indxg2p( ix, descx( mb_ ), myrow, descx( rsrc_ ),
414  $ nprow )
415  iroffa = mod( ia-1, desca( mb_ ) )
416  iroffaf = mod( iaf-1, descaf( mb_ ) )
417  icoffa = mod( ja-1, desca( nb_ ) )
418  iroffb = mod( ib-1, descb( mb_ ) )
419  iroffx = mod( ix-1, descx( mb_ ) )
420  CALL infog2l( ia, ja, desca, nprow, npcol, myrow,
421  $ mycol, iia, jja, iarow, iacol )
422  np = numroc( n+iroffa, desca( mb_ ), myrow, iarow, nprow )
423  IF( myrow.EQ.iarow )
424  $ np = np-iroffa
425  nq = numroc( n+icoffa, desca( nb_ ), mycol, iacol, npcol )
426  IF( mycol.EQ.iacol )
427  $ nq = nq-icoffa
428  lwmin = 3*desca( lld_ )
429  lrwmin = max( 2*nq, np )
430  nofact = lsame( fact, 'N' )
431  equil = lsame( fact, 'E' )
432  IF( nofact .OR. equil ) THEN
433  equed = 'N'
434  rcequ = .false.
435  ELSE
436  rcequ = lsame( equed, 'Y' )
437  smlnum = pslamch( ictxt, 'Safe minimum' )
438  bignum = one / smlnum
439  END IF
440  IF( .NOT.nofact .AND. .NOT.equil .AND.
441  $ .NOT.lsame( fact, 'F' ) ) THEN
442  info = -1
443  ELSE IF( .NOT.lsame( uplo, 'U' ) .AND.
444  $ .NOT.lsame( uplo, 'L' ) ) THEN
445  info = -2
446  ELSE IF( iroffa.NE.0 ) THEN
447  info = -6
448  ELSE IF( icoffa.NE.0 .OR. iroffa.NE.icoffa ) THEN
449  info = -7
450  ELSE IF( desca( mb_ ).NE.desca( nb_ ) ) THEN
451  info = -(800+nb_)
452  ELSE IF( iafrow.NE.iarow ) THEN
453  info = -10
454  ELSE IF( iroffaf.NE.0 ) THEN
455  info = -10
456  ELSE IF( ictxt.NE.descaf( ctxt_ ) ) THEN
457  info = -(1200+ctxt_)
458  ELSE IF( lsame( fact, 'F' ) .AND. .NOT.
459  $ ( rcequ .OR. lsame( equed, 'N' ) ) ) THEN
460  info = -13
461  ELSE
462  IF( rcequ ) THEN
463 *
464  smin = bignum
465  smax = zero
466  DO 10 j = iia, iia + np - 1
467  smin = min( smin, sr( j ) )
468  smax = max( smax, sr( j ) )
469  10 CONTINUE
470  CALL sgamn2d( ictxt, 'Columnwise', ' ', 1, 1, smin,
471  $ 1, idumm, idumm, -1, -1, mycol )
472  CALL sgamx2d( ictxt, 'Columnwise', ' ', 1, 1, smax,
473  $ 1, idumm, idumm, -1, -1, mycol )
474  IF( smin.LE.zero ) THEN
475  info = -14
476  ELSE IF( n.GT.0 ) THEN
477  scond = max( smin, smlnum ) / min( smax, bignum )
478  ELSE
479  scond = one
480  END IF
481  END IF
482  END IF
483  END IF
484 *
485  work( 1 ) = real( lwmin )
486  rwork( 1 ) = real( lrwmin )
487  lquery = ( lwork.EQ.-1 .OR. lrwork.EQ.-1 )
488  IF( info.EQ.0 ) THEN
489  IF( ibrow.NE.iarow ) THEN
490  info = -18
491  ELSE IF( ixrow.NE.ibrow ) THEN
492  info = -22
493  ELSE IF( descb( mb_ ).NE.desca( nb_ ) ) THEN
494  info = -(2000+nb_)
495  ELSE IF( ictxt.NE.descb( ctxt_ ) ) THEN
496  info = -(2000+ctxt_)
497  ELSE IF( ictxt.NE.descx( ctxt_ ) ) THEN
498  info = -(2400+ctxt_)
499  ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
500  info = -28
501  ELSE IF( lrwork.LT.lrwmin .AND. .NOT.lquery ) THEN
502  info = -30
503  END IF
504  idum1( 1 ) = ichar( fact )
505  idum2( 1 ) = 1
506  idum1( 2 ) = ichar( uplo )
507  idum2( 2 ) = 2
508  IF( lsame( fact, 'F' ) ) THEN
509  idum1( 3 ) = ichar( equed )
510  idum2( 3 ) = 13
511  IF( lwork.EQ.-1 ) THEN
512  idum1( 4 ) = -1
513  ELSE
514  idum1( 4 ) = 1
515  END IF
516  idum2( 4 ) = 28
517  IF( lrwork.EQ.-1 ) THEN
518  idum1( 5 ) = -1
519  ELSE
520  idum1( 5 ) = 1
521  END IF
522  idum2( 5 ) = 30
523  CALL pchk2mat( n, 3, n, 3, ia, ja, desca, 8, n, 3, nrhs,
524  $ 4, ib, jb, descb, 19, 5, idum1, idum2,
525  $ info )
526  ELSE
527  IF( lwork.EQ.-1 ) THEN
528  idum1( 3 ) = -1
529  ELSE
530  idum1( 3 ) = 1
531  END IF
532  idum2( 3 ) = 28
533  IF( lrwork.EQ.-1 ) THEN
534  idum1( 4 ) = -1
535  ELSE
536  idum1( 4 ) = 1
537  END IF
538  idum2( 4 ) = 30
539  CALL pchk2mat( n, 3, n, 3, ia, ja, desca, 8, n, 3, nrhs,
540  $ 4, ib, jb, descb, 19, 4, idum1, idum2,
541  $ info )
542  END IF
543  END IF
544  END IF
545 *
546  IF( info.NE.0 ) THEN
547  CALL pxerbla( ictxt, 'PCPOSVX', -info )
548  RETURN
549  ELSE IF( lquery ) THEN
550  RETURN
551  END IF
552 *
553  IF( equil ) THEN
554 *
555 * Compute row and column scalings to equilibrate the matrix A.
556 *
557  CALL pcpoequ( n, a, ia, ja, desca, sr, sc, scond, amax,
558  $ infequ )
559 *
560  IF( infequ.EQ.0 ) THEN
561 *
562 * Equilibrate the matrix.
563 *
564  CALL pclaqsy( uplo, n, a, ia, ja, desca, sr, sc, scond,
565  $ amax, equed )
566  rcequ = lsame( equed, 'Y' )
567  END IF
568  END IF
569 *
570 * Scale the right-hand side.
571 *
572  CALL infog2l( ib, jb, descb, nprow, npcol, myrow, mycol, iib,
573  $ jjb, ibrow, ibcol )
574  ldb = descb( lld_ )
575  iroff = mod( ib-1, descb( mb_ ) )
576  icoff = mod( jb-1, descb( nb_ ) )
577  np = numroc( n+iroff, descb( mb_ ), myrow, ibrow, nprow )
578  nrhsq = numroc( nrhs+icoff, descb( nb_ ), mycol, ibcol, npcol )
579  IF( myrow.EQ.ibrow ) np = np-iroff
580  IF( mycol.EQ.ibcol ) nrhsq = nrhsq-icoff
581 *
582  IF( rcequ ) THEN
583  DO 30 j = jjb, jjb+nrhsq-1
584  DO 20 i = iib, iib+np-1
585  b( i + ( j-1 )*ldb ) = sr( i )*b( i + ( j-1 )*ldb )
586  20 CONTINUE
587  30 CONTINUE
588  END IF
589 *
590  IF( nofact .OR. equil ) THEN
591 *
592 * Compute the Cholesky factorization A = U'*U or A = L*L'.
593 *
594  CALL pclacpy( 'Full', n, n, a, ia, ja, desca, af, iaf, jaf,
595  $ descaf )
596  CALL pcpotrf( uplo, n, af, iaf, jaf, descaf, info )
597 *
598 * Return if INFO is non-zero.
599 *
600  IF( info.NE.0 ) THEN
601  IF( info.GT.0 )
602  $ rcond = zero
603  RETURN
604  END IF
605  END IF
606 *
607 * Compute the norm of the matrix A.
608 *
609  anorm = pclanhe( '1', uplo, n, a, ia, ja, desca, rwork )
610 *
611 * Compute the reciprocal of the condition number of A.
612 *
613  CALL pcpocon( uplo, n, af, iaf, jaf, descaf, anorm, rcond, work,
614  $ lwork, rwork, lrwork, info )
615 *
616 * Return if the matrix is singular to working precision.
617 *
618  IF( rcond.LT.pslamch( ictxt, 'Epsilon' ) ) THEN
619  info = ia + n
620  RETURN
621  END IF
622 *
623 * Compute the solution matrix X.
624 *
625  CALL pclacpy( 'Full', n, nrhs, b, ib, jb, descb, x, ix, jx,
626  $ descx )
627  CALL pcpotrs( uplo, n, nrhs, af, iaf, jaf, descaf, x, ix, jx,
628  $ descx, info )
629 *
630 * Use iterative refinement to improve the computed solution and
631 * compute error bounds and backward error estimates for it.
632 *
633  CALL pcporfs( uplo, n, nrhs, a, ia, ja, desca, af, iaf, jaf,
634  $ descaf, b, ib, jb, descb, x, ix, jx, descx, ferr,
635  $ berr, work, lwork, rwork, lrwork, info )
636 *
637 * Transform the solution matrix X to a solution of the original
638 * system.
639 *
640  CALL infog2l( ix, jx, descx, nprow, npcol, myrow, mycol, iix,
641  $ jjx, ixrow, ixcol )
642  ldx = descx( lld_ )
643  iroff = mod( ix-1, descx( mb_ ) )
644  icoff = mod( jx-1, descx( nb_ ) )
645  np = numroc( n+iroff, descx( mb_ ), myrow, ixrow, nprow )
646  nrhsq = numroc( nrhs+icoff, descx( nb_ ), mycol, ixcol, npcol )
647  IF( myrow.EQ.ibrow ) np = np-iroff
648  IF( mycol.EQ.ibcol ) nrhsq = nrhsq-icoff
649 *
650  IF( rcequ ) THEN
651  DO 50 j = jjx, jjx+nrhsq-1
652  DO 40 i = iix, iix+np-1
653  x( i + ( j-1 )*ldx ) = sr( i )*x( i + ( j-1 )*ldx )
654  40 CONTINUE
655  50 CONTINUE
656  DO 60 j = jjx, jjx+nrhsq-1
657  ferr( j ) = ferr( j ) / scond
658  60 CONTINUE
659  END IF
660 *
661  work( 1 ) = real( lwmin )
662  rwork( 1 ) = real( lrwmin )
663  RETURN
664 *
665 * End of PCPOSVX
666 *
667  END
max
#define max(A, B)
Definition: pcgemr.c:180
infog2l
subroutine infog2l(GRINDX, GCINDX, DESC, NPROW, NPCOL, MYROW, MYCOL, LRINDX, LCINDX, RSRC, CSRC)
Definition: infog2l.f:3
pcporfs
subroutine pcporfs(UPLO, N, NRHS, A, IA, JA, DESCA, AF, IAF, JAF, DESCAF, B, IB, JB, DESCB, X, IX, JX, DESCX, FERR, BERR, WORK, LWORK, RWORK, LRWORK, INFO)
Definition: pcporfs.f:4
pcpoequ
subroutine pcpoequ(N, A, IA, JA, DESCA, SR, SC, SCOND, AMAX, INFO)
Definition: pcpoequ.f:3
pcpotrf
subroutine pcpotrf(UPLO, N, A, IA, JA, DESCA, INFO)
Definition: pcpotrf.f:2
pcpocon
subroutine pcpocon(UPLO, N, A, IA, JA, DESCA, ANORM, RCOND, WORK, LWORK, RWORK, LRWORK, INFO)
Definition: pcpocon.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
pcposvx
subroutine pcposvx(FACT, UPLO, N, NRHS, A, IA, JA, DESCA, AF, IAF, JAF, DESCAF, EQUED, SR, SC, B, IB, JB, DESCB, X, IX, JX, DESCX, RCOND, FERR, BERR, WORK, LWORK, RWORK, LRWORK, INFO)
Definition: pcposvx.f:5
pclacpy
subroutine pclacpy(UPLO, M, N, A, IA, JA, DESCA, B, IB, JB, DESCB)
Definition: pclacpy.f:3
pclaqsy
subroutine pclaqsy(UPLO, N, A, IA, JA, DESCA, SR, SC, SCOND, AMAX, EQUED)
Definition: pclaqsy.f:3
chk1mat
subroutine chk1mat(MA, MAPOS0, NA, NAPOS0, IA, JA, DESCA, DESCAPOS0, INFO)
Definition: chk1mat.f:3
pxerbla
subroutine pxerbla(ICTXT, SRNAME, INFO)
Definition: pxerbla.f:2
min
#define min(A, B)
Definition: pcgemr.c:181
pcpotrs
subroutine pcpotrs(UPLO, N, NRHS, A, IA, JA, DESCA, B, IB, JB, DESCB, INFO)
Definition: pcpotrs.f:3