LAPACK  3.10.0 LAPACK: Linear Algebra PACKage

◆ clasyf_rook()

 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.

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.```
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 182 of file clasyf_rook.f.

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