158
159
160
161
162
163 IMPLICIT NONE
164
165
166 CHARACTER UPLO
167 INTEGER N, LDA, LTB, LWORK, INFO
168
169
170 INTEGER IPIV( * ), IPIV2( * )
171 COMPLEX*16 A( LDA, * ), TB( * ), WORK( * )
172
173
174
175
176 COMPLEX*16 CZERO, CONE
177 parameter( czero = ( 0.0d+0, 0.0d+0 ),
178 $ cone = ( 1.0d+0, 0.0d+0 ) )
179
180
181 LOGICAL UPPER, TQUERY, WQUERY
182 INTEGER I, J, K, I1, I2, TD
183 INTEGER LDTB, NB, KB, JB, NT, IINFO
184 COMPLEX*16 PIV
185
186
187 LOGICAL LSAME
188 INTEGER ILAENV
190
191
195
196
197 INTRINSIC min, max
198
199
200
201
202
203 info = 0
204 upper =
lsame( uplo,
'U' )
205 wquery = ( lwork.EQ.-1 )
206 tquery = ( ltb.EQ.-1 )
207 IF( .NOT.upper .AND. .NOT.
lsame( uplo,
'L' ) )
THEN
208 info = -1
209 ELSE IF( n.LT.0 ) THEN
210 info = -2
211 ELSE IF( lda.LT.max( 1, n ) ) THEN
212 info = -4
213 ELSE IF ( ltb .LT. 4*n .AND. .NOT.tquery ) THEN
214 info = -6
215 ELSE IF ( lwork .LT. n .AND. .NOT.wquery ) THEN
216 info = -10
217 END IF
218
219 IF( info.NE.0 ) THEN
220 CALL xerbla(
'ZSYTRF_AA_2STAGE', -info )
221 RETURN
222 END IF
223
224
225
226 nb =
ilaenv( 1,
'ZSYTRF_AA_2STAGE', uplo, n, -1, -1, -1 )
227 IF( info.EQ.0 ) THEN
228 IF( tquery ) THEN
229 tb( 1 ) = (3*nb+1)*n
230 END IF
231 IF( wquery ) THEN
232 work( 1 ) = n*nb
233 END IF
234 END IF
235 IF( tquery .OR. wquery ) THEN
236 RETURN
237 END IF
238
239
240
241 IF ( n.EQ.0 ) THEN
242 RETURN
243 ENDIF
244
245
246
247 ldtb = ltb/n
248 IF( ldtb .LT. 3*nb+1 ) THEN
249 nb = (ldtb-1)/3
250 END IF
251 IF( lwork .LT. nb*n ) THEN
252 nb = lwork/n
253 END IF
254
255
256
257 nt = (n+nb-1)/nb
258 td = 2*nb
259 kb = min(nb, n)
260
261
262
263 DO j = 1, kb
264 ipiv( j ) = j
265 END DO
266
267
268
269 tb( 1 ) = nb
270
271 IF( upper ) THEN
272
273
274
275
276
277 DO j = 0, nt-1
278
279
280
281 kb = min(nb, n-j*nb)
282 DO i = 1, j-1
283 IF( i.EQ.1 ) THEN
284
285 IF( i .EQ. (j-1) ) THEN
286 jb = nb+kb
287 ELSE
288 jb = 2*nb
289 END IF
290 CALL zgemm(
'NoTranspose',
'NoTranspose',
291 $ nb, kb, jb,
292 $ cone, tb( td+1 + (i*nb)*ldtb ), ldtb-1,
293 $ a( (i-1)*nb+1, j*nb+1 ), lda,
294 $ czero, work( i*nb+1 ), n )
295 ELSE
296
297 IF( i .EQ. (j-1) ) THEN
298 jb = 2*nb+kb
299 ELSE
300 jb = 3*nb
301 END IF
302 CALL zgemm(
'NoTranspose',
'NoTranspose',
303 $ nb, kb, jb,
304 $ cone, tb( td+nb+1 + ((i-1)*nb)*ldtb ),
305 $ ldtb-1,
306 $ a( (i-2)*nb+1, j*nb+1 ), lda,
307 $ czero, work( i*nb+1 ), n )
308 END IF
309 END DO
310
311
312
313 CALL zlacpy(
'Upper', kb, kb, a( j*nb+1, j*nb+1 ), lda,
314 $ tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
315 IF( j.GT.1 ) THEN
316
317 CALL zgemm(
'Transpose',
'NoTranspose',
318 $ kb, kb, (j-1)*nb,
319 $ -cone, a( 1, j*nb+1 ), lda,
320 $ work( nb+1 ), n,
321 $ cone, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
322
323 CALL zgemm(
'Transpose',
'NoTranspose',
324 $ kb, nb, kb,
325 $ cone, a( (j-1)*nb+1, j*nb+1 ), lda,
326 $ tb( td+nb+1 + ((j-1)*nb)*ldtb ), ldtb-1,
327 $ czero, work( 1 ), n )
328 CALL zgemm(
'NoTranspose',
'NoTranspose',
329 $ kb, kb, nb,
330 $ -cone, work( 1 ), n,
331 $ a( (j-2)*nb+1, j*nb+1 ), lda,
332 $ cone, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
333 END IF
334
335
336
337 DO i = 1, kb
338 DO k = i+1, kb
339 tb( td+(k-i)+1 + (j*nb+i-1)*ldtb )
340 $ = tb( td-(k-(i+1)) + (j*nb+k-1)*ldtb )
341 END DO
342 END DO
343 IF( j.GT.0 ) THEN
344
345
346
347 CALL ztrsm(
'L',
'U',
'T',
'N', kb, kb, cone,
348 $ a( (j-1)*nb+1, j*nb+1 ), lda,
349 $ tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
350 CALL ztrsm(
'R',
'U',
'N',
'N', kb, kb, cone,
351 $ a( (j-1)*nb+1, j*nb+1 ), lda,
352 $ tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
353 END IF
354
355 IF( j.LT.nt-1 ) THEN
356 IF( j.GT.0 ) THEN
357
358
359
360 IF( j.EQ.1 ) THEN
361 CALL zgemm(
'NoTranspose',
'NoTranspose',
362 $ kb, kb, kb,
363 $ cone, tb( td+1 + (j*nb)*ldtb ), ldtb-1,
364 $ a( (j-1)*nb+1, j*nb+1 ), lda,
365 $ czero, work( j*nb+1 ), n )
366 ELSE
367 CALL zgemm(
'NoTranspose',
'NoTranspose',
368 $ kb, kb, nb+kb,
369 $ cone, tb( td+nb+1 + ((j-1)*nb)*ldtb ),
370 $ ldtb-1,
371 $ a( (j-2)*nb+1, j*nb+1 ), lda,
372 $ czero, work( j*nb+1 ), n )
373 END IF
374
375
376
377 CALL zgemm(
'Transpose',
'NoTranspose',
378 $ nb, n-(j+1)*nb, j*nb,
379 $ -cone, work( nb+1 ), n,
380 $ a( 1, (j+1)*nb+1 ), lda,
381 $ cone, a( j*nb+1, (j+1)*nb+1 ), lda )
382 END IF
383
384
385
386 DO k = 1, nb
387 CALL zcopy( n-(j+1)*nb,
388 $ a( j*nb+k, (j+1)*nb+1 ), lda,
389 $ work( 1+(k-1)*n ), 1 )
390 END DO
391
392
393
394 CALL zgetrf( n-(j+1)*nb, nb,
395 $ work, n,
396 $ ipiv( (j+1)*nb+1 ), iinfo )
397
398
399
400
401
402
403 DO k = 1, nb
404 CALL zcopy( n-(j+1)*nb,
405 $ work( 1+(k-1)*n ), 1,
406 $ a( j*nb+k, (j+1)*nb+1 ), lda )
407 END DO
408
409
410
411 kb = min(nb, n-(j+1)*nb)
412 CALL zlaset(
'Full', kb, nb, czero, czero,
413 $ tb( td+nb+1 + (j*nb)*ldtb), ldtb-1 )
414 CALL zlacpy(
'Upper', kb, nb,
415 $ work, n,
416 $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
417 IF( j.GT.0 ) THEN
418 CALL ztrsm(
'R',
'U',
'N',
'U', kb, nb, cone,
419 $ a( (j-1)*nb+1, j*nb+1 ), lda,
420 $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
421 END IF
422
423
424
425
426 DO k = 1, nb
427 DO i = 1, kb
428 tb( td-nb+k-i+1 + (j*nb+nb+i-1)*ldtb )
429 $ = tb( td+nb+i-k+1 + (j*nb+k-1)*ldtb )
430 END DO
431 END DO
432 CALL zlaset(
'Lower', kb, nb, czero, cone,
433 $ a( j*nb+1, (j+1)*nb+1), lda )
434
435
436
437 DO k = 1, kb
438
439 ipiv( (j+1)*nb+k ) = ipiv( (j+1)*nb+k ) + (j+1)*nb
440
441 i1 = (j+1)*nb+k
442 i2 = ipiv( (j+1)*nb+k )
443 IF( i1.NE.i2 ) THEN
444
445 CALL zswap( k-1, a( (j+1)*nb+1, i1 ), 1,
446 $ a( (j+1)*nb+1, i2 ), 1 )
447
448 IF( i2.GT.(i1+1) )
449 $
CALL zswap( i2-i1-1, a( i1, i1+1 ), lda,
450 $ a( i1+1, i2 ), 1 )
451
452 IF( i2.LT.n )
453 $
CALL zswap( n-i2, a( i1, i2+1 ), lda,
454 $ a( i2, i2+1 ), lda )
455
456 piv = a( i1, i1 )
457 a( i1, i1 ) = a( i2, i2 )
458 a( i2, i2 ) = piv
459
460 IF( j.GT.0 ) THEN
461 CALL zswap( j*nb, a( 1, i1 ), 1,
462 $ a( 1, i2 ), 1 )
463 END IF
464 ENDIF
465 END DO
466 END IF
467 END DO
468 ELSE
469
470
471
472
473
474 DO j = 0, nt-1
475
476
477
478 kb = min(nb, n-j*nb)
479 DO i = 1, j-1
480 IF( i.EQ.1 ) THEN
481
482 IF( i .EQ. (j-1) ) THEN
483 jb = nb+kb
484 ELSE
485 jb = 2*nb
486 END IF
487 CALL zgemm(
'NoTranspose',
'Transpose',
488 $ nb, kb, jb,
489 $ cone, tb( td+1 + (i*nb)*ldtb ), ldtb-1,
490 $ a( j*nb+1, (i-1)*nb+1 ), lda,
491 $ czero, work( i*nb+1 ), n )
492 ELSE
493
494 IF( i .EQ. (j-1) ) THEN
495 jb = 2*nb+kb
496 ELSE
497 jb = 3*nb
498 END IF
499 CALL zgemm(
'NoTranspose',
'Transpose',
500 $ nb, kb, jb,
501 $ cone, tb( td+nb+1 + ((i-1)*nb)*ldtb ),
502 $ ldtb-1,
503 $ a( j*nb+1, (i-2)*nb+1 ), lda,
504 $ czero, work( i*nb+1 ), n )
505 END IF
506 END DO
507
508
509
510 CALL zlacpy(
'Lower', kb, kb, a( j*nb+1, j*nb+1 ), lda,
511 $ tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
512 IF( j.GT.1 ) THEN
513
514 CALL zgemm(
'NoTranspose',
'NoTranspose',
515 $ kb, kb, (j-1)*nb,
516 $ -cone, a( j*nb+1, 1 ), lda,
517 $ work( nb+1 ), n,
518 $ cone, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
519
520 CALL zgemm(
'NoTranspose',
'NoTranspose',
521 $ kb, nb, kb,
522 $ cone, a( j*nb+1, (j-1)*nb+1 ), lda,
523 $ tb( td+nb+1 + ((j-1)*nb)*ldtb ), ldtb-1,
524 $ czero, work( 1 ), n )
525 CALL zgemm(
'NoTranspose',
'Transpose',
526 $ kb, kb, nb,
527 $ -cone, work( 1 ), n,
528 $ a( j*nb+1, (j-2)*nb+1 ), lda,
529 $ cone, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
530 END IF
531
532
533
534 DO i = 1, kb
535 DO k = i+1, kb
536 tb( td-(k-(i+1)) + (j*nb+k-1)*ldtb )
537 $ = tb( td+(k-i)+1 + (j*nb+i-1)*ldtb )
538 END DO
539 END DO
540 IF( j.GT.0 ) THEN
541
542
543
544 CALL ztrsm(
'L',
'L',
'N',
'N', kb, kb, cone,
545 $ a( j*nb+1, (j-1)*nb+1 ), lda,
546 $ tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
547 CALL ztrsm(
'R',
'L',
'T',
'N', kb, kb, cone,
548 $ a( j*nb+1, (j-1)*nb+1 ), lda,
549 $ tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
550 END IF
551
552
553
554 DO i = 1, kb
555 DO k = i+1, kb
556 tb( td-(k-(i+1)) + (j*nb+k-1)*ldtb )
557 $ = tb( td+(k-i)+1 + (j*nb+i-1)*ldtb )
558 END DO
559 END DO
560
561 IF( j.LT.nt-1 ) THEN
562 IF( j.GT.0 ) THEN
563
564
565
566 IF( j.EQ.1 ) THEN
567 CALL zgemm(
'NoTranspose',
'Transpose',
568 $ kb, kb, kb,
569 $ cone, tb( td+1 + (j*nb)*ldtb ), ldtb-1,
570 $ a( j*nb+1, (j-1)*nb+1 ), lda,
571 $ czero, work( j*nb+1 ), n )
572 ELSE
573 CALL zgemm(
'NoTranspose',
'Transpose',
574 $ kb, kb, nb+kb,
575 $ cone, tb( td+nb+1 + ((j-1)*nb)*ldtb ),
576 $ ldtb-1,
577 $ a( j*nb+1, (j-2)*nb+1 ), lda,
578 $ czero, work( j*nb+1 ), n )
579 END IF
580
581
582
583 CALL zgemm(
'NoTranspose',
'NoTranspose',
584 $ n-(j+1)*nb, nb, j*nb,
585 $ -cone, a( (j+1)*nb+1, 1 ), lda,
586 $ work( nb+1 ), n,
587 $ cone, a( (j+1)*nb+1, j*nb+1 ), lda )
588 END IF
589
590
591
592 CALL zgetrf( n-(j+1)*nb, nb,
593 $ a( (j+1)*nb+1, j*nb+1 ), lda,
594 $ ipiv( (j+1)*nb+1 ), iinfo )
595
596
597
598
599
600
601 kb = min(nb, n-(j+1)*nb)
602 CALL zlaset(
'Full', kb, nb, czero, czero,
603 $ tb( td+nb+1 + (j*nb)*ldtb), ldtb-1 )
604 CALL zlacpy(
'Upper', kb, nb,
605 $ a( (j+1)*nb+1, j*nb+1 ), lda,
606 $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
607 IF( j.GT.0 ) THEN
608 CALL ztrsm(
'R',
'L',
'T',
'U', kb, nb, cone,
609 $ a( j*nb+1, (j-1)*nb+1 ), lda,
610 $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
611 END IF
612
613
614
615
616 DO k = 1, nb
617 DO i = 1, kb
618 tb( td-nb+k-i+1 + (j*nb+nb+i-1)*ldtb ) =
619 $ tb( td+nb+i-k+1 + (j*nb+k-1)*ldtb )
620 END DO
621 END DO
622 CALL zlaset(
'Upper', kb, nb, czero, cone,
623 $ a( (j+1)*nb+1, j*nb+1 ), lda )
624
625
626
627 DO k = 1, kb
628
629 ipiv( (j+1)*nb+k ) = ipiv( (j+1)*nb+k ) + (j+1)*nb
630
631 i1 = (j+1)*nb+k
632 i2 = ipiv( (j+1)*nb+k )
633 IF( i1.NE.i2 ) THEN
634
635 CALL zswap( k-1, a( i1, (j+1)*nb+1 ), lda,
636 $ a( i2, (j+1)*nb+1 ), lda )
637
638 IF( i2.GT.(i1+1) )
639 $
CALL zswap( i2-i1-1, a( i1+1, i1 ), 1,
640 $ a( i2, i1+1 ), lda )
641
642 IF( i2.LT.n )
643 $
CALL zswap( n-i2, a( i2+1, i1 ), 1,
644 $ a( i2+1, i2 ), 1 )
645
646 piv = a( i1, i1 )
647 a( i1, i1 ) = a( i2, i2 )
648 a( i2, i2 ) = piv
649
650 IF( j.GT.0 ) THEN
651 CALL zswap( j*nb, a( i1, 1 ), lda,
652 $ a( i2, 1 ), lda )
653 END IF
654 ENDIF
655 END DO
656
657
658
659
660
661 END IF
662 END DO
663 END IF
664
665
666 CALL zgbtrf( n, n, nb, nb, tb, ldtb, ipiv2, info )
667
668 RETURN
669
670
671
subroutine xerbla(srname, info)
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
subroutine zgbtrf(m, n, kl, ku, ab, ldab, ipiv, info)
ZGBTRF
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
subroutine zgetrf(m, n, a, lda, ipiv, info)
ZGETRF
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine zlaset(uplo, m, n, alpha, beta, a, lda)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
subroutine zlaswp(n, a, lda, k1, k2, ipiv, incx)
ZLASWP performs a series of row interchanges on a general rectangular matrix.
logical function lsame(ca, cb)
LSAME
subroutine zswap(n, zx, incx, zy, incy)
ZSWAP
subroutine ztrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
ZTRSM