001:       SUBROUTINE DLASD0( N, SQRE, D, E, U, LDU, VT, LDVT, SMLSIZ, IWORK,
002:      $                   WORK, INFO )
003: *
004: *  -- LAPACK auxiliary routine (version 3.2) --
005: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       INTEGER            INFO, LDU, LDVT, N, SMLSIZ, SQRE
010: *     ..
011: *     .. Array Arguments ..
012:       INTEGER            IWORK( * )
013:       DOUBLE PRECISION   D( * ), E( * ), U( LDU, * ), VT( LDVT, * ),
014:      $                   WORK( * )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  Using a divide and conquer approach, DLASD0 computes the singular
021: *  value decomposition (SVD) of a real upper bidiagonal N-by-M
022: *  matrix B with diagonal D and offdiagonal E, where M = N + SQRE.
023: *  The algorithm computes orthogonal matrices U and VT such that
024: *  B = U * S * VT. The singular values S are overwritten on D.
025: *
026: *  A related subroutine, DLASDA, computes only the singular values,
027: *  and optionally, the singular vectors in compact form.
028: *
029: *  Arguments
030: *  =========
031: *
032: *  N      (input) INTEGER
033: *         On entry, the row dimension of the upper bidiagonal matrix.
034: *         This is also the dimension of the main diagonal array D.
035: *
036: *  SQRE   (input) INTEGER
037: *         Specifies the column dimension of the bidiagonal matrix.
038: *         = 0: The bidiagonal matrix has column dimension M = N;
039: *         = 1: The bidiagonal matrix has column dimension M = N+1;
040: *
041: *  D      (input/output) DOUBLE PRECISION array, dimension (N)
042: *         On entry D contains the main diagonal of the bidiagonal
043: *         matrix.
044: *         On exit D, if INFO = 0, contains its singular values.
045: *
046: *  E      (input) DOUBLE PRECISION array, dimension (M-1)
047: *         Contains the subdiagonal entries of the bidiagonal matrix.
048: *         On exit, E has been destroyed.
049: *
050: *  U      (output) DOUBLE PRECISION array, dimension at least (LDQ, N)
051: *         On exit, U contains the left singular vectors.
052: *
053: *  LDU    (input) INTEGER
054: *         On entry, leading dimension of U.
055: *
056: *  VT     (output) DOUBLE PRECISION array, dimension at least (LDVT, M)
057: *         On exit, VT' contains the right singular vectors.
058: *
059: *  LDVT   (input) INTEGER
060: *         On entry, leading dimension of VT.
061: *
062: *  SMLSIZ (input) INTEGER
063: *         On entry, maximum size of the subproblems at the
064: *         bottom of the computation tree.
065: *
066: *  IWORK  (workspace) INTEGER work array.
067: *         Dimension must be at least (8 * N)
068: *
069: *  WORK   (workspace) DOUBLE PRECISION work array.
070: *         Dimension must be at least (3 * M**2 + 2 * M)
071: *
072: *  INFO   (output) INTEGER
073: *          = 0:  successful exit.
074: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
075: *          > 0:  if INFO = 1, an singular value did not converge
076: *
077: *  Further Details
078: *  ===============
079: *
080: *  Based on contributions by
081: *     Ming Gu and Huan Ren, Computer Science Division, University of
082: *     California at Berkeley, USA
083: *
084: *  =====================================================================
085: *
086: *     .. Local Scalars ..
087:       INTEGER            I, I1, IC, IDXQ, IDXQC, IM1, INODE, ITEMP, IWK,
088:      $                   J, LF, LL, LVL, M, NCC, ND, NDB1, NDIML, NDIMR,
089:      $                   NL, NLF, NLP1, NLVL, NR, NRF, NRP1, SQREI
090:       DOUBLE PRECISION   ALPHA, BETA
091: *     ..
092: *     .. External Subroutines ..
093:       EXTERNAL           DLASD1, DLASDQ, DLASDT, XERBLA
094: *     ..
095: *     .. Executable Statements ..
096: *
097: *     Test the input parameters.
098: *
099:       INFO = 0
100: *
101:       IF( N.LT.0 ) THEN
102:          INFO = -1
103:       ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
104:          INFO = -2
105:       END IF
106: *
107:       M = N + SQRE
108: *
109:       IF( LDU.LT.N ) THEN
110:          INFO = -6
111:       ELSE IF( LDVT.LT.M ) THEN
112:          INFO = -8
113:       ELSE IF( SMLSIZ.LT.3 ) THEN
114:          INFO = -9
115:       END IF
116:       IF( INFO.NE.0 ) THEN
117:          CALL XERBLA( 'DLASD0', -INFO )
118:          RETURN
119:       END IF
120: *
121: *     If the input matrix is too small, call DLASDQ to find the SVD.
122: *
123:       IF( N.LE.SMLSIZ ) THEN
124:          CALL DLASDQ( 'U', SQRE, N, M, N, 0, D, E, VT, LDVT, U, LDU, U,
125:      $                LDU, WORK, INFO )
126:          RETURN
127:       END IF
128: *
129: *     Set up the computation tree.
130: *
131:       INODE = 1
132:       NDIML = INODE + N
133:       NDIMR = NDIML + N
134:       IDXQ = NDIMR + N
135:       IWK = IDXQ + N
136:       CALL DLASDT( N, NLVL, ND, IWORK( INODE ), IWORK( NDIML ),
137:      $             IWORK( NDIMR ), SMLSIZ )
138: *
139: *     For the nodes on bottom level of the tree, solve
140: *     their subproblems by DLASDQ.
141: *
142:       NDB1 = ( ND+1 ) / 2
143:       NCC = 0
144:       DO 30 I = NDB1, ND
145: *
146: *     IC : center row of each node
147: *     NL : number of rows of left  subproblem
148: *     NR : number of rows of right subproblem
149: *     NLF: starting row of the left   subproblem
150: *     NRF: starting row of the right  subproblem
151: *
152:          I1 = I - 1
153:          IC = IWORK( INODE+I1 )
154:          NL = IWORK( NDIML+I1 )
155:          NLP1 = NL + 1
156:          NR = IWORK( NDIMR+I1 )
157:          NRP1 = NR + 1
158:          NLF = IC - NL
159:          NRF = IC + 1
160:          SQREI = 1
161:          CALL DLASDQ( 'U', SQREI, NL, NLP1, NL, NCC, D( NLF ), E( NLF ),
162:      $                VT( NLF, NLF ), LDVT, U( NLF, NLF ), LDU,
163:      $                U( NLF, NLF ), LDU, WORK, INFO )
164:          IF( INFO.NE.0 ) THEN
165:             RETURN
166:          END IF
167:          ITEMP = IDXQ + NLF - 2
168:          DO 10 J = 1, NL
169:             IWORK( ITEMP+J ) = J
170:    10    CONTINUE
171:          IF( I.EQ.ND ) THEN
172:             SQREI = SQRE
173:          ELSE
174:             SQREI = 1
175:          END IF
176:          NRP1 = NR + SQREI
177:          CALL DLASDQ( 'U', SQREI, NR, NRP1, NR, NCC, D( NRF ), E( NRF ),
178:      $                VT( NRF, NRF ), LDVT, U( NRF, NRF ), LDU,
179:      $                U( NRF, NRF ), LDU, WORK, INFO )
180:          IF( INFO.NE.0 ) THEN
181:             RETURN
182:          END IF
183:          ITEMP = IDXQ + IC
184:          DO 20 J = 1, NR
185:             IWORK( ITEMP+J-1 ) = J
186:    20    CONTINUE
187:    30 CONTINUE
188: *
189: *     Now conquer each subproblem bottom-up.
190: *
191:       DO 50 LVL = NLVL, 1, -1
192: *
193: *        Find the first node LF and last node LL on the
194: *        current level LVL.
195: *
196:          IF( LVL.EQ.1 ) THEN
197:             LF = 1
198:             LL = 1
199:          ELSE
200:             LF = 2**( LVL-1 )
201:             LL = 2*LF - 1
202:          END IF
203:          DO 40 I = LF, LL
204:             IM1 = I - 1
205:             IC = IWORK( INODE+IM1 )
206:             NL = IWORK( NDIML+IM1 )
207:             NR = IWORK( NDIMR+IM1 )
208:             NLF = IC - NL
209:             IF( ( SQRE.EQ.0 ) .AND. ( I.EQ.LL ) ) THEN
210:                SQREI = SQRE
211:             ELSE
212:                SQREI = 1
213:             END IF
214:             IDXQC = IDXQ + NLF - 1
215:             ALPHA = D( IC )
216:             BETA = E( IC )
217:             CALL DLASD1( NL, NR, SQREI, D( NLF ), ALPHA, BETA,
218:      $                   U( NLF, NLF ), LDU, VT( NLF, NLF ), LDVT,
219:      $                   IWORK( IDXQC ), IWORK( IWK ), WORK, INFO )
220:             IF( INFO.NE.0 ) THEN
221:                RETURN
222:             END IF
223:    40    CONTINUE
224:    50 CONTINUE
225: *
226:       RETURN
227: *
228: *     End of DLASD0
229: *
230:       END
231: