ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pdlaed1.f
Go to the documentation of this file.
1  SUBROUTINE pdlaed1( N, N1, D, ID, Q, IQ, JQ, DESCQ, RHO, WORK,
2  $ IWORK, INFO )
3 *
4 * -- ScaLAPACK auxiliary routine (version 1.7) --
5 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
6 * and University of California, Berkeley.
7 * December 31, 1998
8 *
9 * .. Scalar Arguments ..
10  INTEGER ID, INFO, IQ, JQ, N, N1
11  DOUBLE PRECISION RHO
12 * ..
13 * .. Array Arguments ..
14  INTEGER DESCQ( * ), IWORK( * )
15  DOUBLE PRECISION D( * ), Q( * ), WORK( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * PDLAED1 computes the updated eigensystem of a diagonal
22 * matrix after modification by a rank-one symmetric matrix,
23 * in parallel.
24 *
25 * T = Q(in) ( D(in) + RHO * Z*Z' ) Q'(in) = Q(out) * D(out) * Q'(out)
26 *
27 * where Z = Q'u, u is a vector of length N with ones in the
28 * N1 and N1 + 1 th elements and zeros elsewhere.
29 *
30 * The eigenvectors of the original matrix are stored in Q, and the
31 * eigenvalues are in D. The algorithm consists of three stages:
32 *
33 * The first stage consists of deflating the size of the problem
34 * when there are multiple eigenvalues or if there is a zero in
35 * the Z vector. For each such occurence the dimension of the
36 * secular equation problem is reduced by one. This stage is
37 * performed by the routine PDLAED2.
38 *
39 * The second stage consists of calculating the updated
40 * eigenvalues. This is done by finding the roots of the secular
41 * equation via the routine SLAED4 (as called by PDLAED3).
42 * This routine also calculates the eigenvectors of the current
43 * problem.
44 *
45 * The final stage consists of computing the updated eigenvectors
46 * directly using the updated eigenvalues. The eigenvectors for
47 * the current problem are multiplied with the eigenvectors from
48 * the overall problem.
49 *
50 * Arguments
51 * =========
52 *
53 * N (global input) INTEGER
54 * The order of the tridiagonal matrix T. N >= 0.
55 *
56 *
57 * N1 (input) INTEGER
58 * The location of the last eigenvalue in the leading
59 * sub-matrix.
60 * min(1,N) <= N1 <= N.
61 *
62 * D (global input/output) DOUBLE PRECISION array, dimension (N)
63 * On entry,the eigenvalues of the rank-1-perturbed matrix.
64 * On exit, the eigenvalues of the repaired matrix.
65 *
66 * ID (global input) INTEGER
67 * Q's global row/col index, which points to the beginning
68 * of the submatrix which is to be operated on.
69 *
70 * Q (local output) DOUBLE PRECISION array,
71 * global dimension (N, N),
72 * local dimension ( LLD_Q, LOCc(JQ+N-1))
73 * Q contains the orthonormal eigenvectors of the symmetric
74 * tridiagonal matrix.
75 *
76 * IQ (global input) INTEGER
77 * Q's global row index, which points to the beginning of the
78 * submatrix which is to be operated on.
79 *
80 * JQ (global input) INTEGER
81 * Q's global column index, which points to the beginning of
82 * the submatrix which is to be operated on.
83 *
84 * DESCQ (global and local input) INTEGER array of dimension DLEN_.
85 * The array descriptor for the distributed matrix Z.
86 *
87 * RHO (input) DOUBLE PRECISION
88 * The subdiagonal entry used to create the rank-1 modification.
89 *
90 * WORK (local workspace/output) DOUBLE PRECISION array,
91 * dimension 6*N + 2*NP*NQ
92 *
93 * IWORK (local workspace/output) INTEGER array,
94 * dimension 7*N + 8*NPCOL + 2
95 *
96 * INFO (global output) INTEGER
97 * = 0: successful exit
98 * < 0: If the i-th argument is an array and the j-entry had
99 * an illegal value, then INFO = -(i*100+j), if the i-th
100 * argument is a scalar and had an illegal value, then
101 * INFO = -i.
102 * > 0: The algorithm failed to compute the ith eigenvalue.
103 *
104 * =====================================================================
105 *
106 * .. Parameters ..
107 *
108  INTEGER BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_,
109  $ mb_, nb_, rsrc_, csrc_, lld_
110  parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
111  $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
112  $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
113  DOUBLE PRECISION ZERO, ONE
114  parameter( zero = 0.0d+0, one = 1.0d+0 )
115 * ..
116 * .. Local Scalars ..
117  INTEGER COL, COLTYP, IBUF, ICTOT, ICTXT, IDLMDA, IIQ,
118  $ indcol, indrow, indx, indxc, indxp, indxr, inq,
119  $ ipq, ipq2, ipsm, ipu, ipwork, iq1, iq2, iqcol,
120  $ iqq, iqrow, iw, iz, j, jc, jj2c, jjc, jjq, jnq,
121  $ k, ldq, ldq2, ldu, mycol, myrow, nb, nn, nn1,
122  $ nn2, np, npcol, nprow, nq
123 * ..
124 * .. Local Arrays ..
125  INTEGER DESCQ2( DLEN_ ), DESCU( DLEN_ )
126 * ..
127 * .. External Functions ..
128  INTEGER NUMROC
129  EXTERNAL numroc
130 * ..
131 * .. External Subroutines ..
132  EXTERNAL blacs_gridinfo, dcopy, descinit, infog1l,
133  $ infog2l, pdgemm, pdlaed2, pdlaed3, pdlaedz,
134  $ pdlaset, pxerbla
135 * ..
136 * .. Intrinsic Functions ..
137  INTRINSIC max, min
138 * ..
139 * .. Executable Statements ..
140 *
141 * This is just to keep ftnchek and toolpack/1 happy
142  IF( block_cyclic_2d*csrc_*ctxt_*dlen_*dtype_*lld_*mb_*m_*nb_*n_*
143  $ rsrc_.LT.0 )RETURN
144 *
145 *
146 * Test the input parameters.
147 *
148  CALL blacs_gridinfo( descq( ctxt_ ), nprow, npcol, myrow, mycol )
149  info = 0
150  IF( nprow.EQ.-1 ) THEN
151  info = -( 600+ctxt_ )
152  ELSE IF( n.LT.0 ) THEN
153  info = -1
154  ELSE IF( id.GT.descq( n_ ) ) THEN
155  info = -4
156  ELSE IF( n1.GE.n ) THEN
157  info = -2
158  END IF
159  IF( info.NE.0 ) THEN
160  CALL pxerbla( descq( ctxt_ ), 'PDLAED1', -info )
161  RETURN
162  END IF
163 *
164 * Quick return if possible
165 *
166  IF( n.EQ.0 )
167  $ RETURN
168 *
169 * The following values are integer pointers which indicate
170 * the portion of the workspace used by a particular array
171 * in PDLAED2 and PDLAED3.
172 *
173  ictxt = descq( ctxt_ )
174  nb = descq( nb_ )
175  ldq = descq( lld_ )
176 *
177  CALL infog2l( iq-1+id, jq-1+id, descq, nprow, npcol, myrow, mycol,
178  $ iiq, jjq, iqrow, iqcol )
179 *
180  np = numroc( n, descq( mb_ ), myrow, iqrow, nprow )
181  nq = numroc( n, descq( nb_ ), mycol, iqcol, npcol )
182 *
183  ldq2 = max( np, 1 )
184  ldu = ldq2
185 *
186  iz = 1
187  idlmda = iz + n
188  iw = idlmda + n
189  ipq2 = iw + n
190  ipu = ipq2 + ldq2*nq
191  ibuf = ipu + ldu*nq
192 * (IBUF est de taille 3*N au maximum)
193 *
194  ictot = 1
195  ipsm = ictot + npcol*4
196  indx = ipsm + npcol*4
197  indxc = indx + n
198  indxp = indxc + n
199  indcol = indxp + n
200  coltyp = indcol + n
201  indrow = coltyp + n
202  indxr = indrow + n
203 *
204  CALL descinit( descq2, n, n, nb, nb, iqrow, iqcol, ictxt, ldq2,
205  $ info )
206  CALL descinit( descu, n, n, nb, nb, iqrow, iqcol, ictxt, ldu,
207  $ info )
208 *
209 * Form the z-vector which consists of the last row of Q_1 and the
210 * first row of Q_2.
211 *
212  ipwork = idlmda
213  CALL pdlaedz( n, n1, id, q, iq, jq, ldq, descq, work( iz ),
214  $ work( ipwork ) )
215 *
216 * Deflate eigenvalues.
217 *
218  ipq = iiq + ( jjq-1 )*ldq
219  CALL pdlaed2( ictxt, k, n, n1, nb, d, iqrow, iqcol, q( ipq ), ldq,
220  $ rho, work( iz ), work( iw ), work( idlmda ),
221  $ work( ipq2 ), ldq2, work( ibuf ), iwork( ictot ),
222  $ iwork( ipsm ), npcol, iwork( indx ), iwork( indxc ),
223  $ iwork( indxp ), iwork( indcol ), iwork( coltyp ),
224  $ nn, nn1, nn2, iq1, iq2 )
225 *
226 *
227 * Solve Secular Equation.
228 *
229  IF( k.NE.0 ) THEN
230  CALL pdlaset( 'A', n, n, zero, one, work( ipu ), 1, 1, descu )
231  CALL pdlaed3( ictxt, k, n, nb, d, iqrow, iqcol, rho,
232  $ work( idlmda ), work( iw ), work( iz ),
233  $ work( ipu ), ldq2, work( ibuf ), iwork( indx ),
234  $ iwork( indcol ), iwork( indrow ), iwork( indxr ),
235  $ iwork( indxc ), iwork( ictot ), npcol, info )
236 *
237 * Compute the updated eigenvectors.
238 *
239  iqq = min( iq1, iq2 )
240  IF( nn1.GT.0 ) THEN
241  inq = iq - 1 + id
242  jnq = jq - 1 + id + iqq - 1
243  CALL pdgemm( 'N', 'N', n1, nn, nn1, one, work( ipq2 ), 1,
244  $ iq1, descq2, work( ipu ), iq1, iqq, descu,
245  $ zero, q, inq, jnq, descq )
246  END IF
247  IF( nn2.GT.0 ) THEN
248  inq = iq - 1 + id + n1
249  jnq = jq - 1 + id + iqq - 1
250  CALL pdgemm( 'N', 'N', n-n1, nn, nn2, one, work( ipq2 ),
251  $ n1+1, iq2, descq2, work( ipu ), iq2, iqq,
252  $ descu, zero, q, inq, jnq, descq )
253  END IF
254 *
255  DO 10 j = k + 1, n
256  jc = iwork( indx+j-1 )
257  CALL infog1l( jq-1+jc, nb, npcol, mycol, iqcol, jjc, col )
258  CALL infog1l( jc, nb, npcol, mycol, iqcol, jj2c, col )
259  IF( mycol.EQ.col ) THEN
260  iq2 = ipq2 + ( jj2c-1 )*ldq2
261  inq = ipq + ( jjc-1 )*ldq
262  CALL dcopy( np, work( iq2 ), 1, q( inq ), 1 )
263  END IF
264  10 CONTINUE
265  END IF
266 *
267  20 CONTINUE
268  RETURN
269 *
270 * End of PDLAED1
271 *
272  END
max
#define max(A, B)
Definition: pcgemr.c:180
pdlaed1
subroutine pdlaed1(N, N1, D, ID, Q, IQ, JQ, DESCQ, RHO, WORK, IWORK, INFO)
Definition: pdlaed1.f:3
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
pdlaed2
subroutine pdlaed2(ICTXT, K, N, N1, NB, D, DROW, DCOL, Q, LDQ, RHO, Z, W, DLAMDA, Q2, LDQ2, QBUF, CTOT, PSM, NPCOL, INDX, INDXC, INDXP, INDCOL, COLTYP, NN, NN1, NN2, IB1, IB2)
Definition: pdlaed2.f:5
descinit
subroutine descinit(DESC, M, N, MB, NB, IRSRC, ICSRC, ICTXT, LLD, INFO)
Definition: descinit.f:3
pdlaset
subroutine pdlaset(UPLO, M, N, ALPHA, BETA, A, IA, JA, DESCA)
Definition: pdblastst.f:6862
pxerbla
subroutine pxerbla(ICTXT, SRNAME, INFO)
Definition: pxerbla.f:2
pdlaedz
subroutine pdlaedz(N, N1, ID, Q, IQ, JQ, LDQ, DESCQ, Z, WORK)
Definition: pdlaedz.f:2
pdlaed3
subroutine pdlaed3(ICTXT, K, N, NB, D, DROW, DCOL, RHO, DLAMDA, W, Z, U, LDU, BUF, INDX, INDCOL, INDROW, INDXR, INDXC, CTOT, NPCOL, INFO)
Definition: pdlaed3.f:4
min
#define min(A, B)
Definition: pcgemr.c:181