Routine: PDLARFG()  File: SRC\pdlarfg.f

 
 
# lines: 281
  # code: 281
  # comment: 0
  # blank:0
# Variables:39
# Callers:11
# Callings:0
# Words:106
# Keywords:64
 

 

 
..
     .. Array Arguments ..
     ..
  Purpose
  =======
  PDLARFG generates a real elementary reflector H of order n, such
  that
     H * sub( X ) = H * ( x(iax,jax) ) = ( alpha ),   H' * H = I.
                        (      x     )   (   0   )
  where alpha is a scalar, and sub( X ) is an (N-1)-element real
  distributed vector X(IX:IX+N-2,JX) if INCX = 1 and X(IX,JX:JX+N-2) if
  INCX = DESCX(M_).  H is represented in the form
        H = I - tau * ( 1 ) * ( 1 v' ) ,
                      ( v )
  where tau is a real scalar and v is a real (N-1)-element
  vector.
  If the elements of sub( X ) are all zero, then tau = 0 and H is
  taken to be the unit matrix.
  Otherwise  1 <= tau <= 2.
  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
  Because vectors may be viewed as a subclass of matrices, a
  distributed vector is considered to be a distributed matrix.
  Arguments
  =========
  N       (global input) INTEGER
          The global order of the elementary reflector. N >= 0.
  ALPHA   (local output) DOUBLE PRECISION
          On exit, alpha is computed in the process scope having the
          vector sub( X ).
  IAX     (global input) INTEGER
          The global row index in X of X(IAX,JAX).
  JAX     (global input) INTEGER
          The global column index in X of X(IAX,JAX).
  X       (local input/local output) DOUBLE PRECISION, pointer into the
          local memory to an array of dimension (LLD_X,*). This array
          contains the local pieces of the distributed vector sub( X ).
          Before entry, the incremented array sub( X ) must contain
          the vector x. On exit, it is overwritten with the vector v.
  IX      (global input) INTEGER
          The row index in the global array X indicating the first
          row of sub( X ).
  JX      (global input) INTEGER
          The column index in the global array X indicating the
          first column of sub( X ).
  DESCX   (global and local input) INTEGER array of dimension DLEN_.
          The array descriptor for the distributed matrix X.
  INCX    (global input) INTEGER
          The global increment for the elements of X. Only two values
          of INCX are supported in this version, namely 1 and M_X.
          INCX must not be zero.
  TAU     (local output) DOUBLE PRECISION array, dimension  LOCc(JX)
          if INCX = 1, and LOCr(IX) otherwise. This array contains the
          Householder scalars related to the Householder vectors.
          TAU is tied to the distributed matrix X.
  =====================================================================
     .. Parameters ..
 
Display dynamic version Find AutoScroll Reload FontSize: - + Hide Comments Hide Blanks Frame FullScreen MailPrint

  
001        SUBROUTINE PDLARFG( N , ALPHA , IAX , JAX , X , IX , JX , DESCX , INCX ,
002       $TAU )
003  
004  *     -- ScaLAPACK auxiliary routine(version 1.7) --
005  *     University of Tennessee , Knoxville , Oak Ridge National Laboratory ,
006  *     and University of California , Berkeley.
007  *     May 1 , 1997
008  
009  *     .. Scalar Arguments ..
010        INTEGER IAX , INCX , IX , JAX , JX , N
011        DOUBLE PRECISION ALPHA
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 ONE , ZERO
018        PARAMETER( ONE = 1.0D + 0 , ZERO = 0.0D + 0 )
019  *     ..
020  *     .. Local Scalars ..
021        INTEGER ICTXT , IIAX , INDXTAU , IXCOL , IXROW , J , JJAX ,
022       $KNT , MYCOL , MYROW , NPCOL , NPROW
023        DOUBLE PRECISION BETA , RSAFMN , SAFMIN , XNORM
024  *     ..
025  *     .. External Subroutines ..
026        EXTERNAL BLACS_GRIDINFO , DGEBR2D , DGEBS2D , PDSCAL ,
027       $INFOG2L , PDNRM2
028  *     ..
029  *     .. External Functions ..
030        DOUBLE PRECISION DLAMCH , DLAPY2
031        EXTERNAL DLAMCH , DLAPY2
032  *     ..
033  *     .. Intrinsic Functions ..
034        INTRINSIC ABS , SIGN
035  *     ..
036  *     .. Executable Statements ..
037  
038  *     Get grid parameters.
039  
040        ICTXT = DESCX( CTXT_ )
041        CALL BLACS_GRIDINFO( ICTXT , NPROW , NPCOL , MYROW , MYCOL )
042  
043        IF( INCX.EQ.DESCX( M_ ) ) THEN
044  
045  *         sub( X ) is distributed across a process row.
046  
047            CALL INFOG2L( IX , JAX , DESCX , NPROW , NPCOL , MYROW , MYCOL ,
048       $    IIAX , JJAX , IXROW , IXCOL )
049  
050            IF( MYROW.NE.IXROW )
051       $        RETURN
052  
053  *             Broadcast X(IAX , JAX) across the process row.
054  
055                IF( MYCOL.EQ.IXCOL ) THEN
056                    J = IIAX + (JJAX - 1)*DESCX( LLD_ )
057                    CALL DGEBS2D( ICTXT , 'Rowwise' , ' ' , 1 , 1 , X( J ) , 1 )
058                    ALPHA = X( J )
059                ELSE
060                    CALL DGEBR2D( ICTXT , 'Rowwise' , ' ' , 1 , 1 , ALPHA , 1 ,
061       $            MYROW , IXCOL )
062                END IF
063  
064                INDXTAU = IIAX
065  
066            ELSE
067  
068  *             sub( X ) is distributed across a process column.
069  
070                CALL INFOG2L( IAX , JX , DESCX , NPROW , NPCOL , MYROW , MYCOL ,
071       $        IIAX , JJAX , IXROW , IXCOL )
072  
073                IF( MYCOL.NE.IXCOL )
074       $            RETURN
075  
076  *                 Broadcast X(IAX , JAX) across the process column.
077  
078                    IF( MYROW.EQ.IXROW ) THEN
079                        J = IIAX + (JJAX - 1)*DESCX( LLD_ )
080                        CALL DGEBS2D( ICTXT , 'Columnwise' , ' ' , 1 , 1 , X( J ) , 1 )
081                        ALPHA = X( J )
082                    ELSE
083                        CALL DGEBR2D( ICTXT , 'Columnwise' , ' ' , 1 , 1 , ALPHA , 1 ,
084       $                IXROW , MYCOL )
085                    END IF
086  
087                    INDXTAU = JJAX
088  
089                END IF
090  
091                IF( N.LE.0 ) THEN
092                    TAU( INDXTAU ) = ZERO
093                    RETURN
094                END IF
095  
096                CALL PDNRM2( N - 1 , XNORM , X , IX , JX , DESCX , INCX )
097  
098                IF( XNORM.EQ.ZERO ) THEN
099  
100  *                 H = I
101  
102                    TAU( INDXTAU ) = ZERO
103  
104                ELSE
105  
106  *                 General case
107  
108                    BETA = - SIGN( DLAPY2( ALPHA , XNORM ) , ALPHA )
109                    SAFMIN = DLAMCH( 'S' )
110                    RSAFMN = ONE / SAFMIN
111                    IF( ABS( BETA ).LT.SAFMIN ) THEN
112  
113  *                     XNORM , BETA may be inaccurate ; scale X and recompute them
114  
115                        KNT = 0
116     10 CONTINUE
117        KNT = KNT + 1
118        CALL PDSCAL( N - 1 , RSAFMN , X , IX , JX , DESCX , INCX )
119        BETA = BETA*RSAFMN
120        ALPHA = ALPHA*RSAFMN
121        IF( ABS( BETA ).LT.SAFMIN )
122       $    GO TO 10
123  
124  *         New BETA is at most 1 , at least SAFMIN
125  
126            CALL PDNRM2( N - 1 , XNORM , X , IX , JX , DESCX , INCX )
127            BETA = - SIGN( DLAPY2( ALPHA , XNORM ) , ALPHA )
128            TAU( INDXTAU ) =( BETA - ALPHA ) / BETA
129            CALL PDSCAL( N - 1 , ONE / (ALPHA - BETA) , X , IX , JX , DESCX , INCX )
130  
131  *         If ALPHA is subnormal , it may lose relative accuracy
132  
133            ALPHA = BETA
134            DO 20 J = 1 , KNT
135                ALPHA = ALPHA*SAFMIN
136     20     CONTINUE
137        ELSE
138            TAU( INDXTAU ) =( BETA - ALPHA ) / BETA
139            CALL PDSCAL( N - 1 , ONE / (ALPHA - BETA) , X , IX , JX , DESCX , INCX )
140            ALPHA = BETA
141        END IF
142        END IF
143  
144        RETURN
145  
146  *     End of PDLARFG
147  
148        END