1*> \brief \b DSYTRF
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DSYTRF + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsytrf.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsytrf.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsytrf.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE DSYTRF( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO )
22*
23*       .. Scalar Arguments ..
24*       CHARACTER          UPLO
25*       INTEGER            INFO, LDA, LWORK, N
26*       ..
27*       .. Array Arguments ..
28*       INTEGER            IPIV( * )
29*       DOUBLE PRECISION   A( LDA, * ), WORK( * )
30*       ..
31*
32*
33*> \par Purpose:
34*  =============
35*>
36*> \verbatim
37*>
38*> DSYTRF computes the factorization of a real symmetric matrix A using
39*> the Bunch-Kaufman diagonal pivoting method.  The form of the
40*> factorization is
41*>
42*>    A = U*D*U**T  or  A = L*D*L**T
43*>
44*> where U (or L) is a product of permutation and unit upper (lower)
45*> triangular matrices, and D is symmetric and block diagonal with
46*> 1-by-1 and 2-by-2 diagonal blocks.
47*>
48*> This is the blocked version of the algorithm, calling Level 3 BLAS.
49*> \endverbatim
50*
51*  Arguments:
52*  ==========
53*
54*> \param[in] UPLO
55*> \verbatim
56*>          UPLO is CHARACTER*1
57*>          = 'U':  Upper triangle of A is stored;
58*>          = 'L':  Lower triangle of A is stored.
59*> \endverbatim
60*>
61*> \param[in] N
62*> \verbatim
63*>          N is INTEGER
64*>          The order of the matrix A.  N >= 0.
65*> \endverbatim
66*>
67*> \param[in,out] A
68*> \verbatim
69*>          A is DOUBLE PRECISION array, dimension (LDA,N)
70*>          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
71*>          N-by-N upper triangular part of A contains the upper
72*>          triangular part of the matrix A, and the strictly lower
73*>          triangular part of A is not referenced.  If UPLO = 'L', the
74*>          leading N-by-N lower triangular part of A contains the lower
75*>          triangular part of the matrix A, and the strictly upper
76*>          triangular part of A is not referenced.
77*>
78*>          On exit, the block diagonal matrix D and the multipliers used
79*>          to obtain the factor U or L (see below for further details).
80*> \endverbatim
81*>
82*> \param[in] LDA
83*> \verbatim
84*>          LDA is INTEGER
85*>          The leading dimension of the array A.  LDA >= max(1,N).
86*> \endverbatim
87*>
88*> \param[out] IPIV
89*> \verbatim
90*>          IPIV is INTEGER array, dimension (N)
91*>          Details of the interchanges and the block structure of D.
92*>          If IPIV(k) > 0, then rows and columns k and IPIV(k) were
93*>          interchanged and D(k,k) is a 1-by-1 diagonal block.
94*>          If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
95*>          columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
96*>          is a 2-by-2 diagonal block.  If UPLO = 'L' and IPIV(k) =
97*>          IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
98*>          interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
99*> \endverbatim
100*>
101*> \param[out] WORK
102*> \verbatim
103*>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
104*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
105*> \endverbatim
106*>
107*> \param[in] LWORK
108*> \verbatim
109*>          LWORK is INTEGER
110*>          The length of WORK.  LWORK >=1.  For best performance
111*>          LWORK >= N*NB, where NB is the block size returned by ILAENV.
112*>
113*>          If LWORK = -1, then a workspace query is assumed; the routine
114*>          only calculates the optimal size of the WORK array, returns
115*>          this value as the first entry of the WORK array, and no error
116*>          message related to LWORK is issued by XERBLA.
117*> \endverbatim
118*>
119*> \param[out] INFO
120*> \verbatim
121*>          INFO is INTEGER
122*>          = 0:  successful exit
123*>          < 0:  if INFO = -i, the i-th argument had an illegal value
124*>          > 0:  if INFO = i, D(i,i) is exactly zero.  The factorization
125*>                has been completed, but the block diagonal matrix D is
126*>                exactly singular, and division by zero will occur if it
127*>                is used to solve a system of equations.
128*> \endverbatim
129*
130*  Authors:
131*  ========
132*
133*> \author Univ. of Tennessee
134*> \author Univ. of California Berkeley
135*> \author Univ. of Colorado Denver
136*> \author NAG Ltd.
137*
138*> \date November 2011
139*
140*> \ingroup doubleSYcomputational
141*
142*> \par Further Details:
143*  =====================
144*>
145*> \verbatim
146*>
147*>  If UPLO = 'U', then A = U*D*U**T, where
148*>     U = P(n)*U(n)* ... *P(k)U(k)* ...,
149*>  i.e., U is a product of terms P(k)*U(k), where k decreases from n to
150*>  1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
151*>  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
152*>  defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
153*>  that if the diagonal block D(k) is of order s (s = 1 or 2), then
154*>
155*>             (   I    v    0   )   k-s
156*>     U(k) =  (   0    I    0   )   s
157*>             (   0    0    I   )   n-k
158*>                k-s   s   n-k
159*>
160*>  If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
161*>  If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
162*>  and A(k,k), and v overwrites A(1:k-2,k-1:k).
163*>
164*>  If UPLO = 'L', then A = L*D*L**T, where
165*>     L = P(1)*L(1)* ... *P(k)*L(k)* ...,
166*>  i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
167*>  n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
168*>  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
169*>  defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
170*>  that if the diagonal block D(k) is of order s (s = 1 or 2), then
171*>
172*>             (   I    0     0   )  k-1
173*>     L(k) =  (   0    I     0   )  s
174*>             (   0    v     I   )  n-k-s+1
175*>                k-1   s  n-k-s+1
176*>
177*>  If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
178*>  If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
179*>  and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
180*> \endverbatim
181*>
182*  =====================================================================
183      SUBROUTINE DSYTRF( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO )
184*
185*  -- LAPACK computational routine (version 3.4.0) --
186*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
187*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
188*     November 2011
189*
190*     .. Scalar Arguments ..
191      CHARACTER          UPLO
192      INTEGER            INFO, LDA, LWORK, N
193*     ..
194*     .. Array Arguments ..
195      INTEGER            IPIV( * )
196      DOUBLE PRECISION   A( LDA, * ), WORK( * )
197*     ..
198*
199*  =====================================================================
200*
201*     .. Local Scalars ..
202      LOGICAL            LQUERY, UPPER
203      INTEGER            IINFO, IWS, J, K, KB, LDWORK, LWKOPT, NB, NBMIN
204*     ..
205*     .. External Functions ..
206      LOGICAL            LSAME
207      INTEGER            ILAENV
208      EXTERNAL           LSAME, ILAENV
209*     ..
210*     .. External Subroutines ..
211      EXTERNAL           DLASYF, DSYTF2, XERBLA
212*     ..
213*     .. Intrinsic Functions ..
214      INTRINSIC          MAX
215*     ..
216*     .. Executable Statements ..
217*
218*     Test the input parameters.
219*
220      INFO = 0
221      UPPER = LSAME( UPLO, 'U' )
222      LQUERY = ( LWORK.EQ.-1 )
223      IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
224         INFO = -1
225      ELSE IF( N.LT.0 ) THEN
226         INFO = -2
227      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
228         INFO = -4
229      ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN
230         INFO = -7
231      END IF
232*
233      IF( INFO.EQ.0 ) THEN
234*
235*        Determine the block size
236*
237         NB = ILAENV( 1, 'DSYTRF', UPLO, N, -1, -1, -1 )
238         LWKOPT = N*NB
239         WORK( 1 ) = LWKOPT
240      END IF
241*
242      IF( INFO.NE.0 ) THEN
243         CALL XERBLA( 'DSYTRF', -INFO )
244         RETURN
245      ELSE IF( LQUERY ) THEN
246         RETURN
247      END IF
248*
249      NBMIN = 2
250      LDWORK = N
251      IF( NB.GT.1 .AND. NB.LT.N ) THEN
252         IWS = LDWORK*NB
253         IF( LWORK.LT.IWS ) THEN
254            NB = MAX( LWORK / LDWORK, 1 )
255            NBMIN = MAX( 2, ILAENV( 2, 'DSYTRF', UPLO, N, -1, -1, -1 ) )
256         END IF
257      ELSE
258         IWS = 1
259      END IF
260      IF( NB.LT.NBMIN )
261     $   NB = N
262*
263      IF( UPPER ) THEN
264*
265*        Factorize A as U*D*U**T using the upper triangle of A
266*
267*        K is the main loop index, decreasing from N to 1 in steps of
268*        KB, where KB is the number of columns factorized by DLASYF;
269*        KB is either NB or NB-1, or K for the last block
270*
271         K = N
272   10    CONTINUE
273*
274*        If K < 1, exit from loop
275*
276         IF( K.LT.1 )
277     $      GO TO 40
278*
279         IF( K.GT.NB ) THEN
280*
281*           Factorize columns k-kb+1:k of A and use blocked code to
282*           update columns 1:k-kb
283*
284            CALL DLASYF( UPLO, K, NB, KB, A, LDA, IPIV, WORK, LDWORK,
285     $                   IINFO )
286         ELSE
287*
288*           Use unblocked code to factorize columns 1:k of A
289*
290            CALL DSYTF2( UPLO, K, A, LDA, IPIV, IINFO )
291            KB = K
292         END IF
293*
294*        Set INFO on the first occurrence of a zero pivot
295*
296         IF( INFO.EQ.0 .AND. IINFO.GT.0 )
297     $      INFO = IINFO
298*
299*        Decrease K and return to the start of the main loop
300*
301         K = K - KB
302         GO TO 10
303*
304      ELSE
305*
306*        Factorize A as L*D*L**T using the lower triangle of A
307*
308*        K is the main loop index, increasing from 1 to N in steps of
309*        KB, where KB is the number of columns factorized by DLASYF;
310*        KB is either NB or NB-1, or N-K+1 for the last block
311*
312         K = 1
313   20    CONTINUE
314*
315*        If K > N, exit from loop
316*
317         IF( K.GT.N )
318     $      GO TO 40
319*
320         IF( K.LE.N-NB ) THEN
321*
322*           Factorize columns k:k+kb-1 of A and use blocked code to
323*           update columns k+kb:n
324*
325            CALL DLASYF( UPLO, N-K+1, NB, KB, A( K, K ), LDA, IPIV( K ),
326     $                   WORK, LDWORK, IINFO )
327         ELSE
328*
329*           Use unblocked code to factorize columns k:n of A
330*
331            CALL DSYTF2( UPLO, N-K+1, A( K, K ), LDA, IPIV( K ), IINFO )
332            KB = N - K + 1
333         END IF
334*
335*        Set INFO on the first occurrence of a zero pivot
336*
337         IF( INFO.EQ.0 .AND. IINFO.GT.0 )
338     $      INFO = IINFO + K - 1
339*
340*        Adjust IPIV
341*
342         DO 30 J = K, K + KB - 1
343            IF( IPIV( J ).GT.0 ) THEN
344               IPIV( J ) = IPIV( J ) + K - 1
345            ELSE
346               IPIV( J ) = IPIV( J ) - K + 1
347            END IF
348   30    CONTINUE
349*
350*        Increase K and return to the start of the main loop
351*
352         K = K + KB
353         GO TO 20
354*
355      END IF
356*
357   40 CONTINUE
358      WORK( 1 ) = LWKOPT
359      RETURN
360*
361*     End of DSYTRF
362*
363      END
364