142 SUBROUTINE cpstrf( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO )
151 INTEGER info, lda, n, rank
164 parameter( one = 1.0e+0, zero = 0.0e+0 )
166 parameter( cone = ( 1.0e+0, 0.0e+0 ) )
170 REAL ajj, sstop, stemp
171 INTEGER i, itemp,
j, jb, k, nb, pvt
185 INTRINSIC conjg, max, min,
REAL, sqrt, maxloc
192 upper =
lsame( uplo,
'U' )
193 IF( .NOT.upper .AND. .NOT.
lsame( uplo,
'L' ) )
THEN
195 ELSE IF( n.LT.0 )
THEN
197 ELSE IF( lda.LT.max( 1, n ) )
THEN
201 CALL
xerbla(
'CPSTRF', -info )
212 nb =
ilaenv( 1,
'CPOTRF', uplo, n, -1, -1, -1 )
213 IF( nb.LE.1 .OR. nb.GE.n )
THEN
217 CALL
cpstf2( uplo, n, a( 1, 1 ), lda, piv, rank, tol, work,
232 work( i ) =
REAL( A( I, I ) )
234 pvt = maxloc( work( 1:n ), 1 )
235 ajj =
REAL( A( PVT, PVT ) )
236 IF( ajj.EQ.zero.OR.
sisnan( ajj ) )
THEN
244 IF( tol.LT.zero )
THEN
245 sstop = n *
slamch(
'Epsilon' ) * ajj
259 jb = min( nb, n-k+1 )
268 DO 150
j = k, k + jb - 1
277 work( i ) = work( i ) +
278 $
REAL( CONJG( A( J-1, I ) )*
281 work( n+i ) =
REAL( A( I, I ) ) - work( i )
286 itemp = maxloc( work( (n+
j):(2*n) ), 1 )
289 IF( ajj.LE.sstop.OR.
sisnan( ajj ) )
THEN
299 a( pvt, pvt ) = a(
j,
j )
300 CALL
cswap(
j-1, a( 1,
j ), 1, a( 1, pvt ), 1 )
302 $ CALL
cswap( n-pvt, a(
j, pvt+1 ), lda,
303 $ a( pvt, pvt+1 ), lda )
304 DO 140 i =
j + 1, pvt - 1
305 ctemp = conjg( a(
j, i ) )
306 a(
j, i ) = conjg( a( i, pvt ) )
309 a(
j, pvt ) = conjg( a(
j, pvt ) )
314 work(
j ) = work( pvt )
317 piv( pvt ) = piv(
j )
328 CALL
cgemv(
'Trans',
j-k, n-
j, -cone, a( k,
j+1 ),
329 $ lda, a( k,
j ), 1, cone, a(
j,
j+1 ),
332 CALL
csscal( n-
j, one / ajj, a(
j,
j+1 ), lda )
340 CALL
cherk(
'Upper',
'Conj Trans', n-
j+1, jb, -one,
341 $ a( k,
j ), lda, one, a(
j,
j ), lda )
354 jb = min( nb, n-k+1 )
363 DO 200
j = k, k + jb - 1
372 work( i ) = work( i ) +
373 $
REAL( CONJG( A( I, J-1 ) )*
376 work( n+i ) =
REAL( A( I, I ) ) - work( i )
381 itemp = maxloc( work( (n+
j):(2*n) ), 1 )
384 IF( ajj.LE.sstop.OR.
sisnan( ajj ) )
THEN
394 a( pvt, pvt ) = a(
j,
j )
395 CALL
cswap(
j-1, a(
j, 1 ), lda, a( pvt, 1 ), lda )
397 $ CALL
cswap( n-pvt, a( pvt+1,
j ), 1,
398 $ a( pvt+1, pvt ), 1 )
399 DO 190 i =
j + 1, pvt - 1
400 ctemp = conjg( a( i,
j ) )
401 a( i,
j ) = conjg( a( pvt, i ) )
404 a( pvt,
j ) = conjg( a( pvt,
j ) )
409 work(
j ) = work( pvt )
412 piv( pvt ) = piv(
j )
423 CALL
cgemv(
'No Trans', n-
j,
j-k, -cone,
424 $ a(
j+1, k ), lda, a(
j, k ), lda, cone,
427 CALL
csscal( n-
j, one / ajj, a(
j+1,
j ), 1 )
435 CALL
cherk(
'Lower',
'No Trans', n-
j+1, jb, -one,
436 $ a(
j, k ), lda, one, a(
j,
j ), lda )