LAPACK  3.8.0 LAPACK: Linear Algebra PACKage

## ◆ zlasyf_rook()

 subroutine zlasyf_rook ( character UPLO, integer N, integer NB, integer KB, complex*16, dimension( lda, * ) A, integer LDA, integer, dimension( * ) IPIV, complex*16, dimension( ldw, * ) W, integer LDW, integer INFO )

ZLASYF_ROOK computes a partial factorization of a complex symmetric matrix using the bounded Bunch-Kaufman ("rook") diagonal pivoting method.

Purpose:
``` ZLASYF_ROOK computes a partial factorization of a complex symmetric
matrix A using the bounded Bunch-Kaufman ("rook") diagonal
pivoting method. The partial factorization has the form:

A  =  ( I  U12 ) ( A11  0  ) (  I       0    )  if UPLO = 'U', or:
( 0  U22 ) (  0   D  ) ( U12**T U22**T )

A  =  ( L11  0 ) (  D   0  ) ( L11**T L21**T )  if UPLO = 'L'
( L21  I ) (  0  A22 ) (  0       I    )

where the order of D is at most NB. The actual order is returned in
the argument KB, and is either NB or NB-1, or N if N <= NB.

ZLASYF_ROOK is an auxiliary routine called by ZSYTRF_ROOK. It uses
blocked code (calling Level 3 BLAS) to update the submatrix
A11 (if UPLO = 'U') or A22 (if UPLO = 'L').```
Parameters
 [in] UPLO ``` UPLO is CHARACTER*1 Specifies whether the upper or lower triangular part of the symmetric matrix A is stored: = 'U': Upper triangular = 'L': Lower triangular``` [in] N ``` N is INTEGER The order of the matrix A. N >= 0.``` [in] NB ``` NB is INTEGER The maximum number of columns of the matrix A that should be factored. NB should be at least 2 to allow for 2-by-2 pivot blocks.``` [out] KB ``` KB is INTEGER The number of columns of A that were actually factored. KB is either NB-1 or NB, or N if N <= NB.``` [in,out] A ``` A is COMPLEX*16 array, dimension (LDA,N) On entry, the symmetric matrix A. If UPLO = 'U', the leading n-by-n upper triangular part of A contains the upper triangular part of the matrix A, and the strictly lower triangular part of A is not referenced. If UPLO = 'L', the leading n-by-n lower triangular part of A contains the lower triangular part of the matrix A, and the strictly upper triangular part of A is not referenced. On exit, A contains details of the partial factorization.``` [in] LDA ``` LDA is INTEGER The leading dimension of the array A. LDA >= max(1,N).``` [out] IPIV ``` IPIV is INTEGER array, dimension (N) Details of the interchanges and the block structure of D. If UPLO = 'U': Only the last KB elements of IPIV are set. If IPIV(k) > 0, then rows and columns k and IPIV(k) were interchanged and D(k,k) is a 1-by-1 diagonal block. If IPIV(k) < 0 and IPIV(k-1) < 0, then rows and columns k and -IPIV(k) were interchanged and rows and columns k-1 and -IPIV(k-1) were inerchaged, D(k-1:k,k-1:k) is a 2-by-2 diagonal block. If UPLO = 'L': Only the first KB elements of IPIV are set. If IPIV(k) > 0, then rows and columns k and IPIV(k) were interchanged and D(k,k) is a 1-by-1 diagonal block. If IPIV(k) < 0 and IPIV(k+1) < 0, then rows and columns k and -IPIV(k) were interchanged and rows and columns k+1 and -IPIV(k+1) were inerchaged, D(k:k+1,k:k+1) is a 2-by-2 diagonal block.``` [out] W ` W is COMPLEX*16 array, dimension (LDW,NB)` [in] LDW ``` LDW is INTEGER The leading dimension of the array W. LDW >= max(1,N).``` [out] INFO ``` INFO is INTEGER = 0: successful exit > 0: if INFO = k, D(k,k) is exactly zero. The factorization has been completed, but the block diagonal matrix D is exactly singular.```
Date
November 2013
Contributors:
```  November 2013,  Igor Kozachenko,
Computer Science Division,
University of California, Berkeley

September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
School of Mathematics,
University of Manchester```

Definition at line 186 of file zlasyf_rook.f.

186 *
187 * -- LAPACK computational routine (version 3.5.0) --
188 * -- LAPACK is a software package provided by Univ. of Tennessee, --
189 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
190 * November 2013
191 *
192 * .. Scalar Arguments ..
193  CHARACTER uplo
194  INTEGER info, kb, lda, ldw, n, nb
195 * ..
196 * .. Array Arguments ..
197  INTEGER ipiv( * )
198  COMPLEX*16 a( lda, * ), w( ldw, * )
199 * ..
200 *
201 * =====================================================================
202 *
203 * .. Parameters ..
204  DOUBLE PRECISION zero, one
205  parameter( zero = 0.0d+0, one = 1.0d+0 )
206  DOUBLE PRECISION eight, sevten
207  parameter( eight = 8.0d+0, sevten = 17.0d+0 )
208  COMPLEX*16 cone, czero
209  parameter( cone = ( 1.0d+0, 0.0d+0 ),
210  \$ czero = ( 0.0d+0, 0.0d+0 ) )
211 * ..
212 * .. Local Scalars ..
213  LOGICAL done
214  INTEGER imax, itemp, j, jb, jj, jmax, jp1, jp2, k, kk,
215  \$ kw, kkw, kp, kstep, p, ii
216  DOUBLE PRECISION absakk, alpha, colmax, rowmax, dtemp, sfmin
217  COMPLEX*16 d11, d12, d21, d22, r1, t, z
218 * ..
219 * .. External Functions ..
220  LOGICAL lsame
221  INTEGER izamax
222  DOUBLE PRECISION dlamch
223  EXTERNAL lsame, izamax, dlamch
224 * ..
225 * .. External Subroutines ..
226  EXTERNAL zcopy, zgemm, zgemv, zscal, zswap
227 * ..
228 * .. Intrinsic Functions ..
229  INTRINSIC abs, max, min, sqrt, dimag, dble
230 * ..
231 * .. Statement Functions ..
232  DOUBLE PRECISION cabs1
233 * ..
234 * .. Statement Function definitions ..
235  cabs1( z ) = abs( dble( z ) ) + abs( dimag( z ) )
236 * ..
237 * .. Executable Statements ..
238 *
239  info = 0
240 *
241 * Initialize ALPHA for use in choosing pivot block size.
242 *
243  alpha = ( one+sqrt( sevten ) ) / eight
244 *
245 * Compute machine safe minimum
246 *
247  sfmin = dlamch( 'S' )
248 *
249  IF( lsame( uplo, 'U' ) ) THEN
250 *
251 * Factorize the trailing columns of A using the upper triangle
252 * of A and working backwards, and compute the matrix W = U12*D
253 * for use in updating A11
254 *
255 * K is the main loop index, decreasing from N in steps of 1 or 2
256 *
257  k = n
258  10 CONTINUE
259 *
260 * KW is the column of W which corresponds to column K of A
261 *
262  kw = nb + k - n
263 *
264 * Exit from loop
265 *
266  IF( ( k.LE.n-nb+1 .AND. nb.LT.n ) .OR. k.LT.1 )
267  \$ GO TO 30
268 *
269  kstep = 1
270  p = k
271 *
272 * Copy column K of A to column KW of W and update it
273 *
274  CALL zcopy( k, a( 1, k ), 1, w( 1, kw ), 1 )
275  IF( k.LT.n )
276  \$ CALL zgemv( 'No transpose', k, n-k, -cone, a( 1, k+1 ),
277  \$ lda, w( k, kw+1 ), ldw, cone, w( 1, kw ), 1 )
278 *
279 * Determine rows and columns to be interchanged and whether
280 * a 1-by-1 or 2-by-2 pivot block will be used
281 *
282  absakk = cabs1( w( k, kw ) )
283 *
284 * IMAX is the row-index of the largest off-diagonal element in
285 * column K, and COLMAX is its absolute value.
286 * Determine both COLMAX and IMAX.
287 *
288  IF( k.GT.1 ) THEN
289  imax = izamax( k-1, w( 1, kw ), 1 )
290  colmax = cabs1( w( imax, kw ) )
291  ELSE
292  colmax = zero
293  END IF
294 *
295  IF( max( absakk, colmax ).EQ.zero ) THEN
296 *
297 * Column K is zero or underflow: set INFO and continue
298 *
299  IF( info.EQ.0 )
300  \$ info = k
301  kp = k
302  CALL zcopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
303  ELSE
304 *
305 * ============================================================
306 *
307 * Test for interchange
308 *
309 * Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
310 * (used to handle NaN and Inf)
311 *
312  IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
313 *
314 * no interchange, use 1-by-1 pivot block
315 *
316  kp = k
317 *
318  ELSE
319 *
320  done = .false.
321 *
322 * Loop until pivot found
323 *
324  12 CONTINUE
325 *
326 * Begin pivot search loop body
327 *
328 *
329 * Copy column IMAX to column KW-1 of W and update it
330 *
331  CALL zcopy( imax, a( 1, imax ), 1, w( 1, kw-1 ), 1 )
332  CALL zcopy( k-imax, a( imax, imax+1 ), lda,
333  \$ w( imax+1, kw-1 ), 1 )
334 *
335  IF( k.LT.n )
336  \$ CALL zgemv( 'No transpose', k, n-k, -cone,
337  \$ a( 1, k+1 ), lda, w( imax, kw+1 ), ldw,
338  \$ cone, w( 1, kw-1 ), 1 )
339 *
340 * JMAX is the column-index of the largest off-diagonal
341 * element in row IMAX, and ROWMAX is its absolute value.
342 * Determine both ROWMAX and JMAX.
343 *
344  IF( imax.NE.k ) THEN
345  jmax = imax + izamax( k-imax, w( imax+1, kw-1 ),
346  \$ 1 )
347  rowmax = cabs1( w( jmax, kw-1 ) )
348  ELSE
349  rowmax = zero
350  END IF
351 *
352  IF( imax.GT.1 ) THEN
353  itemp = izamax( imax-1, w( 1, kw-1 ), 1 )
354  dtemp = cabs1( w( itemp, kw-1 ) )
355  IF( dtemp.GT.rowmax ) THEN
356  rowmax = dtemp
357  jmax = itemp
358  END IF
359  END IF
360 *
361 * Equivalent to testing for
362 * CABS1( W( IMAX, KW-1 ) ).GE.ALPHA*ROWMAX
363 * (used to handle NaN and Inf)
364 *
365  IF( .NOT.(cabs1( w( imax, kw-1 ) ).LT.alpha*rowmax ) )
366  \$ THEN
367 *
368 * interchange rows and columns K and IMAX,
369 * use 1-by-1 pivot block
370 *
371  kp = imax
372 *
373 * copy column KW-1 of W to column KW of W
374 *
375  CALL zcopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
376 *
377  done = .true.
378 *
379 * Equivalent to testing for ROWMAX.EQ.COLMAX,
380 * (used to handle NaN and Inf)
381 *
382  ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
383  \$ THEN
384 *
385 * interchange rows and columns K-1 and IMAX,
386 * use 2-by-2 pivot block
387 *
388  kp = imax
389  kstep = 2
390  done = .true.
391  ELSE
392 *
394 *
395  p = imax
396  colmax = rowmax
397  imax = jmax
398 *
399 * Copy updated JMAXth (next IMAXth) column to Kth of W
400 *
401  CALL zcopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
402 *
403  END IF
404 *
405 * End pivot search loop body
406 *
407  IF( .NOT. done ) GOTO 12
408 *
409  END IF
410 *
411 * ============================================================
412 *
413  kk = k - kstep + 1
414 *
415 * KKW is the column of W which corresponds to column KK of A
416 *
417  kkw = nb + kk - n
418 *
419  IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
420 *
421 * Copy non-updated column K to column P
422 *
423  CALL zcopy( k-p, a( p+1, k ), 1, a( p, p+1 ), lda )
424  CALL zcopy( p, a( 1, k ), 1, a( 1, p ), 1 )
425 *
426 * Interchange rows K and P in last N-K+1 columns of A
427 * and last N-K+2 columns of W
428 *
429  CALL zswap( n-k+1, a( k, k ), lda, a( p, k ), lda )
430  CALL zswap( n-kk+1, w( k, kkw ), ldw, w( p, kkw ), ldw )
431  END IF
432 *
433 * Updated column KP is already stored in column KKW of W
434 *
435  IF( kp.NE.kk ) THEN
436 *
437 * Copy non-updated column KK to column KP
438 *
439  a( kp, k ) = a( kk, k )
440  CALL zcopy( k-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
441  \$ lda )
442  CALL zcopy( kp, a( 1, kk ), 1, a( 1, kp ), 1 )
443 *
444 * Interchange rows KK and KP in last N-KK+1 columns
445 * of A and W
446 *
447  CALL zswap( n-kk+1, a( kk, kk ), lda, a( kp, kk ), lda )
448  CALL zswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),
449  \$ ldw )
450  END IF
451 *
452  IF( kstep.EQ.1 ) THEN
453 *
454 * 1-by-1 pivot block D(k): column KW of W now holds
455 *
456 * W(k) = U(k)*D(k)
457 *
458 * where U(k) is the k-th column of U
459 *
460 * Store U(k) in column k of A
461 *
462  CALL zcopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
463  IF( k.GT.1 ) THEN
464  IF( cabs1( a( k, k ) ).GE.sfmin ) THEN
465  r1 = cone / a( k, k )
466  CALL zscal( k-1, r1, a( 1, k ), 1 )
467  ELSE IF( a( k, k ).NE.czero ) THEN
468  DO 14 ii = 1, k - 1
469  a( ii, k ) = a( ii, k ) / a( k, k )
470  14 CONTINUE
471  END IF
472  END IF
473 *
474  ELSE
475 *
476 * 2-by-2 pivot block D(k): columns KW and KW-1 of W now
477 * hold
478 *
479 * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
480 *
481 * where U(k) and U(k-1) are the k-th and (k-1)-th columns
482 * of U
483 *
484  IF( k.GT.2 ) THEN
485 *
486 * Store U(k) and U(k-1) in columns k and k-1 of A
487 *
488  d12 = w( k-1, kw )
489  d11 = w( k, kw ) / d12
490  d22 = w( k-1, kw-1 ) / d12
491  t = cone / ( d11*d22-cone )
492  DO 20 j = 1, k - 2
493  a( j, k-1 ) = t*( (d11*w( j, kw-1 )-w( j, kw ) ) /
494  \$ d12 )
495  a( j, k ) = t*( ( d22*w( j, kw )-w( j, kw-1 ) ) /
496  \$ d12 )
497  20 CONTINUE
498  END IF
499 *
500 * Copy D(k) to A
501 *
502  a( k-1, k-1 ) = w( k-1, kw-1 )
503  a( k-1, k ) = w( k-1, kw )
504  a( k, k ) = w( k, kw )
505  END IF
506  END IF
507 *
508 * Store details of the interchanges in IPIV
509 *
510  IF( kstep.EQ.1 ) THEN
511  ipiv( k ) = kp
512  ELSE
513  ipiv( k ) = -p
514  ipiv( k-1 ) = -kp
515  END IF
516 *
517 * Decrease K and return to the start of the main loop
518 *
519  k = k - kstep
520  GO TO 10
521 *
522  30 CONTINUE
523 *
524 * Update the upper triangle of A11 (= A(1:k,1:k)) as
525 *
526 * A11 := A11 - U12*D*U12**T = A11 - U12*W**T
527 *
528 * computing blocks of NB columns at a time
529 *
530  DO 50 j = ( ( k-1 ) / nb )*nb + 1, 1, -nb
531  jb = min( nb, k-j+1 )
532 *
533 * Update the upper triangle of the diagonal block
534 *
535  DO 40 jj = j, j + jb - 1
536  CALL zgemv( 'No transpose', jj-j+1, n-k, -cone,
537  \$ a( j, k+1 ), lda, w( jj, kw+1 ), ldw, cone,
538  \$ a( j, jj ), 1 )
539  40 CONTINUE
540 *
541 * Update the rectangular superdiagonal block
542 *
543  IF( j.GE.2 )
544  \$ CALL zgemm( 'No transpose', 'Transpose', j-1, jb,
545  \$ n-k, -cone, a( 1, k+1 ), lda, w( j, kw+1 ), ldw,
546  \$ cone, a( 1, j ), lda )
547  50 CONTINUE
548 *
549 * Put U12 in standard form by partially undoing the interchanges
550 * in columns k+1:n
551 *
552  j = k + 1
553  60 CONTINUE
554 *
555  kstep = 1
556  jp1 = 1
557  jj = j
558  jp2 = ipiv( j )
559  IF( jp2.LT.0 ) THEN
560  jp2 = -jp2
561  j = j + 1
562  jp1 = -ipiv( j )
563  kstep = 2
564  END IF
565 *
566  j = j + 1
567  IF( jp2.NE.jj .AND. j.LE.n )
568  \$ CALL zswap( n-j+1, a( jp2, j ), lda, a( jj, j ), lda )
569  jj = j - 1
570  IF( jp1.NE.jj .AND. kstep.EQ.2 )
571  \$ CALL zswap( n-j+1, a( jp1, j ), lda, a( jj, j ), lda )
572  IF( j.LE.n )
573  \$ GO TO 60
574 *
575 * Set KB to the number of columns factorized
576 *
577  kb = n - k
578 *
579  ELSE
580 *
581 * Factorize the leading columns of A using the lower triangle
582 * of A and working forwards, and compute the matrix W = L21*D
583 * for use in updating A22
584 *
585 * K is the main loop index, increasing from 1 in steps of 1 or 2
586 *
587  k = 1
588  70 CONTINUE
589 *
590 * Exit from loop
591 *
592  IF( ( k.GE.nb .AND. nb.LT.n ) .OR. k.GT.n )
593  \$ GO TO 90
594 *
595  kstep = 1
596  p = k
597 *
598 * Copy column K of A to column K of W and update it
599 *
600  CALL zcopy( n-k+1, a( k, k ), 1, w( k, k ), 1 )
601  IF( k.GT.1 )
602  \$ CALL zgemv( 'No transpose', n-k+1, k-1, -cone, a( k, 1 ),
603  \$ lda, w( k, 1 ), ldw, cone, w( k, k ), 1 )
604 *
605 * Determine rows and columns to be interchanged and whether
606 * a 1-by-1 or 2-by-2 pivot block will be used
607 *
608  absakk = cabs1( w( k, k ) )
609 *
610 * IMAX is the row-index of the largest off-diagonal element in
611 * column K, and COLMAX is its absolute value.
612 * Determine both COLMAX and IMAX.
613 *
614  IF( k.LT.n ) THEN
615  imax = k + izamax( n-k, w( k+1, k ), 1 )
616  colmax = cabs1( w( imax, k ) )
617  ELSE
618  colmax = zero
619  END IF
620 *
621  IF( max( absakk, colmax ).EQ.zero ) THEN
622 *
623 * Column K is zero or underflow: set INFO and continue
624 *
625  IF( info.EQ.0 )
626  \$ info = k
627  kp = k
628  CALL zcopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
629  ELSE
630 *
631 * ============================================================
632 *
633 * Test for interchange
634 *
635 * Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
636 * (used to handle NaN and Inf)
637 *
638  IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
639 *
640 * no interchange, use 1-by-1 pivot block
641 *
642  kp = k
643 *
644  ELSE
645 *
646  done = .false.
647 *
648 * Loop until pivot found
649 *
650  72 CONTINUE
651 *
652 * Begin pivot search loop body
653 *
654 *
655 * Copy column IMAX to column K+1 of W and update it
656 *
657  CALL zcopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1)
658  CALL zcopy( n-imax+1, a( imax, imax ), 1,
659  \$ w( imax, k+1 ), 1 )
660  IF( k.GT.1 )
661  \$ CALL zgemv( 'No transpose', n-k+1, k-1, -cone,
662  \$ a( k, 1 ), lda, w( imax, 1 ), ldw,
663  \$ cone, w( k, k+1 ), 1 )
664 *
665 * JMAX is the column-index of the largest off-diagonal
666 * element in row IMAX, and ROWMAX is its absolute value.
667 * Determine both ROWMAX and JMAX.
668 *
669  IF( imax.NE.k ) THEN
670  jmax = k - 1 + izamax( imax-k, w( k, k+1 ), 1 )
671  rowmax = cabs1( w( jmax, k+1 ) )
672  ELSE
673  rowmax = zero
674  END IF
675 *
676  IF( imax.LT.n ) THEN
677  itemp = imax + izamax( n-imax, w( imax+1, k+1 ), 1)
678  dtemp = cabs1( w( itemp, k+1 ) )
679  IF( dtemp.GT.rowmax ) THEN
680  rowmax = dtemp
681  jmax = itemp
682  END IF
683  END IF
684 *
685 * Equivalent to testing for
686 * CABS1( W( IMAX, K+1 ) ).GE.ALPHA*ROWMAX
687 * (used to handle NaN and Inf)
688 *
689  IF( .NOT.( cabs1( w( imax, k+1 ) ).LT.alpha*rowmax ) )
690  \$ THEN
691 *
692 * interchange rows and columns K and IMAX,
693 * use 1-by-1 pivot block
694 *
695  kp = imax
696 *
697 * copy column K+1 of W to column K of W
698 *
699  CALL zcopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
700 *
701  done = .true.
702 *
703 * Equivalent to testing for ROWMAX.EQ.COLMAX,
704 * (used to handle NaN and Inf)
705 *
706  ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
707  \$ THEN
708 *
709 * interchange rows and columns K+1 and IMAX,
710 * use 2-by-2 pivot block
711 *
712  kp = imax
713  kstep = 2
714  done = .true.
715  ELSE
716 *
718 *
719  p = imax
720  colmax = rowmax
721  imax = jmax
722 *
723 * Copy updated JMAXth (next IMAXth) column to Kth of W
724 *
725  CALL zcopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
726 *
727  END IF
728 *
729 * End pivot search loop body
730 *
731  IF( .NOT. done ) GOTO 72
732 *
733  END IF
734 *
735 * ============================================================
736 *
737  kk = k + kstep - 1
738 *
739  IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
740 *
741 * Copy non-updated column K to column P
742 *
743  CALL zcopy( p-k, a( k, k ), 1, a( p, k ), lda )
744  CALL zcopy( n-p+1, a( p, k ), 1, a( p, p ), 1 )
745 *
746 * Interchange rows K and P in first K columns of A
747 * and first K+1 columns of W
748 *
749  CALL zswap( k, a( k, 1 ), lda, a( p, 1 ), lda )
750  CALL zswap( kk, w( k, 1 ), ldw, w( p, 1 ), ldw )
751  END IF
752 *
753 * Updated column KP is already stored in column KK of W
754 *
755  IF( kp.NE.kk ) THEN
756 *
757 * Copy non-updated column KK to column KP
758 *
759  a( kp, k ) = a( kk, k )
760  CALL zcopy( kp-k-1, a( k+1, kk ), 1, a( kp, k+1 ), lda )
761  CALL zcopy( n-kp+1, a( kp, kk ), 1, a( kp, kp ), 1 )
762 *
763 * Interchange rows KK and KP in first KK columns of A and W
764 *
765  CALL zswap( kk, a( kk, 1 ), lda, a( kp, 1 ), lda )
766  CALL zswap( kk, w( kk, 1 ), ldw, w( kp, 1 ), ldw )
767  END IF
768 *
769  IF( kstep.EQ.1 ) THEN
770 *
771 * 1-by-1 pivot block D(k): column k of W now holds
772 *
773 * W(k) = L(k)*D(k)
774 *
775 * where L(k) is the k-th column of L
776 *
777 * Store L(k) in column k of A
778 *
779  CALL zcopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
780  IF( k.LT.n ) THEN
781  IF( cabs1( a( k, k ) ).GE.sfmin ) THEN
782  r1 = cone / a( k, k )
783  CALL zscal( n-k, r1, a( k+1, k ), 1 )
784  ELSE IF( a( k, k ).NE.czero ) THEN
785  DO 74 ii = k + 1, n
786  a( ii, k ) = a( ii, k ) / a( k, k )
787  74 CONTINUE
788  END IF
789  END IF
790 *
791  ELSE
792 *
793 * 2-by-2 pivot block D(k): columns k and k+1 of W now hold
794 *
795 * ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
796 *
797 * where L(k) and L(k+1) are the k-th and (k+1)-th columns
798 * of L
799 *
800  IF( k.LT.n-1 ) THEN
801 *
802 * Store L(k) and L(k+1) in columns k and k+1 of A
803 *
804  d21 = w( k+1, k )
805  d11 = w( k+1, k+1 ) / d21
806  d22 = w( k, k ) / d21
807  t = cone / ( d11*d22-cone )
808  DO 80 j = k + 2, n
809  a( j, k ) = t*( ( d11*w( j, k )-w( j, k+1 ) ) /
810  \$ d21 )
811  a( j, k+1 ) = t*( ( d22*w( j, k+1 )-w( j, k ) ) /
812  \$ d21 )
813  80 CONTINUE
814  END IF
815 *
816 * Copy D(k) to A
817 *
818  a( k, k ) = w( k, k )
819  a( k+1, k ) = w( k+1, k )
820  a( k+1, k+1 ) = w( k+1, k+1 )
821  END IF
822  END IF
823 *
824 * Store details of the interchanges in IPIV
825 *
826  IF( kstep.EQ.1 ) THEN
827  ipiv( k ) = kp
828  ELSE
829  ipiv( k ) = -p
830  ipiv( k+1 ) = -kp
831  END IF
832 *
833 * Increase K and return to the start of the main loop
834 *
835  k = k + kstep
836  GO TO 70
837 *
838  90 CONTINUE
839 *
840 * Update the lower triangle of A22 (= A(k:n,k:n)) as
841 *
842 * A22 := A22 - L21*D*L21**T = A22 - L21*W**T
843 *
844 * computing blocks of NB columns at a time
845 *
846  DO 110 j = k, n, nb
847  jb = min( nb, n-j+1 )
848 *
849 * Update the lower triangle of the diagonal block
850 *
851  DO 100 jj = j, j + jb - 1
852  CALL zgemv( 'No transpose', j+jb-jj, k-1, -cone,
853  \$ a( jj, 1 ), lda, w( jj, 1 ), ldw, cone,
854  \$ a( jj, jj ), 1 )
855  100 CONTINUE
856 *
857 * Update the rectangular subdiagonal block
858 *
859  IF( j+jb.LE.n )
860  \$ CALL zgemm( 'No transpose', 'Transpose', n-j-jb+1, jb,
861  \$ k-1, -cone, a( j+jb, 1 ), lda, w( j, 1 ), ldw,
862  \$ cone, a( j+jb, j ), lda )
863  110 CONTINUE
864 *
865 * Put L21 in standard form by partially undoing the interchanges
866 * in columns 1:k-1
867 *
868  j = k - 1
869  120 CONTINUE
870 *
871  kstep = 1
872  jp1 = 1
873  jj = j
874  jp2 = ipiv( j )
875  IF( jp2.LT.0 ) THEN
876  jp2 = -jp2
877  j = j - 1
878  jp1 = -ipiv( j )
879  kstep = 2
880  END IF
881 *
882  j = j - 1
883  IF( jp2.NE.jj .AND. j.GE.1 )
884  \$ CALL zswap( j, a( jp2, 1 ), lda, a( jj, 1 ), lda )
885  jj = j + 1
886  IF( jp1.NE.jj .AND. kstep.EQ.2 )
887  \$ CALL zswap( j, a( jp1, 1 ), lda, a( jj, 1 ), lda )
888  IF( j.GE.1 )
889  \$ GO TO 120
890 *
891 * Set KB to the number of columns factorized
892 *
893  kb = k - 1
894 *
895  END IF
896  RETURN
897 *
898 * End of ZLASYF_ROOK
899 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:83
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
Definition: zgemv.f:160
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:83
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189
integer function izamax(N, ZX, INCX)
IZAMAX
Definition: izamax.f:73
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:80
Here is the call graph for this function:
Here is the caller graph for this function: