Routine: PZLASMSUB()  File: SRC\pzlasmsub.f

 
 
# lines: 377
  # code: 377
  # comment: 0
  # blank:0
# Variables:59
# Callers:0
# Callings:1
# Words:225
# Keywords:136
 

 

 
..
     .. Array Arguments ..
     ..
  Purpose
  =======
  PZLASMSUB looks for a small subdiagonal element from the bottom
     of the matrix that it can safely set to zero.
  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
  =========
  A       (global input) COMPLEX*16 array, dimension (DESCA(LLD_),*)
          On entry, the Hessenberg matrix whose tridiagonal part is
          being scanned.
          Unchanged on exit.
  DESCA   (global and local input) INTEGER array of dimension DLEN_.
          The array descriptor for the distributed matrix A.
  I       (global input) INTEGER
          The global location of the bottom of the unreduced
          submatrix of A.
          Unchanged on exit.
  L       (global input) INTEGER
          The global location of the top of the unreduced submatrix
          of A.
          Unchanged on exit.
  K       (global output) INTEGER
          On exit, this yields the bottom portion of the unreduced
          submatrix.  This will satisfy: L <= M  <= I-1.
  SMLNUM  (global input) DOUBLE PRECISION
          On entry, a "small number" for the given matrix.
          Unchanged on exit.
  BUF     (local output) COMPLEX*16 array of size LWORK.
  LWORK   (global input) INTEGER
          On exit, LWORK is the size of the work buffer.
          This must be at least 2*Ceil( Ceil( (I-L)/HBL ) /
                                        LCM(NPROW,NPCOL) )
          Here LCM is least common multiple, and NPROWxNPCOL is the
          logical grid size.
   Notes:
     This routine does a global maximum and must be called by all
     processes.
     This code is basically a parallelization of the following snip
     of LAPACK code from ZLAHQR:
        Look for a single small subdiagonal element.
        DO 20 K = I, L + 1, -1
           TST1 = CABS1( H( K-1, K-1 ) ) + CABS1( H( K, K ) )
           IF( TST1.EQ.ZERO )
    $         TST1 = ZLANHS( '1', I-L+1, H( L, L ), LDH, WORK )
           IF( CABS1( H( K, K-1 ) ).LE.MAX( ULP*TST1, SMLNUM ) )
    $         GO TO 30
  20    CONTINUE
  30    CONTINUE
  Further Details
  ===============
  Implemented by:  M. Fahey, May 28, 1999
  =====================================================================
     .. Parameters ..
 
Display dynamic version Find AutoScroll Reload FontSize: - + Hide Comments Hide Blanks Frame FullScreen MailPrint

  
001        SUBROUTINE PZLASMSUB( A , DESCA , I , L , K , SMLNUM , BUF , LWORK )
002  
003  *     -- ScaLAPACK routine(version 1.7) --
004  *     University of Tennessee , Knoxville , Oak Ridge National Laboratory ,
005  *     and University of California , Berkeley.
006  *     July 31 , 2001
007  
008  *     .. Scalar Arguments ..
009        INTEGER I , K , L , LWORK
010        DOUBLE PRECISION SMLNUM
011        INTEGER BLOCK_CYCLIC_2D , CSRC_ , CTXT_ , DLEN_ , DTYPE_ ,
012       $LLD_ , MB_ , M_ , NB_ , N_ , RSRC_
013        PARAMETER( BLOCK_CYCLIC_2D = 1 , DLEN_ = 9 , DTYPE_ = 1 ,
014       $CTXT_ = 2 , M_ = 3 , N_ = 4 , MB_ = 5 , NB_ = 6 ,
015       $RSRC_ = 7 , CSRC_ = 8 , LLD_ = 9 )
016        DOUBLE PRECISION ZERO
017        PARAMETER( ZERO = 0.0D + 0 )
018  *     ..
019  *     .. Local Scalars ..
020        INTEGER CONTXT , DOWN , HBL , IBUF1 , IBUF2 , ICOL1 , ICOL2 ,
021       $II , III , IRCV1 , IRCV2 , IROW1 , IROW2 , ISRC ,
022       $ISTR1 , ISTR2 , ITMP1 , ITMP2 , JJ , JJJ , JSRC , LDA ,
023       $LEFT , MODKM1 , MYCOL , MYROW , NPCOL , NPROW , NUM ,
024       $RIGHT , UP
025        DOUBLE PRECISION TST1 , ULP
026        COMPLEX*16 CDUM , H10 , H11 , H22
027  *     ..
028  *     .. External Functions ..
029        INTEGER ILCM , NUMROC
030        DOUBLE PRECISION PDLAMCH
031        EXTERNAL ILCM , NUMROC , PDLAMCH
032  *     ..
033  *     .. External Subroutines ..
034        EXTERNAL BLACS_GRIDINFO , IGAMX2D , INFOG1L , INFOG2L ,
035       $ZGERV2D , ZGESD2D
036  *     ..
037  *     .. Intrinsic Functions ..
038        INTRINSIC ABS , DBLE , DIMAG , MAX , MOD
039  *     ..
040  *     .. Statement Functions ..
041        DOUBLE PRECISION CABS1
042  *     ..
043  *     .. Statement Function definitions ..
044        CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
045  *     ..
046  *     .. Executable Statements ..
047  
048        HBL = DESCA( MB_ )
049        CONTXT = DESCA( CTXT_ )
050        LDA = DESCA( LLD_ )
051        ULP = PDLAMCH( CONTXT , 'PRECISION' )
052        CALL BLACS_GRIDINFO( CONTXT , NPROW , NPCOL , MYROW , MYCOL )
053        LEFT = MOD( MYCOL + NPCOL - 1 , NPCOL )
054        RIGHT = MOD( MYCOL + 1 , NPCOL )
055        UP = MOD( MYROW + NPROW - 1 , NPROW )
056        DOWN = MOD( MYROW + 1 , NPROW )
057        NUM = NPROW*NPCOL
058  
059  *     BUFFER1 STARTS AT BUF(ISTR1 + 1) AND WILL CONTAINS IBUF1 ELEMENTS
060  *     BUFFER2 STARTS AT BUF(ISTR2 + 1) AND WILL CONTAINS IBUF2 ELEMENTS
061  
062        ISTR1 = 0
063        ISTR2 =(( I - L ) / HBL )
064        IF( ISTR2*HBL.LT.( I - L ) )
065       $    ISTR2 = ISTR2 + 1
066            II = ISTR2 / ILCM( NPROW , NPCOL )
067            IF( II*ILCM( NPROW , NPCOL ).LT.ISTR2 ) THEN
068                ISTR2 = II + 1
069            ELSE
070                ISTR2 = II
071            END IF
072            IF( LWORK.LT.2*ISTR2 ) THEN
073  
074  *             Error !
075  
076                RETURN
077            END IF
078            CALL INFOG2L( I , I , DESCA , NPROW , NPCOL , MYROW , MYCOL , IROW1 ,
079       $    ICOL1 , II , JJ )
080            MODKM1 = MOD( I - 1 + HBL , HBL )
081  
082  *         COPY OUR RELEVANT PIECES OF TRIADIAGONAL THAT WE OWE INTO
083  *         2 BUFFERS TO SEND TO WHOMEVER OWNS H(K , K) AS K MOVES DIAGONALLY
084  *         UP THE TRIDIAGONAL
085  
086            IBUF1 = 0
087            IBUF2 = 0
088            IRCV1 = 0
089            IRCV2 = 0
090            DO 10 K = I , L + 1 , - 1
091                IF(( MODKM1.EQ.0 ) .AND.( DOWN.EQ.II ) .AND.
092       $( RIGHT.EQ.JJ ) ) THEN
093  
094  *             WE MUST PACK H(K - 1 , K - 1) AND SEND IT DIAGONAL DOWN
095  
096                IF(( DOWN.NE.MYROW ) .OR.( RIGHT.NE.MYCOL ) ) THEN
097                    CALL INFOG2L( K - 1 , K - 1 , DESCA , NPROW , NPCOL , MYROW ,
098       $            MYCOL , IROW1 , ICOL1 , ISRC , JSRC )
099                    IBUF1 = IBUF1 + 1
100                    BUF( ISTR1 + IBUF1 ) = A(( ICOL1 - 1 )*LDA + IROW1 )
101                END IF
102            END IF
103            IF(( MODKM1.EQ.0 ) .AND.( MYROW.EQ.II ) .AND.
104       $( RIGHT.EQ.JJ ) ) THEN
105  
106  *         WE MUST PACK H(K , K - 1) AND SEND IT RIGHT
107  
108            IF( NPCOL.GT.1 ) THEN
109                CALL INFOG2L( K , K - 1 , DESCA , NPROW , NPCOL , MYROW , MYCOL ,
110       $        IROW1 , ICOL1 , ISRC , JSRC )
111                IBUF2 = IBUF2 + 1
112                BUF( ISTR2 + IBUF2 ) = A(( ICOL1 - 1 )*LDA + IROW1 )
113            END IF
114        END IF
115  
116  *     ADD UP THE RECEIVES
117  
118        IF(( MYROW.EQ.II ) .AND.( MYCOL.EQ.JJ ) ) THEN
119            IF(( MODKM1.EQ.0 ) .AND.(( NPROW.GT.1 ) .OR.( NPCOL.GT.
120       $    1 ) ) ) THEN
121  
122  *         WE MUST RECEIVE H(K - 1 , K - 1) FROM DIAGONAL UP
123  
124            IRCV1 = IRCV1 + 1
125        END IF
126        IF(( MODKM1.EQ.0 ) .AND.( NPCOL.GT.1 ) ) THEN
127  
128  *         WE MUST RECEIVE H(K , K - 1) FROM LEFT
129  
130            IRCV2 = IRCV2 + 1
131        END IF
132        END IF
133  
134  *     POSSIBLY CHANGE OWNERS(OCCURS ONLY WHEN MOD(K - 1 , HBL) = 0)
135  
136        IF( MODKM1.EQ.0 ) THEN
137            II = II - 1
138            JJ = JJ - 1
139            IF( II.LT.0 )
140       $        II = NPROW - 1
141                IF( JJ.LT.0 )
142       $            JJ = NPCOL - 1
143                END IF
144                MODKM1 = MODKM1 - 1
145                IF( MODKM1.LT.0 )
146       $            MODKM1 = HBL - 1
147     10 CONTINUE
148  
149  *     SEND DATA ON TO THE APPROPRIATE NODE IF THERE IS ANY DATA TO SEND
150  
151        IF( IBUF1.GT.0 ) THEN
152            CALL ZGESD2D( CONTXT , IBUF1 , 1 , BUF( ISTR1 + 1 ) , IBUF1 , DOWN ,
153       $    RIGHT )
154        END IF
155        IF( IBUF2.GT.0 ) THEN
156            CALL ZGESD2D( CONTXT , IBUF2 , 1 , BUF( ISTR2 + 1 ) , IBUF2 , MYROW ,
157       $    RIGHT )
158        END IF
159  
160  *     RECEIVE APPROPRIATE DATA IF THERE IS ANY
161  
162        IF( IRCV1.GT.0 ) THEN
163            CALL ZGERV2D( CONTXT , IRCV1 , 1 , BUF( ISTR1 + 1 ) , IRCV1 , UP ,
164       $    LEFT )
165        END IF
166        IF( IRCV2.GT.0 ) THEN
167            CALL ZGERV2D( CONTXT , IRCV2 , 1 , BUF( ISTR2 + 1 ) , IRCV2 , MYROW ,
168       $    LEFT )
169        END IF
170  
171  *     START MAIN LOOP
172  
173        IBUF1 = 0
174        IBUF2 = 0
175        CALL INFOG2L( I , I , DESCA , NPROW , NPCOL , MYROW , MYCOL , IROW1 ,
176       $ICOL1 , II , JJ )
177        MODKM1 = MOD( I - 1 + HBL , HBL )
178  
179  *     LOOK FOR A SINGLE SMALL SUBDIAGONAL ELEMENT.
180  
181  *     Start loop for subdiagonal search
182  
183        DO 40 K = I , L + 1 , - 1
184            IF(( MYROW.EQ.II ) .AND.( MYCOL.EQ.JJ ) ) THEN
185                IF( MODKM1.EQ.0 ) THEN
186  
187  *                 Grab information from WORK array
188  
189                    IF( NUM.GT.1 ) THEN
190                        IBUF1 = IBUF1 + 1
191                        H11 = BUF( ISTR1 + IBUF1 )
192                    ELSE
193                        H11 = A(( ICOL1 - 2 )*LDA + IROW1 - 1 )
194                    END IF
195                    IF( NPCOL.GT.1 ) THEN
196                        IBUF2 = IBUF2 + 1
197                        H10 = BUF( ISTR2 + IBUF2 )
198                    ELSE
199                        H10 = A(( ICOL1 - 2 )*LDA + IROW1 )
200                    END IF
201                ELSE
202  
203  *                 Information is local
204  
205                    H11 = A(( ICOL1 - 2 )*LDA + IROW1 - 1 )
206                    H10 = A(( ICOL1 - 2 )*LDA + IROW1 )
207                END IF
208                H22 = A(( ICOL1 - 1 )*LDA + IROW1 )
209                TST1 = CABS1( H11 ) + CABS1( H22 )
210                IF( TST1.EQ.ZERO ) THEN
211  
212  *                 FIND SOME NORM OF THE LOCAL H(L : I , L : I)
213  
214                    CALL INFOG1L( L , HBL , NPROW , MYROW , 0 , ITMP1 , III )
215                    IROW2 = NUMROC( I , HBL , MYROW , 0 , NPROW )
216                    CALL INFOG1L( L , HBL , NPCOL , MYCOL , 0 , ITMP2 , III )
217                    ICOL2 = NUMROC( I , HBL , MYCOL , 0 , NPCOL )
218                    DO 30 III = ITMP1 , IROW2
219                        DO 20 JJJ = ITMP2 , ICOL2
220                            TST1 = TST1 + CABS1( A(( JJJ - 1 )*LDA + III ) )
221     20                 CONTINUE
222     30             CONTINUE
223                END IF
224                IF( CABS1( H10 ).LE.MAX( ULP*TST1 , SMLNUM ) )
225       $            GO TO 50
226                    IROW1 = IROW1 - 1
227                    ICOL1 = ICOL1 - 1
228                END IF
229                MODKM1 = MODKM1 - 1
230                IF( MODKM1.LT.0 )
231       $            MODKM1 = HBL - 1
232                    IF(( MODKM1.EQ.HBL - 1 ) .AND.( K.GT.2 ) ) THEN
233                        II = MOD( II + NPROW - 1 , NPROW )
234                        JJ = MOD( JJ + NPCOL - 1 , NPCOL )
235                        CALL INFOG2L( K - 1 , K - 1 , DESCA , NPROW , NPCOL , MYROW , MYCOL ,
236       $                IROW1 , ICOL1 , ITMP1 , ITMP2 )
237                    END IF
238     40 CONTINUE
239     50 CONTINUE
240        CALL IGAMX2D( CONTXT , 'ALL' , ' ' , 1 , 1 , K , 1 , ITMP1 , ITMP2 , - 1 ,
241       $- 1 , - 1 )
242        RETURN
243  
244  *     End of PZLASMSUB
245  
246        END