SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
pclatrz.f
Go to the documentation of this file.
1 SUBROUTINE pclatrz( M, N, L, A, IA, JA, DESCA, TAU, WORK )
2*
3* -- ScaLAPACK routine (version 1.7) --
4* University of Tennessee, Knoxville, Oak Ridge National Laboratory,
5* and University of California, Berkeley.
6* December 31, 1998
7*
8* .. Scalar Arguments ..
9 INTEGER IA, JA, L, M, N
10* ..
11* .. Array Arguments ..
12 INTEGER DESCA( * )
13 COMPLEX A( * ), TAU( * ), WORK( * )
14* ..
15*
16* Purpose
17* =======
18*
19* PCLATRZ reduces the M-by-N ( M<=N ) complex upper trapezoidal
20* matrix sub( A ) = [A(IA:IA+M-1,JA:JA+M-1) A(IA:IA+M-1,JA+N-L:JA+N-1)]
21* to upper triangular form by means of unitary transformations.
22*
23* The upper trapezoidal matrix sub( A ) is factored as
24*
25* sub( A ) = ( R 0 ) * Z,
26*
27* where Z is an N-by-N unitary matrix and R is an M-by-M upper
28* triangular matrix.
29*
30* Notes
31* =====
32*
33* Each global data object is described by an associated description
34* vector. This vector stores the information required to establish
35* the mapping between an object element and its corresponding process
36* and memory location.
37*
38* Let A be a generic term for any 2D block cyclicly distributed array.
39* Such a global array has an associated description vector DESCA.
40* In the following comments, the character _ should be read as
41* "of the global array".
42*
43* NOTATION STORED IN EXPLANATION
44* --------------- -------------- --------------------------------------
45* DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
46* DTYPE_A = 1.
47* CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
48* the BLACS process grid A is distribu-
49* ted over. The context itself is glo-
50* bal, but the handle (the integer
51* value) may vary.
52* M_A (global) DESCA( M_ ) The number of rows in the global
53* array A.
54* N_A (global) DESCA( N_ ) The number of columns in the global
55* array A.
56* MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
57* the rows of the array.
58* NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
59* the columns of the array.
60* RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
61* row of the array A is distributed.
62* CSRC_A (global) DESCA( CSRC_ ) The process column over which the
63* first column of the array A is
64* distributed.
65* LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
66* array. LLD_A >= MAX(1,LOCr(M_A)).
67*
68* Let K be the number of rows or columns of a distributed matrix,
69* and assume that its process grid has dimension p x q.
70* LOCr( K ) denotes the number of elements of K that a process
71* would receive if K were distributed over the p processes of its
72* process column.
73* Similarly, LOCc( K ) denotes the number of elements of K that a
74* process would receive if K were distributed over the q processes of
75* its process row.
76* The values of LOCr() and LOCc() may be determined via a call to the
77* ScaLAPACK tool function, NUMROC:
78* LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
79* LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
80* An upper bound for these quantities may be computed by:
81* LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
82* LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
83*
84* Arguments
85* =========
86*
87* M (global input) INTEGER
88* The number of rows to be operated on, i.e. the number of rows
89* of the distributed submatrix sub( A ). M >= 0.
90*
91* N (global input) INTEGER
92* The number of columns to be operated on, i.e. the number of
93* columns of the distributed submatrix sub( A ). N >= 0.
94*
95* L (global input) INTEGER
96* The columns of the distributed submatrix sub( A ) containing
97* the meaningful part of the Householder reflectors. L > 0.
98*
99* A (local input/local output) COMPLEX pointer into the
100* local memory to an array of dimension (LLD_A, LOCc(JA+N-1)).
101* On entry, the local pieces of the M-by-N distributed matrix
102* sub( A ) which is to be factored. On exit, the leading M-by-M
103* upper triangular part of sub( A ) contains the upper trian-
104* gular matrix R, and elements N-L+1 to N of the first M rows
105* of sub( A ), with the array TAU, represent the unitary matrix
106* Z as a product of M elementary reflectors.
107*
108* IA (global input) INTEGER
109* The row index in the global array A indicating the first
110* row of sub( A ).
111*
112* JA (global input) INTEGER
113* The column index in the global array A indicating the
114* first column of sub( A ).
115*
116* DESCA (global and local input) INTEGER array of dimension DLEN_.
117* The array descriptor for the distributed matrix A.
118*
119* TAU (local output) COMPLEX, array, dimension LOCr(IA+M-1)
120* This array contains the scalar factors of the elementary
121* reflectors. TAU is tied to the distributed matrix A.
122*
123* WORK (local workspace) COMPLEX array, dimension (LWORK)
124* LWORK >= Nq0 + MAX( 1, Mp0 ), where
125*
126* IROFF = MOD( IA-1, MB_A ), ICOFF = MOD( JA-1, NB_A ),
127* IAROW = INDXG2P( IA, MB_A, MYROW, RSRC_A, NPROW ),
128* IACOL = INDXG2P( JA, NB_A, MYCOL, CSRC_A, NPCOL ),
129* Mp0 = NUMROC( M+IROFF, MB_A, MYROW, IAROW, NPROW ),
130* Nq0 = NUMROC( N+ICOFF, NB_A, MYCOL, IACOL, NPCOL ),
131*
132* and NUMROC, INDXG2P are ScaLAPACK tool functions;
133* MYROW, MYCOL, NPROW and NPCOL can be determined by calling
134* the subroutine BLACS_GRIDINFO.
135*
136* Further Details
137* ===============
138*
139* The factorization is obtained by Householder's method. The kth
140* transformation matrix, Z( k ), whose conjugate transpose is used to
141* introduce zeros into the (m - k + 1)th row of sub( A ), is given in
142* the form
143*
144* Z( k ) = ( I 0 ),
145* ( 0 T( k ) )
146*
147* where
148*
149* T( k ) = I - tau*u( k )*u( k )', u( k ) = ( 1 ),
150* ( 0 )
151* ( z( k ) )
152*
153* tau is a scalar and z( k ) is an ( n - m ) element vector.
154* tau and z( k ) are chosen to annihilate the elements of the kth row
155* of sub( A ).
156*
157* The scalar tau is returned in the kth element of TAU and the vector
158* u( k ) in the kth row of sub( A ), such that the elements of z( k )
159* are in a( k, m + 1 ), ..., a( k, n ). The elements of R are returned
160* in the upper triangular part of sub( A ).
161*
162* Z is given by
163*
164* Z = Z( 1 ) * Z( 2 ) * ... * Z( m ).
165*
166* =====================================================================
167*
168* .. Parameters ..
169 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
170 $ LLD_, MB_, M_, NB_, N_, RSRC_
171 parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
172 $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
173 $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
174 COMPLEX ONE, ZERO
175 parameter( one = ( 1.0e+0, 0.0e+0 ),
176 $ zero = ( 0.0e+0, 0.0e+0 ) )
177* ..
178* .. Local Scalars ..
179 INTEGER I, IAROW, ICTXT, II, J, J1, MP, MYCOL, MYROW,
180 $ NPCOL, NPROW
181 COMPLEX AII
182* ..
183* .. Local Arrays ..
184 INTEGER DESCTAU( DLEN_ )
185* ..
186* .. External Subroutines ..
187 EXTERNAL descset, infog1l, pcelset, pclacgv,
188 $ pclarfg, pclarz
189* ..
190* .. External Functions ..
191 INTEGER NUMROC
192 EXTERNAL numroc
193* ..
194* .. Intrinsic Functions ..
195 INTRINSIC conjg, max
196* ..
197* .. Executable Statements ..
198*
199* Quick return if possible
200*
201 IF( m.EQ.0 .OR. n.EQ.0 )
202 $ RETURN
203*
204* Get grid parameters
205*
206 ictxt = desca( ctxt_ )
207 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
208*
209 mp = numroc( ia+m-1, desca( mb_ ), myrow, desca( rsrc_ ),
210 $ nprow )
211*
212 CALL descset( desctau, desca( m_ ), 1, desca( mb_ ), 1,
213 $ desca( rsrc_ ), mycol, ictxt, max( 1, mp ) )
214*
215 IF( m.EQ.n ) THEN
216*
217 CALL infog1l( ia, desca( mb_ ), nprow, myrow, desca( rsrc_ ),
218 $ ii, iarow )
219 DO 10 i = ii, mp
220 tau( i ) = zero
221 10 CONTINUE
222*
223 ELSE
224*
225 aii = zero
226*
227 j1 = ja + n - l
228 DO 20 i = ia+m-1, ia, -1
229 j = ja + i - ia
230*
231* Generate elementary reflector H(i) to annihilate
232* [ A(i, j) A(i,j1:ja+n-1) ]
233*
234 CALL pclacgv( 1, a, i, j, desca, desca( m_ ) )
235 CALL pclacgv( l, a, i, j1, desca, desca( m_ ) )
236 CALL pclarfg( l+1, aii, i, j, a, i, j1, desca, desca( m_ ),
237 $ tau )
238*
239* Apply H(i) to A(ia:i-1,j:ja+n-1) from the right
240*
241 CALL pclarz( 'Right', i-ia, ja+n-j, l, a, i, j1, desca,
242 $ desca( m_ ), tau, a, ia, j, desca, work )
243 CALL pcelset( a, i, j, desca, conjg( aii ) )
244*
245 20 CONTINUE
246*
247 CALL pclacgv( m, tau, ia, 1, desctau, 1 )
248*
249 END IF
250*
251 RETURN
252*
253* End of PCLATRZ
254*
255 END
subroutine descset(desc, m, n, mb, nb, irsrc, icsrc, ictxt, lld)
Definition descset.f:3
subroutine infog1l(gindx, nb, nprocs, myroc, isrcproc, lindx, rocsrc)
Definition infog1l.f:3
subroutine pcelset(a, ia, ja, desca, alpha)
Definition pcelset.f:2
#define max(A, B)
Definition pcgemr.c:180
subroutine pclacgv(n, x, ix, jx, descx, incx)
Definition pclacgv.f:2
subroutine pclarfg(n, alpha, iax, jax, x, ix, jx, descx, incx, tau)
Definition pclarfg.f:3
subroutine pclarz(side, m, n, l, v, iv, jv, descv, incv, tau, c, ic, jc, descc, work)
Definition pclarz.f:3
subroutine pclatrz(m, n, l, a, ia, ja, desca, tau, work)
Definition pclatrz.f:2