ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pchegs2.f
Go to the documentation of this file.
1 *
2 *
3  SUBROUTINE pchegs2( IBTYPE, UPLO, N, A, IA, JA, DESCA, B, IB, JB,
4  $ DESCB, INFO )
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 A( * ), B( * )
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 * PCHEGS2 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 * PCPOTRF.
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 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 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 * PCPOTRF.
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  REAL ONE, HALF
169  parameter( one = 1.0e+0, half = 0.5e+0 )
170  COMPLEX CONE
171  parameter( cone = ( 1.0e+0, 0.0e+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  REAL AKK, BKK
180  COMPLEX CT
181 * ..
182 * .. External Subroutines ..
183  EXTERNAL blacs_exit, blacs_gridinfo, caxpy, cher2,
184  $ chk1mat, clacgv, csscal, ctrmv, ctrsv, infog2l,
185  $ pxerbla
186 * ..
187 * .. Intrinsic Functions ..
188  INTRINSIC mod, real
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, 'PCHEGS2', -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 = real( a( ioffa-lda ) )
290  bkk = real( b( ioffb-ldb ) )
291  akk = akk / bkk**2
292  a( ioffa-lda ) = akk
293  IF( k.LT.n ) THEN
294  CALL csscal( n-k, one / bkk, a( ioffa ), lda )
295  ct = -half*akk
296  CALL clacgv( n-k, a( ioffa ), lda )
297  CALL clacgv( n-k, b( ioffb ), ldb )
298  CALL caxpy( n-k, ct, b( ioffb ), ldb, a( ioffa ),
299  $ lda )
300  CALL cher2( uplo, n-k, -cone, a( ioffa ), lda,
301  $ b( ioffb ), ldb, a( ioffa+1 ), lda )
302  CALL caxpy( n-k, ct, b( ioffb ), ldb, a( ioffa ),
303  $ lda )
304  CALL clacgv( n-k, b( ioffb ), ldb )
305  CALL ctrsv( uplo, 'Conjugate transpose', 'Non-unit',
306  $ n-k, b( ioffb+1 ), ldb, a( ioffa ), lda )
307  CALL clacgv( 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 = real( a( ioffa-1 ) )
331  bkk = real( b( ioffb-1 ) )
332  akk = akk / bkk**2
333  a( ioffa-1 ) = akk
334 *
335  IF( k.LT.n ) THEN
336  CALL csscal( n-k, one / bkk, a( ioffa ), 1 )
337  ct = -half*akk
338  CALL caxpy( n-k, ct, b( ioffb ), 1, a( ioffa ), 1 )
339  CALL cher2( uplo, n-k, -cone, a( ioffa ), 1,
340  $ b( ioffb ), 1, a( ioffa+lda ), lda )
341  CALL caxpy( n-k, ct, b( ioffb ), 1, a( ioffa ), 1 )
342  CALL ctrsv( 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 = real( a( ioffa+k-1 ) )
370  bkk = real( b( ioffb+k-1 ) )
371  CALL ctrmv( uplo, 'No transpose', 'Non-unit', k-1,
372  $ b( iib+( jjb-1 )*ldb ), ldb, a( ioffa ), 1 )
373  ct = half*akk
374  CALL caxpy( k-1, ct, b( ioffb ), 1, a( ioffa ), 1 )
375  CALL cher2( uplo, k-1, cone, a( ioffa ), 1, b( ioffb ),
376  $ 1, a( iia+( jja-1 )*lda ), lda )
377  CALL caxpy( k-1, ct, b( ioffb ), 1, a( ioffa ), 1 )
378  CALL csscal( 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 = real( a( ioffa+( k-1 )*lda ) )
401  bkk = real( b( ioffb+( k-1 )*ldb ) )
402  CALL clacgv( k-1, a( ioffa ), lda )
403  CALL ctrmv( uplo, 'Conjugate transpose', 'Non-unit', k-1,
404  $ b( iib+( jjb-1 )*ldb ), ldb, a( ioffa ),
405  $ lda )
406  ct = half*akk
407  CALL clacgv( k-1, b( ioffb ), ldb )
408  CALL caxpy( k-1, ct, b( ioffb ), ldb, a( ioffa ), lda )
409  CALL cher2( uplo, k-1, cone, a( ioffa ), lda, b( ioffb ),
410  $ ldb, a( iia+( jja-1 )*lda ), lda )
411  CALL caxpy( k-1, ct, b( ioffb ), ldb, a( ioffa ), lda )
412  CALL clacgv( k-1, b( ioffb ), ldb )
413  CALL csscal( k-1, bkk, a( ioffa ), lda )
414  CALL clacgv( 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 PCHEGS2
432 *
433  END
pchegs2
subroutine pchegs2(IBTYPE, UPLO, N, A, IA, JA, DESCA, B, IB, JB, DESCB, INFO)
Definition: pchegs2.f:5
infog2l
subroutine infog2l(GRINDX, GCINDX, DESC, NPROW, NPCOL, MYROW, MYCOL, LRINDX, LCINDX, RSRC, CSRC)
Definition: infog2l.f:3
chk1mat
subroutine chk1mat(MA, MAPOS0, NA, NAPOS0, IA, JA, DESCA, DESCAPOS0, INFO)
Definition: chk1mat.f:3
pxerbla
subroutine pxerbla(ICTXT, SRNAME, INFO)
Definition: pxerbla.f:2