LAPACK 3.3.0

zgttrf.f

Go to the documentation of this file.
00001       SUBROUTINE ZGTTRF( N, DL, D, DU, DU2, IPIV, INFO )
00002 *
00003 *  -- LAPACK routine (version 3.2) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       INTEGER            INFO, N
00010 *     ..
00011 *     .. Array Arguments ..
00012       INTEGER            IPIV( * )
00013       COMPLEX*16         D( * ), DL( * ), DU( * ), DU2( * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  ZGTTRF computes an LU factorization of a complex tridiagonal matrix A
00020 *  using elimination with partial pivoting and row interchanges.
00021 *
00022 *  The factorization has the form
00023 *     A = L * U
00024 *  where L is a product of permutation and unit lower bidiagonal
00025 *  matrices and U is upper triangular with nonzeros in only the main
00026 *  diagonal and first two superdiagonals.
00027 *
00028 *  Arguments
00029 *  =========
00030 *
00031 *  N       (input) INTEGER
00032 *          The order of the matrix A.
00033 *
00034 *  DL      (input/output) COMPLEX*16 array, dimension (N-1)
00035 *          On entry, DL must contain the (n-1) sub-diagonal elements of
00036 *          A.
00037 *
00038 *          On exit, DL is overwritten by the (n-1) multipliers that
00039 *          define the matrix L from the LU factorization of A.
00040 *
00041 *  D       (input/output) COMPLEX*16 array, dimension (N)
00042 *          On entry, D must contain the diagonal elements of A.
00043 *
00044 *          On exit, D is overwritten by the n diagonal elements of the
00045 *          upper triangular matrix U from the LU factorization of A.
00046 *
00047 *  DU      (input/output) COMPLEX*16 array, dimension (N-1)
00048 *          On entry, DU must contain the (n-1) super-diagonal elements
00049 *          of A.
00050 *
00051 *          On exit, DU is overwritten by the (n-1) elements of the first
00052 *          super-diagonal of U.
00053 *
00054 *  DU2     (output) COMPLEX*16 array, dimension (N-2)
00055 *          On exit, DU2 is overwritten by the (n-2) elements of the
00056 *          second super-diagonal of U.
00057 *
00058 *  IPIV    (output) INTEGER array, dimension (N)
00059 *          The pivot indices; for 1 <= i <= n, row i of the matrix was
00060 *          interchanged with row IPIV(i).  IPIV(i) will always be either
00061 *          i or i+1; IPIV(i) = i indicates a row interchange was not
00062 *          required.
00063 *
00064 *  INFO    (output) INTEGER
00065 *          = 0:  successful exit
00066 *          < 0:  if INFO = -k, the k-th argument had an illegal value
00067 *          > 0:  if INFO = k, U(k,k) is exactly zero. The factorization
00068 *                has been completed, but the factor U is exactly
00069 *                singular, and division by zero will occur if it is used
00070 *                to solve a system of equations.
00071 *
00072 *  =====================================================================
00073 *
00074 *     .. Parameters ..
00075       DOUBLE PRECISION   ZERO
00076       PARAMETER          ( ZERO = 0.0D+0 )
00077 *     ..
00078 *     .. Local Scalars ..
00079       INTEGER            I
00080       COMPLEX*16         FACT, TEMP, ZDUM
00081 *     ..
00082 *     .. External Subroutines ..
00083       EXTERNAL           XERBLA
00084 *     ..
00085 *     .. Intrinsic Functions ..
00086       INTRINSIC          ABS, DBLE, DIMAG
00087 *     ..
00088 *     .. Statement Functions ..
00089       DOUBLE PRECISION   CABS1
00090 *     ..
00091 *     .. Statement Function definitions ..
00092       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
00093 *     ..
00094 *     .. Executable Statements ..
00095 *
00096       INFO = 0
00097       IF( N.LT.0 ) THEN
00098          INFO = -1
00099          CALL XERBLA( 'ZGTTRF', -INFO )
00100          RETURN
00101       END IF
00102 *
00103 *     Quick return if possible
00104 *
00105       IF( N.EQ.0 )
00106      $   RETURN
00107 *
00108 *     Initialize IPIV(i) = i and DU2(i) = 0
00109 *
00110       DO 10 I = 1, N
00111          IPIV( I ) = I
00112    10 CONTINUE
00113       DO 20 I = 1, N - 2
00114          DU2( I ) = ZERO
00115    20 CONTINUE
00116 *
00117       DO 30 I = 1, N - 2
00118          IF( CABS1( D( I ) ).GE.CABS1( DL( I ) ) ) THEN
00119 *
00120 *           No row interchange required, eliminate DL(I)
00121 *
00122             IF( CABS1( D( I ) ).NE.ZERO ) THEN
00123                FACT = DL( I ) / D( I )
00124                DL( I ) = FACT
00125                D( I+1 ) = D( I+1 ) - FACT*DU( I )
00126             END IF
00127          ELSE
00128 *
00129 *           Interchange rows I and I+1, eliminate DL(I)
00130 *
00131             FACT = D( I ) / DL( I )
00132             D( I ) = DL( I )
00133             DL( I ) = FACT
00134             TEMP = DU( I )
00135             DU( I ) = D( I+1 )
00136             D( I+1 ) = TEMP - FACT*D( I+1 )
00137             DU2( I ) = DU( I+1 )
00138             DU( I+1 ) = -FACT*DU( I+1 )
00139             IPIV( I ) = I + 1
00140          END IF
00141    30 CONTINUE
00142       IF( N.GT.1 ) THEN
00143          I = N - 1
00144          IF( CABS1( D( I ) ).GE.CABS1( DL( I ) ) ) THEN
00145             IF( CABS1( D( I ) ).NE.ZERO ) THEN
00146                FACT = DL( I ) / D( I )
00147                DL( I ) = FACT
00148                D( I+1 ) = D( I+1 ) - FACT*DU( I )
00149             END IF
00150          ELSE
00151             FACT = D( I ) / DL( I )
00152             D( I ) = DL( I )
00153             DL( I ) = FACT
00154             TEMP = DU( I )
00155             DU( I ) = D( I+1 )
00156             D( I+1 ) = TEMP - FACT*D( I+1 )
00157             IPIV( I ) = I + 1
00158          END IF
00159       END IF
00160 *
00161 *     Check for a zero on the diagonal of U.
00162 *
00163       DO 40 I = 1, N
00164          IF( CABS1( D( I ) ).EQ.ZERO ) THEN
00165             INFO = I
00166             GO TO 50
00167          END IF
00168    40 CONTINUE
00169    50 CONTINUE
00170 *
00171       RETURN
00172 *
00173 *     End of ZGTTRF
00174 *
00175       END
 All Files Functions