ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
psposvx.f
Go to the documentation of this file.
1  SUBROUTINE psposvx( 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, 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, UPLO
13  INTEGER IA, IAF, IB, INFO, IX, JA, JAF, JB, JX, LIWORK,
14  $ LWORK, N, NRHS
15  REAL RCOND
16 * ..
17 * .. Array Arguments ..
18  INTEGER DESCA( * ), DESCAF( * ), DESCB( * ),
19  $ DESCX( * ), IWORK( * )
20  REAL A( * ), AF( * ),
21  $ b( * ), berr( * ), ferr( * ),
22  $ sc( * ), sr( * ), work( * ), x( * )
23 * ..
24 *
25 * Purpose
26 * =======
27 *
28 * PSPOSVX uses the Cholesky factorization A = U**T*U or A = L*L**T to
29 * compute the solution to a real 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) REAL pointer into
157 * the local memory to an array of local dimension
158 * ( LLD_A, LOCc(JA+N-1) ).
159 * On entry, the symmetric 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) REAL 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) REAL 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) REAL 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) REAL 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) REAL 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) REAL 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( PSPOCON( LWORK ), PSPORFS( 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 * IWORK (local workspace/local output) INTEGER array,
321 * dimension (LIWORK)
322 * On exit, IWORK(1) returns the minimal and optimal LIWORK.
323 *
324 * LIWORK (local or global input) INTEGER
325 * The dimension of the array IWORK.
326 * LIWORK is local input and must be at least
327 * LIWORK = DESCA( LLD_ )
328 * LIWORK = LOCr(N_A).
329 *
330 * If LIWORK = -1, then LIWORK is global input and a workspace
331 * query is assumed; the routine only calculates the minimum
332 * and optimal size for all work arrays. Each of these
333 * values is returned in the first entry of the corresponding
334 * work array, and no error message is issued by PXERBLA.
335 *
336 *
337 * INFO (global output) INTEGER
338 * = 0: successful exit
339 * < 0: if INFO = -i, the i-th argument had an illegal value
340 * > 0: if INFO = i, and i is
341 * <= N: if INFO = i, the leading minor of order i of A
342 * is not positive definite, so the factorization
343 * could not be completed, and the solution and error
344 * bounds could not be computed.
345 * = N+1: RCOND is less than machine precision. The
346 * factorization has been completed, but the matrix
347 * is singular to working precision, and the solution
348 * and error bounds have not been computed.
349 *
350 * =====================================================================
351 *
352 * .. Parameters ..
353  INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
354  $ LLD_, MB_, M_, NB_, N_, RSRC_
355  PARAMETER ( BLOCK_CYCLIC_2D = 1, dlen_ = 9, dtype_ = 1,
356  $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
357  $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
358  REAL ONE, ZERO
359  PARAMETER ( ONE = 1.0e+0, zero = 0.0e+0 )
360 * ..
361 * .. Local Scalars ..
362  LOGICAL EQUIL, LQUERY, NOFACT, RCEQU
363  INTEGER I, IACOL, IAROW, IAFROW, IBROW, IBCOL, ICOFF,
364  $ ICOFFA, ICTXT, IDUMM, IIA, IIB, IIX, INFEQU,
365  $ iroff, iroffa, iroffaf, iroffb, iroffx, ixcol,
366  $ ixrow, j, jja, jjb, jjx, ldb, ldx, liwmin,
367  $ lwmin, mycol, myrow, np, npcol, nprow, nrhsq,
368  $ nq
369  REAL AMAX, ANORM, BIGNUM, SCOND, SMAX, SMIN, SMLNUM
370 * ..
371 * .. Local Arrays ..
372  INTEGER IDUM1( 5 ), IDUM2( 5 )
373 * ..
374 * .. External Subroutines ..
375  EXTERNAL blacs_gridinfo, chk1mat, pchk2mat, infog2l,
376  $ pspocon, pspoequ, psporfs,
377  $ pspotrf, pspotrs,
378  $ pslacpy, pslaqsy, pxerbla,
379  $ sgamn2d, sgamx2d
380 * ..
381 * .. External Functions ..
382  LOGICAL LSAME
383  INTEGER INDXG2P, NUMROC
384  REAL PSLANSY, PSLAMCH
385  EXTERNAL indxg2p, lsame, numroc, pslansy, pslamch
386 * ..
387 * .. Intrinsic Functions ..
388  INTRINSIC ichar, max, min, mod
389 * ..
390 * .. Executable Statements ..
391 *
392 * Get grid parameters
393 *
394  ictxt = desca( ctxt_ )
395  CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
396 *
397 * Test the input parameters
398 *
399  info = 0
400  IF( nprow.EQ.-1 ) THEN
401  info = -(800+ctxt_)
402  ELSE
403  CALL chk1mat( n, 3, n, 3, ia, ja, desca, 8, info )
404  IF( lsame( fact, 'F' ) )
405  $ CALL chk1mat( n, 3, n, 3, iaf, jaf, descaf, 12, info )
406  CALL chk1mat( n, 3, nrhs, 4, ib, jb, descb, 20, info )
407  IF( info.EQ.0 ) THEN
408  iarow = indxg2p( ia, desca( mb_ ), myrow, desca( rsrc_ ),
409  $ nprow )
410  iafrow = indxg2p( iaf, descaf( mb_ ), myrow,
411  $ descaf( rsrc_ ), nprow )
412  ibrow = indxg2p( ib, descb( mb_ ), myrow, descb( rsrc_ ),
413  $ nprow )
414  ixrow = indxg2p( ix, descx( mb_ ), myrow, descx( rsrc_ ),
415  $ nprow )
416  iroffa = mod( ia-1, desca( mb_ ) )
417  iroffaf = mod( iaf-1, descaf( mb_ ) )
418  icoffa = mod( ja-1, desca( nb_ ) )
419  iroffb = mod( ib-1, descb( mb_ ) )
420  iroffx = mod( ix-1, descx( mb_ ) )
421  CALL infog2l( ia, ja, desca, nprow, npcol, myrow,
422  $ mycol, iia, jja, iarow, iacol )
423  np = numroc( n+iroffa, desca( mb_ ), myrow, iarow, nprow )
424  IF( myrow.EQ.iarow )
425  $ np = np-iroffa
426  nq = numroc( n+icoffa, desca( nb_ ), mycol, iacol, npcol )
427  IF( mycol.EQ.iacol )
428  $ nq = nq-icoffa
429  lwmin = 3*desca( lld_ )
430  liwmin = np
431  nofact = lsame( fact, 'N' )
432  equil = lsame( fact, 'E' )
433  IF( nofact .OR. equil ) THEN
434  equed = 'N'
435  rcequ = .false.
436  ELSE
437  rcequ = lsame( equed, 'Y' )
438  smlnum = pslamch( ictxt, 'Safe minimum' )
439  bignum = one / smlnum
440  END IF
441  IF( .NOT.nofact .AND. .NOT.equil .AND.
442  $ .NOT.lsame( fact, 'F' ) ) THEN
443  info = -1
444  ELSE IF( .NOT.lsame( uplo, 'U' ) .AND.
445  $ .NOT.lsame( uplo, 'L' ) ) THEN
446  info = -2
447  ELSE IF( iroffa.NE.0 ) THEN
448  info = -6
449  ELSE IF( icoffa.NE.0 .OR. iroffa.NE.icoffa ) THEN
450  info = -7
451  ELSE IF( desca( mb_ ).NE.desca( nb_ ) ) THEN
452  info = -(800+nb_)
453  ELSE IF( iafrow.NE.iarow ) THEN
454  info = -10
455  ELSE IF( iroffaf.NE.0 ) THEN
456  info = -10
457  ELSE IF( ictxt.NE.descaf( ctxt_ ) ) THEN
458  info = -(1200+ctxt_)
459  ELSE IF( lsame( fact, 'F' ) .AND. .NOT.
460  $ ( rcequ .OR. lsame( equed, 'N' ) ) ) THEN
461  info = -13
462  ELSE
463  IF( rcequ ) THEN
464 *
465  smin = bignum
466  smax = zero
467  DO 10 j = iia, iia + np - 1
468  smin = min( smin, sr( j ) )
469  smax = max( smax, sr( j ) )
470  10 CONTINUE
471  CALL sgamn2d( ictxt, 'Columnwise', ' ', 1, 1, smin,
472  $ 1, idumm, idumm, -1, -1, mycol )
473  CALL sgamx2d( ictxt, 'Columnwise', ' ', 1, 1, smax,
474  $ 1, idumm, idumm, -1, -1, mycol )
475  IF( smin.LE.zero ) THEN
476  info = -14
477  ELSE IF( n.GT.0 ) THEN
478  scond = max( smin, smlnum ) / min( smax, bignum )
479  ELSE
480  scond = one
481  END IF
482  END IF
483  END IF
484  END IF
485 *
486  work( 1 ) = real( lwmin )
487  iwork( 1 ) = liwmin
488  lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
489  IF( info.EQ.0 ) THEN
490  IF( ibrow.NE.iarow ) THEN
491  info = -18
492  ELSE IF( ixrow.NE.ibrow ) THEN
493  info = -22
494  ELSE IF( descb( mb_ ).NE.desca( nb_ ) ) THEN
495  info = -(2000+nb_)
496  ELSE IF( ictxt.NE.descb( ctxt_ ) ) THEN
497  info = -(2000+ctxt_)
498  ELSE IF( ictxt.NE.descx( ctxt_ ) ) THEN
499  info = -(2400+ctxt_)
500  ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
501  info = -28
502  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
503  info = -30
504  END IF
505  idum1( 1 ) = ichar( fact )
506  idum2( 1 ) = 1
507  idum1( 2 ) = ichar( uplo )
508  idum2( 2 ) = 2
509  IF( lsame( fact, 'F' ) ) THEN
510  idum1( 3 ) = ichar( equed )
511  idum2( 3 ) = 13
512  IF( lwork.EQ.-1 ) THEN
513  idum1( 4 ) = -1
514  ELSE
515  idum1( 4 ) = 1
516  END IF
517  idum2( 4 ) = 28
518  IF( liwork.EQ.-1 ) THEN
519  idum1( 5 ) = -1
520  ELSE
521  idum1( 5 ) = 1
522  END IF
523  idum2( 5 ) = 30
524  CALL pchk2mat( n, 3, n, 3, ia, ja, desca, 8, n, 3, nrhs,
525  $ 4, ib, jb, descb, 19, 5, idum1, idum2,
526  $ info )
527  ELSE
528  IF( lwork.EQ.-1 ) THEN
529  idum1( 3 ) = -1
530  ELSE
531  idum1( 3 ) = 1
532  END IF
533  idum2( 3 ) = 28
534  IF( liwork.EQ.-1 ) THEN
535  idum1( 4 ) = -1
536  ELSE
537  idum1( 4 ) = 1
538  END IF
539  idum2( 4 ) = 30
540  CALL pchk2mat( n, 3, n, 3, ia, ja, desca, 8, n, 3, nrhs,
541  $ 4, ib, jb, descb, 19, 4, idum1, idum2,
542  $ info )
543  END IF
544  END IF
545  END IF
546 *
547  IF( info.NE.0 ) THEN
548  CALL pxerbla( ictxt, 'PSPOSVX', -info )
549  RETURN
550  ELSE IF( lquery ) THEN
551  RETURN
552  END IF
553 *
554  IF( equil ) THEN
555 *
556 * Compute row and column scalings to equilibrate the matrix A.
557 *
558  CALL pspoequ( n, a, ia, ja, desca, sr, sc, scond, amax,
559  $ infequ )
560 *
561  IF( infequ.EQ.0 ) THEN
562 *
563 * Equilibrate the matrix.
564 *
565  CALL pslaqsy( uplo, n, a, ia, ja, desca, sr, sc, scond,
566  $ amax, equed )
567  rcequ = lsame( equed, 'Y' )
568  END IF
569  END IF
570 *
571 * Scale the right-hand side.
572 *
573  CALL infog2l( ib, jb, descb, nprow, npcol, myrow, mycol, iib,
574  $ jjb, ibrow, ibcol )
575  ldb = descb( lld_ )
576  iroff = mod( ib-1, descb( mb_ ) )
577  icoff = mod( jb-1, descb( nb_ ) )
578  np = numroc( n+iroff, descb( mb_ ), myrow, ibrow, nprow )
579  nrhsq = numroc( nrhs+icoff, descb( nb_ ), mycol, ibcol, npcol )
580  IF( myrow.EQ.ibrow ) np = np-iroff
581  IF( mycol.EQ.ibcol ) nrhsq = nrhsq-icoff
582 *
583  IF( rcequ ) THEN
584  DO 30 j = jjb, jjb+nrhsq-1
585  DO 20 i = iib, iib+np-1
586  b( i + ( j-1 )*ldb ) = sr( i )*b( i + ( j-1 )*ldb )
587  20 CONTINUE
588  30 CONTINUE
589  END IF
590 *
591  IF( nofact .OR. equil ) THEN
592 *
593 * Compute the Cholesky factorization A = U'*U or A = L*L'.
594 *
595  CALL pslacpy( 'Full', n, n, a, ia, ja, desca, af, iaf, jaf,
596  $ descaf )
597  CALL pspotrf( uplo, n, af, iaf, jaf, descaf, info )
598 *
599 * Return if INFO is non-zero.
600 *
601  IF( info.NE.0 ) THEN
602  IF( info.GT.0 )
603  $ rcond = zero
604  RETURN
605  END IF
606  END IF
607 *
608 * Compute the norm of the matrix A.
609 *
610  anorm = pslansy( '1', uplo, n, a, ia, ja, desca, work )
611 *
612 * Compute the reciprocal of the condition number of A.
613 *
614  CALL pspocon( uplo, n, af, iaf, jaf, descaf, anorm, rcond, work,
615  $ lwork, iwork, liwork, info )
616 *
617 * Return if the matrix is singular to working precision.
618 *
619  IF( rcond.LT.pslamch( ictxt, 'Epsilon' ) ) THEN
620  info = ia + n
621  RETURN
622  END IF
623 *
624 * Compute the solution matrix X.
625 *
626  CALL pslacpy( 'Full', n, nrhs, b, ib, jb, descb, x, ix, jx,
627  $ descx )
628  CALL pspotrs( uplo, n, nrhs, af, iaf, jaf, descaf, x, ix, jx,
629  $ descx, info )
630 *
631 * Use iterative refinement to improve the computed solution and
632 * compute error bounds and backward error estimates for it.
633 *
634  CALL psporfs( uplo, n, nrhs, a, ia, ja, desca, af, iaf, jaf,
635  $ descaf, b, ib, jb, descb, x, ix, jx, descx, ferr,
636  $ berr, work, lwork, iwork, liwork, info )
637 *
638 * Transform the solution matrix X to a solution of the original
639 * system.
640 *
641  CALL infog2l( ix, jx, descx, nprow, npcol, myrow, mycol, iix,
642  $ jjx, ixrow, ixcol )
643  ldx = descx( lld_ )
644  iroff = mod( ix-1, descx( mb_ ) )
645  icoff = mod( jx-1, descx( nb_ ) )
646  np = numroc( n+iroff, descx( mb_ ), myrow, ixrow, nprow )
647  nrhsq = numroc( nrhs+icoff, descx( nb_ ), mycol, ixcol, npcol )
648  IF( myrow.EQ.ibrow ) np = np-iroff
649  IF( mycol.EQ.ibcol ) nrhsq = nrhsq-icoff
650 *
651  IF( rcequ ) THEN
652  DO 50 j = jjx, jjx+nrhsq-1
653  DO 40 i = iix, iix+np-1
654  x( i + ( j-1 )*ldx ) = sr( i )*x( i + ( j-1 )*ldx )
655  40 CONTINUE
656  50 CONTINUE
657  DO 60 j = jjx, jjx+nrhsq-1
658  ferr( j ) = ferr( j ) / scond
659  60 CONTINUE
660  END IF
661 *
662  work( 1 ) = real( lwmin )
663  iwork( 1 ) = liwmin
664  RETURN
665 *
666 * End of PSPOSVX
667 *
668  END
pspoequ
subroutine pspoequ(N, A, IA, JA, DESCA, SR, SC, SCOND, AMAX, INFO)
Definition: pspoequ.f:3
max
#define max(A, B)
Definition: pcgemr.c:180
infog2l
subroutine infog2l(GRINDX, GCINDX, DESC, NPROW, NPCOL, MYROW, MYCOL, LRINDX, LCINDX, RSRC, CSRC)
Definition: infog2l.f:3
psposvx
subroutine psposvx(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, IWORK, LIWORK, INFO)
Definition: psposvx.f:5
pspocon
subroutine pspocon(UPLO, N, A, IA, JA, DESCA, ANORM, RCOND, WORK, LWORK, IWORK, LIWORK, INFO)
Definition: pspocon.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
pslacpy
subroutine pslacpy(UPLO, M, N, A, IA, JA, DESCA, B, IB, JB, DESCB)
Definition: pslacpy.f:3
pspotrs
subroutine pspotrs(UPLO, N, NRHS, A, IA, JA, DESCA, B, IB, JB, DESCB, INFO)
Definition: pspotrs.f:3
chk1mat
subroutine chk1mat(MA, MAPOS0, NA, NAPOS0, IA, JA, DESCA, DESCAPOS0, INFO)
Definition: chk1mat.f:3
pslaqsy
subroutine pslaqsy(UPLO, N, A, IA, JA, DESCA, SR, SC, SCOND, AMAX, EQUED)
Definition: pslaqsy.f:3
pspotrf
subroutine pspotrf(UPLO, N, A, IA, JA, DESCA, INFO)
Definition: pspotrf.f:2
pxerbla
subroutine pxerbla(ICTXT, SRNAME, INFO)
Definition: pxerbla.f:2
min
#define min(A, B)
Definition: pcgemr.c:181
psporfs
subroutine psporfs(UPLO, N, NRHS, A, IA, JA, DESCA, AF, IAF, JAF, DESCAF, B, IB, JB, DESCB, X, IX, JX, DESCX, FERR, BERR, WORK, LWORK, IWORK, LIWORK, INFO)
Definition: psporfs.f:4