LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dsytf2()

subroutine dsytf2 ( character  UPLO,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
integer  INFO 
)

DSYTF2 computes the factorization of a real symmetric indefinite matrix, using the diagonal pivoting method (unblocked algorithm).

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

Purpose:
 DSYTF2 computes the factorization of a real symmetric matrix A using
 the Bunch-Kaufman diagonal pivoting method:

    A = U*D*U**T  or  A = L*D*L**T

 where U (or L) is a product of permutation and unit upper (lower)
 triangular matrices, U**T is the transpose of U, and D is symmetric and
 block diagonal with 1-by-1 and 2-by-2 diagonal blocks.

 This is the unblocked version of the algorithm, calling Level 2 BLAS.
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,out]A
          A is DOUBLE PRECISION 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, the block diagonal matrix D and the multipliers used
          to obtain the factor U or L (see below for further details).
[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':
             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) = IPIV(k-1) < 0, then rows and columns
             k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
             is a 2-by-2 diagonal block.

          If UPLO = 'L':
             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) = IPIV(k+1) < 0, then rows and columns
             k+1 and -IPIV(k) were interchanged and D(k:k+1,k:k+1)
             is a 2-by-2 diagonal block.
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -k, the k-th argument had an illegal value
          > 0: if INFO = k, D(k,k) is exactly zero.  The factorization
               has been completed, but the block diagonal matrix D is
               exactly singular, and division by zero will occur if it
               is used to solve a system of equations.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
  If UPLO = 'U', then A = U*D*U**T, where
     U = P(n)*U(n)* ... *P(k)U(k)* ...,
  i.e., U is a product of terms P(k)*U(k), where k decreases from n to
  1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
  defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
  that if the diagonal block D(k) is of order s (s = 1 or 2), then

             (   I    v    0   )   k-s
     U(k) =  (   0    I    0   )   s
             (   0    0    I   )   n-k
                k-s   s   n-k

  If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
  If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
  and A(k,k), and v overwrites A(1:k-2,k-1:k).

  If UPLO = 'L', then A = L*D*L**T, where
     L = P(1)*L(1)* ... *P(k)*L(k)* ...,
  i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
  n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
  defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
  that if the diagonal block D(k) is of order s (s = 1 or 2), then

             (   I    0     0   )  k-1
     L(k) =  (   0    I     0   )  s
             (   0    v     I   )  n-k-s+1
                k-1   s  n-k-s+1

  If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
  If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
  and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
Contributors:
  09-29-06 - patch from
    Bobby Cheng, MathWorks

    Replace l.204 and l.372
         IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
    by
         IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN

  01-01-96 - Based on modifications by
    J. Lewis, Boeing Computer Services Company
    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
  1-96 - Based on modifications by J. Lewis, Boeing Computer Services
         Company

Definition at line 193 of file dsytf2.f.

194*
195* -- LAPACK computational routine --
196* -- LAPACK is a software package provided by Univ. of Tennessee, --
197* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
198*
199* .. Scalar Arguments ..
200 CHARACTER UPLO
201 INTEGER INFO, LDA, N
202* ..
203* .. Array Arguments ..
204 INTEGER IPIV( * )
205 DOUBLE PRECISION A( LDA, * )
206* ..
207*
208* =====================================================================
209*
210* .. Parameters ..
211 DOUBLE PRECISION ZERO, ONE
212 parameter( zero = 0.0d+0, one = 1.0d+0 )
213 DOUBLE PRECISION EIGHT, SEVTEN
214 parameter( eight = 8.0d+0, sevten = 17.0d+0 )
215* ..
216* .. Local Scalars ..
217 LOGICAL UPPER
218 INTEGER I, IMAX, J, JMAX, K, KK, KP, KSTEP
219 DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D11, D12, D21, D22, R1,
220 $ ROWMAX, T, WK, WKM1, WKP1
221* ..
222* .. External Functions ..
223 LOGICAL LSAME, DISNAN
224 INTEGER IDAMAX
225 EXTERNAL lsame, idamax, disnan
226* ..
227* .. External Subroutines ..
228 EXTERNAL dscal, dswap, dsyr, xerbla
229* ..
230* .. Intrinsic Functions ..
231 INTRINSIC abs, max, sqrt
232* ..
233* .. Executable Statements ..
234*
235* Test the input parameters.
236*
237 info = 0
238 upper = lsame( uplo, 'U' )
239 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
240 info = -1
241 ELSE IF( n.LT.0 ) THEN
242 info = -2
243 ELSE IF( lda.LT.max( 1, n ) ) THEN
244 info = -4
245 END IF
246 IF( info.NE.0 ) THEN
247 CALL xerbla( 'DSYTF2', -info )
248 RETURN
249 END IF
250*
251* Initialize ALPHA for use in choosing pivot block size.
252*
253 alpha = ( one+sqrt( sevten ) ) / eight
254*
255 IF( upper ) THEN
256*
257* Factorize A as U*D*U**T using the upper triangle of A
258*
259* K is the main loop index, decreasing from N to 1 in steps of
260* 1 or 2
261*
262 k = n
263 10 CONTINUE
264*
265* If K < 1, exit from loop
266*
267 IF( k.LT.1 )
268 $ GO TO 70
269 kstep = 1
270*
271* Determine rows and columns to be interchanged and whether
272* a 1-by-1 or 2-by-2 pivot block will be used
273*
274 absakk = abs( a( k, k ) )
275*
276* IMAX is the row-index of the largest off-diagonal element in
277* column K, and COLMAX is its absolute value.
278* Determine both COLMAX and IMAX.
279*
280 IF( k.GT.1 ) THEN
281 imax = idamax( k-1, a( 1, k ), 1 )
282 colmax = abs( a( imax, k ) )
283 ELSE
284 colmax = zero
285 END IF
286*
287 IF( (max( absakk, colmax ).EQ.zero) .OR. disnan(absakk) ) THEN
288*
289* Column K is zero or underflow, or contains a NaN:
290* set INFO and continue
291*
292 IF( info.EQ.0 )
293 $ info = k
294 kp = k
295 ELSE
296 IF( absakk.GE.alpha*colmax ) THEN
297*
298* no interchange, use 1-by-1 pivot block
299*
300 kp = k
301 ELSE
302*
303* JMAX is the column-index of the largest off-diagonal
304* element in row IMAX, and ROWMAX is its absolute value
305*
306 jmax = imax + idamax( k-imax, a( imax, imax+1 ), lda )
307 rowmax = abs( a( imax, jmax ) )
308 IF( imax.GT.1 ) THEN
309 jmax = idamax( imax-1, a( 1, imax ), 1 )
310 rowmax = max( rowmax, abs( a( jmax, imax ) ) )
311 END IF
312*
313 IF( absakk.GE.alpha*colmax*( colmax / rowmax ) ) THEN
314*
315* no interchange, use 1-by-1 pivot block
316*
317 kp = k
318 ELSE IF( abs( a( imax, imax ) ).GE.alpha*rowmax ) THEN
319*
320* interchange rows and columns K and IMAX, use 1-by-1
321* pivot block
322*
323 kp = imax
324 ELSE
325*
326* interchange rows and columns K-1 and IMAX, use 2-by-2
327* pivot block
328*
329 kp = imax
330 kstep = 2
331 END IF
332 END IF
333*
334 kk = k - kstep + 1
335 IF( kp.NE.kk ) THEN
336*
337* Interchange rows and columns KK and KP in the leading
338* submatrix A(1:k,1:k)
339*
340 CALL dswap( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
341 CALL dswap( kk-kp-1, a( kp+1, kk ), 1, a( kp, kp+1 ),
342 $ lda )
343 t = a( kk, kk )
344 a( kk, kk ) = a( kp, kp )
345 a( kp, kp ) = t
346 IF( kstep.EQ.2 ) THEN
347 t = a( k-1, k )
348 a( k-1, k ) = a( kp, k )
349 a( kp, k ) = t
350 END IF
351 END IF
352*
353* Update the leading submatrix
354*
355 IF( kstep.EQ.1 ) THEN
356*
357* 1-by-1 pivot block D(k): column k now holds
358*
359* W(k) = U(k)*D(k)
360*
361* where U(k) is the k-th column of U
362*
363* Perform a rank-1 update of A(1:k-1,1:k-1) as
364*
365* A := A - U(k)*D(k)*U(k)**T = A - W(k)*1/D(k)*W(k)**T
366*
367 r1 = one / a( k, k )
368 CALL dsyr( uplo, k-1, -r1, a( 1, k ), 1, a, lda )
369*
370* Store U(k) in column k
371*
372 CALL dscal( k-1, r1, a( 1, k ), 1 )
373 ELSE
374*
375* 2-by-2 pivot block D(k): columns k and k-1 now hold
376*
377* ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
378*
379* where U(k) and U(k-1) are the k-th and (k-1)-th columns
380* of U
381*
382* Perform a rank-2 update of A(1:k-2,1:k-2) as
383*
384* A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**T
385* = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )**T
386*
387 IF( k.GT.2 ) THEN
388*
389 d12 = a( k-1, k )
390 d22 = a( k-1, k-1 ) / d12
391 d11 = a( k, k ) / d12
392 t = one / ( d11*d22-one )
393 d12 = t / d12
394*
395 DO 30 j = k - 2, 1, -1
396 wkm1 = d12*( d11*a( j, k-1 )-a( j, k ) )
397 wk = d12*( d22*a( j, k )-a( j, k-1 ) )
398 DO 20 i = j, 1, -1
399 a( i, j ) = a( i, j ) - a( i, k )*wk -
400 $ a( i, k-1 )*wkm1
401 20 CONTINUE
402 a( j, k ) = wk
403 a( j, k-1 ) = wkm1
404 30 CONTINUE
405*
406 END IF
407*
408 END IF
409 END IF
410*
411* Store details of the interchanges in IPIV
412*
413 IF( kstep.EQ.1 ) THEN
414 ipiv( k ) = kp
415 ELSE
416 ipiv( k ) = -kp
417 ipiv( k-1 ) = -kp
418 END IF
419*
420* Decrease K and return to the start of the main loop
421*
422 k = k - kstep
423 GO TO 10
424*
425 ELSE
426*
427* Factorize A as L*D*L**T using the lower triangle of A
428*
429* K is the main loop index, increasing from 1 to N in steps of
430* 1 or 2
431*
432 k = 1
433 40 CONTINUE
434*
435* If K > N, exit from loop
436*
437 IF( k.GT.n )
438 $ GO TO 70
439 kstep = 1
440*
441* Determine rows and columns to be interchanged and whether
442* a 1-by-1 or 2-by-2 pivot block will be used
443*
444 absakk = abs( a( k, k ) )
445*
446* IMAX is the row-index of the largest off-diagonal element in
447* column K, and COLMAX is its absolute value.
448* Determine both COLMAX and IMAX.
449*
450 IF( k.LT.n ) THEN
451 imax = k + idamax( n-k, a( k+1, k ), 1 )
452 colmax = abs( a( imax, k ) )
453 ELSE
454 colmax = zero
455 END IF
456*
457 IF( (max( absakk, colmax ).EQ.zero) .OR. disnan(absakk) ) THEN
458*
459* Column K is zero or underflow, or contains a NaN:
460* set INFO and continue
461*
462 IF( info.EQ.0 )
463 $ info = k
464 kp = k
465 ELSE
466 IF( absakk.GE.alpha*colmax ) THEN
467*
468* no interchange, use 1-by-1 pivot block
469*
470 kp = k
471 ELSE
472*
473* JMAX is the column-index of the largest off-diagonal
474* element in row IMAX, and ROWMAX is its absolute value
475*
476 jmax = k - 1 + idamax( imax-k, a( imax, k ), lda )
477 rowmax = abs( a( imax, jmax ) )
478 IF( imax.LT.n ) THEN
479 jmax = imax + idamax( n-imax, a( imax+1, imax ), 1 )
480 rowmax = max( rowmax, abs( a( jmax, imax ) ) )
481 END IF
482*
483 IF( absakk.GE.alpha*colmax*( colmax / rowmax ) ) THEN
484*
485* no interchange, use 1-by-1 pivot block
486*
487 kp = k
488 ELSE IF( abs( a( imax, imax ) ).GE.alpha*rowmax ) THEN
489*
490* interchange rows and columns K and IMAX, use 1-by-1
491* pivot block
492*
493 kp = imax
494 ELSE
495*
496* interchange rows and columns K+1 and IMAX, use 2-by-2
497* pivot block
498*
499 kp = imax
500 kstep = 2
501 END IF
502 END IF
503*
504 kk = k + kstep - 1
505 IF( kp.NE.kk ) THEN
506*
507* Interchange rows and columns KK and KP in the trailing
508* submatrix A(k:n,k:n)
509*
510 IF( kp.LT.n )
511 $ CALL dswap( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ), 1 )
512 CALL dswap( kp-kk-1, a( kk+1, kk ), 1, a( kp, kk+1 ),
513 $ lda )
514 t = a( kk, kk )
515 a( kk, kk ) = a( kp, kp )
516 a( kp, kp ) = t
517 IF( kstep.EQ.2 ) THEN
518 t = a( k+1, k )
519 a( k+1, k ) = a( kp, k )
520 a( kp, k ) = t
521 END IF
522 END IF
523*
524* Update the trailing submatrix
525*
526 IF( kstep.EQ.1 ) THEN
527*
528* 1-by-1 pivot block D(k): column k now holds
529*
530* W(k) = L(k)*D(k)
531*
532* where L(k) is the k-th column of L
533*
534 IF( k.LT.n ) THEN
535*
536* Perform a rank-1 update of A(k+1:n,k+1:n) as
537*
538* A := A - L(k)*D(k)*L(k)**T = A - W(k)*(1/D(k))*W(k)**T
539*
540 d11 = one / a( k, k )
541 CALL dsyr( uplo, n-k, -d11, a( k+1, k ), 1,
542 $ a( k+1, k+1 ), lda )
543*
544* Store L(k) in column K
545*
546 CALL dscal( n-k, d11, a( k+1, k ), 1 )
547 END IF
548 ELSE
549*
550* 2-by-2 pivot block D(k)
551*
552 IF( k.LT.n-1 ) THEN
553*
554* Perform a rank-2 update of A(k+2:n,k+2:n) as
555*
556* A := A - ( (A(k) A(k+1))*D(k)**(-1) ) * (A(k) A(k+1))**T
557*
558* where L(k) and L(k+1) are the k-th and (k+1)-th
559* columns of L
560*
561 d21 = a( k+1, k )
562 d11 = a( k+1, k+1 ) / d21
563 d22 = a( k, k ) / d21
564 t = one / ( d11*d22-one )
565 d21 = t / d21
566*
567 DO 60 j = k + 2, n
568*
569 wk = d21*( d11*a( j, k )-a( j, k+1 ) )
570 wkp1 = d21*( d22*a( j, k+1 )-a( j, k ) )
571*
572 DO 50 i = j, n
573 a( i, j ) = a( i, j ) - a( i, k )*wk -
574 $ a( i, k+1 )*wkp1
575 50 CONTINUE
576*
577 a( j, k ) = wk
578 a( j, k+1 ) = wkp1
579*
580 60 CONTINUE
581 END IF
582 END IF
583 END IF
584*
585* Store details of the interchanges in IPIV
586*
587 IF( kstep.EQ.1 ) THEN
588 ipiv( k ) = kp
589 ELSE
590 ipiv( k ) = -kp
591 ipiv( k+1 ) = -kp
592 END IF
593*
594* Increase K and return to the start of the main loop
595*
596 k = k + kstep
597 GO TO 40
598*
599 END IF
600*
601 70 CONTINUE
602*
603 RETURN
604*
605* End of DSYTF2
606*
logical function disnan(DIN)
DISNAN tests input for NaN.
Definition: disnan.f:59
integer function idamax(N, DX, INCX)
IDAMAX
Definition: idamax.f:71
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.f:82
subroutine dsyr(UPLO, N, ALPHA, X, INCX, A, LDA)
DSYR
Definition: dsyr.f:132
Here is the call graph for this function:
Here is the caller graph for this function: