Routine: PCUNGQL()  File: SRC\pcungql.f

 
 
# lines: 300
  # code: 300
  # comment: 0
  # blank:0
# Variables:43
# Callers:0
# Callings:4
# Words:122
# Keywords:68
 

 

..
     .. Array Arguments ..
     ..
  Purpose
  =======
  PCUNGQL generates an M-by-N complex distributed matrix Q denoting
  A(IA:IA+M-1,JA:JA+N-1) with orthonormal columns, which is defined as
  the last N columns of a product of K elementary reflectors of order M
        Q  =  H(k) . . . H(2) H(1)
  as returned by PCGEQLF.
  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
  =========
  M       (global input) INTEGER
          The number of rows to be operated on i.e the number of rows
          of the distributed submatrix Q. M >= 0.
  N       (global input) INTEGER
          The number of columns to be operated on i.e the number of
          columns of the distributed submatrix Q. M >= N >= 0.
  K       (global input) INTEGER
          The number of elementary reflectors whose product defines the
          matrix Q. N >= K >= 0.
  A       (local input/local output) COMPLEX pointer into the
          local memory to an array of dimension (LLD_A,LOCc(JA+N-1)).
          On entry, the j-th column must contain the vector which
          defines the elementary reflector H(j), JA+N-K <= j <= JA+N-1,
          as returned by PCGEQLF in the K columns of its distributed
          matrix argument A(IA:*,JA+N-K:JA+N-1). On exit, this array
          contains the local pieces of the M-by-N distributed matrix Q.
  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.
  TAU     (local input) COMPLEX, array, dimension LOCc(JA+N-1)
          This array contains the scalar factors TAU(j) of the
          elementary reflectors H(j) as returned by PCGEQLF.
          TAU is tied to the distributed matrix A.
  WORK    (local workspace/local output) COMPLEX 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 >= NB_A * ( NqA0 + MpA0 + NB_A ), where
          IROFFA = MOD( IA-1, MB_A ), ICOFFA = MOD( JA-1, NB_A ),
          IAROW = INDXG2P( IA, MB_A, MYROW, RSRC_A, NPROW ),
          IACOL = INDXG2P( JA, NB_A, MYCOL, CSRC_A, NPCOL ),
          MpA0 = NUMROC( M+IROFFA, MB_A, MYROW, IAROW, NPROW ),
          NqA0 = NUMROC( N+ICOFFA, NB_A, MYCOL, IACOL, NPCOL ),
          INDXG2P and NUMROC are ScaLAPACK tool functions;
          MYROW, MYCOL, NPROW and NPCOL can be determined by calling
          the subroutine BLACS_GRIDINFO.
          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.
  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 PCUNGQL( M , N , K , A , IA , JA , DESCA , TAU , WORK , LWORK ,
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 , K , LWORK , M , N
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        COMPLEX ZERO
017        PARAMETER( ZERO =( 0.0E + 0 , 0.0E + 0 ) )
018  *     ..
019  *     .. Local Scalars ..
020        LOGICAL LQUERY
021        CHARACTER COLBTOP , ROWBTOP
022        INTEGER IACOL , IAROW , ICTXT , IINFO , IPW , J , JB , JN ,
023       $LWMIN , MPA0 , MYCOL , MYROW , NPCOL , NPROW , NQA0
024  *     ..
025  *     .. Local Arrays ..
026        INTEGER IDUM1( 2 ) , IDUM2( 2 )
027  *     ..
028  *     .. External Subroutines ..
029        EXTERNAL BLACS_GRIDINFO , CHK1MAT , PCHK1MAT , PCLARFB ,
030       $PCLARFT , PCLASET , PCUNG2L , PB_TOPGET ,
031       $PB_TOPSET , PXERBLA
032  *     ..
033  *     .. External Functions ..
034        INTEGER ICEIL , INDXG2P , NUMROC
035        EXTERNAL ICEIL , INDXG2P , NUMROC
036  *     ..
037  *     .. Intrinsic Functions ..
038        INTRINSIC CMPLX , MIN , MOD , REAL
039  *     ..
040  *     .. Executable Statements ..
041  
042  *     Get grid parameters
043  
044        ICTXT = DESCA( CTXT_ )
045        CALL BLACS_GRIDINFO( ICTXT , NPROW , NPCOL , MYROW , MYCOL )
046  
047  *     Test the input parameters
048  
049        INFO = 0
050        IF( NPROW.EQ. - 1 ) THEN
051            INFO = - (700 + CTXT_)
052        ELSE
053            CALL CHK1MAT( M , 1 , N , 2 , IA , JA , DESCA , 7 , INFO )
054            IF( INFO.EQ.0 ) THEN
055                IAROW = INDXG2P( IA , DESCA( MB_ ) , MYROW , DESCA( RSRC_ ) ,
056       $        NPROW )
057                IACOL = INDXG2P( JA , DESCA( NB_ ) , MYCOL , DESCA( CSRC_ ) ,
058       $        NPCOL )
059                MPA0 = NUMROC( M + MOD( IA - 1 , DESCA( MB_ ) ) , DESCA( MB_ ) ,
060       $        MYROW , IAROW , NPROW )
061                NQA0 = NUMROC( N + MOD( JA - 1 , DESCA( NB_ ) ) , DESCA( NB_ ) ,
062       $        MYCOL , IACOL , NPCOL )
063                LWMIN = DESCA( NB_ ) * ( MPA0 + NQA0 + DESCA( NB_ ) )
064  
065                WORK( 1 ) = CMPLX( REAL( LWMIN ) )
066                LQUERY =( LWORK.EQ. - 1 )
067                IF( N.GT.M ) THEN
068                    INFO = - 2
069                ELSE IF( K.LT.0 .OR. K.GT.N ) THEN
070                    INFO = - 3
071                ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
072                    INFO = - 10
073                END IF
074            END IF
075            IDUM1( 1 ) = K
076            IDUM2( 1 ) = 3
077            IF( LWORK.EQ. - 1 ) THEN
078                IDUM1( 2 ) = - 1
079            ELSE
080                IDUM1( 2 ) = 1
081            END IF
082            IDUM2( 2 ) = 10
083            CALL PCHK1MAT( M , 1 , N , 2 , IA , JA , DESCA , 7 , 2 , IDUM1 , IDUM2 ,
084       $    INFO )
085        END IF
086        IF( INFO.NE.0 ) THEN
087            CALL PXERBLA( ICTXT , 'PCUNGQL' , - INFO )
088            RETURN
089        ELSE IF( LQUERY ) THEN
090            RETURN
091        END IF
092  
093  *     Quick return if possible
094  
095        IF( N.LE.0 )
096       $    RETURN
097  
098            IPW = DESCA( NB_ )*DESCA( NB_ ) + 1
099            JN = MIN( ICEIL( JA + N - K , DESCA( NB_ ) )*DESCA( NB_ ) , JA + N - 1 )
100            CALL PB_TOPGET( ICTXT , 'Broadcast' , 'Rowwise' , ROWBTOP )
101            CALL PB_TOPGET( ICTXT , 'Broadcast' , 'Columnwise' , COLBTOP )
102            CALL PB_TOPSET( ICTXT , 'Broadcast' , 'Rowwise' , 'I - ring' )
103            CALL PB_TOPSET( ICTXT , 'Broadcast' , 'Columnwise' , ' ' )
104  
105  *         Set A(ia + m - n + jn - ja + 1 : ia - m + 1 , ja : jn) to zero.
106  
107            CALL PCLASET ( 'All' , N - JN + JA - 1 , JN - JA + 1 , ZERO , ZERO , A ,
108       $    IA + M - N + JN - JA + 1 , JA , DESCA )
109  
110  *         Use unblocked code for the first or only block.
111  
112            CALL PCUNG2L ( M - N + JN - JA + 1 , JN - JA + 1 , JN - JA - N + K + 1 , A , IA , JA , DESCA ,
113       $    TAU , WORK , LWORK , IINFO )
114  
115  *         Use blocked code
116  
117            DO 10 J = JN + 1 , JA + N - 1 , DESCA( NB_ )
118                JB = MIN( JA + N - J , DESCA( NB_ ) )
119  
120  *             Form the triangular factor of the block reflector
121  *             H = H(j + jb - 1) . . . H(j + 1) H(j)
122  
123                CALL PCLARFT ( 'Backward' , 'Columnwise' , M - N + J + JB - JA , JB ,
124       $        A , IA , J , DESCA , TAU , WORK , WORK( IPW ) )
125  
126  *             Apply H to A(ia : ia + m - n + j + jb - ja - 1 , ja : j - 1) from the left
127  
128                CALL PCLARFB ( 'Left' , 'No transpose' , 'Backward' ,
129       $        'Columnwise' , M - N + J + JB - JA , J - JA , JB , A , IA ,
130       $        J , DESCA , WORK , A , IA , JA , DESCA , WORK( IPW ) )
131  
132  *             Apply H to rows ia : m - k + i + ib - 1 of current block
133  
134                CALL PCUNG2L ( M - N + J + JB - JA , JB , JB , A , IA , J , DESCA , TAU , WORK ,
135       $        LWORK , IINFO )
136  
137  *             Set rows ia + m - n + j + jb - ja : ia + m - 1 , j : j + jb - 1 of current block to
138  *             zero
139  
140                CALL PCLASET ( 'All' , N - J - JB + JA , JB , ZERO , ZERO , A ,
141       $        IA + M - N + J + JB - JA , J , DESCA )
142  
143     10     CONTINUE
144  
145            CALL PB_TOPGET( ICTXT , 'Broadcast' , 'Rowwise' , ROWBTOP )
146            CALL PB_TOPGET( ICTXT , 'Broadcast' , 'Columnwise' , COLBTOP )
147  
148            WORK( 1 ) = CMPLX( REAL( LWMIN ) )
149  
150            RETURN
151  
152  *         End of PCUNGQL
153  
154        END