SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ pzhegs2()

subroutine pzhegs2 ( integer  ibtype,
character  uplo,
integer  n,
complex*16, dimension( * )  a,
integer  ia,
integer  ja,
integer, dimension( * )  desca,
complex*16, dimension( * )  b,
integer  ib,
integer  jb,
integer, dimension( * )  descb,
integer  info 
)

Definition at line 3 of file pzhegs2.f.

5*
6* -- ScaLAPACK routine (version 1.7) --
7* University of Tennessee, Knoxville, Oak Ridge National Laboratory,
8* and University of California, Berkeley.
9* May 1, 1997
10*
11* .. Scalar Arguments ..
12 CHARACTER UPLO
13 INTEGER IA, IB, IBTYPE, INFO, JA, JB, N
14* ..
15* .. Array Arguments ..
16 INTEGER DESCA( * ), DESCB( * )
17 COMPLEX*16 A( * ), B( * )
18* ..
19*
20* Purpose
21* =======
22*
23* PZHEGS2 reduces a complex Hermitian-definite generalized eigenproblem
24* to standard form.
25*
26* In the following sub( A ) denotes A( IA:IA+N-1, JA:JA+N-1 ) and
27* sub( B ) denotes B( IB:IB+N-1, JB:JB+N-1 ).
28*
29* If IBTYPE = 1, the problem is sub( A )*x = lambda*sub( B )*x,
30* and sub( A ) is overwritten by inv(U**H)*sub( A )*inv(U) or
31* inv(L)*sub( A )*inv(L**H)
32*
33* If IBTYPE = 2 or 3, the problem is sub( A )*sub( B )*x = lambda*x or
34* sub( B )*sub( A )*x = lambda*x, and sub( A ) is overwritten by
35* U*sub( A )*U**H or L**H*sub( A )*L.
36*
37* sub( B ) must have been previously factorized as U**H*U or L*L**H by
38* PZPOTRF.
39*
40* Notes
41* =====
42*
43* Each global data object is described by an associated description
44* vector. This vector stores the information required to establish
45* the mapping between an object element and its corresponding process
46* and memory location.
47*
48* Let A be a generic term for any 2D block cyclicly distributed array.
49* Such a global array has an associated description vector DESCA.
50* In the following comments, the character _ should be read as
51* "of the global array".
52*
53* NOTATION STORED IN EXPLANATION
54* --------------- -------------- --------------------------------------
55* DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
56* DTYPE_A = 1.
57* CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
58* the BLACS process grid A is distribu-
59* ted over. The context itself is glo-
60* bal, but the handle (the integer
61* value) may vary.
62* M_A (global) DESCA( M_ ) The number of rows in the global
63* array A.
64* N_A (global) DESCA( N_ ) The number of columns in the global
65* array A.
66* MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
67* the rows of the array.
68* NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
69* the columns of the array.
70* RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
71* row of the array A is distributed.
72* CSRC_A (global) DESCA( CSRC_ ) The process column over which the
73* first column of the array A is
74* distributed.
75* LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
76* array. LLD_A >= MAX(1,LOCr(M_A)).
77*
78* Let K be the number of rows or columns of a distributed matrix,
79* and assume that its process grid has dimension p x q.
80* LOCr( K ) denotes the number of elements of K that a process
81* would receive if K were distributed over the p processes of its
82* process column.
83* Similarly, LOCc( K ) denotes the number of elements of K that a
84* process would receive if K were distributed over the q processes of
85* its process row.
86* The values of LOCr() and LOCc() may be determined via a call to the
87* ScaLAPACK tool function, NUMROC:
88* LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
89* LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
90* An upper bound for these quantities may be computed by:
91* LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
92* LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
93*
94* Arguments
95* =========
96*
97* IBTYPE (global input) INTEGER
98* = 1: compute inv(U**H)*sub( A )*inv(U) or
99* inv(L)*sub( A )*inv(L**H);
100* = 2 or 3: compute U*sub( A )*U**H or L**H*sub( A )*L.
101*
102* UPLO (global input) CHARACTER
103* = 'U': Upper triangle of sub( A ) is stored and sub( B ) is
104* factored as U**H*U;
105* = 'L': Lower triangle of sub( A ) is stored and sub( B ) is
106* factored as L*L**H.
107*
108* N (global input) INTEGER
109* The order of the matrices sub( A ) and sub( B ). N >= 0.
110*
111* A (local input/local output) COMPLEX*16 pointer into the
112* local memory to an array of dimension (LLD_A, LOCc(JA+N-1)).
113* On entry, this array contains the local pieces of the
114* N-by-N Hermitian distributed matrix sub( A ). If UPLO = 'U',
115* the leading N-by-N upper triangular part of sub( A ) contains
116* the upper triangular part of the matrix, and its strictly
117* lower triangular part is not referenced. If UPLO = 'L', the
118* leading N-by-N lower triangular part of sub( A ) contains
119* the lower triangular part of the matrix, and its strictly
120* upper triangular part is not referenced.
121*
122* On exit, if INFO = 0, the transformed matrix, stored in the
123* same format as sub( A ).
124*
125* IA (global input) INTEGER
126* A's global row index, which points to the beginning of the
127* submatrix which is to be operated on.
128*
129* JA (global input) INTEGER
130* A's global column index, which points to the beginning of
131* the submatrix which is to be operated on.
132*
133* DESCA (global and local input) INTEGER array of dimension DLEN_.
134* The array descriptor for the distributed matrix A.
135*
136* B (local input) COMPLEX*16 pointer into the local memory
137* to an array of dimension (LLD_B, LOCc(JB+N-1)). On entry,
138* this array contains the local pieces of the triangular factor
139* from the Cholesky factorization of sub( B ), as returned by
140* PZPOTRF.
141*
142* IB (global input) INTEGER
143* B's global row index, which points to the beginning of the
144* submatrix which is to be operated on.
145*
146* JB (global input) INTEGER
147* B's global column index, which points to the beginning of
148* the submatrix which is to be operated on.
149*
150* DESCB (global and local input) INTEGER array of dimension DLEN_.
151* The array descriptor for the distributed matrix B.
152*
153* INFO (global output) INTEGER
154* = 0: successful exit
155* < 0: If the i-th argument is an array and the j-entry had
156* an illegal value, then INFO = -(i*100+j), if the i-th
157* argument is a scalar and had an illegal value, then
158* INFO = -i.
159*
160* =====================================================================
161*
162* .. Parameters ..
163 INTEGER BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_,
164 $ MB_, NB_, RSRC_, CSRC_, LLD_
165 parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
166 $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
167 $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
168 DOUBLE PRECISION ONE, HALF
169 parameter( one = 1.0d+0, half = 0.5d+0 )
170 COMPLEX*16 CONE
171 parameter( cone = ( 1.0d+0, 0.0d+0 ) )
172* ..
173* .. Local Scalars ..
174 LOGICAL UPPER
175 INTEGER IACOL, IAROW, IBCOL, IBROW, ICOFFA, ICOFFB,
176 $ ICTXT, IIA, IIB, IOFFA, IOFFB, IROFFA, IROFFB,
177 $ JJA, JJB, K, LDA, LDB, MYCOL, MYROW, NPCOL,
178 $ NPROW
179 DOUBLE PRECISION AKK, BKK
180 COMPLEX*16 CT
181* ..
182* .. External Subroutines ..
183 EXTERNAL blacs_exit, blacs_gridinfo, chk1mat, infog2l,
184 $ pxerbla, zaxpy, zdscal, zher2, zlacgv, ztrmv,
185 $ ztrsv
186* ..
187* .. Intrinsic Functions ..
188 INTRINSIC dble, mod
189* ..
190* .. External Functions ..
191 LOGICAL LSAME
192 INTEGER INDXG2P
193 EXTERNAL lsame, indxg2p
194* ..
195* .. Executable Statements ..
196* This is just to keep ftnchek happy
197 IF( block_cyclic_2d*csrc_*ctxt_*dlen_*dtype_*lld_*mb_*m_*nb_*n_*
198 $ rsrc_.LT.0 )RETURN
199*
200* Get grid parameters
201*
202 ictxt = desca( ctxt_ )
203 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
204*
205* Test the input parameters.
206*
207 info = 0
208 IF( nprow.EQ.-1 ) THEN
209 info = -( 700+ctxt_ )
210 ELSE
211 upper = lsame( uplo, 'U' )
212 CALL chk1mat( n, 3, n, 3, ia, ja, desca, 7, info )
213 CALL chk1mat( n, 3, n, 3, ib, jb, descb, 11, info )
214 IF( info.EQ.0 ) THEN
215 iarow = indxg2p( ia, desca( mb_ ), myrow, desca( rsrc_ ),
216 $ nprow )
217 ibrow = indxg2p( ib, descb( mb_ ), myrow, descb( rsrc_ ),
218 $ nprow )
219 iacol = indxg2p( ja, desca( nb_ ), mycol, desca( csrc_ ),
220 $ npcol )
221 ibcol = indxg2p( jb, descb( nb_ ), mycol, descb( csrc_ ),
222 $ npcol )
223 iroffa = mod( ia-1, desca( mb_ ) )
224 icoffa = mod( ja-1, desca( nb_ ) )
225 iroffb = mod( ib-1, descb( mb_ ) )
226 icoffb = mod( jb-1, descb( nb_ ) )
227 IF( ibtype.LT.1 .OR. ibtype.GT.3 ) THEN
228 info = -1
229 ELSE IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
230 info = -2
231 ELSE IF( n.LT.0 ) THEN
232 info = -3
233 ELSE IF( n+icoffa.GT.desca( nb_ ) ) THEN
234 info = -3
235 ELSE IF( iroffa.NE.0 ) THEN
236 info = -5
237 ELSE IF( icoffa.NE.0 ) THEN
238 info = -6
239 ELSE IF( desca( mb_ ).NE.desca( nb_ ) ) THEN
240 info = -( 700+nb_ )
241 ELSE IF( iroffb.NE.0 .OR. ibrow.NE.iarow ) THEN
242 info = -9
243 ELSE IF( icoffb.NE.0 .OR. ibcol.NE.iacol ) THEN
244 info = -10
245 ELSE IF( descb( mb_ ).NE.desca( mb_ ) ) THEN
246 info = -( 1100+mb_ )
247 ELSE IF( descb( nb_ ).NE.desca( nb_ ) ) THEN
248 info = -( 1100+nb_ )
249 ELSE IF( ictxt.NE.descb( ctxt_ ) ) THEN
250 info = -( 1100+ctxt_ )
251 END IF
252 END IF
253 END IF
254*
255 IF( info.NE.0 ) THEN
256 CALL pxerbla( ictxt, 'PZHEGS2', -info )
257 CALL blacs_exit( ictxt )
258 RETURN
259 END IF
260*
261* Quick return if possible
262*
263 IF( n.EQ.0 .OR. ( myrow.NE.iarow .OR. mycol.NE.iacol ) )
264 $ RETURN
265*
266* Compute local information
267*
268 lda = desca( lld_ )
269 ldb = descb( lld_ )
270 CALL infog2l( ia, ja, desca, nprow, npcol, myrow, mycol, iia, jja,
271 $ iarow, iacol )
272 CALL infog2l( ia, ja, desca, nprow, npcol, myrow, mycol, iib, jjb,
273 $ ibrow, ibcol )
274*
275 IF( ibtype.EQ.1 ) THEN
276*
277 IF( upper ) THEN
278*
279 ioffa = iia + jja*lda
280 ioffb = iib + jjb*ldb
281*
282* Compute inv(U')*sub( A )*inv(U)
283*
284 DO 10 k = 1, n
285*
286* Update the upper triangle of
287* A(ia+k-1:ia+n-a,ia+k-1:ia+n-1)
288*
289 akk = dble( a( ioffa-lda ) )
290 bkk = dble( b( ioffb-ldb ) )
291 akk = akk / bkk**2
292 a( ioffa-lda ) = akk
293 IF( k.LT.n ) THEN
294 CALL zdscal( n-k, one / bkk, a( ioffa ), lda )
295 ct = -half*akk
296 CALL zlacgv( n-k, a( ioffa ), lda )
297 CALL zlacgv( n-k, b( ioffb ), ldb )
298 CALL zaxpy( n-k, ct, b( ioffb ), ldb, a( ioffa ),
299 $ lda )
300 CALL zher2( uplo, n-k, -cone, a( ioffa ), lda,
301 $ b( ioffb ), ldb, a( ioffa+1 ), lda )
302 CALL zaxpy( n-k, ct, b( ioffb ), ldb, a( ioffa ),
303 $ lda )
304 CALL zlacgv( n-k, b( ioffb ), ldb )
305 CALL ztrsv( uplo, 'Conjugate transpose', 'Non-unit',
306 $ n-k, b( ioffb+1 ), ldb, a( ioffa ), lda )
307 CALL zlacgv( n-k, a( ioffa ), lda )
308 END IF
309*
310* A( IOFFA ) -> A( K, K+1 )
311* B( IOFFB ) -> B( K, K+1 )
312*
313 ioffa = ioffa + lda + 1
314 ioffb = ioffb + ldb + 1
315*
316 10 CONTINUE
317*
318 ELSE
319*
320 ioffa = iia + 1 + ( jja-1 )*lda
321 ioffb = iib + 1 + ( jjb-1 )*ldb
322*
323* Compute inv(L)*sub( A )*inv(L')
324*
325 DO 20 k = 1, n
326*
327* Update the lower triangle of
328* A(ia+k-1:ia+n-a,ia+k-1:ia+n-1)
329*
330 akk = dble( a( ioffa-1 ) )
331 bkk = dble( b( ioffb-1 ) )
332 akk = akk / bkk**2
333 a( ioffa-1 ) = akk
334*
335 IF( k.LT.n ) THEN
336 CALL zdscal( n-k, one / bkk, a( ioffa ), 1 )
337 ct = -half*akk
338 CALL zaxpy( n-k, ct, b( ioffb ), 1, a( ioffa ), 1 )
339 CALL zher2( uplo, n-k, -cone, a( ioffa ), 1,
340 $ b( ioffb ), 1, a( ioffa+lda ), lda )
341 CALL zaxpy( n-k, ct, b( ioffb ), 1, a( ioffa ), 1 )
342 CALL ztrsv( uplo, 'No transpose', 'Non-unit', n-k,
343 $ b( ioffb+ldb ), ldb, a( ioffa ), 1 )
344 END IF
345*
346* A( IOFFA ) -> A( K+1, K )
347* B( IOFFB ) -> B( K+1, K )
348*
349 ioffa = ioffa + lda + 1
350 ioffb = ioffb + ldb + 1
351*
352 20 CONTINUE
353*
354 END IF
355*
356 ELSE
357*
358 IF( upper ) THEN
359*
360 ioffa = iia + ( jja-1 )*lda
361 ioffb = iib + ( jjb-1 )*ldb
362*
363* Compute U*sub( A )*U'
364*
365 DO 30 k = 1, n
366*
367* Update the upper triangle of A(ia:ia+k-1,ja:ja+k-1)
368*
369 akk = dble( a( ioffa+k-1 ) )
370 bkk = dble( b( ioffb+k-1 ) )
371 CALL ztrmv( uplo, 'No transpose', 'Non-unit', k-1,
372 $ b( iib+( jjb-1 )*ldb ), ldb, a( ioffa ), 1 )
373 ct = half*akk
374 CALL zaxpy( k-1, ct, b( ioffb ), 1, a( ioffa ), 1 )
375 CALL zher2( uplo, k-1, cone, a( ioffa ), 1, b( ioffb ),
376 $ 1, a( iia+( jja-1 )*lda ), lda )
377 CALL zaxpy( k-1, ct, b( ioffb ), 1, a( ioffa ), 1 )
378 CALL zdscal( k-1, bkk, a( ioffa ), 1 )
379 a( ioffa+k-1 ) = akk*bkk**2
380*
381* A( IOFFA ) -> A( 1, K )
382* B( IOFFB ) -> B( 1, K )
383*
384 ioffa = ioffa + lda
385 ioffb = ioffb + ldb
386*
387 30 CONTINUE
388*
389 ELSE
390*
391 ioffa = iia + ( jja-1 )*lda
392 ioffb = iib + ( jjb-1 )*ldb
393*
394* Compute L'*sub( A )*L
395*
396 DO 40 k = 1, n
397*
398* Update the lower triangle of A(ia:ia+k-1,ja:ja+k-1)
399*
400 akk = dble( a( ioffa+( k-1 )*lda ) )
401 bkk = dble( b( ioffb+( k-1 )*ldb ) )
402 CALL zlacgv( k-1, a( ioffa ), lda )
403 CALL ztrmv( uplo, 'Conjugate transpose', 'Non-unit', k-1,
404 $ b( iib+( jjb-1 )*ldb ), ldb, a( ioffa ),
405 $ lda )
406 ct = half*akk
407 CALL zlacgv( k-1, b( ioffb ), ldb )
408 CALL zaxpy( k-1, ct, b( ioffb ), ldb, a( ioffa ), lda )
409 CALL zher2( uplo, k-1, cone, a( ioffa ), lda, b( ioffb ),
410 $ ldb, a( iia+( jja-1 )*lda ), lda )
411 CALL zaxpy( k-1, ct, b( ioffb ), ldb, a( ioffa ), lda )
412 CALL zlacgv( k-1, b( ioffb ), ldb )
413 CALL zdscal( k-1, bkk, a( ioffa ), lda )
414 CALL zlacgv( k-1, a( ioffa ), lda )
415 a( ioffa+( k-1 )*lda ) = akk*bkk**2
416*
417* A( IOFFA ) -> A( K, 1 )
418* B( IOFFB ) -> B( K, 1 )
419*
420 ioffa = ioffa + 1
421 ioffb = ioffb + 1
422*
423 40 CONTINUE
424*
425 END IF
426*
427 END IF
428*
429 RETURN
430*
431* End of PZHEGS2
432*
chk1mat
subroutine chk1mat(ma, mapos0, na, napos0, ia, ja, desca, descapos0, info)
Definition chk1mat.f:3
indxg2p
integer function indxg2p(indxglob, nb, iproc, isrcproc, nprocs)
Definition indxg2p.f:2
infog2l
subroutine infog2l(grindx, gcindx, desc, nprow, npcol, myrow, mycol, lrindx, lcindx, rsrc, csrc)
Definition infog2l.f:3
pxerbla
subroutine pxerbla(ictxt, srname, info)
Definition pxerbla.f:2
lsame
logical function lsame(ca, cb)
Definition tools.f:1724
Here is the call graph for this function:
Here is the caller graph for this function: