ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pdsygs2.f
Go to the documentation of this file.
1 *
2 *
3  SUBROUTINE pdsygs2( 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  DOUBLE PRECISION A( * ), B( * )
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 * PDSYGS2 reduces a real symmetric-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**T)*sub( A )*inv(U) or
31 * inv(L)*sub( A )*inv(L**T)
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**T or L**T*sub( A )*L.
36 *
37 * sub( B ) must have been previously factorized as U**T*U or L*L**T by
38 * PDPOTRF.
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**T)*sub( A )*inv(U) or
99 * inv(L)*sub( A )*inv(L**T);
100 * = 2 or 3: compute U*sub( A )*U**T or L**T*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**T*U;
105 * = 'L': Lower triangle of sub( A ) is stored and sub( B ) is
106 * factored as L*L**T.
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) DOUBLE PRECISION 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 symmetric 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) DOUBLE PRECISION 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 * PDPOTRF.
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 * ..
171 * .. Local Scalars ..
172  LOGICAL UPPER
173  INTEGER IACOL, IAROW, IBCOL, IBROW, ICOFFA, ICOFFB,
174  $ ictxt, iia, iib, ioffa, ioffb, iroffa, iroffb,
175  $ jja, jjb, k, lda, ldb, mycol, myrow, npcol,
176  $ nprow
177  DOUBLE PRECISION AKK, BKK, CT
178 * ..
179 * .. External Subroutines ..
180  EXTERNAL blacs_exit, blacs_gridinfo, chk1mat, daxpy,
181  $ dscal, dsyr2, dtrmv, dtrsv, infog2l, pxerbla
182 * ..
183 * .. Intrinsic Functions ..
184  INTRINSIC mod
185 * ..
186 * .. External Functions ..
187  LOGICAL LSAME
188  INTEGER INDXG2P
189  EXTERNAL lsame, indxg2p
190 * ..
191 * .. Executable Statements ..
192 * This is just to keep ftnchek happy
193  IF( block_cyclic_2d*csrc_*ctxt_*dlen_*dtype_*lld_*mb_*m_*nb_*n_*
194  $ rsrc_.LT.0 )RETURN
195 *
196 * Get grid parameters
197 *
198  ictxt = desca( ctxt_ )
199  CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
200 *
201 * Test the input parameters.
202 *
203  info = 0
204  IF( nprow.EQ.-1 ) THEN
205  info = -( 700+ctxt_ )
206  ELSE
207  upper = lsame( uplo, 'U' )
208  CALL chk1mat( n, 3, n, 3, ia, ja, desca, 7, info )
209  CALL chk1mat( n, 3, n, 3, ib, jb, descb, 11, info )
210  IF( info.EQ.0 ) THEN
211  iarow = indxg2p( ia, desca( mb_ ), myrow, desca( rsrc_ ),
212  $ nprow )
213  ibrow = indxg2p( ib, descb( mb_ ), myrow, descb( rsrc_ ),
214  $ nprow )
215  iacol = indxg2p( ja, desca( nb_ ), mycol, desca( csrc_ ),
216  $ npcol )
217  ibcol = indxg2p( jb, descb( nb_ ), mycol, descb( csrc_ ),
218  $ npcol )
219  iroffa = mod( ia-1, desca( mb_ ) )
220  icoffa = mod( ja-1, desca( nb_ ) )
221  iroffb = mod( ib-1, descb( mb_ ) )
222  icoffb = mod( jb-1, descb( nb_ ) )
223  IF( ibtype.LT.1 .OR. ibtype.GT.3 ) THEN
224  info = -1
225  ELSE IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
226  info = -2
227  ELSE IF( n.LT.0 ) THEN
228  info = -3
229  ELSE IF( n+icoffa.GT.desca( nb_ ) ) THEN
230  info = -3
231  ELSE IF( iroffa.NE.0 ) THEN
232  info = -5
233  ELSE IF( icoffa.NE.0 ) THEN
234  info = -6
235  ELSE IF( desca( mb_ ).NE.desca( nb_ ) ) THEN
236  info = -( 700+nb_ )
237  ELSE IF( iroffb.NE.0 .OR. ibrow.NE.iarow ) THEN
238  info = -9
239  ELSE IF( icoffb.NE.0 .OR. ibcol.NE.iacol ) THEN
240  info = -10
241  ELSE IF( descb( mb_ ).NE.desca( mb_ ) ) THEN
242  info = -( 1100+mb_ )
243  ELSE IF( descb( nb_ ).NE.desca( nb_ ) ) THEN
244  info = -( 1100+nb_ )
245  ELSE IF( ictxt.NE.descb( ctxt_ ) ) THEN
246  info = -( 1100+ctxt_ )
247  END IF
248  END IF
249  END IF
250 *
251  IF( info.NE.0 ) THEN
252  CALL pxerbla( ictxt, 'PDSYGS2', -info )
253  CALL blacs_exit( ictxt )
254  RETURN
255  END IF
256 *
257 * Quick return if possible
258 *
259  IF( n.EQ.0 .OR. ( myrow.NE.iarow .OR. mycol.NE.iacol ) )
260  $ RETURN
261 *
262 * Compute local information
263 *
264  lda = desca( lld_ )
265  ldb = descb( lld_ )
266  CALL infog2l( ia, ja, desca, nprow, npcol, myrow, mycol, iia, jja,
267  $ iarow, iacol )
268  CALL infog2l( ia, ja, desca, nprow, npcol, myrow, mycol, iib, jjb,
269  $ ibrow, ibcol )
270 *
271  IF( ibtype.EQ.1 ) THEN
272 *
273  IF( upper ) THEN
274 *
275  ioffa = iia + jja*lda
276  ioffb = iib + jjb*ldb
277 *
278 * Compute inv(U')*sub( A )*inv(U)
279 *
280  DO 10 k = 1, n
281 *
282 * Update the upper triangle of
283 * A(ia+k-1:ia+n-a,ia+k-1:ia+n-1)
284 *
285  akk = a( ioffa-lda )
286  bkk = b( ioffb-ldb )
287  akk = akk / bkk**2
288  a( ioffa-lda ) = akk
289  IF( k.LT.n ) THEN
290  CALL dscal( n-k, one / bkk, a( ioffa ), lda )
291  ct = -half*akk
292  CALL daxpy( n-k, ct, b( ioffb ), ldb, a( ioffa ),
293  $ lda )
294  CALL dsyr2( uplo, n-k, -one, a( ioffa ), lda,
295  $ b( ioffb ), ldb, a( ioffa+1 ), lda )
296  CALL daxpy( n-k, ct, b( ioffb ), ldb, a( ioffa ),
297  $ lda )
298  CALL dtrsv( uplo, 'Transpose', 'Non-unit', n-k,
299  $ b( ioffb+1 ), ldb, a( ioffa ), lda )
300  END IF
301 *
302 * A( IOFFA ) -> A( K, K+1 )
303 * B( IOFFB ) -> B( K, K+1 )
304 *
305  ioffa = ioffa + lda + 1
306  ioffb = ioffb + ldb + 1
307 *
308  10 CONTINUE
309 *
310  ELSE
311 *
312  ioffa = iia + 1 + ( jja-1 )*lda
313  ioffb = iib + 1 + ( jjb-1 )*ldb
314 *
315 * Compute inv(L)*sub( A )*inv(L')
316 *
317  DO 20 k = 1, n
318 *
319 * Update the lower triangle of
320 * A(ia+k-1:ia+n-a,ia+k-1:ia+n-1)
321 *
322  akk = a( ioffa-1 )
323  bkk = b( ioffb-1 )
324  akk = akk / bkk**2
325  a( ioffa-1 ) = akk
326 *
327  IF( k.LT.n ) THEN
328  CALL dscal( n-k, one / bkk, a( ioffa ), 1 )
329  ct = -half*akk
330  CALL daxpy( n-k, ct, b( ioffb ), 1, a( ioffa ), 1 )
331  CALL dsyr2( uplo, n-k, -one, a( ioffa ), 1,
332  $ b( ioffb ), 1, a( ioffa+lda ), lda )
333  CALL daxpy( n-k, ct, b( ioffb ), 1, a( ioffa ), 1 )
334  CALL dtrsv( uplo, 'No transpose', 'Non-unit', n-k,
335  $ b( ioffb+ldb ), ldb, a( ioffa ), 1 )
336  END IF
337 *
338 * A( IOFFA ) -> A( K+1, K )
339 * B( IOFFB ) -> B( K+1, K )
340 *
341  ioffa = ioffa + lda + 1
342  ioffb = ioffb + ldb + 1
343 *
344  20 CONTINUE
345 *
346  END IF
347 *
348  ELSE
349 *
350  IF( upper ) THEN
351 *
352  ioffa = iia + ( jja-1 )*lda
353  ioffb = iib + ( jjb-1 )*ldb
354 *
355 * Compute U*sub( A )*U'
356 *
357  DO 30 k = 1, n
358 *
359 * Update the upper triangle of A(ia:ia+k-1,ja:ja+k-1)
360 *
361  akk = a( ioffa+k-1 )
362  bkk = b( ioffb+k-1 )
363  CALL dtrmv( uplo, 'No transpose', 'Non-unit', k-1,
364  $ b( iib+( jjb-1 )*ldb ), ldb, a( ioffa ), 1 )
365  ct = half*akk
366  CALL daxpy( k-1, ct, b( ioffb ), 1, a( ioffa ), 1 )
367  CALL dsyr2( uplo, k-1, one, a( ioffa ), 1, b( ioffb ), 1,
368  $ a( iia+( jja-1 )*lda ), lda )
369  CALL daxpy( k-1, ct, b( ioffb ), 1, a( ioffa ), 1 )
370  CALL dscal( k-1, bkk, a( ioffa ), 1 )
371  a( ioffa+k-1 ) = akk*bkk**2
372 *
373 * A( IOFFA ) -> A( 1, K )
374 * B( IOFFB ) -> B( 1, K )
375 *
376  ioffa = ioffa + lda
377  ioffb = ioffb + ldb
378 *
379  30 CONTINUE
380 *
381  ELSE
382 *
383  ioffa = iia + ( jja-1 )*lda
384  ioffb = iib + ( jjb-1 )*ldb
385 *
386 * Compute L'*sub( A )*L
387 *
388  DO 40 k = 1, n
389 *
390 * Update the lower triangle of A(ia:ia+k-1,ja:ja+k-1)
391 *
392  akk = a( ioffa+( k-1 )*lda )
393  bkk = b( ioffb+( k-1 )*ldb )
394  CALL dtrmv( uplo, 'Transpose', 'Non-unit', k-1,
395  $ b( iib+( jjb-1 )*ldb ), ldb, a( ioffa ),
396  $ lda )
397  ct = half*akk
398  CALL daxpy( k-1, ct, b( ioffb ), ldb, a( ioffa ), lda )
399  CALL dsyr2( uplo, k-1, one, a( ioffa ), lda, b( ioffb ),
400  $ ldb, a( iia+( jja-1 )*lda ), lda )
401  CALL daxpy( k-1, ct, b( ioffb ), ldb, a( ioffa ), lda )
402  CALL dscal( k-1, bkk, a( ioffa ), lda )
403  a( ioffa+( k-1 )*lda ) = akk*bkk**2
404 *
405 * A( IOFFA ) -> A( K, 1 )
406 * B( IOFFB ) -> B( K, 1 )
407 *
408  ioffa = ioffa + 1
409  ioffb = ioffb + 1
410 *
411  40 CONTINUE
412 *
413  END IF
414 *
415  END IF
416 *
417  RETURN
418 *
419 * End of PDSYGS2
420 *
421  END
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
pdsygs2
subroutine pdsygs2(IBTYPE, UPLO, N, A, IA, JA, DESCA, B, IB, JB, DESCB, INFO)
Definition: pdsygs2.f:5