Man Linux: Main Page and Category List

NAME

       PZGESVX - use the LU factorization to compute the solution to a complex
       system   of   linear   equations     A(IA:IA+N-1,JA:JA+N-1)   *   X   =
       B(IB:IB+N-1,JB:JB+NRHS-1),

SYNOPSIS

       SUBROUTINE PZGESVX( 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 )

           CHARACTER       EQUED, FACT, TRANS

           INTEGER         IA, IAF, IB, INFO, IX, JA,  JAF,  JB,  JX,  LRWORK,
                           LWORK, N, NRHS

           DOUBLE          PRECISION RCOND

           INTEGER         DESCA(  *  ),  DESCAF( * ), DESCB( * ), DESCX( * ),
                           IPIV( * )

           DOUBLE          PRECISION BERR( * ), C( * ), FERR( *  ),  R(  *  ),
                           RWORK( * )

           COMPLEX*16      A( * ), AF( * ), B( * ), WORK( * ), X( * )

PURPOSE

       PZGESVX  uses the LU factorization to compute the solution to a complex
       system of linear equations

       where  A(IA:IA+N-1,JA:JA+N-1)  is  an   N-by-N   matrix   and   X   and
       B(IB:IB+N-1,JB:JB+NRHS-1) are N-by-NRHS matrices.

       Error  bounds  on  the  solution  and  a  condition  estimate  are also
       provided.

       Notes
       =====

       Each global data object  is  described  by  an  associated  description
       vector.   This  vector stores the information required to establish the
       mapping between an object element and  its  corresponding  process  and
       memory location.

       Let  A  be  a generic term for any 2D block cyclicly distributed array.
       Such a global array has an associated description vector DESCA.  In the
       following  comments,  the  character _ should be read as "of the global
       array".

       NOTATION        STORED IN      EXPLANATION
       ---------------  --------------  --------------------------------------
       DTYPE_A(global) DESCA( DTYPE_ )The descriptor type.  In this case,
                                      DTYPE_A = 1.
       CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
                                      the BLACS process grid A is distribu-
                                      ted over. The context itself is glo-
                                      bal, but the handle (the integer
                                      value) may vary.
       M_A    (global) DESCA( M_ )    The number of rows in the global
                                      array A.
       N_A    (global) DESCA( N_ )    The number of columns in the global
                                      array A.
       MB_A   (global) DESCA( MB_ )   The blocking factor used to distribute
                                      the rows of the array.
       NB_A   (global) DESCA( NB_ )   The blocking factor used to distribute
                                      the columns of the array.
       RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
                                      row  of  the  array  A  is  distributed.
       CSRC_A (global) DESCA( CSRC_ ) The process column over which the
                                      first column of the array A is
                                      distributed.
       LLD_A  (local)  DESCA( LLD_ )  The leading dimension of the local
                                      array.  LLD_A >= MAX(1,LOCr(M_A)).

       Let K be the number of rows or columns of  a  distributed  matrix,  and
       assume that its process grid has dimension p x q.
       LOCr(  K  )  denotes  the  number of elements of K that a process would
       receive if K were distributed over  the  p  processes  of  its  process
       column.
       Similarly, LOCc( K ) denotes the number of elements of K that a process
       would receive if K were distributed over the q processes of its process
       row.
       The  values  of  LOCr()  and LOCc() may be determined via a call to the
       ScaLAPACK tool function, NUMROC:
               LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
               LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).  An  upper
       bound for these quantities may be computed by:
               LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
               LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A

DESCRIPTION

       In  the  following  description,  A  denotes  A(IA:IA+N-1,JA:JA+N-1), B
       denotes B(IB:IB+N-1,JB:JB+NRHS-1) and X denotes
       X(IX:IX+N-1,JX:JX+NRHS-1).

       The following steps are performed:

       1. If FACT = ’E’, real scaling factors are computed to equilibrate
          the system:
             TRANS = ’N’:  diag(R)*A*diag(C)     *inv(diag(C))*X = diag(R)*B
             TRANS = ’T’: (diag(R)*A*diag(C))**T *inv(diag(R))*X = diag(C)*B
             TRANS = ’C’: (diag(R)*A*diag(C))**H *inv(diag(R))*X = diag(C)*B
          Whether or not the system will be equilibrated depends on the
          scaling of the matrix A, but if equilibration is used, A is
          overwritten by diag(R)*A*diag(C) and B by diag(R)*B (if TRANS=’N’)
          or diag(C)*B (if TRANS = ’T’ or ’C’).

       2. If FACT = ’N’ or ’E’, the LU decomposition is used to factor the
          matrix A (after equilibration if FACT = ’E’) as
             A = P * L * U,
          where P is a permutation matrix, L is a unit lower triangular
          matrix, and U is upper triangular.

       3. The factored form of A is used to estimate the condition number
          of the matrix A.  If the reciprocal of the condition number is
          less than machine precision, steps 4-6 are skipped.

       4. The system of equations is solved for X using the factored form
          of A.

       5. Iterative refinement is applied to improve the computed solution
          matrix and calculate error bounds and backward error estimates
          for it.

       6. If FACT = ’E’ and equilibration was used, the matrix X is
          premultiplied by diag(C) (if TRANS = ’N’) or diag(R) (if
          TRANS = ’T’ or ’C’) so that it solves the original system
          before equilibration.

ARGUMENTS

       FACT    (global input) CHARACTER
               Specifies whether or  not  the  factored  form  of  the  matrix
               A(IA:IA+N-1,JA:JA+N-1) is supplied on entry, and if not,
               whether    the    matrix   A(IA:IA+N-1,JA:JA+N-1)   should   be
               equilibrated  before  it  is  factored.   =  ’F’:   On   entry,
               AF(IAF:IAF+N-1,JAF:JAF+N-1) and IPIV con-
               tain  the factored form of A(IA:IA+N-1,JA:JA+N-1).  If EQUED is
               not   ’N’,   the   matrix   A(IA:IA+N-1,JA:JA+N-1)   has   been
               equilibrated   with   scaling   factors   given  by  R  and  C.
               A(IA:IA+N-1,JA:JA+N-1), AF(IAF:IAF+N-1,JAF:JAF+N-1),  and  IPIV
               are  not  modified.   = ’N’:  The matrix A(IA:IA+N-1,JA:JA+N-1)
               will be copied to
               AF(IAF:IAF+N-1,JAF:JAF+N-1) and factored.
               = ’E’:   The  matrix  A(IA:IA+N-1,JA:JA+N-1)  will  be  equili-
               brated if necessary, then copied to AF(IAF:IAF+N-1,JAF:JAF+N-1)
               and factored.

       TRANS   (global input) CHARACTER
               Specifies the form of the system of equations:
               = ’N’:  A(IA:IA+N-1,JA:JA+N-1) * X(IX:IX+N-1,JX:JX+NRHS-1)
               = B(IB:IB+N-1,JB:JB+NRHS-1)     (No transpose)
               = ’T’:  A(IA:IA+N-1,JA:JA+N-1)**T * X(IX:IX+N-1,JX:JX+NRHS-1)
               = B(IB:IB+N-1,JB:JB+NRHS-1)  (Transpose)
               = ’C’:  A(IA:IA+N-1,JA:JA+N-1)**H * X(IX:IX+N-1,JX:JX+NRHS-1)
               = B(IB:IB+N-1,JB:JB+NRHS-1)  (Conjugate transpose)

       N       (global input) INTEGER
               The number of rows and columns to  be  operated  on,  i.e.  the
               order  of  the distributed submatrix A(IA:IA+N-1,JA:JA+N-1).  N
               >= 0.

       NRHS    (global input) INTEGER
               The number of right-hand sides, i.e., the number of columns  of
               the distributed submatrices B(IB:IB+N-1,JB:JB+NRHS-1) and
               X(IX:IX+N-1,JX:JX+NRHS-1).  NRHS >= 0.

       A       (local input/local output) COMPLEX*16 pointer into
               the   local   memory   to   an   array   of   local   dimension
               (LLD_A,LOCc(JA+N-1)).    On   entry,    the    N-by-N    matrix
               A(IA:IA+N-1,JA:JA+N-1).  If FACT = ’F’ and EQUED is not ’N’,
               then A(IA:IA+N-1,JA:JA+N-1) must have been equilibrated by
               the  scaling  factors in R and/or C.  A(IA:IA+N-1,JA:JA+N-1) is
               not modified if FACT = ’F’ or  ’N’, or if FACT = ’E’ and  EQUED
               = ’N’ on exit.

               On exit, if EQUED .ne. ’N’, A(IA:IA+N-1,JA:JA+N-1) is scaled as
               follows:
               EQUED = ’R’:  A(IA:IA+N-1,JA:JA+N-1) :=
               diag(R) * A(IA:IA+N-1,JA:JA+N-1)
               EQUED = ’C’:  A(IA:IA+N-1,JA:JA+N-1) :=
               A(IA:IA+N-1,JA:JA+N-1) * diag(C)
               EQUED = ’B’:  A(IA:IA+N-1,JA:JA+N-1) :=
               diag(R) * A(IA:IA+N-1,JA:JA+N-1) * diag(C).

       IA      (global input) INTEGER
               The row index in the global array A indicating the first row of
               sub( A ).

       JA      (global input) INTEGER
               The  column  index  in  the global array A indicating the first
               column of sub( A ).

       DESCA   (global and local input) INTEGER array of dimension DLEN_.
               The array descriptor for the distributed matrix A.

       AF      (local input or local output) COMPLEX*16 pointer
               into  the  local  memory  to  an  array  of   local   dimension
               (LLD_AF,LOCc(JA+N-1)).      If     FACT     =     ’F’,     then
               AF(IAF:IAF+N-1,JAF:JAF+N-1) is an input argument and  on  entry
               contains   the   factors   L   and  U  from  the  factorization
               A(IA:IA+N-1,JA:JA+N-1) = P*L*U  as  computed  by  PZGETRF.   If
               EQUED   .ne.   ’N’,  then  AF  is  the  factored  form  of  the
               equilibrated matrix A(IA:IA+N-1,JA:JA+N-1).

               If FACT = ’N’, then AF(IAF:IAF+N-1,JAF:JAF+N-1)  is  an  output
               argument  and  on  exit  returns  the  factors L and U from the
               factorization A(IA:IA+N-1,JA:JA+N-1) = P*L*U of the original
               matrix A(IA:IA+N-1,JA:JA+N-1).

               If FACT = ’E’, then AF(IAF:IAF+N-1,JAF:JAF+N-1)  is  an  output
               argument  and  on  exit  returns  the  factors L and U from the
               factorization A(IA:IA+N-1,JA:JA+N-1) = P*L*U of the equili-
               brated matrix A(IA:IA+N-1,JA:JA+N-1) (see the description of
               A(IA:IA+N-1,JA:JA+N-1)  for  the  form  of   the   equilibrated
               matrix).

       IAF     (global input) INTEGER
               The  row  index in the global array AF indicating the first row
               of sub( AF ).

       JAF     (global input) INTEGER
               The column index in the global array AF  indicating  the  first
               column of sub( AF ).

       DESCAF  (global and local input) INTEGER array of dimension DLEN_.
               The array descriptor for the distributed matrix AF.

       IPIV    (local input or local output) INTEGER array, dimension
               LOCr(M_A)+MB_A. If FACT = ’F’, then IPIV is an input argu- ment
               and  on  entry  contains  the  pivot  indices  from  the   fac-
               torization   A(IA:IA+N-1,JA:JA+N-1)  =  P*L*U  as  computed  by
               PZGETRF; IPIV(i) -> The global row  local  row  i  was  swapped
               with.  This array must be aligned with A( IA:IA+N-1, * ).

               If  FACT  =  ’N’,  then  IPIV is an output argument and on exit
               contains   the   pivot   indices   from    the    factorization
               A(IA:IA+N-1,JA:JA+N-1) = P*L*U of the original matrix
               A(IA:IA+N-1,JA:JA+N-1).

               If  FACT  =  ’E’,  then  IPIV is an output argument and on exit
               contains   the   pivot   indices   from    the    factorization
               A(IA:IA+N-1,JA:JA+N-1) = P*L*U of the equilibrated matrix
               A(IA:IA+N-1,JA:JA+N-1).

       EQUED   (global input or global output) CHARACTER
               Specifies  the form of equilibration that was done.  = ’N’:  No
               equilibration (always true if FACT = ’N’).
               = ’R’:  Row  equilibration,  i.e.,  A(IA:IA+N-1,JA:JA+N-1)  has
               been  premultiplied  by diag(R).  = ’C’:  Column equilibration,
               i.e.,  A(IA:IA+N-1,JA:JA+N-1)  has   been   postmultiplied   by
               diag(C).  = ’B’:  Both row and column equilibration, i.e.,
               A(IA:IA+N-1,JA:JA+N-1) has been replaced by
               diag(R)  * A(IA:IA+N-1,JA:JA+N-1) * diag(C).  EQUED is an input
               variable if FACT = ’F’; otherwise, it is an output variable.

       R       (local input or local output) DOUBLE PRECISION array,
               dimension   LOCr(M_A).    The    row    scale    factors    for
               A(IA:IA+N-1,JA:JA+N-1).
               If  EQUED = ’R’ or ’B’, A(IA:IA+N-1,JA:JA+N-1) is multiplied on
               the left by diag(R); if EQUED=’N’ or ’C’, R is not acces-  sed.
               R is an input variable if FACT = ’F’; otherwise, R is an output
               variable.  If FACT = ’F’ and EQUED = ’R’ or ’B’,  each  element
               of  R  must  be  positive.   R  is  replicated in every process
               column, and is aligned with the distributed matrix A.

       C       (local input or local output) DOUBLE PRECISION array,
               dimension   LOCc(N_A).    The   column   scale   factors    for
               A(IA:IA+N-1,JA:JA+N-1).
               If  EQUED = ’C’ or ’B’, A(IA:IA+N-1,JA:JA+N-1) is multiplied on
               the right by diag(C); if EQUED = ’N’ or ’R’, C is not accessed.
               C is an input variable if FACT = ’F’; otherwise, C is an output
               variable.  If FACT = ’F’ and EQUED = ’C’ or C is replicated  in
               every  process  row, and is aligned with the distributed matrix
               A.

       B       (local input/local output) COMPLEX*16 pointer
               into  the  local  memory  to  an  array  of   local   dimension
               (LLD_B,LOCc(JB+NRHS-1)  ).   On entry, the N-by-NRHS right-hand
               side matrix B(IB:IB+N-1,JB:JB+NRHS-1). On exit, if
               EQUED = ’N’,  B(IB:IB+N-1,JB:JB+NRHS-1)  is  not  modified;  if
               TRANS  =  ’N’  and  EQUED  =  ’R’  or  ’B’, B is overwritten by
               diag(R)*B(IB:IB+N-1,JB:JB+NRHS-1); if TRANS = ’T’ or ’C’
               and EQUED = ’C’ or ’B’, B(IB:IB+N-1,JB:JB+NRHS-1) is over-
               written by diag(C)*B(IB:IB+N-1,JB:JB+NRHS-1).

       IB      (global input) INTEGER
               The row index in the global array B indicating the first row of
               sub( B ).

       JB      (global input) INTEGER
               The  column  index  in  the global array B indicating the first
               column of sub( B ).

       DESCB   (global and local input) INTEGER array of dimension DLEN_.
               The array descriptor for the distributed matrix B.

       X       (local input/local output) COMPLEX*16 pointer
               into the local memory to an array of  local  dimension  (LLD_X,
               LOCc(JX+NRHS-1)).   If  INFO = 0, the N-by-NRHS solution matrix
               X(IX:IX+N-1,JX:JX+NRHS-1) to the original
               system of equations.  Note that A(IA:IA+N-1,JA:JA+N-1) and
               B(IB:IB+N-1,JB:JB+NRHS-1) are modified on exit  if  EQUED  .ne.
               ’N’,   and   the   solution   to  the  equilibrated  system  is
               inv(diag(C))*X(IX:IX+N-1,JX:JX+NRHS-1) if TRANS = ’N’ and EQUED
               =  ’C’  or  ’B’,  or  inv(diag(R))*X(IX:IX+N-1,JX:JX+NRHS-1) if
               TRANS = ’T’ or ’C’ and EQUED = ’R’ or ’B’.

       IX      (global input) INTEGER
               The row index in the global array X indicating the first row of
               sub( X ).

       JX      (global input) INTEGER
               The  column  index  in  the global array X indicating the first
               column of sub( X ).

       DESCX   (global and local input) INTEGER array of dimension DLEN_.
               The array descriptor for the distributed matrix X.

       RCOND   (global output) DOUBLE PRECISION
               The estimate of the reciprocal condition number of  the  matrix
               A(IA:IA+N-1,JA:JA+N-1) after equilibration (if done).  If RCOND
               is less than the machine precision (in particular, if  RCOND  =
               0),   the  matrix  is  singular  to  working  precision.   This
               condition is indicated by a return code of INFO > 0.

       FERR    (local output) DOUBLE PRECISION array, dimension LOCc(N_B)
               The estimated forward error bounds  for  each  solution  vector
               X(j)    (the    j-th    column    of    the   solution   matrix
               X(IX:IX+N-1,JX:JX+NRHS-1).  If  XTRUE  is  the  true  solution,
               FERR(j)  bounds  the  magnitude of the largest entry in (X(j) -
               XTRUE) divided by the magnitude of the largest entry  in  X(j).
               The  estimate  is as reliable as the estimate for RCOND, and is
               almost always a slight overestimate of the true error.  FERR is
               replicated  in  every  process  row,  and  is  aligned with the
               matrices B and X.

       BERR    (local output) DOUBLE PRECISION array, dimension LOCc(N_B).
               The componentwise relative  backward  error  of  each  solution
               vector X(j) (i.e., the smallest relative change in any entry of
               A(IA:IA+N-1,JA:JA+N-1) or
               B(IB:IB+N-1,JB:JB+NRHS-1) that makes X(j) an  exact  solution).
               BERR  is  replicated  in every process row, and is aligned with
               the matrices B and X.

       WORK    (local workspace/local output) COMPLEX*16 array,
               dimension (LWORK) On exit,  WORK(1)  returns  the  minimal  and
               optimal LWORK.

       LWORK   (local or global input) INTEGER
               The dimension of the array WORK.  LWORK is local input and must
               be at least LWORK = MAX( PZGECON( LWORK ), PZGERFS( LWORK ) ) +
               LOCr( N_A ).

               If LWORK = -1, then LWORK is global input and a workspace query
               is assumed; the routine only calculates the minimum and optimal
               size  for  all work arrays. Each of these values is returned in
               the first entry of the corresponding work array, and  no  error
               message is issued by PXERBLA.

       RWORK   (local workspace/local output) DOUBLE PRECISION array,
               dimension  (LRWORK)  On  exit, RWORK(1) returns the minimal and
               optimal LRWORK.

       LRWORK  (local or global input) INTEGER
               The dimension of the array RWORK.  LRWORK is  local  input  and
               must be at least LRWORK = 2*LOCc(N_A).

               If  LRWORK  =  -1,  then LRWORK is global input and a workspace
               query is assumed; the routine only calculates the  minimum  and
               optimal  size  for  all  work  arrays.  Each of these values is
               returned in the first entry of the  corresponding  work  array,
               and no error message is issued by PXERBLA.

       INFO    (global output) INTEGER
               = 0:  successful exit
               < 0:  if INFO = -i, the i-th argument had an illegal value
               > 0:  if INFO = i, and i is
               <= N:  U(IA+I-1,IA+I-1) is exactly zero.  The factorization has
               been completed, but the factor U is exactly  singular,  so  the
               solution  and error bounds could not be computed.  = N+1: RCOND
               is less than machine precision.   The  factorization  has  been
               completed, but the matrix is singular to working precision, and
               the solution and error bounds have not been computed.