Routine: PZGECON()  File: SRC\pzgecon.f

 
 
# lines: 422
  # code: 422
  # comment: 0
  # blank:0
# Variables:76
# Callers:1
# Callings:4
# Words:242
# Keywords:147
 

 

..
     .. Array Arguments ..
     ..
  Purpose
  =======
  PZGECON estimates the reciprocal of the condition number of a general
  distributed complex matrix A(IA:IA+N-1,JA:JA+N-1), in either the
  1-norm or the infinity-norm, using the LU factorization computed by
  PZGETRF.
  An estimate is obtained for norm(inv(A(IA:IA+N-1,JA:JA+N-1))), and
  the reciprocal of the condition number is computed as
             RCOND = 1 / ( norm( A(IA:IA+N-1,JA:JA+N-1)      ) *
                           norm( inv(A(IA:IA+N-1,JA:JA+N-1)) ) ).
  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
  Arguments
  =========
  NORM    (global input) CHARACTER
          Specifies whether the 1-norm condition number or the
          infinity-norm condition number is required:
          = '1' or 'O':  1-norm
          = 'I':         Infinity-norm
  N       (global input) INTEGER
          The order of the distributed matrix A(IA:IA+N-1,JA:JA+N-1).
          N >= 0.
  A       (local input) COMPLEX*16 pointer into the local memory
          to an array of dimension ( LLD_A, LOCc(JA+N-1) ). On entry,
          this array contains the local pieces of the factors L and U
          from the factorization A(IA:IA+N-1,JA:JA+N-1) = P*L*U; the
          unit diagonal elements of L are not stored.
  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.
  ANORM   (global input) DOUBLE PRECISION
          If NORM = '1' or 'O', the 1-norm of the original distributed
          matrix A(IA:IA+N-1,JA:JA+N-1).
          If NORM = 'I', the infinity-norm of the original distributed
          matrix A(IA:IA+N-1,JA:JA+N-1).
  RCOND   (global output) DOUBLE PRECISION
          The reciprocal of the condition number of the distributed
          matrix A(IA:IA+N-1,JA:JA+N-1), computed as
             RCOND = 1 / ( norm( A(IA:IA+N-1,JA:JA+N-1)      ) *
                           norm( inv(A(IA:IA+N-1,JA:JA+N-1)) ) ).
  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 >= 2*LOCr(N+MOD(IA-1,MB_A)) +
          MAX( 2, MAX(NB_A*CEIL(NPROW-1,NPCOL),LOCc(N+MOD(JA-1,NB_A)) +
          NB_A*CEIL(NPCOL-1,NPROW)) ).
          LOCr and LOCc values can be computed using the ScaLAPACK
          tool function NUMROC; NPROW and NPCOL can be determined by
          calling the subroutine BLACS_GRIDINFO.
          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+MOD(JA-1,NB_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 the i-th argument is an array and the j-entry had
                an illegal value, then INFO = -(i*100+j), if the i-th
                argument is a scalar and had an illegal value, then
                INFO = -i.
  =====================================================================
     .. Parameters ..

 
Display dynamic version Find AutoScroll Reload FontSize: - + Hide Comments Hide Blanks Frame FullScreen MailPrint

 
001        SUBROUTINE PZGECON( NORM , N , A , IA , JA , DESCA , ANORM , RCOND , WORK ,
002       $LWORK , RWORK , LRWORK , INFO )
003  
004  *     -- ScaLAPACK routine(version 1.7) --
005  *     University of Tennessee , Knoxville , Oak Ridge National Laboratory ,
006  *     and University of California , Berkeley.
007  *     May 25 , 2001
008  
009  *     .. Scalar Arguments ..
010        CHARACTER NORM
011        INTEGER IA , INFO , JA , LRWORK , LWORK , N
012        DOUBLE PRECISION ANORM , RCOND
013        INTEGER BLOCK_CYCLIC_2D , CSRC_ , CTXT_ , DLEN_ , DTYPE_ ,
014       $LLD_ , MB_ , M_ , NB_ , N_ , RSRC_
015        PARAMETER( BLOCK_CYCLIC_2D = 1 , DLEN_ = 9 , DTYPE_ = 1 ,
016       $CTXT_ = 2 , M_ = 3 , N_ = 4 , MB_ = 5 , NB_ = 6 ,
017       $RSRC_ = 7 , CSRC_ = 8 , LLD_ = 9 )
018        DOUBLE PRECISION ONE , ZERO
019        PARAMETER( ONE = 1.0D + 0 , ZERO = 0.0D + 0 )
020  *     ..
021  *     .. Local Scalars ..
022        LOGICAL LQUERY , ONENRM
023        CHARACTER CBTOP , COLCTOP , NORMIN , ROWCTOP
024        INTEGER IACOL , IAROW , ICOFF , ICTXT , IIA , IPNL , IPNU ,
025       $IPV , IPW , IPX , IROFF , IV , IX , IXX , JJA , JV , JX ,
026       $KASE , KASE1 , LRWMIN , LWMIN , MYCOL , MYROW , NP ,
027       $NPCOL , NPMOD , NPROW , NQ , NQMOD
028        DOUBLE PRECISION AINVNM , SCALE , SL , SMLNUM , SU
029        COMPLEX*16 WMAX , ZDUM
030  *     ..
031  *     .. Local Arrays ..
032        INTEGER DESCV( DLEN_ ) , DESCX( DLEN_ ) , IDUM1( 3 ) ,
033       $IDUM2( 3 )
034  *     ..
035  *     .. External Subroutines ..
036        EXTERNAL BLACS_GRIDINFO , CHK1MAT , DESCSET , INFOG2L ,
037       $PCHK1MAT , PB_TOPGET , PB_TOPSET , PXERBLA ,
038       $PZAMAX , PZLATRS , PZLACON , PZDRSCL ,
039       $ZGEBR2D , ZGEBS2D
040  *     ..
041  *     .. External Functions ..
042        LOGICAL LSAME
043        INTEGER ICEIL , INDXG2P , NUMROC
044        DOUBLE PRECISION PDLAMCH
045        EXTERNAL ICEIL , INDXG2P , LSAME , NUMROC , PDLAMCH
046  *     ..
047  *     .. Intrinsic Functions ..
048        INTRINSIC ABS , DBLE , DIMAG , ICHAR , MAX , MOD
049  *     ..
050  *     .. Statement Functions ..
051        DOUBLE PRECISION CABS1
052  *     ..
053  *     .. Statement Function definitions ..
054        CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
055  *     ..
056  *     .. Executable Statements ..
057  
058  *     Get grid parameters
059  
060        ICTXT = DESCA( CTXT_ )
061        CALL BLACS_GRIDINFO( ICTXT , NPROW , NPCOL , MYROW , MYCOL )
062  
063  *     Test the input parameters
064  
065        INFO = 0
066        IF( NPROW.EQ. - 1 ) THEN
067            INFO = - ( 600 + CTXT_ )
068        ELSE
069            CALL CHK1MAT( N , 2 , N , 2 , IA , JA , DESCA , 6 , INFO )
070            IF( INFO.EQ.0 ) THEN
071                ONENRM = NORM.EQ.'1' .OR. LSAME( NORM , 'O' )
072                IAROW = INDXG2P( IA , DESCA( MB_ ) , MYROW , DESCA( RSRC_ ) ,
073       $        NPROW )
074                IACOL = INDXG2P( JA , DESCA( NB_ ) , MYCOL , DESCA( CSRC_ ) ,
075       $        NPCOL )
076                NPMOD = NUMROC( N + MOD( IA - 1 , DESCA( MB_ ) ) , DESCA( MB_ ) ,
077       $        MYROW , IAROW , NPROW )
078                NQMOD = NUMROC( N + MOD( JA - 1 , DESCA( NB_ ) ) , DESCA( NB_ ) ,
079       $        MYCOL , IACOL , NPCOL )
080                LWMIN = 2*NPMOD +
081       $        MAX( 2 , MAX( DESCA( NB_ )*
082       $        MAX( 1 , ICEIL( NPROW - 1 , NPCOL ) ) , NQMOD +
083       $        DESCA( NB_ )*
084       $        MAX( 1 , ICEIL( NPCOL - 1 , NPROW ) ) ) )
085                WORK( 1 ) = DBLE( LWMIN )
086                LRWMIN = 2*NQMOD
087                RWORK( 1 ) = DBLE( LRWMIN )
088                LQUERY =( LWORK.EQ. - 1 .OR. LRWORK.EQ. - 1 )
089  
090                IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM , 'I' ) ) THEN
091                    INFO = - 1
092                ELSE IF( ANORM.LT.ZERO ) THEN
093                    INFO = - 7
094                ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
095                    INFO = - 10
096                ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN
097                    INFO = - 12
098                END IF
099            END IF
100  
101            IF( ONENRM ) THEN
102                IDUM1( 1 ) = ICHAR( '1' )
103            ELSE
104                IDUM1( 1 ) = ICHAR( 'I' )
105            END IF
106            IDUM2( 1 ) = 1
107            IF( LWORK.EQ. - 1 ) THEN
108                IDUM1( 2 ) = - 1
109            ELSE
110                IDUM1( 2 ) = 1
111            END IF
112            IDUM2( 2 ) = 10
113            IF( LRWORK.EQ. - 1 ) THEN
114                IDUM1( 3 ) = - 1
115            ELSE
116                IDUM1( 3 ) = 1
117            END IF
118            IDUM2( 3 ) = 12
119            CALL PCHK1MAT( N , 2 , N , 2 , IA , JA , DESCA , 6 , 3 , IDUM1 , IDUM2 ,
120       $    INFO )
121        END IF
122  
123        IF( INFO.NE.0 ) THEN
124            CALL PXERBLA( ICTXT , 'PZGECON' , - INFO )
125            RETURN
126        ELSE IF( LQUERY ) THEN
127            RETURN
128        END IF
129  
130  *     Quick return if possible
131  
132        RCOND = ZERO
133        IF( N.EQ.0 ) THEN
134            RCOND = ONE
135            RETURN
136        ELSE IF( ANORM.EQ.ZERO ) THEN
137            RETURN
138        ELSE IF( N.EQ.1 ) THEN
139            RCOND = ONE
140            RETURN
141        END IF
142  
143        CALL PB_TOPGET( ICTXT , 'Combine' , 'Columnwise' , COLCTOP )
144        CALL PB_TOPGET( ICTXT , 'Combine' , 'Rowwise' , ROWCTOP )
145        CALL PB_TOPSET( ICTXT , 'Combine' , 'Columnwise' , '1 - tree' )
146        CALL PB_TOPSET( ICTXT , 'Combine' , 'Rowwise' , '1 - tree' )
147  
148        SMLNUM = PDLAMCH( ICTXT , 'Safe minimum' )
149        IROFF = MOD( IA - 1 , DESCA( MB_ ) )
150        ICOFF = MOD( JA - 1 , DESCA( NB_ ) )
151        CALL INFOG2L( IA , JA , DESCA , NPROW , NPCOL , MYROW , MYCOL , IIA , JJA ,
152       $IAROW , IACOL )
153        NP = NUMROC( N + IROFF , DESCA( MB_ ) , MYROW , IAROW , NPROW )
154        NQ = NUMROC( N + ICOFF , DESCA( NB_ ) , MYCOL , IACOL , NPCOL )
155        IV = IROFF + 1
156        IX = IV
157        JV = ICOFF + 1
158        JX = JV
159  
160        IPX = 1
161        IPV = IPX + NP
162        IPW = IPV + NP
163        IPNL = 1
164        IPNU = IPNL + NQ
165  
166        CALL DESCSET( DESCV , N + IROFF , 1 , DESCA( MB_ ) , 1 , IAROW , MYCOL ,
167       $ICTXT , MAX( 1 , NP ) )
168        CALL DESCSET( DESCX , N + IROFF , 1 , DESCA( MB_ ) , 1 , IAROW , MYCOL ,
169       $ICTXT , MAX( 1 , NP ) )
170  
171  *     Estimate the norm of inv(A).
172  
173        AINVNM = ZERO
174        NORMIN = 'N'
175        IF( ONENRM ) THEN
176            KASE1 = 1
177        ELSE
178            KASE1 = 2
179        END IF
180        KASE = 0
181  
182     10 CONTINUE
183        CALL PZLACON ( N , WORK( IPV ) , IV , JV , DESCV , WORK( IPX ) , IX , JX ,
184       $DESCX , AINVNM , KASE )
185        IF( KASE.NE.0 ) THEN
186            IF( KASE.EQ.KASE1 ) THEN
187  
188  *             Multiply by inv(L).
189  
190                DESCX( CSRC_ ) = IACOL
191                CALL PZLATRS ( 'Lower' , 'No transpose' , 'Unit' , NORMIN ,
192       $        N , A , IA , JA , DESCA , WORK( IPX ) , IX , JX ,
193       $        DESCX , SL , RWORK( IPNL ) , WORK( IPW ) )
194                DESCX( CSRC_ ) = MYCOL
195  
196  *             Multiply by inv(U).
197  
198                DESCX( CSRC_ ) = IACOL
199                CALL PZLATRS ( 'Upper' , 'No transpose' , 'Non - unit' , NORMIN ,
200       $        N , A , IA , JA , DESCA , WORK( IPX ) , IX , JX ,
201       $        DESCX , SU , RWORK( IPNU ) , WORK( IPW ) )
202                DESCX( CSRC_ ) = MYCOL
203            ELSE
204  
205  *             Multiply by inv(U').
206  
207                DESCX( CSRC_ ) = IACOL
208                CALL PZLATRS ( 'Upper' , 'Conjugate transpose' , 'Non - unit' ,
209       $        NORMIN , N , A , IA , JA , DESCA , WORK( IPX ) , IX ,
210       $        JX , DESCX , SU , RWORK( IPNU ) , WORK( IPW ) )
211                DESCX( CSRC_ ) = MYCOL
212  
213  *             Multiply by inv(L').
214  
215                DESCX( CSRC_ ) = IACOL
216                CALL PZLATRS ( 'Lower' , 'Conjugate transpose' , 'Unit' ,
217       $        NORMIN , N , A , IA , JA , DESCA , WORK( IPX ) ,
218       $        IX , JX , DESCX , SL , RWORK( IPNL ) ,
219       $        WORK( IPW ) )
220                DESCX( CSRC_ ) = MYCOL
221            END IF
222  
223  *         Divide X by 1 / (SL*SU) if doing so will not cause overflow.
224  
225            SCALE = SL*SU
226            NORMIN = 'Y'
227            IF( SCALE.NE.ONE ) THEN
228                CALL PZAMAX( N , WMAX , IXX , WORK( IPX ) , IX , JX , DESCX , 1 )
229                IF( DESCX( M_ ).EQ.1 .AND. N.EQ.1 ) THEN
230                    CALL PB_TOPGET( ICTXT , 'Broadcast' , 'Columnwise' , CBTOP )
231                    IF( MYROW.EQ.IAROW ) THEN
232                        CALL ZGEBS2D( ICTXT , 'Column' , CBTOP , 1 , 1 , WMAX , 1 )
233                    ELSE
234                        CALL ZGEBR2D( ICTXT , 'Column' , CBTOP , 1 , 1 , WMAX , 1 ,
235       $                IAROW , MYCOL )
236                    END IF
237                END IF
238                IF( SCALE.LT.CABS1( WMAX )*SMLNUM .OR. SCALE.EQ.ZERO )
239       $            GO TO 20
240                    CALL PZDRSCL ( N , SCALE , WORK( IPX ) , IX , JX , DESCX , 1 )
241                END IF
242                GO TO 10
243            END IF
244  
245  *         Compute the estimate of the reciprocal condition number.
246  
247            IF( AINVNM.NE.ZERO )
248       $        RCOND =( ONE / AINVNM ) / ANORM
249  
250     20 CONTINUE
251  
252        CALL PB_TOPSET( ICTXT , 'Combine' , 'Columnwise' , COLCTOP )
253        CALL PB_TOPSET( ICTXT , 'Combine' , 'Rowwise' , ROWCTOP )
254  
255        RETURN
256  
257  *     End of PZGECON
258  
259        END