ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
psgehrd.f
Go to the documentation of this file.
1  SUBROUTINE psgehrd( N, ILO, IHI, A, IA, JA, DESCA, TAU, WORK,
2  $ LWORK, INFO )
3 *
4 * -- ScaLAPACK routine (version 1.7) --
5 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
6 * and University of California, Berkeley.
7 * May 25, 2001
8 *
9 * .. Scalar Arguments ..
10  INTEGER IA, IHI, ILO, INFO, JA, LWORK, N
11 * ..
12 * .. Array Arguments ..
13  INTEGER DESCA( * )
14  REAL A( * ), TAU( * ), WORK( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * PSGEHRD reduces a real general distributed matrix sub( A )
21 * to upper Hessenberg form H by an orthogonal similarity transforma-
22 * tion: Q' * sub( A ) * Q = H, where
23 * sub( A ) = A(IA+N-1:IA+N-1,JA+N-1:JA+N-1).
24 *
25 * Notes
26 * =====
27 *
28 * Each global data object is described by an associated description
29 * vector. This vector stores the information required to establish
30 * the mapping between an object element and its corresponding process
31 * and memory location.
32 *
33 * Let A be a generic term for any 2D block cyclicly distributed array.
34 * Such a global array has an associated description vector DESCA.
35 * In the following comments, the character _ should be read as
36 * "of the global array".
37 *
38 * NOTATION STORED IN EXPLANATION
39 * --------------- -------------- --------------------------------------
40 * DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
41 * DTYPE_A = 1.
42 * CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
43 * the BLACS process grid A is distribu-
44 * ted over. The context itself is glo-
45 * bal, but the handle (the integer
46 * value) may vary.
47 * M_A (global) DESCA( M_ ) The number of rows in the global
48 * array A.
49 * N_A (global) DESCA( N_ ) The number of columns in the global
50 * array A.
51 * MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
52 * the rows of the array.
53 * NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
54 * the columns of the array.
55 * RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
56 * row of the array A is distributed.
57 * CSRC_A (global) DESCA( CSRC_ ) The process column over which the
58 * first column of the array A is
59 * distributed.
60 * LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
61 * array. LLD_A >= MAX(1,LOCr(M_A)).
62 *
63 * Let K be the number of rows or columns of a distributed matrix,
64 * and assume that its process grid has dimension p x q.
65 * LOCr( K ) denotes the number of elements of K that a process
66 * would receive if K were distributed over the p processes of its
67 * process column.
68 * Similarly, LOCc( K ) denotes the number of elements of K that a
69 * process would receive if K were distributed over the q processes of
70 * its process row.
71 * The values of LOCr() and LOCc() may be determined via a call to the
72 * ScaLAPACK tool function, NUMROC:
73 * LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
74 * LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
75 * An upper bound for these quantities may be computed by:
76 * LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
77 * LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
78 *
79 * Arguments
80 * =========
81 *
82 * N (global input) INTEGER
83 * The number of rows and columns to be operated on, i.e. the
84 * order of the distributed submatrix sub( A ). N >= 0.
85 *
86 * ILO (global input) INTEGER
87 * IHI (global input) INTEGER
88 * It is assumed that sub( A ) is already upper triangular in
89 * rows IA:IA+ILO-2 and IA+IHI:IA+N-1 and columns JA:JA+ILO-2
90 * and JA+IHI:JA+N-1. See Further Details. If N > 0,
91 * 1 <= ILO <= IHI <= N; otherwise set ILO = 1, IHI = N.
92 *
93 * A (local input/local output) REAL pointer into the
94 * local memory to an array of dimension (LLD_A,LOCc(JA+N-1)).
95 * On entry, this array contains the local pieces of the N-by-N
96 * general distributed matrix sub( A ) to be reduced. On exit,
97 * the upper triangle and the first subdiagonal of sub( A ) are
98 * overwritten with the upper Hessenberg matrix H, and the ele-
99 * ments below the first subdiagonal, with the array TAU, repre-
100 * sent the orthogonal matrix Q as a product of elementary
101 * reflectors. See Further Details.
102 *
103 * IA (global input) INTEGER
104 * The row index in the global array A indicating the first
105 * row of sub( A ).
106 *
107 * JA (global input) INTEGER
108 * The column index in the global array A indicating the
109 * first column of sub( A ).
110 *
111 * DESCA (global and local input) INTEGER array of dimension DLEN_.
112 * The array descriptor for the distributed matrix A.
113 *
114 * TAU (local output) REAL array, dimension LOCc(JA+N-2)
115 * The scalar factors of the elementary reflectors (see Further
116 * Details). Elements JA:JA+ILO-2 and JA+IHI:JA+N-2 of TAU are
117 * set to zero. TAU is tied to the distributed matrix A.
118 *
119 * WORK (local workspace/local output) REAL array,
120 * dimension (LWORK)
121 * On exit, WORK( 1 ) returns the minimal and optimal LWORK.
122 *
123 * LWORK (local or global input) INTEGER
124 * The dimension of the array WORK.
125 * LWORK is local input and must be at least
126 * LWORK >= NB*NB + NB*MAX( IHIP+1, IHLP+INLQ )
127 *
128 * where NB = MB_A = NB_A, IROFFA = MOD( IA-1, NB ),
129 * ICOFFA = MOD( JA-1, NB ), IOFF = MOD( IA+ILO-2, NB ),
130 * IAROW = INDXG2P( IA, NB, MYROW, RSRC_A, NPROW ),
131 * IHIP = NUMROC( IHI+IROFFA, NB, MYROW, IAROW, NPROW ),
132 * ILROW = INDXG2P( IA+ILO-1, NB, MYROW, RSRC_A, NPROW ),
133 * IHLP = NUMROC( IHI-ILO+IOFF+1, NB, MYROW, ILROW, NPROW ),
134 * ILCOL = INDXG2P( JA+ILO-1, NB, MYCOL, CSRC_A, NPCOL ),
135 * INLQ = NUMROC( N-ILO+IOFF+1, NB, MYCOL, ILCOL, NPCOL ),
136 *
137 * INDXG2P and NUMROC are ScaLAPACK tool functions;
138 * MYROW, MYCOL, NPROW and NPCOL can be determined by calling
139 * the subroutine BLACS_GRIDINFO.
140 *
141 * If LWORK = -1, then LWORK is global input and a workspace
142 * query is assumed; the routine only calculates the minimum
143 * and optimal size for all work arrays. Each of these
144 * values is returned in the first entry of the corresponding
145 * work array, and no error message is issued by PXERBLA.
146 *
147 * INFO (global output) INTEGER
148 * = 0: successful exit
149 * < 0: If the i-th argument is an array and the j-entry had
150 * an illegal value, then INFO = -(i*100+j), if the i-th
151 * argument is a scalar and had an illegal value, then
152 * INFO = -i.
153 *
154 * Further Details
155 * ===============
156 *
157 * The matrix Q is represented as a product of (ihi-ilo) elementary
158 * reflectors
159 *
160 * Q = H(ilo) H(ilo+1) . . . H(ihi-1).
161 *
162 * Each H(i) has the form
163 *
164 * H(i) = I - tau * v * v'
165 *
166 * where tau is a real scalar, and v is a real vector with
167 * v(1:I) = 0, v(I+1) = 1 and v(IHI+1:N) = 0; v(I+2:IHI) is stored on
168 * exit in A(IA+ILO+I:IA+IHI-1,JA+ILO+I-2), and tau in TAU(JA+ILO+I-2).
169 *
170 * The contents of A(IA:IA+N-1,JA:JA+N-1) are illustrated by the follow-
171 * ing example, with N = 7, ILO = 2 and IHI = 6:
172 *
173 * on entry on exit
174 *
175 * ( a a a a a a a ) ( a a h h h h a )
176 * ( a a a a a a ) ( a h h h h a )
177 * ( a a a a a a ) ( h h h h h h )
178 * ( a a a a a a ) ( v2 h h h h h )
179 * ( a a a a a a ) ( v2 v3 h h h h )
180 * ( a a a a a a ) ( v2 v3 v4 h h h )
181 * ( a ) ( a )
182 *
183 * where a denotes an element of the original matrix sub( A ), H denotes
184 * a modified element of the upper Hessenberg matrix H, and vi denotes
185 * an element of the vector defining H(JA+ILO+I-2).
186 *
187 * Alignment requirements
188 * ======================
189 *
190 * The distributed submatrix sub( A ) must verify some alignment proper-
191 * ties, namely the following expression should be true:
192 * ( MB_A.EQ.NB_A .AND. IROFFA.EQ.ICOFFA )
193 *
194 * =====================================================================
195 *
196 * .. Parameters ..
197  INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
198  $ lld_, mb_, m_, nb_, n_, rsrc_
199  parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
200  $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
201  $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
202  REAL ONE, ZERO
203  parameter( one = 1.0e+0, zero = 0.0e+0 )
204 * ..
205 * .. Local Scalars ..
206  LOGICAL LQUERY
207  CHARACTER COLCTOP, ROWCTOP
208  INTEGER I, IACOL, IAROW, IB, ICOFFA, ICTXT, IHIP,
209  $ ihlp, iia, iinfo, ilcol, ilrow, imcol, inlq,
210  $ ioff, ipt, ipw, ipy, iroffa, j, jj, jja, jy,
211  $ k, l, lwmin, mycol, myrow, nb, npcol, nprow,
212  $ nq
213  REAL EI
214 * ..
215 * .. Local Arrays ..
216  INTEGER DESCY( DLEN_ ), IDUM1( 3 ), IDUM2( 3 )
217 * ..
218 * .. External Subroutines ..
219  EXTERNAL blacs_gridinfo, chk1mat, descset, infog1l,
220  $ infog2l, pchk1mat, psgemm, psgehd2,
221  $ pslahrd, pslarfb, pb_topget, pb_topset, pxerbla
222 * ..
223 * .. External Functions ..
224  INTEGER INDXG2P, NUMROC
225  EXTERNAL indxg2p, numroc
226 * ..
227 * .. Intrinsic Functions ..
228  INTRINSIC float, max, min, mod
229 * ..
230 * .. Executable Statements ..
231 *
232 * Get grid parameters
233 *
234  ictxt = desca( ctxt_ )
235  CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
236 *
237 * Test the input parameters
238 *
239  info = 0
240  IF( nprow.EQ.-1 ) THEN
241  info = -(700+ctxt_)
242  ELSE
243  CALL chk1mat( n, 1, n, 1, ia, ja, desca, 7, info )
244  IF( info.EQ.0 ) THEN
245  nb = desca( nb_ )
246  iroffa = mod( ia-1, nb )
247  icoffa = mod( ja-1, nb )
248  CALL infog2l( ia, ja, desca, nprow, npcol, myrow, mycol,
249  $ iia, jja, iarow, iacol )
250  ihip = numroc( ihi+iroffa, nb, myrow, iarow, nprow )
251  ioff = mod( ia+ilo-2, nb )
252  ilrow = indxg2p( ia+ilo-1, nb, myrow, desca( rsrc_ ),
253  $ nprow )
254  ihlp = numroc( ihi-ilo+ioff+1, nb, myrow, ilrow, nprow )
255  ilcol = indxg2p( ja+ilo-1, nb, mycol, desca( csrc_ ),
256  $ npcol )
257  inlq = numroc( n-ilo+ioff+1, nb, mycol, ilcol, npcol )
258  lwmin = nb*( nb + max( ihip+1, ihlp+inlq ) )
259 *
260  work( 1 ) = float( lwmin )
261  lquery = ( lwork.EQ.-1 )
262  IF( ilo.LT.1 .OR. ilo.GT.max( 1, n ) ) THEN
263  info = -2
264  ELSE IF( ihi.LT.min( ilo, n ) .OR. ihi.GT.n ) THEN
265  info = -3
266 C ELSE IF( IROFFA.NE.ICOFFA .OR. IROFFA.NE.0 ) THEN
267  ELSE IF( iroffa.NE.icoffa ) THEN
268  info = -6
269  ELSE IF( desca( mb_ ).NE.desca( nb_ ) ) THEN
270  info = -(700+nb_)
271  ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
272  info = -10
273  END IF
274  END IF
275  idum1( 1 ) = ilo
276  idum2( 1 ) = 2
277  idum1( 2 ) = ihi
278  idum2( 2 ) = 3
279  IF( lwork.EQ.-1 ) THEN
280  idum1( 3 ) = -1
281  ELSE
282  idum1( 3 ) = 1
283  END IF
284  idum2( 3 ) = 10
285  CALL pchk1mat( n, 1, n, 1, ia, ja, desca, 7, 3, idum1, idum2,
286  $ info )
287  END IF
288 *
289  IF( info.NE.0 ) THEN
290  CALL pxerbla( ictxt, 'PSGEHRD', -info )
291  RETURN
292  ELSE IF( lquery ) THEN
293  RETURN
294  END IF
295 *
296 * Set elements JA:JA+ILO-2 and JA+JHI-1:JA+N-2 of TAU to zero.
297 *
298  nq = numroc( ja+n-2, nb, mycol, desca( csrc_ ), npcol )
299  CALL infog1l( ja+ilo-2, nb, npcol, mycol, desca( csrc_ ), jj,
300  $ imcol )
301  DO 10 j = jja, min( jj, nq )
302  tau( j ) = zero
303  10 CONTINUE
304 *
305  CALL infog1l( ja+ihi-1, nb, npcol, mycol, desca( csrc_ ), jj,
306  $ imcol )
307  DO 20 j = jj, nq
308  tau( j ) = zero
309  20 CONTINUE
310 *
311 * Quick return if possible
312 *
313  IF( ihi-ilo.LE.0 )
314  $ RETURN
315 *
316  CALL pb_topget( ictxt, 'Combine', 'Columnwise', colctop )
317  CALL pb_topget( ictxt, 'Combine', 'Rowwise', rowctop )
318  CALL pb_topset( ictxt, 'Combine', 'Columnwise', '1-tree' )
319  CALL pb_topset( ictxt, 'Combine', 'Rowwise', '1-tree' )
320 *
321  ipt = 1
322  ipy = ipt + nb * nb
323  ipw = ipy + ihip * nb
324  CALL descset( descy, ihi+iroffa, nb, nb, nb, iarow, ilcol, ictxt,
325  $ max( 1, ihip ) )
326 *
327  k = ilo
328  ib = nb - ioff
329  jy = ioff + 1
330 *
331 * Loop over remaining block of columns
332 *
333  DO 30 l = 1, ihi-ilo+ioff-nb, nb
334  i = ia + k - 1
335  j = ja + k - 1
336 *
337 * Reduce columns j:j+ib-1 to Hessenberg form, returning the
338 * matrices V and T of the block reflector H = I - V*T*V'
339 * which performs the reduction, and also the matrix Y = A*V*T
340 *
341  CALL pslahrd( ihi, k, ib, a, ia, j, desca, tau, work( ipt ),
342  $ work( ipy ), 1, jy, descy, work( ipw ) )
343 *
344 * Apply the block reflector H to A(ia:ia+ihi-1,j+ib:ja+ihi-1)
345 * from the right, computing A := A - Y * V'.
346 * V(i+ib,ib-1) must be set to 1.
347 *
348  CALL pselset2( ei, a, i+ib, j+ib-1, desca, one )
349  CALL psgemm( 'No transpose', 'Transpose', ihi, ihi-k-ib+1, ib,
350  $ -one, work( ipy ), 1, jy, descy, a, i+ib, j,
351  $ desca, one, a, ia, j+ib, desca )
352  CALL pselset( a, i+ib, j+ib-1, desca, ei )
353 *
354 * Apply the block reflector H to A(i+1:ia+ihi-1,j+ib:ja+n-1) from
355 * the left
356 *
357  CALL pslarfb( 'Left', 'Transpose', 'Forward', 'Columnwise',
358  $ ihi-k, n-k-ib+1, ib, a, i+1, j, desca,
359  $ work( ipt ), a, i+1, j+ib, desca, work( ipy ) )
360 *
361  k = k + ib
362  ib = nb
363  jy = 1
364  descy( csrc_ ) = mod( descy( csrc_ ) + 1, npcol )
365 *
366  30 CONTINUE
367 *
368 * Use unblocked code to reduce the rest of the matrix
369 *
370  CALL psgehd2( n, k, ihi, a, ia, ja, desca, tau, work, lwork,
371  $ iinfo )
372 *
373  CALL pb_topset( ictxt, 'Combine', 'Columnwise', colctop )
374  CALL pb_topset( ictxt, 'Combine', 'Rowwise', rowctop )
375 *
376  work( 1 ) = float( lwmin )
377 *
378  RETURN
379 *
380 * End of PSGEHRD
381 *
382  END
psgehrd
subroutine psgehrd(N, ILO, IHI, A, IA, JA, DESCA, TAU, WORK, LWORK, INFO)
Definition: psgehrd.f:3
max
#define max(A, B)
Definition: pcgemr.c:180
infog1l
subroutine infog1l(GINDX, NB, NPROCS, MYROC, ISRCPROC, LINDX, ROCSRC)
Definition: infog1l.f:3
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
pselset2
subroutine pselset2(ALPHA, A, IA, JA, DESCA, BETA)
Definition: pselset2.f:2
pchk1mat
subroutine pchk1mat(MA, MAPOS0, NA, NAPOS0, IA, JA, DESCA, DESCAPOS0, NEXTRA, EX, EXPOS, INFO)
Definition: pchkxmat.f:3
psgehd2
subroutine psgehd2(N, ILO, IHI, A, IA, JA, DESCA, TAU, WORK, LWORK, INFO)
Definition: psgehd2.f:3
descset
subroutine descset(DESC, M, N, MB, NB, IRSRC, ICSRC, ICTXT, LLD)
Definition: descset.f:3
pslahrd
subroutine pslahrd(N, K, NB, A, IA, JA, DESCA, TAU, T, Y, IY, JY, DESCY, WORK)
Definition: pslahrd.f:3
chk1mat
subroutine chk1mat(MA, MAPOS0, NA, NAPOS0, IA, JA, DESCA, DESCAPOS0, INFO)
Definition: chk1mat.f:3
pslarfb
subroutine pslarfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, V, IV, JV, DESCV, T, C, IC, JC, DESCC, WORK)
Definition: pslarfb.f:3
pxerbla
subroutine pxerbla(ICTXT, SRNAME, INFO)
Definition: pxerbla.f:2
min
#define min(A, B)
Definition: pcgemr.c:181