LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine clasyf_rook ( character  UPLO,
integer  N,
integer  NB,
integer  KB,
complex, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
complex, dimension( ldw, * )  W,
integer  LDW,
integer  INFO 
)

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

Download CLASYF_ROOK + dependencies [TGZ] [ZIP] [TXT]

Purpose: 
 CLASYF_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.

 CLASYF_ROOK is an auxiliary routine called by CSYTRF_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 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 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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
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 clasyf_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 a( lda, * ), w( ldw, * )
199 * ..
200 *
201 * =====================================================================
202 *
203 * .. Parameters ..
204  REAL zero, one
205  parameter ( zero = 0.0e+0, one = 1.0e+0 )
206  REAL eight, sevten
207  parameter ( eight = 8.0e+0, sevten = 17.0e+0 )
208  COMPLEX cone, czero
209  parameter ( cone = ( 1.0e+0, 0.0e+0 ),
210  $ czero = ( 0.0e+0, 0.0e+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  REAL absakk, alpha, colmax, rowmax, stemp, sfmin
217  COMPLEX d11, d12, d21, d22, r1, t, z
218 * ..
219 * .. External Functions ..
220  LOGICAL lsame
221  INTEGER icamax
222  REAL slamch
223  EXTERNAL lsame, icamax, slamch
224 * ..
225 * .. External Subroutines ..
226  EXTERNAL ccopy, cgemm, cgemv, cscal, cswap
227 * ..
228 * .. Intrinsic Functions ..
229  INTRINSIC abs, max, min, sqrt, aimag, real
230 * ..
231 * .. Statement Functions ..
232  REAL cabs1
233 * ..
234 * .. Statement Function definitions ..
235  cabs1( z ) = abs( REAL( Z ) ) + abs( aimag( 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 = slamch( '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 ccopy( k, a( 1, k ), 1, w( 1, kw ), 1 )
275  IF( k.LT.n )
276  $ CALL cgemv( '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 = icamax( 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 ccopy( 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 ccopy( imax, a( 1, imax ), 1, w( 1, kw-1 ), 1 )
332  CALL ccopy( k-imax, a( imax, imax+1 ), lda,
333  $ w( imax+1, kw-1 ), 1 )
334 *
335  IF( k.LT.n )
336  $ CALL cgemv( '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 + icamax( 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 = icamax( imax-1, w( 1, kw-1 ), 1 )
354  stemp = cabs1( w( itemp, kw-1 ) )
355  IF( stemp.GT.rowmax ) THEN
356  rowmax = stemp
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 ccopy( 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 *
393 * Pivot not found: set params and repeat
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 ccopy( 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 ccopy( k-p, a( p+1, k ), 1, a( p, p+1 ), lda )
424  CALL ccopy( 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 cswap( n-k+1, a( k, k ), lda, a( p, k ), lda )
430  CALL cswap( 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 ccopy( k-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
441  $ lda )
442  CALL ccopy( 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 cswap( n-kk+1, a( kk, kk ), lda, a( kp, kk ), lda )
448  CALL cswap( 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 ccopy( 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 cscal( 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 cgemv( '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 cgemm( '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 cswap( 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 cswap( 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 ccopy( n-k+1, a( k, k ), 1, w( k, k ), 1 )
601  IF( k.GT.1 )
602  $ CALL cgemv( '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 + icamax( 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 ccopy( 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 ccopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1)
658  CALL ccopy( n-imax+1, a( imax, imax ), 1,
659  $ w( imax, k+1 ), 1 )
660  IF( k.GT.1 )
661  $ CALL cgemv( '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 + icamax( 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 + icamax( n-imax, w( imax+1, k+1 ), 1)
678  stemp = cabs1( w( itemp, k+1 ) )
679  IF( stemp.GT.rowmax ) THEN
680  rowmax = stemp
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 ccopy( 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 *
717 * Pivot not found: set params and repeat
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 ccopy( 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 ccopy( p-k, a( k, k ), 1, a( p, k ), lda )
744  CALL ccopy( 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 cswap( k, a( k, 1 ), lda, a( p, 1 ), lda )
750  CALL cswap( 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 ccopy( kp-k-1, a( k+1, kk ), 1, a( kp, k+1 ), lda )
761  CALL ccopy( 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 cswap( kk, a( kk, 1 ), lda, a( kp, 1 ), lda )
766  CALL cswap( 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 ccopy( 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 cscal( 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 cgemv( '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 cgemm( '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 cswap( j, a( jp2, 1 ), lda, a( jj, 1 ), lda )
885  jj = j + 1
886  IF( jp1.NE.jj .AND. kstep.EQ.2 )
887  $ CALL cswap( 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 CLASYF_ROOK
899 *
cscal
subroutine cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:54
cgemv
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
Definition: cgemv.f:160
icamax
integer function icamax(N, CX, INCX)
ICAMAX
Definition: icamax.f:53
ccopy
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:52
cswap
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:52
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
cgemm
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:189
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the call graph for this function:

Here is the caller graph for this function: