ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pslaed3.f
Go to the documentation of this file.
1  SUBROUTINE pslaed3( ICTXT, K, N, NB, D, DROW, DCOL, RHO, DLAMDA,
2  $ W, Z, U, LDU, BUF, INDX, INDCOL, INDROW,
3  $ INDXR, INDXC, CTOT, NPCOL, INFO )
4 *
5 * -- ScaLAPACK auxiliary routine (version 1.7) --
6 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
7 * and University of California, Berkeley.
8 * December 31, 1998
9 *
10 * .. Scalar Arguments ..
11  INTEGER DCOL, DROW, ICTXT, INFO, K, LDU, N, NB, NPCOL
12  REAL RHO
13 * ..
14 * .. Array Arguments ..
15  INTEGER CTOT( 0: NPCOL-1, 4 ), INDCOL( * ),
16  $ INDROW( * ), INDX( * ), INDXC( * ), INDXR( * )
17  REAL BUF( * ), D( * ), DLAMDA( * ), U( LDU, * ),
18  $ W( * ), Z( * )
19 * ..
20 *
21 * Purpose
22 * =======
23 *
24 * PSLAED3 finds the roots of the secular equation, as defined by the
25 * values in D, W, and RHO, between 1 and K. It makes the
26 * appropriate calls to SLAED4
27 *
28 * This code makes very mild assumptions about floating point
29 * arithmetic. It will work on machines with a guard digit in
30 * add/subtract, or on those binary machines without guard digits
31 * which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
32 * It could conceivably fail on hexadecimal or decimal machines
33 * without guard digits, but we know of none.
34 *
35 * Arguments
36 * =========
37 *
38 * ICTXT (global input) INTEGER
39 * The BLACS context handle, indicating the global context of
40 * the operation on the matrix. The context itself is global.
41 *
42 * K (output) INTEGER
43 * The number of non-deflated eigenvalues, and the order of the
44 * related secular equation. 0 <= K <=N.
45 *
46 * N (input) INTEGER
47 * The dimension of the symmetric tridiagonal matrix. N >= 0.
48 *
49 * NB (global input) INTEGER
50 * The blocking factor used to distribute the columns of the
51 * matrix. NB >= 1.
52 *
53 * D (input/output) REAL array, dimension (N)
54 * On entry, D contains the eigenvalues of the two submatrices to
55 * be combined.
56 * On exit, D contains the trailing (N-K) updated eigenvalues
57 * (those which were deflated) sorted into increasing order.
58 *
59 * DROW (global input) INTEGER
60 * The process row over which the first row of the matrix D is
61 * distributed. 0 <= DROW < NPROW.
62 *
63 * DCOL (global input) INTEGER
64 * The process column over which the first column of the
65 * matrix D is distributed. 0 <= DCOL < NPCOL.
66 *
67 * RHO (global input/output) REAL
68 * On entry, the off-diagonal element associated with the rank-1
69 * cut which originally split the two submatrices which are now
70 * being recombined.
71 * On exit, RHO has been modified to the value required by
72 * PSLAED3.
73 *
74 * DLAMDA (global output) REAL array, dimension (N)
75 * A copy of the first K eigenvalues which will be used by
76 * SLAED4 to form the secular equation.
77 *
78 * W (global output) REAL array, dimension (N)
79 * The first k values of the final deflation-altered z-vector
80 * which will be passed to SLAED4.
81 *
82 * Z (global input) REAL array, dimension (N)
83 * On entry, Z contains the updating vector (the last
84 * row of the first sub-eigenvector matrix and the first row of
85 * the second sub-eigenvector matrix).
86 * On exit, the contents of Z have been destroyed by the updating
87 * process.
88 *
89 * U (global output) REAL array
90 * global dimension (N, N), local dimension (LDU, NQ).
91 * (See PSLAED0 for definition of NQ.)
92 * Q contains the orthonormal eigenvectors of the symmetric
93 * tridiagonal matrix.
94 *
95 * LDU (input) INTEGER
96 * The leading dimension of the array U.
97 *
98 * BUF (workspace) REAL array, dimension 3*N
99 *
100 *
101 * INDX (workspace) INTEGER array, dimension (N)
102 * The permutation used to sort the contents of DLAMDA into
103 * ascending order.
104 *
105 * INDCOL (workspace) INTEGER array, dimension (N)
106 *
107 *
108 * INDROW (workspace) INTEGER array, dimension (N)
109 *
110 *
111 * INDXR (workspace) INTEGER array, dimension (N)
112 *
113 *
114 * INDXC (workspace) INTEGER array, dimension (N)
115 *
116 * CTOT (workspace) INTEGER array, dimension( NPCOL, 4)
117 *
118 * NPCOL (global input) INTEGER
119 * The total number of columns over which the distributed
120 * submatrix is distributed.
121 *
122 * INFO (output) INTEGER
123 * = 0: successful exit.
124 * < 0: if INFO = -i, the i-th argument had an illegal value.
125 * > 0: The algorithm failed to compute the ith eigenvalue.
126 *
127 * =====================================================================
128 *
129 * .. Parameters ..
130  REAL ONE
131  PARAMETER ( ONE = 1.0e0 )
132 * ..
133 * .. Local Scalars ..
134  INTEGER COL, GI, I, IINFO, IIU, IPD, IU, J, JJU, JU,
135  $ KK, KL, KLC, KLR, MYCOL, MYKL, MYKLR, MYROW,
136  $ nprow, pdc, pdr, row
137  REAL AUX, TEMP
138 * ..
139 * .. External Functions ..
140  INTEGER INDXG2L
141  REAL SLAMC3, SNRM2
142  EXTERNAL indxg2l, slamc3, snrm2
143 * ..
144 * .. External Subroutines ..
145  EXTERNAL blacs_gridinfo, scopy, sgebr2d, sgebs2d,
146  $ sgerv2d, sgesd2d, slaed4
147 * ..
148 * .. Intrinsic Functions ..
149  INTRINSIC mod, sign, sqrt
150 * ..
151 * .. Executable Statements ..
152 *
153 * Test the input parameters.
154 *
155  iinfo = 0
156 *
157 * Quick return if possible
158 *
159  IF( k.EQ.0 )
160  $ RETURN
161 *
162  CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
163 *
164  row = drow
165  col = dcol
166  DO 20 i = 1, n, nb
167  DO 10 j = 0, nb - 1
168  indrow( i+j ) = row
169  indcol( i+j ) = col
170  10 CONTINUE
171  row = mod( row+1, nprow )
172  col = mod( col+1, npcol )
173  20 CONTINUE
174 *
175  mykl = ctot( mycol, 1 ) + ctot( mycol, 2 ) + ctot( mycol, 3 )
176  klr = mykl / nprow
177  IF( myrow.EQ.drow ) THEN
178  myklr = klr + mod( mykl, nprow )
179  ELSE
180  myklr = klr
181  END IF
182  pdc = 1
183  col = dcol
184  30 CONTINUE
185  IF( mycol.NE.col ) THEN
186  pdc = pdc + ctot( col, 1 ) + ctot( col, 2 ) + ctot( col, 3 )
187  col = mod( col+1, npcol )
188  GO TO 30
189  END IF
190  pdr = pdc
191  kl = klr + mod( mykl, nprow )
192  row = drow
193  40 CONTINUE
194  IF( myrow.NE.row ) THEN
195  pdr = pdr + kl
196  kl = klr
197  row = mod( row+1, nprow )
198  GO TO 40
199  END IF
200 *
201  DO 50 i = 1, k
202  dlamda( i ) = slamc3( dlamda( i ), dlamda( i ) ) - dlamda( i )
203  z( i ) = one
204  50 CONTINUE
205  IF( myklr.GT.0 ) THEN
206  kk = pdr
207  DO 80 i = 1, myklr
208  CALL slaed4( k, kk, dlamda, w, buf, rho, buf( k+i ), iinfo )
209  IF( iinfo.NE.0 ) THEN
210  info = kk
211  END IF
212 *
213 * ..Compute part of z
214 *
215  DO 60 j = 1, kk - 1
216  z( j ) = z( j )*( buf( j ) /
217  $ ( dlamda( j )-dlamda( kk ) ) )
218  60 CONTINUE
219  z( kk ) = z( kk )*buf( kk )
220  DO 70 j = kk + 1, k
221  z( j ) = z( j )*( buf( j ) /
222  $ ( dlamda( j )-dlamda( kk ) ) )
223  70 CONTINUE
224  kk = kk + 1
225  80 CONTINUE
226 *
227  IF( myrow.NE.drow ) THEN
228  CALL scopy( k, z, 1, buf, 1 )
229  CALL sgesd2d( ictxt, k+myklr, 1, buf, k+myklr, drow, mycol )
230  ELSE
231  ipd = 2*k + 1
232  CALL scopy( myklr, buf( k+1 ), 1, buf( ipd ), 1 )
233  IF( klr.GT.0 ) THEN
234  ipd = myklr + ipd
235  row = mod( drow+1, nprow )
236  DO 100 i = 1, nprow - 1
237  CALL sgerv2d( ictxt, k+klr, 1, buf, k+klr, row,
238  $ mycol )
239  CALL scopy( klr, buf( k+1 ), 1, buf( ipd ), 1 )
240  DO 90 j = 1, k
241  z( j ) = z( j )*buf( j )
242  90 CONTINUE
243  ipd = ipd + klr
244  row = mod( row+1, nprow )
245  100 CONTINUE
246  END IF
247  END IF
248  END IF
249 *
250  IF( myrow.EQ.drow ) THEN
251  IF( mycol.NE.dcol .AND. mykl.NE.0 ) THEN
252  CALL scopy( k, z, 1, buf, 1 )
253  CALL scopy( mykl, buf( 2*k+1 ), 1, buf( k+1 ), 1 )
254  CALL sgesd2d( ictxt, k+mykl, 1, buf, k+mykl, myrow, dcol )
255  ELSE IF( mycol.EQ.dcol ) THEN
256  ipd = 2*k + 1
257  col = dcol
258  kl = mykl
259  DO 120 i = 1, npcol - 1
260  ipd = ipd + kl
261  col = mod( col+1, npcol )
262  kl = ctot( col, 1 ) + ctot( col, 2 ) + ctot( col, 3 )
263  IF( kl.NE.0 ) THEN
264  CALL sgerv2d( ictxt, k+kl, 1, buf, k+kl, myrow, col )
265  CALL scopy( kl, buf( k+1 ), 1, buf( ipd ), 1 )
266  DO 110 j = 1, k
267  z( j ) = z( j )*buf( j )
268  110 CONTINUE
269  END IF
270  120 CONTINUE
271  DO 130 i = 1, k
272  z( i ) = sign( sqrt( -z( i ) ), w( i ) )
273  130 CONTINUE
274 *
275  END IF
276  END IF
277 *
278 * Diffusion
279 *
280  IF( myrow.EQ.drow .AND. mycol.EQ.dcol ) THEN
281  CALL scopy( k, z, 1, buf, 1 )
282  CALL scopy( k, buf( 2*k+1 ), 1, buf( k+1 ), 1 )
283  CALL sgebs2d( ictxt, 'All', ' ', 2*k, 1, buf, 2*k )
284  ELSE
285  CALL sgebr2d( ictxt, 'All', ' ', 2*k, 1, buf, 2*k, drow, dcol )
286  CALL scopy( k, buf, 1, z, 1 )
287  END IF
288 *
289 * Copy of D at the good place
290 *
291  klc = 0
292  klr = 0
293  DO 140 i = 1, k
294  gi = indx( i )
295  d( gi ) = buf( k+i )
296  col = indcol( gi )
297  row = indrow( gi )
298  IF( col.EQ.mycol ) THEN
299  klc = klc + 1
300  indxc( klc ) = i
301  END IF
302  IF( row.EQ.myrow ) THEN
303  klr = klr + 1
304  indxr( klr ) = i
305  END IF
306  140 CONTINUE
307 *
308 * Compute eigenvectors of the modified rank-1 modification.
309 *
310  IF( mykl.NE.0 ) THEN
311  DO 180 j = 1, mykl
312  kk = indxc( j )
313  ju = indx( kk )
314  jju = indxg2l( ju, nb, j, j, npcol )
315  CALL slaed4( k, kk, dlamda, w, buf, rho, aux, iinfo )
316  IF( iinfo.NE.0 ) THEN
317  info = kk
318  END IF
319  IF( k.EQ.1 .OR. k.EQ.2 ) THEN
320  DO 150 i = 1, klr
321  kk = indxr( i )
322  iu = indx( kk )
323  iiu = indxg2l( iu, nb, j, j, nprow )
324  u( iiu, jju ) = buf( kk )
325  150 CONTINUE
326  GO TO 180
327  END IF
328 *
329  DO 160 i = 1, k
330  buf( i ) = z( i ) / buf( i )
331  160 CONTINUE
332  temp = snrm2( k, buf, 1 )
333  DO 170 i = 1, klr
334  kk = indxr( i )
335  iu = indx( kk )
336  iiu = indxg2l( iu, nb, j, j, nprow )
337  u( iiu, jju ) = buf( kk ) / temp
338  170 CONTINUE
339  180 CONTINUE
340  END IF
341 *
342  190 CONTINUE
343  RETURN
344 *
345 * End of PSLAED3
346 *
347  END
pslaed3
subroutine pslaed3(ICTXT, K, N, NB, D, DROW, DCOL, RHO, DLAMDA, W, Z, U, LDU, BUF, INDX, INDCOL, INDROW, INDXR, INDXC, CTOT, NPCOL, INFO)
Definition: pslaed3.f:4