001:       SUBROUTINE SLAED0( ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS,
002:      $                   WORK, IWORK, INFO )
003: *
004: *  -- LAPACK routine (version 3.2) --
005: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
006: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       INTEGER            ICOMPQ, INFO, LDQ, LDQS, N, QSIZ
011: *     ..
012: *     .. Array Arguments ..
013:       INTEGER            IWORK( * )
014:       REAL               D( * ), E( * ), Q( LDQ, * ), QSTORE( LDQS, * ),
015:      $                   WORK( * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  SLAED0 computes all eigenvalues and corresponding eigenvectors of a
022: *  symmetric tridiagonal matrix using the divide and conquer method.
023: *
024: *  Arguments
025: *  =========
026: *
027: *  ICOMPQ  (input) INTEGER
028: *          = 0:  Compute eigenvalues only.
029: *          = 1:  Compute eigenvectors of original dense symmetric matrix
030: *                also.  On entry, Q contains the orthogonal matrix used
031: *                to reduce the original matrix to tridiagonal form.
032: *          = 2:  Compute eigenvalues and eigenvectors of tridiagonal
033: *                matrix.
034: *
035: *  QSIZ   (input) INTEGER
036: *         The dimension of the orthogonal matrix used to reduce
037: *         the full matrix to tridiagonal form.  QSIZ >= N if ICOMPQ = 1.
038: *
039: *  N      (input) INTEGER
040: *         The dimension of the symmetric tridiagonal matrix.  N >= 0.
041: *
042: *  D      (input/output) REAL array, dimension (N)
043: *         On entry, the main diagonal of the tridiagonal matrix.
044: *         On exit, its eigenvalues.
045: *
046: *  E      (input) REAL array, dimension (N-1)
047: *         The off-diagonal elements of the tridiagonal matrix.
048: *         On exit, E has been destroyed.
049: *
050: *  Q      (input/output) REAL array, dimension (LDQ, N)
051: *         On entry, Q must contain an N-by-N orthogonal matrix.
052: *         If ICOMPQ = 0    Q is not referenced.
053: *         If ICOMPQ = 1    On entry, Q is a subset of the columns of the
054: *                          orthogonal matrix used to reduce the full
055: *                          matrix to tridiagonal form corresponding to
056: *                          the subset of the full matrix which is being
057: *                          decomposed at this time.
058: *         If ICOMPQ = 2    On entry, Q will be the identity matrix.
059: *                          On exit, Q contains the eigenvectors of the
060: *                          tridiagonal matrix.
061: *
062: *  LDQ    (input) INTEGER
063: *         The leading dimension of the array Q.  If eigenvectors are
064: *         desired, then  LDQ >= max(1,N).  In any case,  LDQ >= 1.
065: *
066: *  QSTORE (workspace) REAL array, dimension (LDQS, N)
067: *         Referenced only when ICOMPQ = 1.  Used to store parts of
068: *         the eigenvector matrix when the updating matrix multiplies
069: *         take place.
070: *
071: *  LDQS   (input) INTEGER
072: *         The leading dimension of the array QSTORE.  If ICOMPQ = 1,
073: *         then  LDQS >= max(1,N).  In any case,  LDQS >= 1.
074: *
075: *  WORK   (workspace) REAL array,
076: *         If ICOMPQ = 0 or 1, the dimension of WORK must be at least
077: *                     1 + 3*N + 2*N*lg N + 2*N**2
078: *                     ( lg( N ) = smallest integer k
079: *                                 such that 2^k >= N )
080: *         If ICOMPQ = 2, the dimension of WORK must be at least
081: *                     4*N + N**2.
082: *
083: *  IWORK  (workspace) INTEGER array,
084: *         If ICOMPQ = 0 or 1, the dimension of IWORK must be at least
085: *                        6 + 6*N + 5*N*lg N.
086: *                        ( lg( N ) = smallest integer k
087: *                                    such that 2^k >= N )
088: *         If ICOMPQ = 2, the dimension of IWORK must be at least
089: *                        3 + 5*N.
090: *
091: *  INFO   (output) INTEGER
092: *          = 0:  successful exit.
093: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
094: *          > 0:  The algorithm failed to compute an eigenvalue while
095: *                working on the submatrix lying in rows and columns
096: *                INFO/(N+1) through mod(INFO,N+1).
097: *
098: *  Further Details
099: *  ===============
100: *
101: *  Based on contributions by
102: *     Jeff Rutter, Computer Science Division, University of California
103: *     at Berkeley, USA
104: *
105: *  =====================================================================
106: *
107: *     .. Parameters ..
108:       REAL               ZERO, ONE, TWO
109:       PARAMETER          ( ZERO = 0.E0, ONE = 1.E0, TWO = 2.E0 )
110: *     ..
111: *     .. Local Scalars ..
112:       INTEGER            CURLVL, CURPRB, CURR, I, IGIVCL, IGIVNM,
113:      $                   IGIVPT, INDXQ, IPERM, IPRMPT, IQ, IQPTR, IWREM,
114:      $                   J, K, LGN, MATSIZ, MSD2, SMLSIZ, SMM1, SPM1,
115:      $                   SPM2, SUBMAT, SUBPBS, TLVLS
116:       REAL               TEMP
117: *     ..
118: *     .. External Subroutines ..
119:       EXTERNAL           SCOPY, SGEMM, SLACPY, SLAED1, SLAED7, SSTEQR,
120:      $                   XERBLA
121: *     ..
122: *     .. External Functions ..
123:       INTEGER            ILAENV
124:       EXTERNAL           ILAENV
125: *     ..
126: *     .. Intrinsic Functions ..
127:       INTRINSIC          ABS, INT, LOG, MAX, REAL
128: *     ..
129: *     .. Executable Statements ..
130: *
131: *     Test the input parameters.
132: *
133:       INFO = 0
134: *
135:       IF( ICOMPQ.LT.0 .OR. ICOMPQ.GT.2 ) THEN
136:          INFO = -1
137:       ELSE IF( ( ICOMPQ.EQ.1 ) .AND. ( QSIZ.LT.MAX( 0, N ) ) ) THEN
138:          INFO = -2
139:       ELSE IF( N.LT.0 ) THEN
140:          INFO = -3
141:       ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
142:          INFO = -7
143:       ELSE IF( LDQS.LT.MAX( 1, N ) ) THEN
144:          INFO = -9
145:       END IF
146:       IF( INFO.NE.0 ) THEN
147:          CALL XERBLA( 'SLAED0', -INFO )
148:          RETURN
149:       END IF
150: *
151: *     Quick return if possible
152: *
153:       IF( N.EQ.0 )
154:      $   RETURN
155: *
156:       SMLSIZ = ILAENV( 9, 'SLAED0', ' ', 0, 0, 0, 0 )
157: *
158: *     Determine the size and placement of the submatrices, and save in
159: *     the leading elements of IWORK.
160: *
161:       IWORK( 1 ) = N
162:       SUBPBS = 1
163:       TLVLS = 0
164:    10 CONTINUE
165:       IF( IWORK( SUBPBS ).GT.SMLSIZ ) THEN
166:          DO 20 J = SUBPBS, 1, -1
167:             IWORK( 2*J ) = ( IWORK( J )+1 ) / 2
168:             IWORK( 2*J-1 ) = IWORK( J ) / 2
169:    20    CONTINUE
170:          TLVLS = TLVLS + 1
171:          SUBPBS = 2*SUBPBS
172:          GO TO 10
173:       END IF
174:       DO 30 J = 2, SUBPBS
175:          IWORK( J ) = IWORK( J ) + IWORK( J-1 )
176:    30 CONTINUE
177: *
178: *     Divide the matrix into SUBPBS submatrices of size at most SMLSIZ+1
179: *     using rank-1 modifications (cuts).
180: *
181:       SPM1 = SUBPBS - 1
182:       DO 40 I = 1, SPM1
183:          SUBMAT = IWORK( I ) + 1
184:          SMM1 = SUBMAT - 1
185:          D( SMM1 ) = D( SMM1 ) - ABS( E( SMM1 ) )
186:          D( SUBMAT ) = D( SUBMAT ) - ABS( E( SMM1 ) )
187:    40 CONTINUE
188: *
189:       INDXQ = 4*N + 3
190:       IF( ICOMPQ.NE.2 ) THEN
191: *
192: *        Set up workspaces for eigenvalues only/accumulate new vectors
193: *        routine
194: *
195:          TEMP = LOG( REAL( N ) ) / LOG( TWO )
196:          LGN = INT( TEMP )
197:          IF( 2**LGN.LT.N )
198:      $      LGN = LGN + 1
199:          IF( 2**LGN.LT.N )
200:      $      LGN = LGN + 1
201:          IPRMPT = INDXQ + N + 1
202:          IPERM = IPRMPT + N*LGN
203:          IQPTR = IPERM + N*LGN
204:          IGIVPT = IQPTR + N + 2
205:          IGIVCL = IGIVPT + N*LGN
206: *
207:          IGIVNM = 1
208:          IQ = IGIVNM + 2*N*LGN
209:          IWREM = IQ + N**2 + 1
210: *
211: *        Initialize pointers
212: *
213:          DO 50 I = 0, SUBPBS
214:             IWORK( IPRMPT+I ) = 1
215:             IWORK( IGIVPT+I ) = 1
216:    50    CONTINUE
217:          IWORK( IQPTR ) = 1
218:       END IF
219: *
220: *     Solve each submatrix eigenproblem at the bottom of the divide and
221: *     conquer tree.
222: *
223:       CURR = 0
224:       DO 70 I = 0, SPM1
225:          IF( I.EQ.0 ) THEN
226:             SUBMAT = 1
227:             MATSIZ = IWORK( 1 )
228:          ELSE
229:             SUBMAT = IWORK( I ) + 1
230:             MATSIZ = IWORK( I+1 ) - IWORK( I )
231:          END IF
232:          IF( ICOMPQ.EQ.2 ) THEN
233:             CALL SSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ),
234:      $                   Q( SUBMAT, SUBMAT ), LDQ, WORK, INFO )
235:             IF( INFO.NE.0 )
236:      $         GO TO 130
237:          ELSE
238:             CALL SSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ),
239:      $                   WORK( IQ-1+IWORK( IQPTR+CURR ) ), MATSIZ, WORK,
240:      $                   INFO )
241:             IF( INFO.NE.0 )
242:      $         GO TO 130
243:             IF( ICOMPQ.EQ.1 ) THEN
244:                CALL SGEMM( 'N', 'N', QSIZ, MATSIZ, MATSIZ, ONE,
245:      $                     Q( 1, SUBMAT ), LDQ, WORK( IQ-1+IWORK( IQPTR+
246:      $                     CURR ) ), MATSIZ, ZERO, QSTORE( 1, SUBMAT ),
247:      $                     LDQS )
248:             END IF
249:             IWORK( IQPTR+CURR+1 ) = IWORK( IQPTR+CURR ) + MATSIZ**2
250:             CURR = CURR + 1
251:          END IF
252:          K = 1
253:          DO 60 J = SUBMAT, IWORK( I+1 )
254:             IWORK( INDXQ+J ) = K
255:             K = K + 1
256:    60    CONTINUE
257:    70 CONTINUE
258: *
259: *     Successively merge eigensystems of adjacent submatrices
260: *     into eigensystem for the corresponding larger matrix.
261: *
262: *     while ( SUBPBS > 1 )
263: *
264:       CURLVL = 1
265:    80 CONTINUE
266:       IF( SUBPBS.GT.1 ) THEN
267:          SPM2 = SUBPBS - 2
268:          DO 90 I = 0, SPM2, 2
269:             IF( I.EQ.0 ) THEN
270:                SUBMAT = 1
271:                MATSIZ = IWORK( 2 )
272:                MSD2 = IWORK( 1 )
273:                CURPRB = 0
274:             ELSE
275:                SUBMAT = IWORK( I ) + 1
276:                MATSIZ = IWORK( I+2 ) - IWORK( I )
277:                MSD2 = MATSIZ / 2
278:                CURPRB = CURPRB + 1
279:             END IF
280: *
281: *     Merge lower order eigensystems (of size MSD2 and MATSIZ - MSD2)
282: *     into an eigensystem of size MATSIZ.
283: *     SLAED1 is used only for the full eigensystem of a tridiagonal
284: *     matrix.
285: *     SLAED7 handles the cases in which eigenvalues only or eigenvalues
286: *     and eigenvectors of a full symmetric matrix (which was reduced to
287: *     tridiagonal form) are desired.
288: *
289:             IF( ICOMPQ.EQ.2 ) THEN
290:                CALL SLAED1( MATSIZ, D( SUBMAT ), Q( SUBMAT, SUBMAT ),
291:      $                      LDQ, IWORK( INDXQ+SUBMAT ),
292:      $                      E( SUBMAT+MSD2-1 ), MSD2, WORK,
293:      $                      IWORK( SUBPBS+1 ), INFO )
294:             ELSE
295:                CALL SLAED7( ICOMPQ, MATSIZ, QSIZ, TLVLS, CURLVL, CURPRB,
296:      $                      D( SUBMAT ), QSTORE( 1, SUBMAT ), LDQS,
297:      $                      IWORK( INDXQ+SUBMAT ), E( SUBMAT+MSD2-1 ),
298:      $                      MSD2, WORK( IQ ), IWORK( IQPTR ),
299:      $                      IWORK( IPRMPT ), IWORK( IPERM ),
300:      $                      IWORK( IGIVPT ), IWORK( IGIVCL ),
301:      $                      WORK( IGIVNM ), WORK( IWREM ),
302:      $                      IWORK( SUBPBS+1 ), INFO )
303:             END IF
304:             IF( INFO.NE.0 )
305:      $         GO TO 130
306:             IWORK( I / 2+1 ) = IWORK( I+2 )
307:    90    CONTINUE
308:          SUBPBS = SUBPBS / 2
309:          CURLVL = CURLVL + 1
310:          GO TO 80
311:       END IF
312: *
313: *     end while
314: *
315: *     Re-merge the eigenvalues/vectors which were deflated at the final
316: *     merge step.
317: *
318:       IF( ICOMPQ.EQ.1 ) THEN
319:          DO 100 I = 1, N
320:             J = IWORK( INDXQ+I )
321:             WORK( I ) = D( J )
322:             CALL SCOPY( QSIZ, QSTORE( 1, J ), 1, Q( 1, I ), 1 )
323:   100    CONTINUE
324:          CALL SCOPY( N, WORK, 1, D, 1 )
325:       ELSE IF( ICOMPQ.EQ.2 ) THEN
326:          DO 110 I = 1, N
327:             J = IWORK( INDXQ+I )
328:             WORK( I ) = D( J )
329:             CALL SCOPY( N, Q( 1, J ), 1, WORK( N*I+1 ), 1 )
330:   110    CONTINUE
331:          CALL SCOPY( N, WORK, 1, D, 1 )
332:          CALL SLACPY( 'A', N, N, WORK( N+1 ), N, Q, LDQ )
333:       ELSE
334:          DO 120 I = 1, N
335:             J = IWORK( INDXQ+I )
336:             WORK( I ) = D( J )
337:   120    CONTINUE
338:          CALL SCOPY( N, WORK, 1, D, 1 )
339:       END IF
340:       GO TO 140
341: *
342:   130 CONTINUE
343:       INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1
344: *
345:   140 CONTINUE
346:       RETURN
347: *
348: *     End of SLAED0
349: *
350:       END
351: