142
143
144
145
146
147
148 CHARACTER UPLO
149 INTEGER INFO, KD, LDAB, N
150
151
152 DOUBLE PRECISION AB( LDAB, * )
153
154
155
156
157
158 DOUBLE PRECISION ONE, ZERO
159 parameter( one = 1.0d+0, zero = 0.0d+0 )
160 INTEGER NBMAX, LDWORK
161 parameter( nbmax = 32, ldwork = nbmax+1 )
162
163
164 INTEGER I, I2, I3, IB, II, J, JJ, NB
165
166
167 DOUBLE PRECISION WORK( LDWORK, NBMAX )
168
169
170 LOGICAL LSAME
171 INTEGER ILAENV
173
174
176
177
178 INTRINSIC min
179
180
181
182
183
184 info = 0
185 IF( ( .NOT.
lsame( uplo,
'U' ) ) .AND.
186 $ ( .NOT.
lsame( uplo,
'L' ) ) )
THEN
187 info = -1
188 ELSE IF( n.LT.0 ) THEN
189 info = -2
190 ELSE IF( kd.LT.0 ) THEN
191 info = -3
192 ELSE IF( ldab.LT.kd+1 ) THEN
193 info = -5
194 END IF
195 IF( info.NE.0 ) THEN
196 CALL xerbla(
'DPBTRF', -info )
197 RETURN
198 END IF
199
200
201
202 IF( n.EQ.0 )
203 $ RETURN
204
205
206
207 nb =
ilaenv( 1,
'DPBTRF', uplo, n, kd, -1, -1 )
208
209
210
211
212 nb = min( nb, nbmax )
213
214 IF( nb.LE.1 .OR. nb.GT.kd ) THEN
215
216
217
218 CALL dpbtf2( uplo, n, kd, ab, ldab, info )
219 ELSE
220
221
222
223 IF(
lsame( uplo,
'U' ) )
THEN
224
225
226
227
228
229
230
231 DO 20 j = 1, nb
232 DO 10 i = 1, j - 1
233 work( i, j ) = zero
234 10 CONTINUE
235 20 CONTINUE
236
237
238
239 DO 70 i = 1, n, nb
240 ib = min( nb, n-i+1 )
241
242
243
244 CALL dpotf2( uplo, ib, ab( kd+1, i ), ldab-1, ii )
245 IF( ii.NE.0 ) THEN
246 info = i + ii - 1
247 GO TO 150
248 END IF
249 IF( i+ib.LE.n ) THEN
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265 i2 = min( kd-ib, n-i-ib+1 )
266 i3 = min( ib, n-i-kd+1 )
267
268 IF( i2.GT.0 ) THEN
269
270
271
272 CALL dtrsm(
'Left',
'Upper',
'Transpose',
273 $ 'Non-unit', ib, i2, one, ab( kd+1, i ),
274 $ ldab-1, ab( kd+1-ib, i+ib ), ldab-1 )
275
276
277
278 CALL dsyrk(
'Upper',
'Transpose', i2, ib, -one,
279 $ ab( kd+1-ib, i+ib ), ldab-1, one,
280 $ ab( kd+1, i+ib ), ldab-1 )
281 END IF
282
283 IF( i3.GT.0 ) THEN
284
285
286
287 DO 40 jj = 1, i3
288 DO 30 ii = jj, ib
289 work( ii, jj ) = ab( ii-jj+1, jj+i+kd-1 )
290 30 CONTINUE
291 40 CONTINUE
292
293
294
295 CALL dtrsm(
'Left',
'Upper',
'Transpose',
296 $ 'Non-unit', ib, i3, one, ab( kd+1, i ),
297 $ ldab-1, work, ldwork )
298
299
300
301 IF( i2.GT.0 )
302 $
CALL dgemm(
'Transpose',
'No Transpose', i2, i3,
303 $ ib, -one, ab( kd+1-ib, i+ib ),
304 $ ldab-1, work, ldwork, one,
305 $ ab( 1+ib, i+kd ), ldab-1 )
306
307
308
309 CALL dsyrk(
'Upper',
'Transpose', i3, ib, -one,
310 $ work, ldwork, one, ab( kd+1, i+kd ),
311 $ ldab-1 )
312
313
314
315 DO 60 jj = 1, i3
316 DO 50 ii = jj, ib
317 ab( ii-jj+1, jj+i+kd-1 ) = work( ii, jj )
318 50 CONTINUE
319 60 CONTINUE
320 END IF
321 END IF
322 70 CONTINUE
323 ELSE
324
325
326
327
328
329
330
331 DO 90 j = 1, nb
332 DO 80 i = j + 1, nb
333 work( i, j ) = zero
334 80 CONTINUE
335 90 CONTINUE
336
337
338
339 DO 140 i = 1, n, nb
340 ib = min( nb, n-i+1 )
341
342
343
344 CALL dpotf2( uplo, ib, ab( 1, i ), ldab-1, ii )
345 IF( ii.NE.0 ) THEN
346 info = i + ii - 1
347 GO TO 150
348 END IF
349 IF( i+ib.LE.n ) THEN
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365 i2 = min( kd-ib, n-i-ib+1 )
366 i3 = min( ib, n-i-kd+1 )
367
368 IF( i2.GT.0 ) THEN
369
370
371
372 CALL dtrsm(
'Right',
'Lower',
'Transpose',
373 $ 'Non-unit', i2, ib, one, ab( 1, i ),
374 $ ldab-1, ab( 1+ib, i ), ldab-1 )
375
376
377
378 CALL dsyrk(
'Lower',
'No Transpose', i2, ib, -one,
379 $ ab( 1+ib, i ), ldab-1, one,
380 $ ab( 1, i+ib ), ldab-1 )
381 END IF
382
383 IF( i3.GT.0 ) THEN
384
385
386
387 DO 110 jj = 1, ib
388 DO 100 ii = 1, min( jj, i3 )
389 work( ii, jj ) = ab( kd+1-jj+ii, jj+i-1 )
390 100 CONTINUE
391 110 CONTINUE
392
393
394
395 CALL dtrsm(
'Right',
'Lower',
'Transpose',
396 $ 'Non-unit', i3, ib, one, ab( 1, i ),
397 $ ldab-1, work, ldwork )
398
399
400
401 IF( i2.GT.0 )
402 $
CALL dgemm(
'No transpose',
'Transpose', i3, i2,
403 $ ib, -one, work, ldwork,
404 $ ab( 1+ib, i ), ldab-1, one,
405 $ ab( 1+kd-ib, i+ib ), ldab-1 )
406
407
408
409 CALL dsyrk(
'Lower',
'No Transpose', i3, ib, -one,
410 $ work, ldwork, one, ab( 1, i+kd ),
411 $ ldab-1 )
412
413
414
415 DO 130 jj = 1, ib
416 DO 120 ii = 1, min( jj, i3 )
417 ab( kd+1-jj+ii, jj+i-1 ) = work( ii, jj )
418 120 CONTINUE
419 130 CONTINUE
420 END IF
421 END IF
422 140 CONTINUE
423 END IF
424 END IF
425 RETURN
426
427 150 CONTINUE
428 RETURN
429
430
431
subroutine xerbla(srname, info)
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
subroutine dsyrk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
DSYRK
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
logical function lsame(ca, cb)
LSAME
subroutine dpbtf2(uplo, n, kd, ab, ldab, info)
DPBTF2 computes the Cholesky factorization of a symmetric/Hermitian positive definite band matrix (un...
subroutine dpotf2(uplo, n, a, lda, info)
DPOTF2 computes the Cholesky factorization of a symmetric/Hermitian positive definite matrix (unblock...
subroutine dtrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
DTRSM