SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
pslaed0.f
Go to the documentation of this file.
1 SUBROUTINE pslaed0( N, D, E, Q, IQ, JQ, DESCQ, WORK, IWORK, INFO )
2*
3* -- ScaLAPACK auxiliary routine (version 1.7) --
4* University of Tennessee, Knoxville, Oak Ridge National Laboratory,
5* and University of California, Berkeley.
6* December 31, 1998
7*
8* .. Scalar Arguments ..
9 INTEGER INFO, IQ, JQ, N
10* ..
11* .. Array Arguments ..
12 INTEGER DESCQ( * ), IWORK( * )
13 REAL D( * ), E( * ), Q( * ), WORK( * )
14* ..
15*
16* Purpose
17* =======
18*
19* PSLAED0 computes all eigenvalues and corresponding eigenvectors of a
20* symmetric tridiagonal matrix using the divide and conquer method.
21*
22*
23* Arguments
24* =========
25*
26* N (global input) INTEGER
27* The order of the tridiagonal matrix T. N >= 0.
28*
29* D (global input/output) REAL array, dimension (N)
30* On entry, the diagonal elements of the tridiagonal matrix.
31* On exit, if INFO = 0, the eigenvalues in descending order.
32*
33* E (global input/output) REAL array, dimension (N-1)
34* On entry, the subdiagonal elements of the tridiagonal matrix.
35* On exit, E has been destroyed.
36*
37* Q (local output) REAL array,
38* global dimension (N, N),
39* local dimension ( LLD_Q, LOCc(JQ+N-1))
40* Q contains the orthonormal eigenvectors of the symmetric
41* tridiagonal matrix.
42* On output, Q is distributed across the P processes in block
43* cyclic format.
44*
45* IQ (global input) INTEGER
46* Q's global row index, which points to the beginning of the
47* submatrix which is to be operated on.
48*
49* JQ (global input) INTEGER
50* Q's global column index, which points to the beginning of
51* the submatrix which is to be operated on.
52*
53* DESCQ (global and local input) INTEGER array of dimension DLEN_.
54* The array descriptor for the distributed matrix Z.
55*
56*
57* WORK (local workspace ) REAL array, dimension (LWORK)
58* LWORK = 6*N + 2*NP*NQ, with
59* NP = NUMROC( N, MB_Q, MYROW, IQROW, NPROW )
60* NQ = NUMROC( N, NB_Q, MYCOL, IQCOL, NPCOL )
61* IQROW = INDXG2P( IQ, NB_Q, MYROW, RSRC_Q, NPROW )
62* IQCOL = INDXG2P( JQ, MB_Q, MYCOL, CSRC_Q, NPCOL )
63*
64* IWORK (local workspace/output) INTEGER array, dimension (LIWORK)
65* LIWORK = 2 + 7*N + 8*NPCOL
66*
67* INFO (global output) INTEGER
68* = 0: successful exit
69* < 0: If the i-th argument is an array and the j-entry had
70* an illegal value, then INFO = -(i*100+j), if the i-th
71* argument is a scalar and had an illegal value, then
72* INFO = -i.
73* > 0: The algorithm failed to compute the INFO/(N+1) th
74* eigenvalue while working on the submatrix lying in
75* global rows and columns mod(INFO,N+1).
76*
77* =====================================================================
78*
79* .. Parameters ..
80*
81 INTEGER BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_,
82 $ MB_, NB_, RSRC_, CSRC_, LLD_
83 parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
84 $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
85 $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
86* ..
87* .. Local Scalars ..
88 INTEGER I, ID, IDCOL, IDROW, IID, IINFO, IIQ, IM1, IM2,
89 $ IPQ, IQCOL, IQROW, J, JJD, JJQ, LDQ, MATSIZ,
90 $ MYCOL, MYROW, N1, NB, NBL, NBL1, NPCOL, NPROW,
91 $ SUBPBS, TSUBPBS
92* ..
93* .. External Subroutines ..
94 EXTERNAL blacs_gridinfo, infog2l, pslaed1, pxerbla,
95 $ sgebr2d, sgebs2d, sgerv2d, sgesd2d, ssteqr
96* ..
97* .. Intrinsic Functions ..
98 INTRINSIC abs, min
99* ..
100* .. Executable Statements ..
101*
102* This is just to keep ftnchek and toolpack/1 happy
103 IF( block_cyclic_2d*csrc_*ctxt_*dlen_*dtype_*lld_*mb_*m_*nb_*n_*
104 $ rsrc_.LT.0 )RETURN
105*
106* Test the input parameters.
107*
108 CALL blacs_gridinfo( descq( ctxt_ ), nprow, npcol, myrow, mycol )
109 info = 0
110 IF( descq( nb_ ).GT.n .OR. n.LT.2 )
111 $ info = -1
112 IF( info.NE.0 ) THEN
113 CALL pxerbla( descq( ctxt_ ), 'PSLAED0', -info )
114 RETURN
115 END IF
116*
117 nb = descq( nb_ )
118 ldq = descq( lld_ )
119 CALL infog2l( iq, jq, descq, nprow, npcol, myrow, mycol, iiq, jjq,
120 $ iqrow, iqcol )
121*
122* Determine the size and placement of the submatrices, and save in
123* the leading elements of IWORK.
124*
125 tsubpbs = ( n-1 ) / nb + 1
126 iwork( 1 ) = tsubpbs
127 subpbs = 1
128 10 CONTINUE
129 IF( iwork( subpbs ).GT.1 ) THEN
130 DO 20 j = subpbs, 1, -1
131 iwork( 2*j ) = ( iwork( j )+1 ) / 2
132 iwork( 2*j-1 ) = iwork( j ) / 2
133 20 CONTINUE
134 subpbs = 2*subpbs
135 GO TO 10
136 END IF
137 DO 30 j = 2, subpbs
138 iwork( j ) = iwork( j ) + iwork( j-1 )
139 30 CONTINUE
140*
141* Divide the matrix into TSUBPBS submatrices of size at most NB
142* using rank-1 modifications (cuts).
143*
144 DO 40 i = nb + 1, n, nb
145 im1 = i - 1
146 d( im1 ) = d( im1 ) - abs( e( im1 ) )
147 d( i ) = d( i ) - abs( e( im1 ) )
148 40 CONTINUE
149*
150* Solve each submatrix eigenproblem at the bottom of the divide and
151* conquer tree. D is the same on each process.
152*
153 DO 50 id = 1, n, nb
154 CALL infog2l( iq-1+id, jq-1+id, descq, nprow, npcol, myrow,
155 $ mycol, iid, jjd, idrow, idcol )
156 matsiz = min( nb, n-id+1 )
157 IF( myrow.EQ.idrow .AND. mycol.EQ.idcol ) THEN
158 ipq = iid + ( jjd-1 )*ldq
159 CALL ssteqr( 'I', matsiz, d( id ), e( id ), q( ipq ), ldq,
160 $ work, info )
161 IF( info.NE.0 ) THEN
162 CALL pxerbla( descq( ctxt_ ), 'SSTEQR', -info )
163 RETURN
164 END IF
165 IF( myrow.NE.iqrow .OR. mycol.NE.iqcol ) THEN
166 CALL sgesd2d( descq( ctxt_ ), matsiz, 1, d( id ), matsiz,
167 $ iqrow, iqcol )
168 END IF
169 ELSE IF( myrow.EQ.iqrow .AND. mycol.EQ.iqcol ) THEN
170 CALL sgerv2d( descq( ctxt_ ), matsiz, 1, d( id ), matsiz,
171 $ idrow, idcol )
172 END IF
173 50 CONTINUE
174*
175 IF( myrow.EQ.iqrow .AND. mycol.EQ.iqcol ) THEN
176 CALL sgebs2d( descq( ctxt_ ), 'A', ' ', n, 1, d, n )
177 ELSE
178 CALL sgebr2d( descq( ctxt_ ), 'A', ' ', n, 1, d, n, iqrow,
179 $ iqcol )
180 END IF
181*
182* Successively merge eigensystems of adjacent submatrices
183* into eigensystem for the corresponding larger matrix.
184*
185* while ( SUBPBS > 1 )
186*
187 60 CONTINUE
188 IF( subpbs.GT.1 ) THEN
189 im2 = subpbs - 2
190 DO 80 i = 0, im2, 2
191 IF( i.EQ.0 ) THEN
192 nbl = iwork( 2 )
193 nbl1 = iwork( 1 )
194 IF( nbl1.EQ.0 )
195 $ GO TO 70
196 id = 1
197 matsiz = min( n, nbl*nb )
198 n1 = nbl1*nb
199 ELSE
200 nbl = iwork( i+2 ) - iwork( i )
201 nbl1 = nbl / 2
202 IF( nbl1.EQ.0 )
203 $ GO TO 70
204 id = iwork( i )*nb + 1
205 matsiz = min( nb*nbl, n-id+1 )
206 n1 = nbl1*nb
207 END IF
208*
209* Merge lower order eigensystems (of size N1 and MATSIZ - N1)
210* into an eigensystem of size MATSIZ.
211*
212 CALL pslaed1( matsiz, n1, d( id ), id, q, iq, jq, descq,
213 $ e( id+n1-1 ), work, iwork( subpbs+1 ), iinfo )
214 IF( iinfo.NE.0 ) THEN
215 info = iinfo*( n+1 ) + id
216 END IF
217 70 CONTINUE
218 iwork( i / 2+1 ) = iwork( i+2 )
219 80 CONTINUE
220 subpbs = subpbs / 2
221 GO TO 60
222 END IF
223*
224* end while
225*
226 90 CONTINUE
227 RETURN
228*
229* End of PSLAED0
230*
231 END
subroutine infog2l(grindx, gcindx, desc, nprow, npcol, myrow, mycol, lrindx, lcindx, rsrc, csrc)
Definition infog2l.f:3
#define min(A, B)
Definition pcgemr.c:181
subroutine pslaed0(n, d, e, q, iq, jq, descq, work, iwork, info)
Definition pslaed0.f:2
subroutine pslaed1(n, n1, d, id, q, iq, jq, descq, rho, work, iwork, info)
Definition pslaed1.f:3
subroutine pxerbla(ictxt, srname, info)
Definition pxerbla.f:2