Routine: PZTRCON()  File: SRC\pztrcon.f

 
 
# lines: 431
  # code: 431
  # comment: 0
  # blank:0
# Variables:77
# Callers:0
# Callings:5
# Words:257
# Keywords:156
 

 

..
     .. Array Arguments ..
     ..
  Purpose
  =======
  PZTRCON estimates the reciprocal of the condition number of a
  triangular distributed matrix A(IA:IA+N-1,JA:JA+N-1), in either the
  1-norm or the infinity-norm.
  The norm of A(IA:IA+N-1,JA:JA+N-1) is computed and an estimate is
  obtained for norm(inv(A(IA:IA+N-1,JA:JA+N-1))), then 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.
  UPLO    (global input) CHARACTER
          = 'U':  A(IA:IA+N-1,JA:JA+N-1) is upper triangular;
          = 'L':  A(IA:IA+N-1,JA:JA+N-1) is lower triangular.
  DIAG    (global input) CHARACTER
          = 'N':  A(IA:IA+N-1,JA:JA+N-1) is non-unit triangular;
          = 'U':  A(IA:IA+N-1,JA:JA+N-1) is unit triangular.
  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) ). This array
          contains the local pieces of the triangular distributed
          matrix A(IA:IA+N-1,JA:JA+N-1). If UPLO = 'U', the leading
          N-by-N upper triangular part of this distributed matrix con-
          tains the upper triangular matrix, and its strictly lower
          triangular part is not referenced.  If UPLO = 'L', the
          leading N-by-N lower triangular part of this ditributed
          matrix contains the lower triangular matrix, and the strictly
          upper triangular part is not referenced. If DIAG = 'U', the
          diagonal elements of A(IA:IA+N-1,JA:JA+N-1) are also not
          referenced and are assumed to be 1.
  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.
  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(P-1,Q),LOCc(N+MOD(JA-1,NB_A)) +
          NB_A*CEIL(Q-1,P)) ).
          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 >= 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 PZTRCON( NORM , UPLO , DIAG , N , A , IA , JA , DESCA , RCOND ,
002       $WORK , 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 DIAG , NORM , UPLO
011        INTEGER IA , JA , INFO , LRWORK , LWORK , N
012        DOUBLE PRECISION 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 , NOUNIT , ONENRM , UPPER
023        CHARACTER CBTOP , COLCTOP , NORMIN , ROWCTOP
024        INTEGER IACOL , IAROW , ICOFF , ICTXT , IIA , IPN , IPV , IPW ,
025       $IPX , IROFF , IV , IX , IXX , JJA , JV , JX , KASE ,
026       $KASE1 , LRWMIN , LWMIN , MYCOL , MYROW , NP , NPCOL ,
027       $NPMOD , NPROW , NQMOD
028        DOUBLE PRECISION AINVNM , ANORM , SCALE , SMLNUM
029        COMPLEX*16 WMAX , ZDUM
030  *     ..
031  *     .. Local Arrays ..
032        INTEGER DESCV( DLEN_ ) , DESCX( DLEN_ ) , IDUM1( 5 ) ,
033       $IDUM2( 5 )
034  *     ..
035  *     .. External Subroutines ..
036        EXTERNAL BLACS_GRIDINFO , CHK1MAT , DESCSET , INFOG2L ,
037       $PB_TOPGET , PB_TOPSET , PXERBLA , PCHK1MAT ,
038       $PZAMAX , PZLATRS , PZLACON , PZDRSCL ,
039       $ZGEBR2D , ZGEBS2D
040  *     ..
041  *     .. External Functions ..
042        LOGICAL LSAME
043        INTEGER ICEIL , INDXG2P , NUMROC
044        DOUBLE PRECISION PDLAMCH , PZLANTR
045        EXTERNAL ICEIL , INDXG2P , LSAME , NUMROC , PDLAMCH ,
046       $PZLANTR
047  *     ..
048  *     .. Intrinsic Functions ..
049        INTRINSIC ABS , DBLE , DIMAG , ICHAR , MAX , MOD
050  *     ..
051  *     .. Statement Functions ..
052        DOUBLE PRECISION CABS1
053  *     ..
054  *     .. Statement Function definitions ..
055        CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
056  *     ..
057  *     .. Executable Statements ..
058  
059  *     Get grid parameters
060  
061        ICTXT = DESCA( CTXT_ )
062        CALL BLACS_GRIDINFO( ICTXT , NPROW , NPCOL , MYROW , MYCOL )
063  
064  *     Test the input parameters
065  
066        INFO = 0
067        IF( NPROW.EQ. - 1 ) THEN
068            INFO = - ( 800 + CTXT_ )
069        ELSE
070            CALL CHK1MAT( N , 4 , N , 4 , IA , JA , DESCA , 8 , INFO )
071            IF( INFO.EQ.0 ) THEN
072                UPPER = LSAME( UPLO , 'U' )
073                ONENRM = NORM.EQ.'1' .OR. LSAME( NORM , 'O' )
074                NOUNIT = LSAME( DIAG , 'N' )
075                IAROW = INDXG2P( IA , DESCA( MB_ ) , MYROW , DESCA( RSRC_ ) ,
076       $        NPROW )
077                IACOL = INDXG2P( JA , DESCA( NB_ ) , MYCOL , DESCA( CSRC_ ) ,
078       $        NPCOL )
079                NPMOD = NUMROC( N + MOD( IA - 1 , DESCA( MB_ ) ) , DESCA( MB_ ) ,
080       $        MYROW , IAROW , NPROW )
081                NQMOD = NUMROC( N + MOD( JA - 1 , DESCA( NB_ ) ) , DESCA( NB_ ) ,
082       $        MYCOL , IACOL , NPCOL )
083                LWMIN = 2*NPMOD +
084       $        MAX( 2 , MAX( DESCA( NB_ )*
085       $        MAX( 1 , ICEIL( NPROW - 1 , NPCOL ) ) , NQMOD +
086       $        DESCA( NB_ )*
087       $        MAX( 1 , ICEIL( NPCOL - 1 , NPROW ) ) ) )
088                WORK( 1 ) = DBLE( LWMIN )
089                LRWMIN = NQMOD
090                RWORK( 1 ) = DBLE( LRWMIN )
091                LQUERY =( LWORK.EQ. - 1 .OR. LRWORK.EQ. - 1 )
092  
093                IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM , 'I' ) ) THEN
094                    INFO = - 1
095                ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO , 'L' ) ) THEN
096                    INFO = - 2
097                ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG , 'U' ) ) THEN
098                    INFO = - 3
099                ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
100                    INFO = - 11
101                ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN
102                    INFO = - 13
103                END IF
104            END IF
105  
106            IF( ONENRM ) THEN
107                IDUM1( 1 ) = ICHAR( '1' )
108            ELSE
109                IDUM1( 1 ) = ICHAR( 'I' )
110            END IF
111            IDUM2( 1 ) = 1
112            IF( UPPER ) THEN
113                IDUM1( 2 ) = ICHAR( 'U' )
114            ELSE
115                IDUM1( 2 ) = ICHAR( 'L' )
116            END IF
117            IDUM2( 2 ) = 2
118            IF( NOUNIT ) THEN
119                IDUM1( 3 ) = ICHAR( 'N' )
120            ELSE
121                IDUM1( 3 ) = ICHAR( 'U' )
122            END IF
123            IDUM2( 3 ) = 3
124            IF( LWORK.EQ. - 1 ) THEN
125                IDUM1( 4 ) = - 1
126            ELSE
127                IDUM1( 4 ) = 1
128            END IF
129            IDUM2( 4 ) = 11
130            IF( LRWORK.EQ. - 1 ) THEN
131                IDUM1( 5 ) = - 1
132            ELSE
133                IDUM1( 5 ) = 1
134            END IF
135            IDUM2( 5 ) = 13
136            CALL PCHK1MAT( N , 4 , N , 4 , IA , JA , DESCA , 8 , 5 , IDUM1 , IDUM2 ,
137       $    INFO )
138        END IF
139  
140        IF( INFO.NE.0 ) THEN
141            CALL PXERBLA( ICTXT , 'PZTRCON' , - INFO )
142            RETURN
143        ELSE IF( LQUERY ) THEN
144            RETURN
145        END IF
146  
147  *     Quick return if possible
148  
149        IF( N.EQ.0 ) THEN
150            RCOND = ONE
151            RETURN
152        END IF
153  
154        CALL PB_TOPGET( ICTXT , 'Combine' , 'Columnwise' , COLCTOP )
155        CALL PB_TOPGET( ICTXT , 'Combine' , 'Rowwise' , ROWCTOP )
156        CALL PB_TOPSET( ICTXT , 'Combine' , 'Columnwise' , '1 - tree' )
157        CALL PB_TOPSET( ICTXT , 'Combine' , 'Rowwise' , '1 - tree' )
158  
159        RCOND = ZERO
160        SMLNUM = PDLAMCH( ICTXT , 'Safe minimum' )*DBLE( MAX( 1 , N ) )
161        CALL INFOG2L( IA , JA , DESCA , NPROW , NPCOL , MYROW , MYCOL , IIA , JJA ,
162       $IAROW , IACOL )
163        IROFF = MOD( IA - 1 , DESCA( MB_ ) )
164        ICOFF = MOD( JA - 1 , DESCA( NB_ ) )
165        NP = NUMROC( N + IROFF , DESCA( MB_ ) , MYROW , IAROW , NPROW )
166        IV = IROFF + 1
167        IX = IV
168        JV = ICOFF + 1
169        JX = JV
170  
171        IPX = 1
172        IPV = IPX + NP
173        IPW = IPV + NP
174        IPN = 1
175  
176        CALL DESCSET( DESCV , N + IROFF , 1 , DESCA( MB_ ) , 1 , IAROW , MYCOL ,
177       $ICTXT , MAX( 1 , NP ) )
178        CALL DESCSET( DESCX , N + IROFF , 1 , DESCA( MB_ ) , 1 , IAROW , MYCOL ,
179       $ICTXT , MAX( 1 , NP ) )
180  
181  *     Compute the norm of the triangular matrix A.
182  
183        ANORM = PZLANTR( NORM , UPLO , DIAG , N , N , A , IA , JA , DESCA , RWORK )
184  
185  *     Continue only if ANORM > 0.
186  
187        IF( ANORM.GT.ZERO ) THEN
188  
189  *         Estimate the norm of the inverse of A.
190  
191            AINVNM = ZERO
192            NORMIN = 'N'
193            IF( ONENRM ) THEN
194                KASE1 = 1
195            ELSE
196                KASE1 = 2
197            END IF
198            KASE = 0
199     10 CONTINUE
200        CALL PZLACON ( N , WORK( IPV ) , IV , JV , DESCV , WORK( IPX ) ,
201       $IX , JX , DESCX , AINVNM , KASE )
202        IF( KASE.NE.0 ) THEN
203            IF( KASE.EQ.KASE1 ) THEN
204  
205  *             Multiply by inv(A).
206  
207                DESCX( CSRC_ ) = IACOL
208                CALL PZLATRS ( UPLO , 'No transpose' , DIAG , NORMIN ,
209       $        N , A , IA , JA , DESCA , WORK( IPX ) , IX , JX ,
210       $        DESCX , SCALE , RWORK( IPN ) , WORK( IPW ) )
211                DESCX( CSRC_ ) = MYCOL
212            ELSE
213  
214  *             Multiply by inv(A').
215  
216                DESCX( CSRC_ ) = IACOL
217                CALL PZLATRS ( UPLO , 'Conjugate transpose' , DIAG , NORMIN ,
218       $        N , A , IA , JA , DESCA , WORK( IPX ) , IX , JX ,
219       $        DESCX , SCALE , RWORK( IPN ) , WORK( IPW ) )
220                DESCX( CSRC_ ) = MYCOL
221            END IF
222            NORMIN = 'Y'
223  
224  *         Multiply by 1 / SCALE if doing so will not cause overflow.
225  
226            IF( SCALE.NE.ONE ) THEN
227                CALL PZAMAX( N , WMAX , IXX , WORK( IPX ) , IX , JX ,
228       $        DESCX , 1 )
229                IF( DESCX( M_ ).EQ.1 .AND. N.EQ.1 ) THEN
230                    CALL PB_TOPGET( ICTXT , 'Broadcast' , 'Columnwise' ,
231       $            CBTOP )
232                    IF( MYROW.EQ.IAROW ) THEN
233                        CALL ZGEBS2D( ICTXT , 'Column' , CBTOP , 1 , 1 , WMAX ,
234       $                1 )
235                    ELSE
236                        CALL ZGEBR2D( ICTXT , 'Column' , CBTOP , 1 , 1 , WMAX ,
237       $                1 , IAROW , MYCOL )
238                    END IF
239                END IF
240                IF( SCALE.LT.CABS1( WMAX )*SMLNUM .OR. SCALE.EQ.ZERO )
241       $            GO TO 20
242                    CALL PZDRSCL ( N , SCALE , WORK( IPX ) , IX , JX , DESCX , 1 )
243                END IF
244                GO TO 10
245            END IF
246  
247  *         Compute the estimate of the reciprocal condition number.
248  
249            IF( AINVNM.NE.ZERO )
250       $        RCOND =( ONE / ANORM ) / AINVNM
251            END IF
252  
253     20 CONTINUE
254  
255        CALL PB_TOPSET( ICTXT , 'Combine' , 'Columnwise' , COLCTOP )
256        CALL PB_TOPSET( ICTXT , 'Combine' , 'Rowwise' , ROWCTOP )
257  
258        RETURN
259  
260  *     End of PZTRCON
261  
262        END