ScaLAPACK  2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
psgesvx.f
Go to the documentation of this file.
00001       SUBROUTINE PSGESVX( FACT, TRANS, N, NRHS, A, IA, JA, DESCA, AF,
00002      $                    IAF, JAF, DESCAF, IPIV, EQUED, R, C, B, IB,
00003      $                    JB, DESCB, X, IX, JX, DESCX, RCOND, FERR,
00004      $                    BERR, WORK, LWORK, IWORK, LIWORK, INFO )
00005 *
00006 *  -- ScaLAPACK routine (version 1.7) --
00007 *     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
00008 *     and University of California, Berkeley.
00009 *     December 31, 1998
00010 *
00011 *     .. Scalar Arguments ..
00012       CHARACTER          EQUED, FACT, TRANS
00013       INTEGER            IA, IAF, IB, INFO, IX, JA, JAF, JB, JX, LIWORK,
00014      $                   LWORK, N, NRHS
00015       REAL               RCOND
00016 *     ..
00017 *     .. Array Arguments ..
00018       INTEGER            DESCA( * ), DESCAF( * ), DESCB( * ),
00019      $                   DESCX( * ), IPIV( * ), IWORK( * )
00020       REAL               A( * ), AF( * ), B( * ), BERR( * ), C( * ),
00021      $                   FERR( * ), R( * ), WORK( * ), X( * )
00022 *     ..
00023 *
00024 *  Purpose
00025 *  =======
00026 *
00027 *  PSGESVX uses the LU factorization to compute the solution to a real
00028 *  system of linear equations
00029 *
00030 *        A(IA:IA+N-1,JA:JA+N-1) * X = B(IB:IB+N-1,JB:JB+NRHS-1),
00031 *
00032 *  where A(IA:IA+N-1,JA:JA+N-1) is an N-by-N matrix and X and
00033 *  B(IB:IB+N-1,JB:JB+NRHS-1) are N-by-NRHS matrices.
00034 *
00035 *  Error bounds on the solution and a condition estimate are also
00036 *  provided.
00037 *
00038 *  Notes
00039 *  =====
00040 *
00041 *  Each global data object is described by an associated description
00042 *  vector.  This vector stores the information required to establish
00043 *  the mapping between an object element and its corresponding process
00044 *  and memory location.
00045 *
00046 *  Let A be a generic term for any 2D block cyclicly distributed array.
00047 *  Such a global array has an associated description vector DESCA.
00048 *  In the following comments, the character _ should be read as
00049 *  "of the global array".
00050 *
00051 *  NOTATION        STORED IN      EXPLANATION
00052 *  --------------- -------------- --------------------------------------
00053 *  DTYPE_A(global) DESCA( DTYPE_ )The descriptor type.  In this case,
00054 *                                 DTYPE_A = 1.
00055 *  CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
00056 *                                 the BLACS process grid A is distribu-
00057 *                                 ted over. The context itself is glo-
00058 *                                 bal, but the handle (the integer
00059 *                                 value) may vary.
00060 *  M_A    (global) DESCA( M_ )    The number of rows in the global
00061 *                                 array A.
00062 *  N_A    (global) DESCA( N_ )    The number of columns in the global
00063 *                                 array A.
00064 *  MB_A   (global) DESCA( MB_ )   The blocking factor used to distribute
00065 *                                 the rows of the array.
00066 *  NB_A   (global) DESCA( NB_ )   The blocking factor used to distribute
00067 *                                 the columns of the array.
00068 *  RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
00069 *                                 row of the array A is distributed.
00070 *  CSRC_A (global) DESCA( CSRC_ ) The process column over which the
00071 *                                 first column of the array A is
00072 *                                 distributed.
00073 *  LLD_A  (local)  DESCA( LLD_ )  The leading dimension of the local
00074 *                                 array.  LLD_A >= MAX(1,LOCr(M_A)).
00075 *
00076 *  Let K be the number of rows or columns of a distributed matrix,
00077 *  and assume that its process grid has dimension p x q.
00078 *  LOCr( K ) denotes the number of elements of K that a process
00079 *  would receive if K were distributed over the p processes of its
00080 *  process column.
00081 *  Similarly, LOCc( K ) denotes the number of elements of K that a
00082 *  process would receive if K were distributed over the q processes of
00083 *  its process row.
00084 *  The values of LOCr() and LOCc() may be determined via a call to the
00085 *  ScaLAPACK tool function, NUMROC:
00086 *          LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
00087 *          LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
00088 *  An upper bound for these quantities may be computed by:
00089 *          LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
00090 *          LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
00091 *
00092 *  Description
00093 *  ===========
00094 *
00095 *  In the following description, A denotes A(IA:IA+N-1,JA:JA+N-1),
00096 *  B denotes B(IB:IB+N-1,JB:JB+NRHS-1) and X denotes
00097 *  X(IX:IX+N-1,JX:JX+NRHS-1).
00098 *
00099 *  The following steps are performed:
00100 *
00101 *  1. If FACT = 'E', real scaling factors are computed to equilibrate
00102 *     the system:
00103 *        TRANS = 'N':  diag(R)*A*diag(C)     *inv(diag(C))*X = diag(R)*B
00104 *        TRANS = 'T': (diag(R)*A*diag(C))**T *inv(diag(R))*X = diag(C)*B
00105 *        TRANS = 'C': (diag(R)*A*diag(C))**H *inv(diag(R))*X = diag(C)*B
00106 *     Whether or not the system will be equilibrated depends on the
00107 *     scaling of the matrix A, but if equilibration is used, A is
00108 *     overwritten by diag(R)*A*diag(C) and B by diag(R)*B (if TRANS='N')
00109 *     or diag(C)*B (if TRANS = 'T' or 'C').
00110 *
00111 *  2. If FACT = 'N' or 'E', the LU decomposition is used to factor the
00112 *     matrix A (after equilibration if FACT = 'E') as
00113 *        A = P * L * U,
00114 *     where P is a permutation matrix, L is a unit lower triangular
00115 *     matrix, and U is upper triangular.
00116 *
00117 *  3. The factored form of A is used to estimate the condition number
00118 *     of the matrix A.  If the reciprocal of the condition number is
00119 *     less than machine precision, steps 4-6 are skipped.
00120 *
00121 *  4. The system of equations is solved for X using the factored form
00122 *     of A.
00123 *
00124 *  5. Iterative refinement is applied to improve the computed solution
00125 *     matrix and calculate error bounds and backward error estimates
00126 *     for it.
00127 *
00128 *  6. If FACT = 'E' and equilibration was used, the matrix X is
00129 *     premultiplied by diag(C) (if TRANS = 'N') or diag(R) (if
00130 *     TRANS = 'T' or 'C') so that it solves the original system
00131 *     before equilibration.
00132 *
00133 *  Arguments
00134 *  =========
00135 *
00136 *  FACT    (global input) CHARACTER
00137 *          Specifies whether or not the factored form of the matrix
00138 *          A(IA:IA+N-1,JA:JA+N-1) is supplied on entry, and if not,
00139 *          whether the matrix A(IA:IA+N-1,JA:JA+N-1) should be
00140 *          equilibrated before it is factored.
00141 *          = 'F':  On entry, AF(IAF:IAF+N-1,JAF:JAF+N-1) and IPIV con-
00142 *                  tain the factored form of A(IA:IA+N-1,JA:JA+N-1).
00143 *                  If EQUED is not 'N', the matrix
00144 *                  A(IA:IA+N-1,JA:JA+N-1) has been equilibrated with
00145 *                  scaling factors given by R and C.
00146 *                  A(IA:IA+N-1,JA:JA+N-1), AF(IAF:IAF+N-1,JAF:JAF+N-1),
00147 *                  and IPIV are not modified.
00148 *          = 'N':  The matrix A(IA:IA+N-1,JA:JA+N-1) will be copied to
00149 *                  AF(IAF:IAF+N-1,JAF:JAF+N-1) and factored.
00150 *          = 'E':  The matrix A(IA:IA+N-1,JA:JA+N-1) will be equili-
00151 *                  brated if necessary, then copied to
00152 *                  AF(IAF:IAF+N-1,JAF:JAF+N-1) and factored.
00153 *
00154 *  TRANS   (global input) CHARACTER
00155 *          Specifies the form of the system of equations:
00156 *          = 'N':  A(IA:IA+N-1,JA:JA+N-1) * X(IX:IX+N-1,JX:JX+NRHS-1)
00157 *                = B(IB:IB+N-1,JB:JB+NRHS-1)     (No transpose)
00158 *          = 'T':  A(IA:IA+N-1,JA:JA+N-1)**T * X(IX:IX+N-1,JX:JX+NRHS-1)
00159 *                = B(IB:IB+N-1,JB:JB+NRHS-1)  (Transpose)
00160 *          = 'C':  A(IA:IA+N-1,JA:JA+N-1)**H * X(IX:IX+N-1,JX:JX+NRHS-1)
00161 *                = B(IB:IB+N-1,JB:JB+NRHS-1)  (Transpose)
00162 *
00163 *  N       (global input) INTEGER
00164 *          The number of rows and columns to be operated on, i.e. the
00165 *          order of the distributed submatrix A(IA:IA+N-1,JA:JA+N-1).
00166 *          N >= 0.
00167 *
00168 *  NRHS    (global input) INTEGER
00169 *          The number of right-hand sides, i.e., the number of columns
00170 *          of the distributed submatrices B(IB:IB+N-1,JB:JB+NRHS-1) and
00171 *          X(IX:IX+N-1,JX:JX+NRHS-1).  NRHS >= 0.
00172 *
00173 *  A       (local input/local output) REAL pointer into
00174 *          the local memory to an array of local dimension
00175 *          (LLD_A,LOCc(JA+N-1)).  On entry, the N-by-N matrix
00176 *          A(IA:IA+N-1,JA:JA+N-1).  If FACT = 'F' and EQUED is not 'N',
00177 *          then A(IA:IA+N-1,JA:JA+N-1) must have been equilibrated by
00178 *          the scaling factors in R and/or C.  A(IA:IA+N-1,JA:JA+N-1) is
00179 *          not modified if FACT = 'F' or  'N', or if FACT = 'E' and
00180 *          EQUED = 'N' on exit.
00181 *
00182 *          On exit, if EQUED .ne. 'N', A(IA:IA+N-1,JA:JA+N-1) is scaled
00183 *          as follows:
00184 *          EQUED = 'R':  A(IA:IA+N-1,JA:JA+N-1) :=
00185 *                                      diag(R) * A(IA:IA+N-1,JA:JA+N-1)
00186 *          EQUED = 'C':  A(IA:IA+N-1,JA:JA+N-1) :=
00187 *                                      A(IA:IA+N-1,JA:JA+N-1) * diag(C)
00188 *          EQUED = 'B':  A(IA:IA+N-1,JA:JA+N-1) :=
00189 *                            diag(R) * A(IA:IA+N-1,JA:JA+N-1) * diag(C).
00190 *
00191 *  IA      (global input) INTEGER
00192 *          The row index in the global array A indicating the first
00193 *          row of sub( A ).
00194 *
00195 *  JA      (global input) INTEGER
00196 *          The column index in the global array A indicating the
00197 *          first column of sub( A ).
00198 *
00199 *  DESCA   (global and local input) INTEGER array of dimension DLEN_.
00200 *          The array descriptor for the distributed matrix A.
00201 *
00202 *  AF      (local input or local output) REAL pointer
00203 *          into the local memory to an array of local dimension
00204 *          (LLD_AF,LOCc(JA+N-1)).  If FACT = 'F', then
00205 *          AF(IAF:IAF+N-1,JAF:JAF+N-1) is an input argument and on
00206 *          entry contains the factors L and U from the factorization
00207 *          A(IA:IA+N-1,JA:JA+N-1) = P*L*U as computed by PSGETRF.
00208 *          If EQUED .ne. 'N', then AF is the factored form of the
00209 *          equilibrated matrix A(IA:IA+N-1,JA:JA+N-1).
00210 *
00211 *          If FACT = 'N', then AF(IAF:IAF+N-1,JAF:JAF+N-1) is an output
00212 *          argument and on exit returns the factors L and U from the
00213 *          factorization A(IA:IA+N-1,JA:JA+N-1) = P*L*U of the original
00214 *          matrix A(IA:IA+N-1,JA:JA+N-1).
00215 *
00216 *          If FACT = 'E', then AF(IAF:IAF+N-1,JAF:JAF+N-1) is an output
00217 *          argument and on exit returns the factors L and U from the
00218 *          factorization A(IA:IA+N-1,JA:JA+N-1) = P*L*U of the equili-
00219 *          brated matrix A(IA:IA+N-1,JA:JA+N-1) (see the description of
00220 *          A(IA:IA+N-1,JA:JA+N-1) for the form of the equilibrated
00221 *          matrix).
00222 *
00223 *  IAF     (global input) INTEGER
00224 *          The row index in the global array AF indicating the first
00225 *          row of sub( AF ).
00226 *
00227 *  JAF     (global input) INTEGER
00228 *          The column index in the global array AF indicating the
00229 *          first column of sub( AF ).
00230 *
00231 *  DESCAF  (global and local input) INTEGER array of dimension DLEN_.
00232 *          The array descriptor for the distributed matrix AF.
00233 *
00234 *  IPIV    (local input or local output) INTEGER array, dimension
00235 *          LOCr(M_A)+MB_A. If FACT = 'F', then IPIV is an input argu-
00236 *          ment and on entry contains the pivot indices from the fac-
00237 *          torization A(IA:IA+N-1,JA:JA+N-1) = P*L*U as computed by
00238 *          PSGETRF; IPIV(i) -> The global row local row i was
00239 *          swapped with.  This array must be aligned with
00240 *          A( IA:IA+N-1, * ).
00241 *
00242 *          If FACT = 'N', then IPIV is an output argument and on exit
00243 *          contains the pivot indices from the factorization
00244 *          A(IA:IA+N-1,JA:JA+N-1) = P*L*U of the original matrix
00245 *          A(IA:IA+N-1,JA:JA+N-1).
00246 *
00247 *          If FACT = 'E', then IPIV is an output argument and on exit
00248 *          contains the pivot indices from the factorization
00249 *          A(IA:IA+N-1,JA:JA+N-1) = P*L*U of the equilibrated matrix
00250 *          A(IA:IA+N-1,JA:JA+N-1).
00251 *
00252 *  EQUED   (global input or global output) CHARACTER
00253 *          Specifies the form of equilibration that was done.
00254 *          = 'N':  No equilibration (always true if FACT = 'N').
00255 *          = 'R':  Row equilibration, i.e., A(IA:IA+N-1,JA:JA+N-1) has
00256 *                  been premultiplied by diag(R).
00257 *          = 'C':  Column equilibration, i.e., A(IA:IA+N-1,JA:JA+N-1)
00258 *                  has been postmultiplied by diag(C).
00259 *          = 'B':  Both row and column equilibration, i.e.,
00260 *                  A(IA:IA+N-1,JA:JA+N-1) has been replaced by
00261 *                  diag(R) * A(IA:IA+N-1,JA:JA+N-1) * diag(C).
00262 *          EQUED is an input variable if FACT = 'F'; otherwise, it is an
00263 *          output variable.
00264 *
00265 *  R       (local input or local output) REAL array,
00266 *                                                  dimension LOCr(M_A).
00267 *          The row scale factors for A(IA:IA+N-1,JA:JA+N-1).
00268 *          If EQUED = 'R' or 'B', A(IA:IA+N-1,JA:JA+N-1) is multiplied
00269 *          on the left by diag(R); if EQUED='N' or 'C', R is not acces-
00270 *          sed.  R is an input variable if FACT = 'F'; otherwise, R is
00271 *          an output variable.
00272 *          If FACT = 'F' and EQUED = 'R' or 'B', each element of R must
00273 *          be positive.
00274 *          R is replicated in every process column, and is aligned
00275 *          with the distributed matrix A.
00276 *
00277 *  C       (local input or local output) REAL array,
00278 *                                                  dimension LOCc(N_A).
00279 *          The column scale factors for A(IA:IA+N-1,JA:JA+N-1).
00280 *          If EQUED = 'C' or 'B', A(IA:IA+N-1,JA:JA+N-1) is multiplied
00281 *          on the right by diag(C); if EQUED = 'N' or 'R', C is not
00282 *          accessed.  C is an input variable if FACT = 'F'; otherwise,
00283 *          C is an output variable.  If FACT = 'F' and EQUED = 'C' or
00284 *          'B', each element of C must be positive.
00285 *          C is replicated in every process row, and is aligned with
00286 *          the distributed matrix A.
00287 *
00288 *  B       (local input/local output) REAL pointer
00289 *          into the local memory to an array of local dimension
00290 *          (LLD_B,LOCc(JB+NRHS-1) ).  On entry, the N-by-NRHS right-hand
00291 *          side matrix B(IB:IB+N-1,JB:JB+NRHS-1). On exit, if
00292 *          EQUED = 'N', B(IB:IB+N-1,JB:JB+NRHS-1) is not modified; if
00293 *          TRANS = 'N' and EQUED = 'R' or 'B', B is overwritten by
00294 *          diag(R)*B(IB:IB+N-1,JB:JB+NRHS-1); if TRANS = 'T' or 'C'
00295 *          and EQUED = 'C' or 'B', B(IB:IB+N-1,JB:JB+NRHS-1) is over-
00296 *          written by diag(C)*B(IB:IB+N-1,JB:JB+NRHS-1).
00297 *
00298 *  IB      (global input) INTEGER
00299 *          The row index in the global array B indicating the first
00300 *          row of sub( B ).
00301 *
00302 *  JB      (global input) INTEGER
00303 *          The column index in the global array B indicating the
00304 *          first column of sub( B ).
00305 *
00306 *  DESCB   (global and local input) INTEGER array of dimension DLEN_.
00307 *          The array descriptor for the distributed matrix B.
00308 *
00309 *  X       (local input/local output) REAL pointer
00310 *          into the local memory to an array of local dimension
00311 *          (LLD_X, LOCc(JX+NRHS-1)).  If INFO = 0, the N-by-NRHS
00312 *          solution matrix X(IX:IX+N-1,JX:JX+NRHS-1) to the original
00313 *          system of equations.  Note that A(IA:IA+N-1,JA:JA+N-1) and
00314 *          B(IB:IB+N-1,JB:JB+NRHS-1) are modified on exit if
00315 *          EQUED .ne. 'N', and the solution to the equilibrated system
00316 *          is inv(diag(C))*X(IX:IX+N-1,JX:JX+NRHS-1) if TRANS = 'N'
00317 *          and EQUED = 'C' or 'B', or
00318 *          inv(diag(R))*X(IX:IX+N-1,JX:JX+NRHS-1) if TRANS = 'T' or 'C'
00319 *          and EQUED = 'R' or 'B'.
00320 *
00321 *  IX      (global input) INTEGER
00322 *          The row index in the global array X indicating the first
00323 *          row of sub( X ).
00324 *
00325 *  JX      (global input) INTEGER
00326 *          The column index in the global array X indicating the
00327 *          first column of sub( X ).
00328 *
00329 *  DESCX   (global and local input) INTEGER array of dimension DLEN_.
00330 *          The array descriptor for the distributed matrix X.
00331 *
00332 *  RCOND   (global output) REAL
00333 *          The estimate of the reciprocal condition number of the matrix
00334 *          A(IA:IA+N-1,JA:JA+N-1) after equilibration (if done).  If
00335 *          RCOND is less than the machine precision (in particular, if
00336 *          RCOND = 0), the matrix is singular to working precision.
00337 *          This condition is indicated by a return code of INFO > 0.
00338 *
00339 *  FERR    (local output) REAL array, dimension LOCc(N_B)
00340 *          The estimated forward error bounds for each solution vector
00341 *          X(j) (the j-th column of the solution matrix
00342 *          X(IX:IX+N-1,JX:JX+NRHS-1). If XTRUE is the true solution,
00343 *          FERR(j) bounds the magnitude of the largest entry in
00344 *          (X(j) - XTRUE) divided by the magnitude of the largest entry
00345 *          in X(j).  The estimate is as reliable as the estimate for
00346 *          RCOND, and is almost always a slight overestimate of the
00347 *          true error.  FERR is replicated in every process row, and is
00348 *          aligned with the matrices B and X.
00349 *
00350 *  BERR    (local output) REAL array, dimension LOCc(N_B).
00351 *          The componentwise relative backward error of each solution
00352 *          vector X(j) (i.e., the smallest relative change in
00353 *          any entry of A(IA:IA+N-1,JA:JA+N-1) or
00354 *          B(IB:IB+N-1,JB:JB+NRHS-1) that makes X(j) an exact solution).
00355 *          BERR is replicated in every process row, and is aligned
00356 *          with the matrices B and X.
00357 *
00358 *  WORK    (local workspace/local output) REAL array,
00359 *                                                    dimension (LWORK)
00360 *          On exit, WORK(1) returns the minimal and optimal LWORK.
00361 *
00362 *  LWORK   (local or global input) INTEGER
00363 *          The dimension of the array WORK.
00364 *          LWORK is local input and must be at least
00365 *          LWORK = MAX( PSGECON( LWORK ), PSGERFS( LWORK ) )
00366 *                  + LOCr( N_A ).
00367 *
00368 *          If LWORK = -1, then LWORK is global input and a workspace
00369 *          query is assumed; the routine only calculates the minimum
00370 *          and optimal size for all work arrays. Each of these
00371 *          values is returned in the first entry of the corresponding
00372 *          work array, and no error message is issued by PXERBLA.
00373 *
00374 *  IWORK   (local workspace/local output) INTEGER array,
00375 *                                                  dimension (LIWORK)
00376 *          On exit, IWORK(1) returns the minimal and optimal LIWORK.
00377 *
00378 *  LIWORK  (local or global input) INTEGER
00379 *          The dimension of the array IWORK.
00380 *          LIWORK is local input and must be at least
00381 *          LIWORK = LOCr(N_A).
00382 *
00383 *          If LIWORK = -1, then LIWORK is global input and a workspace
00384 *          query is assumed; the routine only calculates the minimum
00385 *          and optimal size for all work arrays. Each of these
00386 *          values is returned in the first entry of the corresponding
00387 *          work array, and no error message is issued by PXERBLA.
00388 *
00389 *
00390 *  INFO    (global output) INTEGER
00391 *          = 0:  successful exit
00392 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00393 *          > 0:  if INFO = i, and i is
00394 *                <= N:  U(IA+I-1,IA+I-1) is exactly zero.  The
00395 *                       factorization has been completed, but the
00396 *                       factor U is exactly singular, so the solution
00397 *                       and error bounds could not be computed.
00398 *                = N+1: RCOND is less than machine precision.  The
00399 *                       factorization has been completed, but the
00400 *                       matrix is singular to working precision, and
00401 *                       the solution and error bounds have not been
00402 *                       computed.
00403 *
00404 *  =====================================================================
00405 *
00406 *     .. Parameters ..
00407       INTEGER            BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
00408      $                   LLD_, MB_, M_, NB_, N_, RSRC_
00409       PARAMETER          ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1,
00410      $                     CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6,
00411      $                     RSRC_ = 7, CSRC_ = 8, LLD_ = 9 )
00412       REAL               ONE, ZERO
00413       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
00414 *     ..
00415 *     .. Local Scalars ..
00416       LOGICAL            COLEQU, EQUIL, LQUERY, NOFACT, NOTRAN, ROWEQU
00417       CHARACTER          NORM
00418       INTEGER            CONWRK, I, IACOL, IAROW, IAFROW, IBROW, IBCOL,
00419      $                   ICOFFA, ICOFFB, ICOFFX, ICTXT, IDUMM,
00420      $                   IIA, IIB, IIX,
00421      $                   INFEQU, IROFFA, IROFFAF, IROFFB,
00422      $                   IROFFX, IXCOL, IXROW, J, JJA, JJB, JJX,
00423      $                   LCM, LCMQ,
00424      $                   LIWMIN, LWMIN, MYCOL, MYROW, NP, NPCOL, NPROW,
00425      $                   NQ, NQB, NRHSQ, RFSWRK
00426       REAL               AMAX, ANORM, BIGNUM, COLCND, RCMAX, RCMIN,
00427      $                   ROWCND, SMLNUM
00428 *     ..
00429 *     .. Local Arrays ..
00430       INTEGER            CDESC( DLEN_ ), IDUM1( 5 ), IDUM2( 5 )
00431 *     ..
00432 *     .. External Subroutines ..
00433       EXTERNAL           BLACS_GRIDINFO, CHK1MAT, DESCSET, PCHK2MAT,
00434      $                   INFOG2L, PSGECON, PSGEEQU, PSGERFS,
00435      $                   PSGETRF, PSGETRS, PSLACPY,
00436      $                   PSLAQGE, PSCOPY, PXERBLA, SGEBR2D,
00437      $                   SGEBS2D, SGAMN2D, SGAMX2D
00438 *     ..
00439 *     .. External Functions ..
00440       LOGICAL            LSAME
00441       INTEGER            ICEIL, ILCM, INDXG2P, NUMROC
00442       REAL               PSLAMCH, PSLANGE
00443       EXTERNAL           ICEIL, ILCM, INDXG2P, LSAME, NUMROC, PSLANGE,
00444      $                   PSLAMCH
00445 *     ..
00446 *     .. Intrinsic Functions ..
00447       INTRINSIC          ICHAR, MAX, MIN, MOD, REAL
00448 *     ..
00449 *     .. Executable Statements ..
00450 *
00451 *     Get grid parameters
00452 *
00453       ICTXT = DESCA( CTXT_ )
00454       CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )
00455 *
00456 *     Test the input parameters
00457 *
00458       INFO = 0
00459       IF( NPROW.EQ.-1 ) THEN
00460          INFO = -(800+CTXT_)
00461       ELSE
00462          CALL CHK1MAT( N, 3, N, 3, IA, JA, DESCA, 8, INFO )
00463          IF( LSAME( FACT, 'F' ) )
00464      $      CALL CHK1MAT( N, 3, N, 3, IAF, JAF, DESCAF, 12, INFO )
00465          CALL CHK1MAT( N, 3, NRHS, 4, IB, JB, DESCB, 20, INFO )
00466          CALL CHK1MAT( N, 3, NRHS, 4, IX, JX, DESCX, 24, INFO )
00467          NOFACT = LSAME( FACT, 'N' )
00468          EQUIL = LSAME( FACT, 'E' )
00469          NOTRAN = LSAME( TRANS, 'N' )
00470          IF( NOFACT .OR. EQUIL ) THEN
00471             EQUED = 'N'
00472             ROWEQU = .FALSE.
00473             COLEQU = .FALSE.
00474          ELSE
00475             ROWEQU = LSAME( EQUED, 'R' ) .OR. LSAME( EQUED, 'B' )
00476             COLEQU = LSAME( EQUED, 'C' ) .OR. LSAME( EQUED, 'B' )
00477             SMLNUM = PSLAMCH( ICTXT, 'Safe minimum' )
00478             BIGNUM = ONE / SMLNUM
00479          END IF
00480          IF( INFO.EQ.0 ) THEN
00481             IAROW = INDXG2P( IA, DESCA( MB_ ), MYROW, DESCA( RSRC_ ),
00482      $                       NPROW )
00483             IAFROW = INDXG2P( IAF, DESCAF( MB_ ), MYROW,
00484      $                        DESCAF( RSRC_ ), NPROW )
00485             IBROW = INDXG2P( IB, DESCB( MB_ ), MYROW, DESCB( RSRC_ ),
00486      $                       NPROW )
00487             IXROW = INDXG2P( IX, DESCX( MB_ ), MYROW, DESCX( RSRC_ ),
00488      $                       NPROW )
00489             IROFFA = MOD( IA-1, DESCA( MB_ ) )
00490             IROFFAF = MOD( IAF-1, DESCAF( MB_ ) )
00491             ICOFFA = MOD( JA-1, DESCA( NB_ ) )
00492             IROFFB = MOD( IB-1, DESCB( MB_ ) )
00493             ICOFFB = MOD( JB-1, DESCB( NB_ ) )
00494             IROFFX = MOD( IX-1, DESCX( MB_ ) )
00495             ICOFFX = MOD( JX-1, DESCX( NB_ ) )
00496             CALL INFOG2L( IA, JA, DESCA, NPROW, NPCOL, MYROW,
00497      $                    MYCOL, IIA, JJA, IAROW, IACOL )
00498             NP = NUMROC( N+IROFFA, DESCA( MB_ ), MYROW, IAROW,
00499      $                   NPROW )
00500             IF( MYROW.EQ.IAROW )
00501      $         NP = NP-IROFFA
00502             NQ = NUMROC( N+ICOFFA, DESCA( NB_ ), MYCOL, IACOL,
00503      $                   NPCOL )
00504             IF( MYCOL.EQ.IACOL )
00505      $         NQ = NQ-ICOFFA
00506             NQB = ICEIL( N+IROFFA, DESCA( NB_ )*NPCOL )
00507             LCM = ILCM( NPROW, NPCOL )
00508             LCMQ = LCM / NPCOL
00509             CONWRK = 2*NP + 2*NQ + MAX( 2, MAX( DESCA( NB_ )*
00510      $               MAX( 1, ICEIL( NPROW-1, NPCOL ) ), NQ +
00511      $               DESCA( NB_ )*
00512      $               MAX( 1, ICEIL( NPCOL-1, NPROW ) ) ) )
00513             RFSWRK = 3*NP
00514             IF( LSAME( TRANS, 'N' ) ) THEN
00515                RFSWRK = RFSWRK + NP + NQ +
00516      $                 ICEIL( NQB, LCMQ )*DESCA( NB_ )
00517             ELSE IF( LSAME( TRANS, 'T' ).OR.LSAME( TRANS, 'C' ) ) THEN
00518                RFSWRK = RFSWRK + NP + NQ
00519             END IF
00520             LWMIN = MAX( CONWRK, RFSWRK )
00521             WORK( 1 ) = REAL( LWMIN )
00522             LIWMIN = NP
00523             IWORK( 1 ) = LIWMIN
00524             IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND.
00525      $          .NOT.LSAME( FACT, 'F' ) ) THEN
00526                INFO = -1
00527             ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND.
00528      $         .NOT. LSAME( TRANS, 'C' ) ) THEN
00529                INFO = -2
00530             ELSE IF( IROFFA.NE.0 ) THEN
00531                INFO = -6
00532             ELSE IF( ICOFFA.NE.0 .OR. IROFFA.NE.ICOFFA ) THEN
00533                INFO = -7
00534             ELSE IF( DESCA( MB_ ).NE.DESCA( NB_ ) ) THEN
00535                INFO = -(800+NB_)
00536             ELSE IF( IAFROW.NE.IAROW ) THEN
00537                INFO = -10
00538             ELSE IF( IROFFAF.NE.0 ) THEN
00539                INFO = -10
00540             ELSE IF( ICTXT.NE.DESCAF( CTXT_ ) ) THEN
00541                INFO = -(1200+CTXT_)
00542             ELSE IF( LSAME( FACT, 'F' ) .AND. .NOT.
00543      $         ( ROWEQU .OR. COLEQU .OR. LSAME( EQUED, 'N' ) ) ) THEN
00544                INFO = -13
00545             ELSE
00546                IF( ROWEQU ) THEN
00547                   RCMIN = BIGNUM
00548                   RCMAX = ZERO
00549                   DO 10 J = IIA, IIA + NP - 1
00550                      RCMIN = MIN( RCMIN, R( J ) )
00551                      RCMAX = MAX( RCMAX, R( J ) )
00552    10             CONTINUE
00553                   CALL SGAMN2D( ICTXT, 'Columnwise', ' ', 1, 1, RCMIN,
00554      $                          1, IDUMM, IDUMM, -1, -1, MYCOL )
00555                   CALL SGAMX2D( ICTXT, 'Columnwise', ' ', 1, 1, RCMAX,
00556      $                          1, IDUMM, IDUMM, -1, -1, MYCOL )
00557                   IF( RCMIN.LE.ZERO ) THEN
00558                      INFO = -14
00559                   ELSE IF( N.GT.0 ) THEN
00560                      ROWCND = MAX( RCMIN, SMLNUM ) /
00561      $                        MIN( RCMAX, BIGNUM )
00562                   ELSE
00563                      ROWCND = ONE
00564                   END IF
00565                END IF
00566                IF( COLEQU .AND. INFO.EQ.0 ) THEN
00567                   RCMIN = BIGNUM
00568                   RCMAX = ZERO
00569                   DO 20 J = JJA, JJA+NQ-1
00570                      RCMIN = MIN( RCMIN, C( J ) )
00571                      RCMAX = MAX( RCMAX, C( J ) )
00572    20             CONTINUE
00573                   CALL SGAMN2D( ICTXT, 'Rowwise', ' ', 1, 1, RCMIN,
00574      $                          1, IDUMM, IDUMM, -1, -1, MYCOL )
00575                   CALL SGAMX2D( ICTXT, 'Rowwise', ' ', 1, 1, RCMAX,
00576      $                          1, IDUMM, IDUMM, -1, -1, MYCOL )
00577                   IF( RCMIN.LE.ZERO ) THEN
00578                      INFO = -15
00579                   ELSE IF( N.GT.0 ) THEN
00580                      COLCND = MAX( RCMIN, SMLNUM ) /
00581      $                        MIN( RCMAX, BIGNUM )
00582                   ELSE
00583                      COLCND = ONE
00584                   END IF
00585                END IF
00586             END IF
00587          END IF
00588 *
00589          WORK( 1 ) = REAL( LWMIN )
00590          IWORK( 1 ) = LIWMIN
00591          LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
00592          IF( INFO.EQ.0 ) THEN
00593             IF( IBROW.NE.IAROW ) THEN
00594                INFO = -18
00595             ELSE IF( IXROW.NE.IBROW ) THEN
00596                INFO = -22
00597             ELSE IF( DESCB( MB_ ).NE.DESCA( NB_ ) ) THEN
00598                INFO = -(2000+NB_)
00599             ELSE IF( ICTXT.NE.DESCB( CTXT_ ) ) THEN
00600                INFO = -(2000+CTXT_)
00601             ELSE IF( DESCX( MB_ ).NE.DESCA( NB_ ) ) THEN
00602                INFO = -(2400+NB_)
00603             ELSE IF( ICTXT.NE.DESCX( CTXT_ ) ) THEN
00604                INFO = -(2400+CTXT_)
00605             ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
00606                INFO = -29
00607             ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
00608                INFO = -31
00609             END IF
00610             IDUM1( 1 ) = ICHAR( FACT )
00611             IDUM2( 1 ) = 1
00612             IDUM1( 2 ) = ICHAR( TRANS )
00613             IDUM2( 2 ) = 2
00614             IF( LSAME( FACT, 'F' ) ) THEN
00615                IDUM1( 3 ) = ICHAR( EQUED )
00616                IDUM2( 3 ) = 14
00617                IF( LWORK.EQ.-1 ) THEN
00618                   IDUM1( 4 ) = -1
00619                ELSE
00620                   IDUM1( 4 ) = 1
00621                END IF
00622                IDUM2( 4 ) = 29
00623                IF( LIWORK.EQ.-1 ) THEN
00624                   IDUM1( 5 ) = -1
00625                ELSE
00626                   IDUM1( 5 ) = 1
00627                END IF
00628                IDUM2( 5 ) = 31
00629                CALL PCHK2MAT( N, 3, N, 3, IA, JA, DESCA, 8, N, 3,
00630      $                        NRHS, 4, IB, JB, DESCB, 20, 5, IDUM1,
00631      $                        IDUM2, INFO )
00632             ELSE
00633                IF( LWORK.EQ.-1 ) THEN
00634                   IDUM1( 3 ) = -1
00635                ELSE
00636                   IDUM1( 3 ) = 1
00637                END IF
00638                IDUM2( 3 ) = 29
00639                IF( LIWORK.EQ.-1 ) THEN
00640                   IDUM1( 4 ) = -1
00641                ELSE
00642                   IDUM1( 4 ) = 1
00643                END IF
00644                IDUM2( 4 ) = 31
00645                CALL PCHK2MAT( N, 3, N, 3, IA, JA, DESCA, 8, N, 3,
00646      $                        NRHS, 4, IB, JB, DESCB, 20, 4, IDUM1,
00647      $                        IDUM2, INFO )
00648             END IF
00649          END IF
00650       END IF
00651 *
00652       IF( INFO.NE.0 ) THEN
00653          CALL PXERBLA( ICTXT, 'PSGESVX', -INFO )
00654          RETURN
00655       ELSE IF( LQUERY ) THEN
00656          RETURN
00657       END IF
00658 *
00659       IF( EQUIL ) THEN
00660 *
00661 *        Compute row and column scalings to equilibrate the matrix A.
00662 *
00663          CALL PSGEEQU( N, N, A, IA, JA, DESCA, R, C, ROWCND, COLCND,
00664      $                 AMAX, INFEQU )
00665          IF( INFEQU.EQ.0 ) THEN
00666 *
00667 *           Equilibrate the matrix.
00668 *
00669             CALL PSLAQGE( N, N, A, IA, JA, DESCA, R, C, ROWCND, COLCND,
00670      $                    AMAX, EQUED )
00671             ROWEQU = LSAME( EQUED, 'R' ) .OR. LSAME( EQUED, 'B' )
00672             COLEQU = LSAME( EQUED, 'C' ) .OR. LSAME( EQUED, 'B' )
00673          END IF
00674       END IF
00675 *
00676 *     Scale the right-hand side.
00677 *
00678       CALL INFOG2L( IB, JB, DESCB, NPROW, NPCOL, MYROW, MYCOL, IIB,
00679      $              JJB, IBROW, IBCOL )
00680       NP = NUMROC( N+IROFFB, DESCB( MB_ ), MYROW, IBROW, NPROW )
00681       NRHSQ = NUMROC( NRHS+ICOFFB, DESCB( NB_ ), MYCOL, IBCOL, NPCOL )
00682       IF( MYROW.EQ.IBROW )
00683      $   NP = NP-IROFFB
00684       IF( MYCOL.EQ.IBCOL )
00685      $   NRHSQ = NRHSQ-ICOFFB
00686 *
00687       IF( NOTRAN ) THEN
00688          IF( ROWEQU ) THEN
00689             DO 40 J = JJB, JJB+NRHSQ-1
00690                DO 30 I = IIB, IIB+NP-1
00691                   B( I+( J-1 )*DESCB( LLD_ ) ) = R( I )*
00692      $                                     B( I+( J-1 )*DESCB( LLD_ ) )
00693    30          CONTINUE
00694    40       CONTINUE
00695          END IF
00696       ELSE IF( COLEQU ) THEN
00697 *
00698 *        Transpose the Column scale factors
00699 *
00700          CALL DESCSET( CDESC, 1, N+ICOFFA, 1, DESCA( NB_ ), MYROW,
00701      $                 IACOL, ICTXT, 1 )
00702          CALL PSCOPY( N, C, 1, JA, CDESC, CDESC( LLD_ ), WORK, IB, JB,
00703      $                DESCB, 1 )
00704          IF( MYCOL.EQ.IBCOL ) THEN
00705             CALL SGEBS2D( ICTXT, 'Rowwise', ' ', NP, 1, WORK( IIB ),
00706      $                    DESCB( LLD_ ) )
00707          ELSE
00708             CALL SGEBR2D( ICTXT, 'Rowwise', ' ', NP, 1, WORK( IIB ),
00709      $                    DESCB( LLD_ ), MYROW, IBCOL )
00710          END IF
00711          DO 60 J = JJB, JJB+NRHSQ-1
00712             DO 50 I = IIB, IIB+NP-1
00713                B( I+( J-1 )*DESCB( LLD_ ) ) = WORK( I )*
00714      $                                     B( I+( J-1 )*DESCB( LLD_ ) )
00715    50       CONTINUE
00716    60    CONTINUE
00717       END IF
00718 *
00719       IF( NOFACT.OR.EQUIL ) THEN
00720 *
00721 *        Compute the LU factorization of A.
00722 *
00723          CALL PSLACPY( 'Full', N, N, A, IA, JA, DESCA, AF, IAF, JAF,
00724      $                 DESCAF )
00725          CALL PSGETRF( N, N, AF, IAF, JAF, DESCAF, IPIV, INFO )
00726 *
00727 *        Return if INFO is non-zero.
00728 *
00729          IF( INFO.NE.0 ) THEN
00730             IF( INFO.GT.0 )
00731      $         RCOND = ZERO
00732             RETURN
00733          END IF
00734       END IF
00735 *
00736 *     Compute the norm of the matrix A.
00737 *
00738       IF( NOTRAN ) THEN
00739          NORM = '1'
00740       ELSE
00741          NORM = 'I'
00742       END IF
00743       ANORM = PSLANGE( NORM, N, N, A, IA, JA, DESCA, WORK )
00744 *
00745 *     Compute the reciprocal of the condition number of A.
00746 *
00747       CALL PSGECON( NORM, N, AF, IAF, JAF, DESCAF, ANORM, RCOND, WORK,
00748      $              LWORK, IWORK, LIWORK, INFO )
00749 *
00750 *     Return if the matrix is singular to working precision.
00751 *
00752       IF( RCOND.LT.PSLAMCH( ICTXT, 'Epsilon' ) ) THEN
00753          INFO = IA + N
00754          RETURN
00755       END IF
00756 *
00757 *     Compute the solution matrix X.
00758 *
00759       CALL PSLACPY( 'Full', N, NRHS, B, IB, JB, DESCB, X, IX, JX,
00760      $              DESCX )
00761       CALL PSGETRS( TRANS, N, NRHS, AF, IAF, JAF, DESCAF, IPIV, X, IX,
00762      $              JX, DESCX, INFO )
00763 *
00764 *     Use iterative refinement to improve the computed solution and
00765 *     compute error bounds and backward error estimates for it.
00766 *
00767       CALL PSGERFS( TRANS, N, NRHS, A, IA, JA, DESCA, AF, IAF, JAF,
00768      $              DESCAF, IPIV, B, IB, JB, DESCB, X, IX, JX, DESCX,
00769      $              FERR, BERR, WORK, LWORK, IWORK, LIWORK, INFO )
00770 *
00771 *     Transform the solution matrix X to a solution of the original
00772 *     system.
00773 *
00774       CALL INFOG2L( IX, JX, DESCX, NPROW, NPCOL, MYROW, MYCOL, IIX,
00775      $              JJX, IXROW, IXCOL )
00776       NP = NUMROC( N+IROFFX, DESCX( MB_ ), MYROW, IXROW, NPROW )
00777       NRHSQ = NUMROC( NRHS+ICOFFX, DESCX( NB_ ), MYCOL, IXCOL, NPCOL )
00778       IF( MYROW.EQ.IBROW )
00779      $   NP = NP-IROFFX
00780       IF( MYCOL.EQ.IBCOL )
00781      $   NRHSQ = NRHSQ-ICOFFX
00782 *
00783       IF( NOTRAN ) THEN
00784          IF( COLEQU ) THEN
00785 *
00786 *           Transpose the column scaling factors
00787 *
00788             CALL DESCSET( CDESC, 1, N+ICOFFA, 1, DESCA( NB_ ), MYROW,
00789      $                    IACOL, ICTXT, 1 )
00790             CALL PSCOPY( N, C, 1, JA, CDESC, CDESC( LLD_ ), WORK, IX,
00791      $                   JX, DESCX, 1 )
00792             IF( MYCOL.EQ.IBCOL ) THEN
00793                 CALL SGEBS2D( ICTXT, 'Rowwise', ' ', NP, 1,
00794      $                        WORK( IIX ), DESCX( LLD_ ) )
00795             ELSE
00796                 CALL SGEBR2D( ICTXT, 'Rowwise', ' ', NP, 1,
00797      $                        WORK( IIX ), DESCX( LLD_ ), MYROW, IBCOL )
00798             END IF
00799 *
00800             DO 80 J = JJX, JJX+NRHSQ-1
00801                DO 70 I = IIX, IIX+NP-1
00802                   X( I+( J-1 )*DESCX( LLD_ ) ) = WORK( I )*
00803      $                                      X( I+( J-1 )*DESCX( LLD_ ) )
00804    70          CONTINUE
00805    80       CONTINUE
00806             DO 90 J = JJX, JJX+NRHSQ-1
00807                FERR( J ) = FERR( J ) / COLCND
00808    90       CONTINUE
00809          END IF
00810       ELSE IF( ROWEQU ) THEN
00811          DO 110 J = JJX, JJX+NRHSQ-1
00812             DO 100 I = IIX, IIX+NP-1
00813                X( I+( J-1 )*DESCX( LLD_ ) ) = R( I )*
00814      $                                     X( I+( J-1 )*DESCX( LLD_ ) )
00815   100       CONTINUE
00816   110    CONTINUE
00817          DO 120 J = JJX, JJX+NRHSQ-1
00818             FERR( J ) = FERR( J ) / ROWCND
00819   120    CONTINUE
00820       END IF
00821 *
00822       WORK( 1 ) = REAL( LWMIN )
00823       IWORK( 1 ) = LIWMIN
00824 *
00825       RETURN
00826 *
00827 *     End of PSGESVX
00828 *
00829       END