Routine: PSTRCON()  File: SRC\pstrcon.f

 
 
# lines: 426
  # code: 426
  # comment: 0
  # blank:0
# Variables:76
# Callers:0
# Callings:5
# Words:246
# Keywords:156
 

 

..
     .. Array Arguments ..
     ..
  Purpose
  =======
  PSTRCON 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) REAL 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) REAL
          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) REAL 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)) + LOCc(N+MOD(JA-1,NB_A))
          + MAX( 2, MAX( NB_A*MAX( 1, CEIL(NPROW-1,NPCOL) ),
                         LOCc(N+MOD(JA-1,NB_A)) +
                         NB_A*MAX( 1, CEIL(NPCOL-1,NPROW) ) ).
          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.
  IWORK   (local workspace/local output) INTEGER array,
                                                   dimension (LIWORK)
          On exit, IWORK(1) returns the minimal and optimal LIWORK.
  LIWORK  (local or global input) INTEGER
          The dimension of the array IWORK.
          LIWORK is local input and must be at least
          LIWORK >= LOCr(N+MOD(IA-1,MB_A)).
          If LIWORK = -1, then LIWORK 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 PSTRCON( NORM , UPLO , DIAG , N , A , IA , JA , DESCA , RCOND ,
002       $WORK , LWORK , IWORK , LIWORK , 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 , LIWORK , LWORK , N
012        REAL 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        REAL ONE , ZERO
019        PARAMETER( ONE = 1.0E + 0 , ZERO = 0.0E + 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 , LIWMIN , LWMIN , MYCOL , MYROW , NP , NPCOL ,
027       $NPMOD , NPROW , NQ , NQMOD
028        REAL AINVNM , ANORM , SCALE , SMLNUM
029        REAL WMAX
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       $PCHK1MAT , PSAMAX , PSLATRS , PSLACON ,
038       $PSRSCL , PB_TOPGET , PB_TOPSET , PXERBLA , SGEBR2D ,
039       $SGEBS2D
040  *     ..
041  *     .. External Functions ..
042        LOGICAL LSAME
043        INTEGER ICEIL , INDXG2P , NUMROC
044        REAL PSLAMCH , PSLANTR
045        EXTERNAL ICEIL , INDXG2P , LSAME , NUMROC , PSLAMCH ,
046       $PSLANTR
047  *     ..
048  *     .. Intrinsic Functions ..
049        INTRINSIC ABS , ICHAR , MAX , MOD , REAL
050  *     ..
051  *     .. Executable Statements ..
052  
053  *     Get grid parameters
054  
055        ICTXT = DESCA( CTXT_ )
056        CALL BLACS_GRIDINFO( ICTXT , NPROW , NPCOL , MYROW , MYCOL )
057  
058  *     Test the input parameters
059  
060        INFO = 0
061        IF( NPROW.EQ. - 1 ) THEN
062            INFO = - ( 800 + CTXT_ )
063        ELSE
064            CALL CHK1MAT( N , 4 , N , 4 , IA , JA , DESCA , 8 , INFO )
065            IF( INFO.EQ.0 ) THEN
066                UPPER = LSAME( UPLO , 'U' )
067                ONENRM = NORM.EQ.'1' .OR. LSAME( NORM , 'O' )
068                NOUNIT = LSAME( DIAG , 'N' )
069                IAROW = INDXG2P( IA , DESCA( MB_ ) , MYROW , DESCA( RSRC_ ) ,
070       $        NPROW )
071                IACOL = INDXG2P( JA , DESCA( NB_ ) , MYCOL , DESCA( CSRC_ ) ,
072       $        NPCOL )
073                NPMOD = NUMROC( N + MOD( IA - 1 , DESCA( MB_ ) ) , DESCA( MB_ ) ,
074       $        MYROW , IAROW , NPROW )
075                NQMOD = NUMROC( N + MOD( JA - 1 , DESCA( NB_ ) ) , DESCA( NB_ ) ,
076       $        MYCOL , IACOL , NPCOL )
077                LWMIN = 2*NPMOD + NQMOD +
078       $        MAX( 2 , MAX( DESCA( NB_ )*
079       $        MAX( 1 , ICEIL( NPROW - 1 , NPCOL ) ) , NQMOD +
080       $        DESCA( NB_ )*
081       $        MAX( 1 , ICEIL( NPCOL - 1 , NPROW ) ) ) )
082                WORK( 1 ) = REAL( LWMIN )
083                LIWMIN = NPMOD
084                IWORK( 1 ) = LIWMIN
085                LQUERY =( LWORK.EQ. - 1 .OR. LIWORK.EQ. - 1 )
086  
087                IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM , 'I' ) ) THEN
088                    INFO = - 1
089                ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO , 'L' ) ) THEN
090                    INFO = - 2
091                ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG , 'U' ) ) THEN
092                    INFO = - 3
093                ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
094                    INFO = - 11
095                ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
096                    INFO = - 13
097                END IF
098            END IF
099  
100            IF( ONENRM ) THEN
101                IDUM1( 1 ) = ICHAR( '1' )
102            ELSE
103                IDUM1( 1 ) = ICHAR( 'I' )
104            END IF
105            IDUM2( 1 ) = 1
106            IF( UPPER ) THEN
107                IDUM1( 2 ) = ICHAR( 'U' )
108            ELSE
109                IDUM1( 2 ) = ICHAR( 'L' )
110            END IF
111            IDUM2( 2 ) = 2
112            IF( NOUNIT ) THEN
113                IDUM1( 3 ) = ICHAR( 'N' )
114            ELSE
115                IDUM1( 3 ) = ICHAR( 'U' )
116            END IF
117            IDUM2( 3 ) = 3
118            IF( LWORK.EQ. - 1 ) THEN
119                IDUM1( 4 ) = - 1
120            ELSE
121                IDUM1( 4 ) = 1
122            END IF
123            IDUM2( 4 ) = 11
124            IF( LIWORK.EQ. - 1 ) THEN
125                IDUM1( 5 ) = - 1
126            ELSE
127                IDUM1( 5 ) = 1
128            END IF
129            IDUM2( 5 ) = 13
130            CALL PCHK1MAT( N , 4 , N , 4 , IA , JA , DESCA , 8 , 5 , IDUM1 , IDUM2 ,
131       $    INFO )
132        END IF
133  
134        IF( INFO.NE.0 ) THEN
135            CALL PXERBLA( ICTXT , 'PSTRCON' , - INFO )
136            RETURN
137        ELSE IF( LQUERY ) THEN
138            RETURN
139        END IF
140  
141  *     Quick return if possible
142  
143        IF( N.EQ.0 ) THEN
144            RCOND = ONE
145            RETURN
146        END IF
147  
148        CALL PB_TOPGET( ICTXT , 'Combine' , 'Columnwise' , COLCTOP )
149        CALL PB_TOPGET( ICTXT , 'Combine' , 'Rowwise' , ROWCTOP )
150        CALL PB_TOPSET( ICTXT , 'Combine' , 'Columnwise' , '1 - tree' )
151        CALL PB_TOPSET( ICTXT , 'Combine' , 'Rowwise' , '1 - tree' )
152  
153        RCOND = ZERO
154        SMLNUM = PSLAMCH( ICTXT , 'Safe minimum' )*REAL( MAX( 1 , N ) )
155        CALL INFOG2L( IA , JA , DESCA , NPROW , NPCOL , MYROW , MYCOL , IIA , JJA ,
156       $IAROW , IACOL )
157        IROFF = MOD( IA - 1 , DESCA( MB_ ) )
158        ICOFF = MOD( JA - 1 , DESCA( NB_ ) )
159        NP = NUMROC( N + IROFF , DESCA( MB_ ) , MYROW , IAROW , NPROW )
160        NQ = NUMROC( N + ICOFF , DESCA( NB_ ) , MYCOL , IACOL , NPCOL )
161        IV = IROFF + 1
162        IX = IV
163        JV = ICOFF + 1
164        JX = JV
165  
166        IPX = 1
167        IPV = IPX + NP
168        IPN = IPV + NP
169        IPW = IPN + NQ
170  
171        CALL DESCSET( DESCV , N + IROFF , 1 , DESCA( MB_ ) , 1 , IAROW , MYCOL ,
172       $ICTXT , MAX( 1 , NP ) )
173        CALL DESCSET( DESCX , N + IROFF , 1 , DESCA( MB_ ) , 1 , IAROW , MYCOL ,
174       $ICTXT , MAX( 1 , NP ) )
175  
176  *     Compute the norm of the triangular matrix A.
177  
178        ANORM = PSLANTR( NORM , UPLO , DIAG , N , N , A , IA , JA , DESCA , WORK )
179  
180  *     Continue only if ANORM > 0.
181  
182        IF( ANORM.GT.ZERO ) THEN
183  
184  *         Estimate the norm of the inverse of A.
185  
186            AINVNM = ZERO
187            NORMIN = 'N'
188            IF( ONENRM ) THEN
189                KASE1 = 1
190            ELSE
191                KASE1 = 2
192            END IF
193            KASE = 0
194     10 CONTINUE
195        CALL PSLACON ( N , WORK( IPV ) , IV , JV , DESCV , WORK( IPX ) ,
196       $IX , JX , DESCX , IWORK , AINVNM , KASE )
197        IF( KASE.NE.0 ) THEN
198            IF( KASE.EQ.KASE1 ) THEN
199  
200  *             Multiply by inv(A).
201  
202                DESCX( CSRC_ ) = IACOL
203                CALL PSLATRS ( UPLO , 'No transpose' , DIAG , NORMIN ,
204       $        N , A , IA , JA , DESCA , WORK( IPX ) , IX , JX ,
205       $        DESCX , SCALE , WORK( IPN ) , WORK( IPW ) )
206                DESCX( CSRC_ ) = MYCOL
207            ELSE
208  
209  *             Multiply by inv(A').
210  
211                DESCX( CSRC_ ) = IACOL
212                CALL PSLATRS ( UPLO , 'Transpose' , DIAG , NORMIN ,
213       $        N , A , IA , JA , DESCA , WORK( IPX ) , IX , JX ,
214       $        DESCX , SCALE , WORK( IPN ) , WORK( IPW ) )
215                DESCX( CSRC_ ) = MYCOL
216            END IF
217            NORMIN = 'Y'
218  
219  *         Multiply by 1 / SCALE if doing so will not cause overflow.
220  
221            IF( SCALE.NE.ONE ) THEN
222                CALL PSAMAX( N , WMAX , IXX , WORK( IPX ) , IX , JX ,
223       $        DESCX , 1 )
224                IF( DESCX( M_ ).EQ.1 .AND. N.EQ.1 ) THEN
225                    CALL PB_TOPGET( ICTXT , 'Broadcast' , 'Columnwise' ,
226       $            CBTOP )
227                    IF( MYROW.EQ.IAROW ) THEN
228                        CALL SGEBS2D( ICTXT , 'Column' , CBTOP , 1 , 1 , WMAX ,
229       $                1 )
230                    ELSE
231                        CALL SGEBR2D( ICTXT , 'Column' , CBTOP , 1 , 1 , WMAX ,
232       $                1 , IAROW , MYCOL )
233                    END IF
234                END IF
235                IF( SCALE.LT.ABS( WMAX )*SMLNUM .OR. SCALE.EQ.ZERO )
236       $            GO TO 20
237                    CALL PSRSCL ( N , SCALE , WORK( IPX ) , IX , JX , DESCX , 1 )
238                END IF
239                GO TO 10
240            END IF
241  
242  *         Compute the estimate of the reciprocal condition number.
243  
244            IF( AINVNM.NE.ZERO )
245       $        RCOND =( ONE / ANORM ) / AINVNM
246            END IF
247  
248     20 CONTINUE
249  
250        CALL PB_TOPSET( ICTXT , 'Combine' , 'Columnwise' , COLCTOP )
251        CALL PB_TOPSET( ICTXT , 'Combine' , 'Rowwise' , ROWCTOP )
252  
253        RETURN
254  
255  *     End of PSTRCON
256  
257        END