LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
csytrf.f
Go to the documentation of this file.
00001 *> \brief \b CSYTRF
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download CSYTRF + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/csytrf.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/csytrf.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/csytrf.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE CSYTRF( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO )
00022 * 
00023 *       .. Scalar Arguments ..
00024 *       CHARACTER          UPLO
00025 *       INTEGER            INFO, LDA, LWORK, N
00026 *       ..
00027 *       .. Array Arguments ..
00028 *       INTEGER            IPIV( * )
00029 *       COMPLEX            A( LDA, * ), WORK( * )
00030 *       ..
00031 *  
00032 *
00033 *> \par Purpose:
00034 *  =============
00035 *>
00036 *> \verbatim
00037 *>
00038 *> CSYTRF computes the factorization of a complex symmetric matrix A
00039 *> using the Bunch-Kaufman diagonal pivoting method.  The form of the
00040 *> factorization is
00041 *>
00042 *>    A = U*D*U**T  or  A = L*D*L**T
00043 *>
00044 *> where U (or L) is a product of permutation and unit upper (lower)
00045 *> triangular matrices, and D is symmetric and block diagonal with
00046 *> with 1-by-1 and 2-by-2 diagonal blocks.
00047 *>
00048 *> This is the blocked version of the algorithm, calling Level 3 BLAS.
00049 *> \endverbatim
00050 *
00051 *  Arguments:
00052 *  ==========
00053 *
00054 *> \param[in] UPLO
00055 *> \verbatim
00056 *>          UPLO is CHARACTER*1
00057 *>          = 'U':  Upper triangle of A is stored;
00058 *>          = 'L':  Lower triangle of A is stored.
00059 *> \endverbatim
00060 *>
00061 *> \param[in] N
00062 *> \verbatim
00063 *>          N is INTEGER
00064 *>          The order of the matrix A.  N >= 0.
00065 *> \endverbatim
00066 *>
00067 *> \param[in,out] A
00068 *> \verbatim
00069 *>          A is COMPLEX array, dimension (LDA,N)
00070 *>          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
00071 *>          N-by-N upper triangular part of A contains the upper
00072 *>          triangular part of the matrix A, and the strictly lower
00073 *>          triangular part of A is not referenced.  If UPLO = 'L', the
00074 *>          leading N-by-N lower triangular part of A contains the lower
00075 *>          triangular part of the matrix A, and the strictly upper
00076 *>          triangular part of A is not referenced.
00077 *>
00078 *>          On exit, the block diagonal matrix D and the multipliers used
00079 *>          to obtain the factor U or L (see below for further details).
00080 *> \endverbatim
00081 *>
00082 *> \param[in] LDA
00083 *> \verbatim
00084 *>          LDA is INTEGER
00085 *>          The leading dimension of the array A.  LDA >= max(1,N).
00086 *> \endverbatim
00087 *>
00088 *> \param[out] IPIV
00089 *> \verbatim
00090 *>          IPIV is INTEGER array, dimension (N)
00091 *>          Details of the interchanges and the block structure of D.
00092 *>          If IPIV(k) > 0, then rows and columns k and IPIV(k) were
00093 *>          interchanged and D(k,k) is a 1-by-1 diagonal block.
00094 *>          If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
00095 *>          columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
00096 *>          is a 2-by-2 diagonal block.  If UPLO = 'L' and IPIV(k) =
00097 *>          IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
00098 *>          interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
00099 *> \endverbatim
00100 *>
00101 *> \param[out] WORK
00102 *> \verbatim
00103 *>          WORK is COMPLEX array, dimension (MAX(1,LWORK))
00104 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00105 *> \endverbatim
00106 *>
00107 *> \param[in] LWORK
00108 *> \verbatim
00109 *>          LWORK is INTEGER
00110 *>          The length of WORK.  LWORK >=1.  For best performance
00111 *>          LWORK >= N*NB, where NB is the block size returned by ILAENV.
00112 *>
00113 *>          If LWORK = -1, then a workspace query is assumed; the routine
00114 *>          only calculates the optimal size of the WORK array, returns
00115 *>          this value as the first entry of the WORK array, and no error
00116 *>          message related to LWORK is issued by XERBLA.
00117 *> \endverbatim
00118 *>
00119 *> \param[out] INFO
00120 *> \verbatim
00121 *>          INFO is INTEGER
00122 *>          = 0:  successful exit
00123 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
00124 *>          > 0:  if INFO = i, D(i,i) is exactly zero.  The factorization
00125 *>                has been completed, but the block diagonal matrix D is
00126 *>                exactly singular, and division by zero will occur if it
00127 *>                is used to solve a system of equations.
00128 *> \endverbatim
00129 *
00130 *  Authors:
00131 *  ========
00132 *
00133 *> \author Univ. of Tennessee 
00134 *> \author Univ. of California Berkeley 
00135 *> \author Univ. of Colorado Denver 
00136 *> \author NAG Ltd. 
00137 *
00138 *> \date November 2011
00139 *
00140 *> \ingroup complexSYcomputational
00141 *
00142 *> \par Further Details:
00143 *  =====================
00144 *>
00145 *> \verbatim
00146 *>
00147 *>  If UPLO = 'U', then A = U*D*U**T, where
00148 *>     U = P(n)*U(n)* ... *P(k)U(k)* ...,
00149 *>  i.e., U is a product of terms P(k)*U(k), where k decreases from n to
00150 *>  1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
00151 *>  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
00152 *>  defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
00153 *>  that if the diagonal block D(k) is of order s (s = 1 or 2), then
00154 *>
00155 *>             (   I    v    0   )   k-s
00156 *>     U(k) =  (   0    I    0   )   s
00157 *>             (   0    0    I   )   n-k
00158 *>                k-s   s   n-k
00159 *>
00160 *>  If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
00161 *>  If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
00162 *>  and A(k,k), and v overwrites A(1:k-2,k-1:k).
00163 *>
00164 *>  If UPLO = 'L', then A = L*D*L**T, where
00165 *>     L = P(1)*L(1)* ... *P(k)*L(k)* ...,
00166 *>  i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
00167 *>  n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
00168 *>  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
00169 *>  defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
00170 *>  that if the diagonal block D(k) is of order s (s = 1 or 2), then
00171 *>
00172 *>             (   I    0     0   )  k-1
00173 *>     L(k) =  (   0    I     0   )  s
00174 *>             (   0    v     I   )  n-k-s+1
00175 *>                k-1   s  n-k-s+1
00176 *>
00177 *>  If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
00178 *>  If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
00179 *>  and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
00180 *> \endverbatim
00181 *>
00182 *  =====================================================================
00183       SUBROUTINE CSYTRF( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO )
00184 *
00185 *  -- LAPACK computational routine (version 3.4.0) --
00186 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00187 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00188 *     November 2011
00189 *
00190 *     .. Scalar Arguments ..
00191       CHARACTER          UPLO
00192       INTEGER            INFO, LDA, LWORK, N
00193 *     ..
00194 *     .. Array Arguments ..
00195       INTEGER            IPIV( * )
00196       COMPLEX            A( LDA, * ), WORK( * )
00197 *     ..
00198 *
00199 *  =====================================================================
00200 *
00201 *     .. Local Scalars ..
00202       LOGICAL            LQUERY, UPPER
00203       INTEGER            IINFO, IWS, J, K, KB, LDWORK, LWKOPT, NB, NBMIN
00204 *     ..
00205 *     .. External Functions ..
00206       LOGICAL            LSAME
00207       INTEGER            ILAENV
00208       EXTERNAL           LSAME, ILAENV
00209 *     ..
00210 *     .. External Subroutines ..
00211       EXTERNAL           CLASYF, CSYTF2, XERBLA
00212 *     ..
00213 *     .. Intrinsic Functions ..
00214       INTRINSIC          MAX
00215 *     ..
00216 *     .. Executable Statements ..
00217 *
00218 *     Test the input parameters.
00219 *
00220       INFO = 0
00221       UPPER = LSAME( UPLO, 'U' )
00222       LQUERY = ( LWORK.EQ.-1 )
00223       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00224          INFO = -1
00225       ELSE IF( N.LT.0 ) THEN
00226          INFO = -2
00227       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00228          INFO = -4
00229       ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN
00230          INFO = -7
00231       END IF
00232 *
00233       IF( INFO.EQ.0 ) THEN
00234 *
00235 *        Determine the block size
00236 *
00237          NB = ILAENV( 1, 'CSYTRF', UPLO, N, -1, -1, -1 )
00238          LWKOPT = N*NB
00239          WORK( 1 ) = LWKOPT
00240       END IF
00241 *
00242       IF( INFO.NE.0 ) THEN
00243          CALL XERBLA( 'CSYTRF', -INFO )
00244          RETURN
00245       ELSE IF( LQUERY ) THEN
00246          RETURN
00247       END IF
00248 *
00249       NBMIN = 2
00250       LDWORK = N
00251       IF( NB.GT.1 .AND. NB.LT.N ) THEN
00252          IWS = LDWORK*NB
00253          IF( LWORK.LT.IWS ) THEN
00254             NB = MAX( LWORK / LDWORK, 1 )
00255             NBMIN = MAX( 2, ILAENV( 2, 'CSYTRF', UPLO, N, -1, -1, -1 ) )
00256          END IF
00257       ELSE
00258          IWS = 1
00259       END IF
00260       IF( NB.LT.NBMIN )
00261      $   NB = N
00262 *
00263       IF( UPPER ) THEN
00264 *
00265 *        Factorize A as U*D*U**T using the upper triangle of A
00266 *
00267 *        K is the main loop index, decreasing from N to 1 in steps of
00268 *        KB, where KB is the number of columns factorized by CLASYF;
00269 *        KB is either NB or NB-1, or K for the last block
00270 *
00271          K = N
00272    10    CONTINUE
00273 *
00274 *        If K < 1, exit from loop
00275 *
00276          IF( K.LT.1 )
00277      $      GO TO 40
00278 *
00279          IF( K.GT.NB ) THEN
00280 *
00281 *           Factorize columns k-kb+1:k of A and use blocked code to
00282 *           update columns 1:k-kb
00283 *
00284             CALL CLASYF( UPLO, K, NB, KB, A, LDA, IPIV, WORK, N, IINFO )
00285          ELSE
00286 *
00287 *           Use unblocked code to factorize columns 1:k of A
00288 *
00289             CALL CSYTF2( UPLO, K, A, LDA, IPIV, IINFO )
00290             KB = K
00291          END IF
00292 *
00293 *        Set INFO on the first occurrence of a zero pivot
00294 *
00295          IF( INFO.EQ.0 .AND. IINFO.GT.0 )
00296      $      INFO = IINFO
00297 *
00298 *        Decrease K and return to the start of the main loop
00299 *
00300          K = K - KB
00301          GO TO 10
00302 *
00303       ELSE
00304 *
00305 *        Factorize A as L*D*L**T using the lower triangle of A
00306 *
00307 *        K is the main loop index, increasing from 1 to N in steps of
00308 *        KB, where KB is the number of columns factorized by CLASYF;
00309 *        KB is either NB or NB-1, or N-K+1 for the last block
00310 *
00311          K = 1
00312    20    CONTINUE
00313 *
00314 *        If K > N, exit from loop
00315 *
00316          IF( K.GT.N )
00317      $      GO TO 40
00318 *
00319          IF( K.LE.N-NB ) THEN
00320 *
00321 *           Factorize columns k:k+kb-1 of A and use blocked code to
00322 *           update columns k+kb:n
00323 *
00324             CALL CLASYF( UPLO, N-K+1, NB, KB, A( K, K ), LDA, IPIV( K ),
00325      $                   WORK, N, IINFO )
00326          ELSE
00327 *
00328 *           Use unblocked code to factorize columns k:n of A
00329 *
00330             CALL CSYTF2( UPLO, N-K+1, A( K, K ), LDA, IPIV( K ), IINFO )
00331             KB = N - K + 1
00332          END IF
00333 *
00334 *        Set INFO on the first occurrence of a zero pivot
00335 *
00336          IF( INFO.EQ.0 .AND. IINFO.GT.0 )
00337      $      INFO = IINFO + K - 1
00338 *
00339 *        Adjust IPIV
00340 *
00341          DO 30 J = K, K + KB - 1
00342             IF( IPIV( J ).GT.0 ) THEN
00343                IPIV( J ) = IPIV( J ) + K - 1
00344             ELSE
00345                IPIV( J ) = IPIV( J ) - K + 1
00346             END IF
00347    30    CONTINUE
00348 *
00349 *        Increase K and return to the start of the main loop
00350 *
00351          K = K + KB
00352          GO TO 20
00353 *
00354       END IF
00355 *
00356    40 CONTINUE
00357       WORK( 1 ) = LWKOPT
00358       RETURN
00359 *
00360 *     End of CSYTRF
00361 *
00362       END
 All Files Functions