Routine: PSLARFT()  File: SRC\pslarft.f

 
 
# lines: 538
  # code: 538
  # comment: 0
  # blank:0
# Variables:47
# Callers:12
# Callings:0
# Words:245
# Keywords:164
 

 

..
     .. Array Arguments ..
     ..
  Purpose
  =======
  PSLARFT forms the triangular factor T of a real block reflector H
  of order n, which is defined as a product of k elementary reflectors.
  If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular;
  If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular.
  If STOREV = 'C', the vector which defines the elementary reflector
  H(i) is stored in the i-th column of the distributed matrix V, and
     H  =  I - V * T * V'
  If STOREV = 'R', the vector which defines the elementary reflector
  H(i) is stored in the i-th row of the distributed matrix V, and
     H  =  I - V' * T * V
  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
  =========
  DIRECT  (global input) CHARACTER*1
          Specifies the order in which the elementary reflectors are
          multiplied to form the block reflector:
          = 'F': H = H(1) H(2) . . . H(k) (Forward)
          = 'B': H = H(k) . . . H(2) H(1) (Backward)
  STOREV  (global input) CHARACTER*1
          Specifies how the vectors which define the elementary
          reflectors are stored (see also Further Details):
          = 'C': columnwise
          = 'R': rowwise
  N       (global input) INTEGER
          The order of the block reflector H. N >= 0.
  K       (global input) INTEGER
          The order of the triangular factor T (= the number of
          elementary reflectors). 1 <= K <= MB_V (= NB_V).
  V       (input/output) REAL pointer into the local memory
          to an array of local dimension (LOCr(IV+N-1),LOCc(JV+K-1))
          if STOREV = 'C', and (LOCr(IV+K-1),LOCc(JV+N-1)) if
          STOREV = 'R'. The distributed matrix V contains the
          Householder vectors. See further details.
  IV      (global input) INTEGER
          The row index in the global array V indicating the first
          row of sub( V ).
  JV      (global input) INTEGER
          The column index in the global array V indicating the
          first column of sub( V ).
  DESCV   (global and local input) INTEGER array of dimension DLEN_.
          The array descriptor for the distributed matrix V.
  TAU     (local input) REAL, array, dimension LOCr(IV+K-1)
          if INCV = M_V, and LOCc(JV+K-1) otherwise. This array
          contains the Householder scalars related to the Householder
          vectors.  TAU is tied to the distributed matrix V.
  T       (local output) REAL array, dimension (NB_V,NB_V)
          if STOREV = 'Col', and (MB_V,MB_V) otherwise. It contains
          the k-by-k triangular factor of the block reflector asso-
          ciated with V. If DIRECT = 'F', T is upper triangular;
          if DIRECT = 'B', T is lower triangular.
  WORK    (local workspace) REAL array,
                                          dimension (K*(K-1)/2)
  Further Details
  ===============
  The shape of the matrix V and the storage of the vectors which define
  the H(i) is best illustrated by the following example with n = 5 and
  k = 3. The elements equal to 1 are not stored; the corresponding
  array elements are modified but restored on exit. The rest of the
  array is not used.
  DIRECT = 'F' and STOREV = 'C':   DIRECT = 'F' and STOREV = 'R':
  V( IV:IV+N-1,    (  1       )    V( IV:IV+K-1,    (  1 v1 v1 v1 v1 )
     JV:JV+K-1 ) = ( v1  1    )       JV:JV+N-1 ) = (     1 v2 v2 v2 )
                   ( v1 v2  1 )                     (        1 v3 v3 )
                   ( v1 v2 v3 )
                   ( v1 v2 v3 )
  DIRECT = 'B' and STOREV = 'C':   DIRECT = 'B' and STOREV = 'R':
  V( IV:IV+N-1,    ( v1 v2 v3 )    V( IV:IV+K-1,    ( v1 v1  1       )
     JV:JV+K-1 ) = ( v1 v2 v3 )       JV:JV+N-1 ) = ( v2 v2 v2  1    )
                   (  1 v2 v3 )                     ( v3 v3 v3 v3  1 )
                   (     1 v3 )
                   (        1 )
  =====================================================================
     .. Parameters ..

 
Display dynamic version Find AutoScroll Reload FontSize: - + Hide Comments Hide Blanks Frame FullScreen MailPrint

 
001        SUBROUTINE PSLARFT( DIRECT , STOREV , N , K , V , IV , JV , DESCV , TAU ,
002       $T , WORK )
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        CHARACTER DIRECT , STOREV
011        INTEGER IV , JV , K , N
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        REAL ONE , ZERO
018        PARAMETER( ONE = 1.0E + 0 , ZERO = 0.0E + 0 )
019  *     ..
020  *     .. Local Scalars ..
021        LOGICAL FORWARD
022        INTEGER ICOFF , ICTXT , II , IIV , IROFF , IVCOL , IVROW ,
023       $ITMP0 , ITMP1 , IW , JJ , JJV , LDV , MICOL , MIROW ,
024       $MYCOL , MYROW , NP , NPCOL , NPROW , NQ
025        REAL VII
026  *     ..
027  *     .. External Subroutines ..
028        EXTERNAL BLACS_GRIDINFO , INFOG2L , SCOPY , SGEMV ,
029       $SGSUM2D , SLASET , STRMV
030  *     ..
031  *     .. External Functions ..
032        LOGICAL LSAME
033        INTEGER INDXG2P , NUMROC
034        EXTERNAL INDXG2P , LSAME , NUMROC
035  *     ..
036  *     .. Intrinsic Functions ..
037        INTRINSIC MOD
038  *     ..
039  *     .. Executable Statements ..
040  
041  *     Quick return if possible
042  
043        IF( N.LE.0 .OR. K.LE.0 )
044       $    RETURN
045  
046            ICTXT = DESCV( CTXT_ )
047            CALL BLACS_GRIDINFO( ICTXT , NPROW , NPCOL , MYROW , MYCOL )
048  
049            FORWARD = LSAME( DIRECT , 'F' )
050            CALL INFOG2L( IV , JV , DESCV , NPROW , NPCOL , MYROW , MYCOL ,
051       $    IIV , JJV , IVROW , IVCOL )
052  
053            IF( LSAME( STOREV , 'C' ) .AND. MYCOL.EQ.IVCOL ) THEN
054  
055                IW = 1
056                LDV = DESCV( LLD_ )
057                IROFF = MOD( IV - 1 , DESCV( MB_ ) )
058  
059                IF( FORWARD ) THEN
060  
061  *                 DIRECT = 'Forward' , STOREV = 'Columnwise'
062  
063                    NP = NUMROC( N + IROFF , DESCV( MB_ ) , MYROW , IVROW , NPROW )
064                    IF( MYROW.EQ.IVROW ) THEN
065                        NP = NP - IROFF
066                        II = IIV + 1
067                    ELSE
068                        II = IIV
069                    END IF
070                    IF( IROFF + 1.EQ.DESCV( MB_ ) ) THEN
071                        MIROW = MOD( IVROW + 1 , NPROW )
072                    ELSE
073                        MIROW = IVROW
074                    END IF
075                    ITMP0 = 0
076  
077                    DO 10 JJ = JJV + 1 , JJV + K - 1
078  
079                        IF( MYROW.EQ.MIROW ) THEN
080                            VII = V( II + (JJ - 1)*LDV )
081                            V( II + (JJ - 1)*LDV ) = ONE
082                        END IF
083  
084  *                     T(1 : i - 1 , i) = - tau( jv + i - 1 ) *
085  *                     V(iv + i - 1 : iv + n - 1 , jv : jv + i - 2)' * V(iv + i - 1 : iv + n - 1 , jv + i - 1)
086  
087                        ITMP0 = ITMP0 + 1
088                        IF( NP - II + IIV.GT.0 ) THEN
089                            CALL SGEMV( 'Transpose' , NP - II + IIV , ITMP0 ,
090       $                    - TAU( JJ ) , V( II + (JJV - 1)*LDV ) , LDV ,
091       $                    V( II + (JJ - 1)*LDV ) , 1 , ZERO ,
092       $                    WORK( IW ) , 1 )
093                        ELSE
094                            CALL SLASET( 'All' , ITMP0 , 1 , ZERO , ZERO , WORK( IW ) ,
095       $                    ITMP0 )
096                        END IF
097  
098                        IW = IW + ITMP0
099                        IF( MYROW.EQ.MIROW ) THEN
100                            V( II + (JJ - 1)*LDV ) = VII
101                            II = II + 1
102                        END IF
103  
104                        IF( MOD( IV + ITMP0 , DESCV( MB_ ) ).EQ.0 )
105       $                    MIROW = MOD( MIROW + 1 , NPROW )
106  
107     10             CONTINUE
108  
109                    CALL SGSUM2D( ICTXT , 'Columnwise' , ' ' , IW - 1 , 1 , WORK , IW - 1 ,
110       $            IVROW , MYCOL )
111  
112                    IF( MYROW.EQ.IVROW ) THEN
113  
114                        IW = 1
115                        ITMP0 = 0
116                        ITMP1 = 1
117  
118                        T( ITMP1 ) = TAU( JJV )
119  
120                        DO 20 JJ = JJV + 1 , JJV + K - 1
121  
122  *                         T(1 : j - 1 , j) = T(1 : j - 1 , 1 : j - 1) * T(1 : j - 1 , j)
123  
124                            ITMP0 = ITMP0 + 1
125                            ITMP1 = ITMP1 + DESCV( NB_ )
126                            CALL SCOPY( ITMP0 , WORK( IW ) , 1 , T( ITMP1 ) , 1 )
127                            IW = IW + ITMP0
128  
129                            CALL STRMV( 'Upper' , 'No transpose' , 'Non - unit' ,
130       $                    ITMP0 , T , DESCV( NB_ ) , T( ITMP1 ) , 1 )
131                            T(ITMP1 + ITMP0) = TAU( JJ )
132  
133     20                 CONTINUE
134  
135                    END IF
136  
137                ELSE
138  
139  *                 DIRECT = 'Backward' , STOREV = 'Columnwise'
140  
141                    NP = NUMROC( N + IROFF - 1 , DESCV( MB_ ) , MYROW , IVROW , NPROW )
142                    IF( MYROW.EQ.IVROW )
143       $                NP = NP - IROFF
144                        MIROW = INDXG2P( IV + N - 2 , DESCV( MB_ ) , MYROW ,
145       $                DESCV( RSRC_ ) , NPROW )
146                        II = IIV + NP - 1
147                        ITMP0 = 0
148  
149                        DO 30 JJ = JJV + K - 2 , JJV , - 1
150  
151                            IF( MYROW.EQ.MIROW ) THEN
152                                VII = V( II + (JJ - 1)*LDV )
153                                V( II + (JJ - 1)*LDV ) = ONE
154                            END IF
155  
156  *                         T(1 : i - 1 , i) = - tau( jv + i - 1 ) *
157  *                         V(iv : iv + n - k + i - 1 , jv + i : jv + k - 1)' * V(iv : iv + n - k + i - 1 , jv + i - 1)
158  
159                            ITMP0 = ITMP0 + 1
160                            IF( II - IIV + 1.GT.0 ) THEN
161                                CALL SGEMV( 'Transpose' , II - IIV + 1 , ITMP0 , - TAU( JJ ) ,
162       $                        V( IIV + JJ*LDV ) , LDV ,
163       $                        V( IIV + (JJ - 1)*LDV ) , 1 , ZERO ,
164       $                        WORK( IW ) , 1 )
165                            ELSE
166                                CALL SLASET( 'All' , ITMP0 , 1 , ZERO , ZERO , WORK( IW ) ,
167       $                        ITMP0 )
168                            END IF
169  
170                            IW = IW + ITMP0
171                            IF( MYROW.EQ.MIROW ) THEN
172                                V( II + (JJ - 1)*LDV ) = VII
173                                II = II - 1
174                            END IF
175  
176                            IF( MOD( IV + N - ITMP0 - 2 , DESCV(MB_) ).EQ.0 )
177       $                        MIROW = MOD( MIROW + NPROW - 1 , NPROW )
178  
179     30                 CONTINUE
180  
181                        CALL SGSUM2D( ICTXT , 'Columnwise' , ' ' , IW - 1 , 1 , WORK , IW - 1 ,
182       $                IVROW , MYCOL )
183  
184                        IF( MYROW.EQ.IVROW ) THEN
185  
186                            IW = 1
187                            ITMP0 = 0
188                            ITMP1 = K + 1 + (K - 1) * DESCV( NB_ )
189  
190                            T( ITMP1 - 1 ) = TAU( JJV + K - 1 )
191  
192                            DO 40 JJ = JJV + K - 2 , JJV , - 1
193  
194  *                             T(j + 1 : k , j) = T(j + 1 : k , j + 1 : k) * T(j + 1 : k , j)
195  
196                                ITMP0 = ITMP0 + 1
197                                ITMP1 = ITMP1 - DESCV( NB_ ) - 1
198                                CALL SCOPY( ITMP0 , WORK( IW ) , 1 , T( ITMP1 ) , 1 )
199                                IW = IW + ITMP0
200  
201                                CALL STRMV( 'Lower' , 'No transpose' , 'Non - unit' ,
202       $                        ITMP0 , T( ITMP1 + DESCV( NB_ ) ) ,
203       $                        DESCV( NB_ ) , T( ITMP1 ) , 1 )
204                                T( ITMP1 - 1 ) = TAU( JJ )
205  
206     40                     CONTINUE
207  
208                        END IF
209  
210                    END IF
211  
212                ELSE IF( LSAME( STOREV , 'R' ) .AND. MYROW.EQ.IVROW ) THEN
213  
214                    IW = 1
215                    LDV = DESCV( LLD_ )
216                    ICOFF = MOD( JV - 1 , DESCV( NB_ ) )
217  
218                    IF( FORWARD ) THEN
219  
220  *                     DIRECT = 'Forward' , STOREV = 'Rowwise'
221  
222                        NQ = NUMROC( N + ICOFF , DESCV( NB_ ) , MYCOL , IVCOL , NPCOL )
223                        IF( MYCOL.EQ.IVCOL ) THEN
224                            NQ = NQ - ICOFF
225                            JJ = JJV + 1
226                        ELSE
227                            JJ = JJV
228                        END IF
229                        IF( ICOFF + 1.EQ.DESCV( NB_ ) ) THEN
230                            MICOL = MOD( IVCOL + 1 , NPCOL )
231                        ELSE
232                            MICOL = IVCOL
233                        END IF
234                        ITMP0 = 0
235  
236                        DO 50 II = IIV + 1 , IIV + K - 1
237  
238                            IF( MYCOL.EQ.MICOL ) THEN
239                                VII = V( II + (JJ - 1)*LDV )
240                                V( II + (JJ - 1)*LDV ) = ONE
241                            END IF
242  
243  *                         T(1 : i - 1 , i) = - tau( iv + i - 1 ) *
244  *                         V(iv + i - 1 , jv + i - 1 : jv + n - 1) * V(iv : iv + i - 2 , jv + i - 1 : jv + n - 1)'
245  
246                            ITMP0 = ITMP0 + 1
247                            IF( NQ - JJ + JJV.GT.0 ) THEN
248                                CALL SGEMV( 'No transpose' , ITMP0 , NQ - JJ + JJV ,
249       $                        - TAU(II) , V( IIV + (JJ - 1)*LDV ) , LDV ,
250       $                        V( II + (JJ - 1)*LDV ) , LDV , ZERO ,
251       $                        WORK( IW ) , 1 )
252                            ELSE
253                                CALL SLASET( 'All' , ITMP0 , 1 , ZERO , ZERO ,
254       $                        WORK( IW ) , ITMP0 )
255                            END IF
256  
257                            IW = IW + ITMP0
258                            IF( MYCOL.EQ.MICOL ) THEN
259                                V( II + (JJ - 1)*LDV ) = VII
260                                JJ = JJ + 1
261                            END IF
262  
263                            IF( MOD( JV + ITMP0 , DESCV( NB_ ) ).EQ.0 )
264       $                        MICOL = MOD( MICOL + 1 , NPCOL )
265  
266     50                 CONTINUE
267  
268                        CALL SGSUM2D( ICTXT , 'Rowwise' , ' ' , IW - 1 , 1 , WORK , IW - 1 ,
269       $                MYROW , IVCOL )
270  
271                        IF( MYCOL.EQ.IVCOL ) THEN
272  
273                            IW = 1
274                            ITMP0 = 0
275                            ITMP1 = 1
276  
277                            T( ITMP1 ) = TAU( IIV )
278  
279                            DO 60 II = IIV + 1 , IIV + K - 1
280  
281  *                             T(1 : i - 1 , i) = T(1 : i - 1 , 1 : i - 1) * T(1 : i - 1 , i)
282  
283                                ITMP0 = ITMP0 + 1
284                                ITMP1 = ITMP1 + DESCV( MB_ )
285                                CALL SCOPY( ITMP0 , WORK( IW ) , 1 , T( ITMP1 ) , 1 )
286                                IW = IW + ITMP0
287  
288                                CALL STRMV( 'Upper' , 'No transpose' , 'Non - unit' ,
289       $                        ITMP0 , T , DESCV( MB_ ) , T( ITMP1 ) , 1 )
290                                T( ITMP1 + ITMP0 ) = TAU( II )
291  
292     60                     CONTINUE
293  
294                        END IF
295  
296                    ELSE
297  
298  *                     DIRECT = 'Backward' , STOREV = 'Rowwise'
299  
300                        NQ = NUMROC( N + ICOFF - 1 , DESCV( NB_ ) , MYCOL , IVCOL , NPCOL )
301                        IF( MYCOL.EQ.IVCOL )
302       $                    NQ = NQ - ICOFF
303                            MICOL = INDXG2P( JV + N - 2 , DESCV( NB_ ) , MYCOL ,
304       $                    DESCV( CSRC_ ) , NPCOL )
305                            JJ = JJV + NQ - 1
306                            ITMP0 = 0
307  
308                            DO 70 II = IIV + K - 2 , IIV , - 1
309  
310                                IF( MYCOL.EQ.MICOL ) THEN
311                                    VII = V( II + (JJ - 1)*LDV )
312                                    V( II + (JJ - 1)*LDV ) = ONE
313                                END IF
314  
315  *                             T(i + 1 : k , i) = - tau( iv + i - 1 ) *
316  *                             V(iv + i : iv + k - 1 , jv : jv + n - k + i - 1)' * V(iv + i - 1 , jv : jv + n - k + i - 1)'
317  
318                                ITMP0 = ITMP0 + 1
319                                IF( JJ - JJV + 1.GT.0 ) THEN
320                                    CALL SGEMV( 'No transpose' , ITMP0 , JJ - JJV + 1 ,
321       $                            - TAU( II ) , V( II + 1 + (JJV - 1)*LDV ) , LDV ,
322       $                            V( II + (JJV - 1)*LDV ) , LDV , ZERO ,
323       $                            WORK( IW ) , 1 )
324                                ELSE
325                                    CALL SLASET( 'All' , ITMP0 , 1 , ZERO , ZERO ,
326       $                            WORK( IW ) , ITMP0 )
327                                END IF
328  
329                                IW = IW + ITMP0
330                                IF( MYCOL.EQ.MICOL ) THEN
331                                    V( II + (JJ - 1)*LDV ) = VII
332                                    JJ = JJ - 1
333                                END IF
334  
335                                IF( MOD( JV + N - ITMP0 - 2 , DESCV( NB_ ) ).EQ.0 )
336       $                            MICOL = MOD( MICOL + NPCOL - 1 , NPCOL )
337  
338     70                     CONTINUE
339  
340                            CALL SGSUM2D( ICTXT , 'Rowwise' , ' ' , IW - 1 , 1 , WORK , IW - 1 ,
341       $                    MYROW , IVCOL )
342  
343                            IF( MYCOL.EQ.IVCOL ) THEN
344  
345                                IW = 1
346                                ITMP0 = 0
347                                ITMP1 = K + 1 + (K - 1) * DESCV( MB_ )
348  
349                                T( ITMP1 - 1 ) = TAU( IIV + K - 1 )
350  
351                                DO 80 II = IIV + K - 2 , IIV , - 1
352  
353  *                                 T(i + 1 : k , i) = T(i + 1 : k , i + 1 : k) * T(i + 1 : k , i)
354  
355                                    ITMP0 = ITMP0 + 1
356                                    ITMP1 = ITMP1 - DESCV( MB_ ) - 1
357                                    CALL SCOPY( ITMP0 , WORK( IW ) , 1 , T( ITMP1 ) , 1 )
358                                    IW = IW + ITMP0
359  
360                                    CALL STRMV( 'Lower' , 'No transpose' , 'Non - unit' ,
361       $                            ITMP0 , T( ITMP1 + DESCV( MB_ ) ) ,
362       $                            DESCV( MB_ ) , T( ITMP1 ) , 1 )
363                                    T( ITMP1 - 1 ) = TAU( II )
364  
365     80                         CONTINUE
366  
367                            END IF
368  
369                        END IF
370  
371                    END IF
372  
373                    RETURN
374  
375  *                 End of PSLARFT
376  
377                END