Routine: PDPOEQU()  File: SRC\pdpoequ.f

 
 
# lines: 357
  # code: 357
  # comment: 0
  # blank:0
# Variables:56
# Callers:1
# Callings:1
# Words:172
# Keywords:102
 

 

 
..
     .. Array Arguments ..
     ..
  Purpose
  =======
  PDPOEQU computes row and column scalings intended to
  equilibrate a distributed symmetric positive definite matrix
  sub( A ) = A(IA:IA+N-1,JA:JA+N-1) and reduce its condition number
  (with respect to the two-norm).  SR and SC contain the scale
  factors, S(i) = 1/sqrt(A(i,i)), chosen so that the scaled distri-
  buted matrix B with elements B(i,j) = S(i)*A(i,j)*S(j) has ones on
  the  diagonal.  This choice of SR and SC puts the condition number
  of B within a factor N of the smallest possible condition number
  over all possible diagonal scalings.
  The scaling factor are stored along process rows in SR and along
  process columns in SC. The duplication of information simplifies
  greatly the application of the factors.
  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
  =========
  N       (global input) INTEGER
          The number of rows and columns to be operated on i.e the
          order of the distributed submatrix sub( A ). N >= 0.
  A       (local input) DOUBLE PRECISION pointer into the local memory
          to an array of local dimension ( LLD_A, LOCc(JA+N-1) ), the
          N-by-N symmetric positive definite distributed matrix
          sub( A ) whose scaling factors are to be computed.  Only the
          diagonal elements of sub( A ) are referenced.
  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.
  SR      (local output) DOUBLE PRECISION array, dimension LOCr(M_A)
          If INFO = 0, SR(IA:IA+N-1) contains the row scale factors
          for sub( A ). SR is aligned with the distributed matrix A,
          and replicated across every process column. SR is tied to the
          distributed matrix A.
  SC      (local output) DOUBLE PRECISION array, dimension LOCc(N_A)
          If INFO = 0, SC(JA:JA+N-1) contains the column scale factors
          for A(IA:IA+M-1,JA:JA+N-1). SC is aligned with the distribu-
          ted matrix A, and replicated down every process row. SC is
          tied to the distributed matrix A.
  SCOND   (global output) DOUBLE PRECISION
          If INFO = 0, SCOND contains the ratio of the smallest SR(i)
          (or SC(j)) to the largest SR(i) (or SC(j)), with
          IA <= i <= IA+N-1 and JA <= j <= JA+N-1. If SCOND >= 0.1
          and AMAX is neither too large nor too small, it is not worth
          scaling by SR (or SC).
  AMAX    (global output) DOUBLE PRECISION
          Absolute value of largest matrix element.  If AMAX is very
          close to overflow or very close to underflow, the matrix
          should be scaled.
  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.
          > 0:  If INFO = K, the K-th diagonal entry of sub( A ) is
                nonpositive.
  =====================================================================
     .. Parameters ..
 
Display dynamic version Find AutoScroll Reload FontSize: - + Hide Comments Hide Blanks Frame FullScreen MailPrint

  
001        SUBROUTINE PDPOEQU( N , A , IA , JA , DESCA , SR , SC , SCOND , AMAX ,
002       $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        INTEGER IA , INFO , JA , N
011        DOUBLE PRECISION AMAX , SCOND
012        INTEGER BLOCK_CYCLIC_2D , CSRC_ , CTXT_ , DLEN_ , DTYPE_ ,
013       $LLD_ , MB_ , M_ , NB_ , N_ , RSRC_
014        PARAMETER( BLOCK_CYCLIC_2D = 1 , DLEN_ = 9 , DTYPE_ = 1 ,
015       $CTXT_ = 2 , M_ = 3 , N_ = 4 , MB_ = 5 , NB_ = 6 ,
016       $RSRC_ = 7 , CSRC_ = 8 , LLD_ = 9 )
017        DOUBLE PRECISION ZERO , ONE
018        PARAMETER( ZERO = 0.0D + 0 , ONE = 1.0D + 0 )
019  *     ..
020  *     .. Local Scalars ..
021        CHARACTER ALLCTOP , COLCTOP , ROWCTOP
022        INTEGER IACOL , IAROW , ICOFF , ICTXT , ICURCOL , ICURROW ,
023       $IDUMM , II , IIA , IOFFA , IOFFD , IROFF , J , JB , JJ ,
024       $JJA , JN , LDA , LL , MYCOL , MYROW , NP , NPCOL ,
025       $NPROW , NQ
026        DOUBLE PRECISION AII , SMIN
027  *     ..
028  *     .. Local Arrays ..
029        INTEGER DESCSC( DLEN_ ) , DESCSR( DLEN_ )
030  *     ..
031  *     .. External Subroutines ..
032        EXTERNAL BLACS_GRIDINFO , CHK1MAT , DESCSET , DGAMN2D ,
033       $DGAMX2D , DGSUM2D , IGAMN2D , INFOG2L ,
034       $PCHK1MAT , PB_TOPGET , PXERBLA
035  *     ..
036  *     .. External Functions ..
037        INTEGER ICEIL , NUMROC
038        DOUBLE PRECISION PDLAMCH
039        EXTERNAL ICEIL , NUMROC , PDLAMCH
040  *     ..
041  *     .. Intrinsic Functions ..
042        INTRINSIC MAX , MIN , MOD , SQRT
043  *     ..
044  *     .. Executable Statements ..
045  
046  *     Get grid parameters
047  
048        ICTXT = DESCA( CTXT_ )
049        CALL BLACS_GRIDINFO( ICTXT , NPROW , NPCOL , MYROW , MYCOL )
050  
051  *     Test the input parameters.
052  
053        INFO = 0
054        IF( NPROW.EQ. - 1 ) THEN
055            INFO = - (500 + CTXT_)
056        ELSE
057            CALL CHK1MAT( N , 1 , N , 1 , IA , JA , DESCA , 5 , INFO )
058            CALL PCHK1MAT( N , 1 , N , 1 , IA , JA , DESCA , 5 , 0 , IDUMM , IDUMM ,
059       $    INFO )
060        END IF
061  
062        IF( INFO.NE.0 ) THEN
063            CALL PXERBLA( ICTXT , 'PDPOEQU' , - INFO )
064            RETURN
065        END IF
066  
067  *     Quick return if possible
068  
069        IF( N.EQ.0 ) THEN
070            SCOND = ONE
071            AMAX = ZERO
072            RETURN
073        END IF
074  
075        CALL PB_TOPGET( ICTXT , 'Combine' , 'All' , ALLCTOP )
076        CALL PB_TOPGET( ICTXT , 'Combine' , 'Rowwise' , ROWCTOP )
077        CALL PB_TOPGET( ICTXT , 'Combine' , 'Columnwise' , COLCTOP )
078  
079  *     Compute some local indexes
080  
081        CALL INFOG2L( IA , JA , DESCA , NPROW , NPCOL , MYROW , MYCOL , IIA , JJA ,
082       $IAROW , IACOL )
083        IROFF = MOD( IA - 1 , DESCA( MB_ ) )
084        ICOFF = MOD( JA - 1 , DESCA( NB_ ) )
085        NP = NUMROC( N + IROFF , DESCA( MB_ ) , MYROW , IAROW , NPROW )
086        NQ = NUMROC( N + ICOFF , DESCA( NB_ ) , MYCOL , IACOL , NPCOL )
087        IF( MYROW.EQ.IAROW )
088       $    NP = NP - IROFF
089            IF( MYCOL.EQ.IACOL )
090       $        NQ = NQ - ICOFF
091                JN = MIN( ICEIL( JA , DESCA( NB_ ) ) * DESCA( NB_ ) , JA + N - 1 )
092                LDA = DESCA( LLD_ )
093  
094  *             Assign descriptors for SR and SC arrays
095  
096                CALL DESCSET( DESCSR , N , 1 , DESCA( MB_ ) , 1 , 0 , 0 , ICTXT ,
097       $        MAX( 1 , NP ) )
098                CALL DESCSET( DESCSC , 1 , N , 1 , DESCA( NB_ ) , 0 , 0 , ICTXT , 1 )
099  
100  *             Initialize the scaling factors to zero.
101  
102                DO 10 II = IIA , IIA + NP - 1
103                    SR( II ) = ZERO
104     10         CONTINUE
105  
106                DO 20 JJ = JJA , JJA + NQ - 1
107                    SC( JJ ) = ZERO
108     20         CONTINUE
109  
110  *             Find the minimum and maximum diagonal elements.
111  *             Handle first block separately.
112  
113                II = IIA
114                JJ = JJA
115                JB = JN - JA + 1
116                SMIN = ONE / PDLAMCH( ICTXT , 'S' )
117                AMAX = ZERO
118  
119                IOFFA = II + (JJ - 1)*LDA
120                IF( MYROW.EQ.IAROW .AND. MYCOL.EQ.IACOL ) THEN
121                    IOFFD = IOFFA
122                    DO 30 LL = 0 , JB - 1
123                        AII = A( IOFFD )
124                        SR( II + LL ) = AII
125                        SC( JJ + LL ) = AII
126                        SMIN = MIN( SMIN , AII )
127                        AMAX = MAX( AMAX , AII )
128                        IF( AII.LE.ZERO .AND. INFO.EQ.0 )
129       $                    INFO = LL + 1
130                            IOFFD = IOFFD + LDA + 1
131     30             CONTINUE
132                END IF
133  
134                IF( MYROW.EQ.IAROW ) THEN
135                    II = II + JB
136                    IOFFA = IOFFA + JB
137                END IF
138                IF( MYCOL.EQ.IACOL ) THEN
139                    JJ = JJ + JB
140                    IOFFA = IOFFA + JB*LDA
141                END IF
142                ICURROW = MOD( IAROW + 1 , NPROW )
143                ICURCOL = MOD( IACOL + 1 , NPCOL )
144  
145  *             Loop over remaining blocks of columns
146  
147                DO 50 J = JN + 1 , JA + N - 1 , DESCA( NB_ )
148                    JB = MIN( N - J + JA , DESCA( NB_ ) )
149  
150                    IF( MYROW.EQ.ICURROW .AND. MYCOL.EQ.ICURCOL ) THEN
151                        IOFFD = IOFFA
152                        DO 40 LL = 0 , JB - 1
153                            AII = A( IOFFD )
154                            SR( II + LL ) = AII
155                            SC( JJ + LL ) = AII
156                            SMIN = MIN( SMIN , AII )
157                            AMAX = MAX( AMAX , AII )
158                            IF( AII.LE.ZERO .AND. INFO.EQ.0 )
159       $                        INFO = J + LL - JA + 1
160                                IOFFD = IOFFD + LDA + 1
161     40                 CONTINUE
162                    END IF
163  
164                    IF( MYROW.EQ.ICURROW ) THEN
165                        II = II + JB
166                        IOFFA = IOFFA + JB
167                    END IF
168                    IF( MYCOL.EQ.ICURCOL ) THEN
169                        JJ = JJ + JB
170                        IOFFA = IOFFA + JB*LDA
171                    END IF
172                    ICURROW = MOD( ICURROW + 1 , NPROW )
173                    ICURCOL = MOD( ICURCOL + 1 , NPCOL )
174  
175     50         CONTINUE
176  
177  *             Compute scaling factors
178  
179                CALL DGSUM2D( ICTXT , 'Columnwise' , COLCTOP , 1 , NQ , SC( JJA ) ,
180       $        1 , - 1 , MYCOL )
181                CALL DGSUM2D( ICTXT , 'Rowwise' , ROWCTOP , NP , 1 , SR( IIA ) ,
182       $        MAX( 1 , NP ) , - 1 , MYCOL )
183  
184                CALL DGAMX2D( ICTXT , 'All' , ALLCTOP , 1 , 1 , AMAX , 1 , IDUMM , IDUMM ,
185       $        - 1 , - 1 , MYCOL )
186                CALL DGAMN2D( ICTXT , 'All' , ALLCTOP , 1 , 1 , SMIN , 1 , IDUMM , IDUMM ,
187       $        - 1 , - 1 , MYCOL )
188  
189                IF( SMIN.LE.ZERO ) THEN
190  
191  *                 Find the first non - positive diagonal element and return.
192  
193                    CALL IGAMN2D( ICTXT , 'All' , ALLCTOP , 1 , 1 , INFO , 1 , II , JJ , - 1 ,
194       $            - 1 , MYCOL )
195                    RETURN
196  
197                ELSE
198  
199  *                 Set the scale factors to the reciprocals
200  *                 of the diagonal elements.
201  
202                    DO 60 II = IIA , IIA + NP - 1
203                        SR( II ) = ONE / SQRT( SR( II ) )
204     60             CONTINUE
205  
206                    DO 70 JJ = JJA , JJA + NQ - 1
207                        SC( JJ ) = ONE / SQRT( SC( JJ ) )
208     70             CONTINUE
209  
210  *                 Compute SCOND = min(S(I)) / max(S(I))
211  
212                    SCOND = SQRT( SMIN ) / SQRT( AMAX )
213  
214                END IF
215  
216                RETURN
217  
218  *             End of PDPOEQU
219  
220            END