ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pclabrd.f
Go to the documentation of this file.
1  SUBROUTINE pclabrd( M, N, NB, A, IA, JA, DESCA, D, E, TAUQ, TAUP,
2  $ X, IX, JX, DESCX, Y, IY, JY, DESCY, WORK )
3 *
4 * -- ScaLAPACK auxiliary routine (version 1.7) --
5 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
6 * and University of California, Berkeley.
7 * May 1, 1997
8 *
9 * .. Scalar Arguments ..
10  INTEGER IA, IX, IY, JA, JX, JY, M, N, NB
11 * ..
12 * .. Array Arguments ..
13  INTEGER DESCA( * ), DESCX( * ), DESCY( * )
14  REAL D( * ), E( * )
15  COMPLEX A( * ), TAUP( * ), TAUQ( * ), X( * ), Y( * ),
16  $ work( * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * PCLABRD reduces the first NB rows and columns of a complex general
23 * M-by-N distributed matrix sub( A ) = A(IA:IA+M-1,JA:JA+N-1) to upper
24 * or lower bidiagonal form by an unitary transformation Q' * A * P, and
25 * returns the matrices X and Y which are needed to apply the transfor-
26 * mation to the unreduced part of sub( A ).
27 *
28 * If M >= N, sub( A ) is reduced to upper bidiagonal form; if M < N, to
29 * lower bidiagonal form.
30 *
31 * This is an auxiliary routine called by PCGEBRD.
32 *
33 * Notes
34 * =====
35 *
36 * Each global data object is described by an associated description
37 * vector. This vector stores the information required to establish
38 * the mapping between an object element and its corresponding process
39 * and memory location.
40 *
41 * Let A be a generic term for any 2D block cyclicly distributed array.
42 * Such a global array has an associated description vector DESCA.
43 * In the following comments, the character _ should be read as
44 * "of the global array".
45 *
46 * NOTATION STORED IN EXPLANATION
47 * --------------- -------------- --------------------------------------
48 * DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
49 * DTYPE_A = 1.
50 * CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
51 * the BLACS process grid A is distribu-
52 * ted over. The context itself is glo-
53 * bal, but the handle (the integer
54 * value) may vary.
55 * M_A (global) DESCA( M_ ) The number of rows in the global
56 * array A.
57 * N_A (global) DESCA( N_ ) The number of columns in the global
58 * array A.
59 * MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
60 * the rows of the array.
61 * NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
62 * the columns of the array.
63 * RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
64 * row of the array A is distributed.
65 * CSRC_A (global) DESCA( CSRC_ ) The process column over which the
66 * first column of the array A is
67 * distributed.
68 * LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
69 * array. LLD_A >= MAX(1,LOCr(M_A)).
70 *
71 * Let K be the number of rows or columns of a distributed matrix,
72 * and assume that its process grid has dimension p x q.
73 * LOCr( K ) denotes the number of elements of K that a process
74 * would receive if K were distributed over the p processes of its
75 * process column.
76 * Similarly, LOCc( K ) denotes the number of elements of K that a
77 * process would receive if K were distributed over the q processes of
78 * its process row.
79 * The values of LOCr() and LOCc() may be determined via a call to the
80 * ScaLAPACK tool function, NUMROC:
81 * LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
82 * LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
83 * An upper bound for these quantities may be computed by:
84 * LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
85 * LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
86 *
87 * Arguments
88 * =========
89 *
90 * M (global input) INTEGER
91 * The number of rows to be operated on, i.e. the number of rows
92 * of the distributed submatrix sub( A ). M >= 0.
93 *
94 * N (global input) INTEGER
95 * The number of columns to be operated on, i.e. the number of
96 * columns of the distributed submatrix sub( A ). N >= 0.
97 *
98 * NB (global input) INTEGER
99 * The number of leading rows and columns of sub( A ) to be
100 * reduced.
101 *
102 * A (local input/local output) COMPLEX pointer into the
103 * local memory to an array of dimension (LLD_A,LOCc(JA+N-1)).
104 * On entry, this array contains the local pieces of the
105 * general distributed matrix sub( A ) to be reduced. On exit,
106 * the first NB rows and columns of the matrix are overwritten;
107 * the rest of the distributed matrix sub( A ) is unchanged.
108 * If m >= n, elements on and below the diagonal in the first NB
109 * columns, with the array TAUQ, represent the unitary
110 * matrix Q as a product of elementary reflectors; and
111 * elements above the diagonal in the first NB rows, with the
112 * array TAUP, represent the unitary matrix P as a product
113 * of elementary reflectors.
114 * If m < n, elements below the diagonal in the first NB
115 * columns, with the array TAUQ, represent the unitary
116 * matrix Q as a product of elementary reflectors, and
117 * elements on and above the diagonal in the first NB rows,
118 * with the array TAUP, represent the unitary matrix P as
119 * a product of elementary reflectors.
120 * See Further Details.
121 *
122 * IA (global input) INTEGER
123 * The row index in the global array A indicating the first
124 * row of sub( A ).
125 *
126 * JA (global input) INTEGER
127 * The column index in the global array A indicating the
128 * first column of sub( A ).
129 *
130 * DESCA (global and local input) INTEGER array of dimension DLEN_.
131 * The array descriptor for the distributed matrix A.
132 *
133 * D (local output) REAL array, dimension
134 * LOCr(IA+MIN(M,N)-1) if M >= N; LOCc(JA+MIN(M,N)-1) otherwise.
135 * The distributed diagonal elements of the bidiagonal matrix
136 * B: D(i) = A(ia+i-1,ja+i-1). D is tied to the distributed
137 * matrix A.
138 *
139 * E (local output) REAL array, dimension
140 * LOCr(IA+MIN(M,N)-1) if M >= N; LOCc(JA+MIN(M,N)-2) otherwise.
141 * The distributed off-diagonal elements of the bidiagonal
142 * distributed matrix B:
143 * if m >= n, E(i) = A(ia+i-1,ja+i) for i = 1,2,...,n-1;
144 * if m < n, E(i) = A(ia+i,ja+i-1) for i = 1,2,...,m-1.
145 * E is tied to the distributed matrix A.
146 *
147 * TAUQ (local output) COMPLEX array dimension
148 * LOCc(JA+MIN(M,N)-1). The scalar factors of the elementary
149 * reflectors which represent the unitary matrix Q. TAUQ is
150 * tied to the distributed matrix A. See Further Details.
151 *
152 * TAUP (local output) COMPLEX array, dimension
153 * LOCr(IA+MIN(M,N)-1). The scalar factors of the elementary
154 * reflectors which represent the unitary matrix P. TAUP is
155 * tied to the distributed matrix A. See Further Details.
156 *
157 * X (local output) COMPLEX pointer into the local memory
158 * to an array of dimension (LLD_X,NB). On exit, the local
159 * pieces of the distributed M-by-NB matrix
160 * X(IX:IX+M-1,JX:JX+NB-1) required to update the unreduced
161 * part of sub( A ).
162 *
163 * IX (global input) INTEGER
164 * The row index in the global array X indicating the first
165 * row of sub( X ).
166 *
167 * JX (global input) INTEGER
168 * The column index in the global array X indicating the
169 * first column of sub( X ).
170 *
171 * DESCX (global and local input) INTEGER array of dimension DLEN_.
172 * The array descriptor for the distributed matrix X.
173 *
174 * Y (local output) COMPLEX pointer into the local memory
175 * to an array of dimension (LLD_Y,NB). On exit, the local
176 * pieces of the distributed N-by-NB matrix
177 * Y(IY:IY+N-1,JY:JY+NB-1) required to update the unreduced
178 * part of sub( A ).
179 *
180 * IY (global input) INTEGER
181 * The row index in the global array Y indicating the first
182 * row of sub( Y ).
183 *
184 * JY (global input) INTEGER
185 * The column index in the global array Y indicating the
186 * first column of sub( Y ).
187 *
188 * DESCY (global and local input) INTEGER array of dimension DLEN_.
189 * The array descriptor for the distributed matrix Y.
190 *
191 * WORK (local workspace) COMPLEX array, dimension (LWORK)
192 * LWORK >= NB_A + NQ, with
193 *
194 * NQ = NUMROC( N+MOD( IA-1, NB_Y ), NB_Y, MYCOL, IACOL, NPCOL )
195 * IACOL = INDXG2P( JA, NB_A, MYCOL, CSRC_A, NPCOL )
196 *
197 * INDXG2P and NUMROC are ScaLAPACK tool functions;
198 * MYROW, MYCOL, NPROW and NPCOL can be determined by calling
199 * the subroutine BLACS_GRIDINFO.
200 *
201 * Further Details
202 * ===============
203 *
204 * The matrices Q and P are represented as products of elementary
205 * reflectors:
206 *
207 * Q = H(1) H(2) . . . H(nb) and P = G(1) G(2) . . . G(nb)
208 *
209 * Each H(i) and G(i) has the form:
210 *
211 * H(i) = I - tauq * v * v' and G(i) = I - taup * u * u'
212 *
213 * where tauq and taup are complex scalars, and v and u are complex
214 * vectors.
215 *
216 * If m >= n, v(1:i-1) = 0, v(i) = 1, and v(i:m) is stored on exit in
217 * A(ia+i-1:ia+m-1,ja+i-1); u(1:i) = 0, u(i+1) = 1, and u(i+1:n) is
218 * stored on exit in A(ia+i-1,ja+i:ja+n-1); tauq is stored in
219 * TAUQ(ja+i-1) and taup in TAUP(ia+i-1).
220 *
221 * If m < n, v(1:i) = 0, v(i+1) = 1, and v(i+1:m) is stored on exit in
222 * A(ia+i+1:ia+m-1,ja+i-1); u(1:i-1) = 0, u(i) = 1, and u(i:n) is
223 * stored on exit in A(ia+i-1,ja+i:ja+n-1); tauq is stored in
224 * TAUQ(ja+i-1) and taup in TAUP(ia+i-1).
225 *
226 * The elements of the vectors v and u together form the m-by-nb matrix
227 * V and the nb-by-n matrix U' which are needed, with X and Y, to apply
228 * the transformation to the unreduced part of the matrix, using a block
229 * update of the form: sub( A ) := sub( A ) - V*Y' - X*U'.
230 *
231 * The contents of sub( A ) on exit are illustrated by the following
232 * examples with nb = 2:
233 *
234 * m = 6 and n = 5 (m > n): m = 5 and n = 6 (m < n):
235 *
236 * ( 1 1 u1 u1 u1 ) ( 1 u1 u1 u1 u1 u1 )
237 * ( v1 1 1 u2 u2 ) ( 1 1 u2 u2 u2 u2 )
238 * ( v1 v2 a a a ) ( v1 1 a a a a )
239 * ( v1 v2 a a a ) ( v1 v2 a a a a )
240 * ( v1 v2 a a a ) ( v1 v2 a a a a )
241 * ( v1 v2 a a a )
242 *
243 * where a denotes an element of the original matrix which is unchanged,
244 * vi denotes an element of the vector defining H(i), and ui an element
245 * of the vector defining G(i).
246 *
247 * =====================================================================
248 *
249 * .. Parameters ..
250  INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
251  $ lld_, mb_, m_, nb_, n_, rsrc_
252  parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
253  $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
254  $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
255  COMPLEX ONE, ZERO
256  parameter( one = ( 1.0e+0, 0.0e+0 ),
257  $ zero = ( 0.0e+0, 0.0e+0 ) )
258 * ..
259 * .. Local Scalars ..
260  INTEGER I, IACOL, IAROW, ICTXT, II, IPY, IW, J, JJ,
261  $ jwy, k, mycol, myrow, npcol, nprow
262  COMPLEX ALPHA, TAU
263  INTEGER DESCD( DLEN_ ), DESCE( DLEN_ ),
264  $ desctp( dlen_ ), desctq( dlen_ ),
265  $ descw( dlen_ ), descwy( dlen_ )
266 * ..
267 * .. External Subroutines ..
268  EXTERNAL blacs_gridinfo, descset, infog2l, pccopy,
269  $ pcelget, pcelset, pcgemv, pclacgv,
270  $ pclarfg, pcscal, pselset
271 * ..
272 * .. Intrinsic Functions ..
273  INTRINSIC cmplx, min, mod
274 * ..
275 * .. Executable Statements ..
276 *
277 * Quick return if possible
278 *
279  IF( m.LE.0 .OR. n.LE.0 )
280  $ RETURN
281 *
282  ictxt = desca( ctxt_ )
283  CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
284  CALL infog2l( ia, ja, desca, nprow, npcol, myrow, mycol, ii, jj,
285  $ iarow, iacol )
286  ipy = desca( mb_ ) + 1
287  iw = mod( ia-1, desca( nb_ ) ) + 1
288  alpha = zero
289 *
290  CALL descset( descwy, 1, n+mod( ia-1, descy( nb_ ) ), 1,
291  $ desca( nb_ ), iarow, iacol, ictxt, 1 )
292  CALL descset( descw, desca( mb_ ), 1, desca( mb_ ), 1, iarow,
293  $ iacol, ictxt, desca( mb_ ) )
294  CALL descset( desctq, 1, ja+min(m,n)-1, 1, desca( nb_ ), iarow,
295  $ desca( csrc_ ), desca( ctxt_ ), 1 )
296  CALL descset( desctp, ia+min(m,n)-1, 1, desca( mb_ ), 1,
297  $ desca( rsrc_ ), iacol, desca( ctxt_ ),
298  $ desca( lld_ ) )
299 *
300  IF( m.GE.n ) THEN
301 *
302 * Reduce to upper bidiagonal form
303 *
304  CALL descset( descd, 1, ja+min(m,n)-1, 1, desca( nb_ ), myrow,
305  $ desca( csrc_ ), desca( ctxt_ ), 1 )
306  CALL descset( desce, ia+min(m,n)-1, 1, desca( mb_ ), 1,
307  $ desca( rsrc_ ), mycol, desca( ctxt_ ),
308  $ desca( lld_ ) )
309  DO 10 k = 1, nb
310  i = ia + k - 1
311  j = ja + k - 1
312  jwy = iw + k
313 *
314 * Update A(i:ia+m-1,j)
315 *
316  IF( k.GT.1 ) THEN
317  CALL pcgemv( 'No transpose', m-k+1, k-1, -one, a, i, ja,
318  $ desca, y, iy, jy+k-1, descy, 1, one, a, i,
319  $ j, desca, 1 )
320  CALL pcgemv( 'No transpose', m-k+1, k-1, -one, x, ix+k-1,
321  $ jx, descx, a, ia, j, desca, 1, one, a, i, j,
322  $ desca, 1 )
323  CALL pcelset( a, i-1, j, desca, alpha )
324  END IF
325 *
326 * Generate reflection Q(i) to annihilate A(i+1:ia+m-1,j)
327 *
328  CALL pclarfg( m-k+1, alpha, i, j, a, i+1, j, desca, 1,
329  $ tauq )
330  CALL pselset( d, 1, j, descd, real( alpha ) )
331  CALL pcelset( a, i, j, desca, one )
332 *
333 * Compute Y(IA+I:IA+N-1,J)
334 *
335  CALL pcgemv( 'Conjugate transpose', m-k+1, n-k, one, a, i,
336  $ j+1, desca, a, i, j, desca, 1, zero,
337  $ work( ipy ), 1, jwy, descwy, descwy( m_ ) )
338  CALL pcgemv( 'Conjugate transpose', m-k+1, k-1, one, a, i,
339  $ ja, desca, a, i, j, desca, 1, zero, work, iw,
340  $ 1, descw, 1 )
341  CALL pcgemv( 'Conjugate transpose', k-1, n-k, -one, y, iy,
342  $ jy+k, descy, work, iw, 1, descw, 1, one,
343  $ work( ipy ), 1, jwy, descwy, descwy( m_ ) )
344  CALL pcgemv( 'Conjugate transpose', m-k+1, k-1, one, x,
345  $ ix+k-1, jx, descx, a, i, j, desca, 1, zero,
346  $ work, iw, 1, descw, 1 )
347  CALL pcgemv( 'Conjugate transpose', k-1, n-k, -one, a, ia,
348  $ j+1, desca, work, iw, 1, descw, 1, one,
349  $ work( ipy ), 1, jwy, descwy, descwy( m_ ) )
350 *
351  CALL pcelget( 'Rowwise', ' ', tau, tauq, 1, j, desctq )
352  CALL pcscal( n-k, tau, work( ipy ), 1, jwy, descwy,
353  $ descwy( m_ ) )
354  CALL pclacgv( n-k, work( ipy ), 1, jwy, descwy,
355  $ descwy( m_ ) )
356  CALL pccopy( n-k, work( ipy ), 1, jwy, descwy, descwy( m_ ),
357  $ y, iy+k-1, jy+k, descy, descy( m_ ) )
358 *
359 * Update A(i,j+1:ja+n-1)
360 *
361  CALL pclacgv( n-k, a, i, j+1, desca, desca( m_ ) )
362  CALL pclacgv( k, a, i, ja, desca, desca( m_ ) )
363  CALL pcgemv( 'Conjugate transpose', k, n-k, -one, y, iy,
364  $ jy+k, descy, a, i, ja, desca, desca( m_ ), one,
365  $ a, i, j+1, desca, desca( m_ ) )
366  CALL pclacgv( k, a, i, ja, desca, desca( m_ ) )
367  CALL pclacgv( k-1, x, ix+k-1, jx, descx, descx( m_ ) )
368  CALL pcgemv( 'Conjugate transpose', k-1, n-k, -one, a, ia,
369  $ j+1, desca, x, ix+k-1, jx, descx, descx( m_ ),
370  $ one, a, i, j+1, desca, desca( m_ ) )
371  CALL pclacgv( k-1, x, ix+k-1, jx, descx, descx( m_ ) )
372  CALL pcelset( a, i, j, desca, cmplx( real( alpha ) ) )
373 *
374 * Generate reflection P(i) to annihilate A(i,j+2:ja+n-1)
375 *
376  CALL pclarfg( n-k, alpha, i, j+1, a, i,
377  $ min( j+2, n+ja-1 ), desca, desca( m_ ), taup )
378  CALL pselset( e, i, 1, desce, real( alpha ) )
379  CALL pcelset( a, i, j+1, desca, one )
380 *
381 * Compute X(I+1:IA+M-1,J)
382 *
383  CALL pcgemv( 'No transpose', m-k, n-k, one, a, i+1, j+1,
384  $ desca, a, i, j+1, desca, desca( m_ ), zero, x,
385  $ ix+k, jx+k-1, descx, 1 )
386  CALL pcgemv( 'No transpose', k, n-k, one, y, iy, jy+k,
387  $ descy, a, i, j+1, desca, desca( m_ ), zero,
388  $ work, iw, 1, descw, 1 )
389  CALL pcgemv( 'No transpose', m-k, k, -one, a, i+1, ja,
390  $ desca, work, iw, 1, descw, 1, one, x, ix+k,
391  $ jx+k-1, descx, 1 )
392  CALL pcgemv( 'No transpose', k-1, n-k, one, a, ia, j+1,
393  $ desca, a, i, j+1, desca, desca( m_ ), zero,
394  $ work, iw, 1, descw, 1 )
395  CALL pcgemv( 'No transpose', m-k, k-1, -one, x, ix+k, jx,
396  $ descx, work, iw, 1, descw, 1, one, x, ix+k,
397  $ jx+k-1, descx, 1 )
398 *
399  CALL pcelget( 'Columnwise', ' ', tau, taup, i, 1, desctp )
400  CALL pcscal( m-k, tau, x, ix+k, jx+k-1, descx, 1 )
401  CALL pclacgv( n-k, a, i, j+1, desca, desca( m_ ) )
402  10 CONTINUE
403 *
404  ELSE
405 *
406 * Reduce to lower bidiagonal form
407 *
408  CALL descset( descd, ia+min(m,n)-1, 1, desca( mb_ ), 1,
409  $ desca( rsrc_ ), mycol, desca( ctxt_ ),
410  $ desca( lld_ ) )
411  CALL descset( desce, 1, ja+min(m,n)-1, 1, desca( nb_ ), myrow,
412  $ desca( csrc_ ), desca( ctxt_ ), 1 )
413  DO 20 k = 1, nb
414  i = ia + k - 1
415  j = ja + k - 1
416  jwy = iw + k
417 *
418 * Update A(i,j:ja+n-1)
419 *
420  CALL pclacgv( n-k+1, a, i, j, desca, desca( m_ ) )
421  IF( k.GT.1 ) THEN
422  CALL pclacgv( k-1, a, i, ja, desca, desca( m_ ) )
423  CALL pcgemv( 'Conjugate transpose', k-1, n-k+1, -one, y,
424  $ iy, jy+k-1, descy, a, i, ja, desca,
425  $ desca( m_ ), one, a, i, j, desca,
426  $ desca( m_ ) )
427  CALL pclacgv( k-1, a, i, ja, desca, desca( m_ ) )
428  CALL pclacgv( k-1, x, ix+k-1, jx, descx, descx( m_ ) )
429  CALL pcgemv( 'Conjugate transpose', k-1, n-k+1, -one, a,
430  $ ia, j, desca, x, ix+k-1, jx, descx,
431  $ descx( m_ ), one, a, i, j, desca,
432  $ desca( m_ ) )
433  CALL pclacgv( k-1, x, ix+k-1, jx, descx, descx( m_ ) )
434  CALL pcelset( a, i, j-1, desca, cmplx( real( alpha ) ) )
435  END IF
436 *
437 * Generate reflection P(i) to annihilate A(i,j+1:ja+n-1)
438 *
439  CALL pclarfg( n-k+1, alpha, i, j, a, i, j+1, desca,
440  $ desca( m_ ), taup )
441  CALL pselset( d, i, 1, descd, real( alpha ) )
442  CALL pcelset( a, i, j, desca, one )
443 *
444 * Compute X(i+1:ia+m-1,j)
445 *
446  CALL pcgemv( 'No transpose', m-k, n-k+1, one, a, i+1, j,
447  $ desca, a, i, j, desca, desca( m_ ), zero, x,
448  $ ix+k, jx+k-1, descx, 1 )
449  CALL pcgemv( 'No transpose', k-1, n-k+1, one, y, iy, jy+k-1,
450  $ descy, a, i, j, desca, desca( m_ ), zero,
451  $ work, iw, 1, descw, 1 )
452  CALL pcgemv( 'No transpose', m-k, k-1, -one, a, i+1, ja,
453  $ desca, work, iw, 1, descw, 1, one, x, ix+k,
454  $ jx+k-1, descx, 1 )
455  CALL pcgemv( 'No transpose', k-1, n-k+1, one, a, ia, j,
456  $ desca, a, i, j, desca, desca( m_ ), zero,
457  $ work, iw, 1, descw, 1 )
458  CALL pcgemv( 'No transpose', m-k, k-1, -one, x, ix+k, jx,
459  $ descx, work, iw, 1, descw, 1, one, x, ix+k,
460  $ jx+k-1, descx, 1 )
461 *
462  CALL pcelget( 'Columnwise', ' ', tau, taup, i, 1, desctp )
463  CALL pcscal( m-k, tau, x, ix+k, jx+k-1, descx, 1 )
464  CALL pclacgv( n-k+1, a, i, j, desca, desca( m_ ) )
465 *
466 * Update A(i+1:ia+m-1,j)
467 *
468  CALL pcgemv( 'No transpose', m-k, k-1, -one, a, i+1, ja,
469  $ desca, y, iy, jy+k-1, descy, 1, one, a, i+1, j,
470  $ desca, 1 )
471  CALL pcgemv( 'No transpose', m-k, k, -one, x, ix+k, jx,
472  $ descx, a, ia, j, desca, 1, one, a, i+1, j,
473  $ desca, 1 )
474  CALL pcelset( a, i, j, desca, alpha )
475 *
476 * Generate reflection Q(i) to annihilate A(i+2:ia+m-1,j)
477 *
478  CALL pclarfg( m-k, alpha, i+1, j, a, min( i+2, m+ia-1 ),
479  $ j, desca, 1, tauq )
480  CALL pselset( e, 1, j, desce, real( alpha ) )
481  CALL pcelset( a, i+1, j, desca, one )
482 *
483 * Compute Y(ia+i:ia+n-1,j)
484 *
485  CALL pcgemv( 'Conjugate transpose', m-k, n-k, one, a, i+1,
486  $ j+1, desca, a, i+1, j, desca, 1, zero,
487  $ work( ipy ), 1, jwy, descwy, descwy( m_ ) )
488  CALL pcgemv( 'Conjugate transpose', m-k, k-1, one, a, i+1,
489  $ ja, desca, a, i+1, j, desca, 1, zero, work, iw,
490  $ 1, descw, 1 )
491  CALL pcgemv( 'Conjugate transpose', k-1, n-k, -one, y, iy,
492  $ jy+k, descy, work, iw, 1, descw, 1, one,
493  $ work( ipy ), 1, jwy, descwy, descwy( m_ ) )
494  CALL pcgemv( 'Conjugate transpose', m-k, k, one, x, ix+k,
495  $ jx, descx, a, i+1, j, desca, 1, zero, work, iw,
496  $ 1, descw, 1 )
497  CALL pcgemv( 'Conjugate transpose', k, n-k, -one, a, ia,
498  $ j+1, desca, work, iw, 1, descw, 1, one,
499  $ work( ipy ), 1, jwy, descwy, descwy( m_ ) )
500 *
501  CALL pcelget( 'Rowwise', ' ', tau, tauq, 1, j, desctq )
502  CALL pcscal( n-k, tau, work( ipy ), 1, jwy, descwy,
503  $ descwy( m_ ) )
504  CALL pclacgv( n-k, work( ipy ), 1, jwy, descwy,
505  $ descwy( m_ ) )
506  CALL pccopy( n-k, work( ipy ), 1, jwy, descwy, descwy( m_ ),
507  $ y, iy+k-1, jy+k, descy, descy( m_ ) )
508  20 CONTINUE
509  END IF
510 *
511  RETURN
512 *
513 * End of PCLABRD
514 *
515  END
cmplx
float cmplx[2]
Definition: pblas.h:132
pclarfg
subroutine pclarfg(N, ALPHA, IAX, JAX, X, IX, JX, DESCX, INCX, TAU)
Definition: pclarfg.f:3
pcelget
subroutine pcelget(SCOPE, TOP, ALPHA, A, IA, JA, DESCA)
Definition: pcelget.f:2
infog2l
subroutine infog2l(GRINDX, GCINDX, DESC, NPROW, NPCOL, MYROW, MYCOL, LRINDX, LCINDX, RSRC, CSRC)
Definition: infog2l.f:3
pselset
subroutine pselset(A, IA, JA, DESCA, ALPHA)
Definition: pselset.f:2
pclabrd
subroutine pclabrd(M, N, NB, A, IA, JA, DESCA, D, E, TAUQ, TAUP, X, IX, JX, DESCX, Y, IY, JY, DESCY, WORK)
Definition: pclabrd.f:3
pcelset
subroutine pcelset(A, IA, JA, DESCA, ALPHA)
Definition: pcelset.f:2
pclacgv
subroutine pclacgv(N, X, IX, JX, DESCX, INCX)
Definition: pclacgv.f:2
descset
subroutine descset(DESC, M, N, MB, NB, IRSRC, ICSRC, ICTXT, LLD)
Definition: descset.f:3
min
#define min(A, B)
Definition: pcgemr.c:181