1: *> \brief \b DSYTRF_RK computes the factorization of a real symmetric indefinite matrix using the bounded Bunch-Kaufman (rook) diagonal pivoting method (BLAS3 blocked algorithm).
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DSYTRF_RK + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsytrf_rk.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsytrf_rk.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsytrf_rk.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DSYTRF_RK( UPLO, N, A, LDA, E, IPIV, WORK, LWORK,
22: * INFO )
23: *
24: * .. Scalar Arguments ..
25: * CHARACTER UPLO
26: * INTEGER INFO, LDA, LWORK, N
27: * ..
28: * .. Array Arguments ..
29: * INTEGER IPIV( * )
30: * DOUBLE PRECISION A( LDA, * ), E ( * ), WORK( * )
31: * ..
32: *
33: *
34: *> \par Purpose:
35: * =============
36: *>
37: *> \verbatim
38: *> DSYTRF_RK computes the factorization of a real symmetric matrix A
39: *> using the bounded Bunch-Kaufman (rook) diagonal pivoting method:
40: *>
41: *> A = P*U*D*(U**T)*(P**T) or A = P*L*D*(L**T)*(P**T),
42: *>
43: *> where U (or L) is unit upper (or lower) triangular matrix,
44: *> U**T (or L**T) is the transpose of U (or L), P is a permutation
45: *> matrix, P**T is the transpose of P, and D is symmetric and block
46: *> diagonal with 1-by-1 and 2-by-2 diagonal blocks.
47: *>
48: *> This is the blocked version of the algorithm, calling Level 3 BLAS.
49: *> For more information see Further Details section.
50: *> \endverbatim
51: *
52: * Arguments:
53: * ==========
54: *
55: *> \param[in] UPLO
56: *> \verbatim
57: *> UPLO is CHARACTER*1
58: *> Specifies whether the upper or lower triangular part of the
59: *> symmetric matrix A is stored:
60: *> = 'U': Upper triangular
61: *> = 'L': Lower triangular
62: *> \endverbatim
63: *>
64: *> \param[in] N
65: *> \verbatim
66: *> N is INTEGER
67: *> The order of the matrix A. N >= 0.
68: *> \endverbatim
69: *>
70: *> \param[in,out] A
71: *> \verbatim
72: *> A is DOUBLE PRECISION array, dimension (LDA,N)
73: *> On entry, the symmetric matrix A.
74: *> If UPLO = 'U': the leading N-by-N upper triangular part
75: *> of A contains the upper triangular part of the matrix A,
76: *> and the strictly lower triangular part of A is not
77: *> referenced.
78: *>
79: *> If UPLO = 'L': the leading N-by-N lower triangular part
80: *> of A contains the lower triangular part of the matrix A,
81: *> and the strictly upper triangular part of A is not
82: *> referenced.
83: *>
84: *> On exit, contains:
85: *> a) ONLY diagonal elements of the symmetric block diagonal
86: *> matrix D on the diagonal of A, i.e. D(k,k) = A(k,k);
87: *> (superdiagonal (or subdiagonal) elements of D
88: *> are stored on exit in array E), and
89: *> b) If UPLO = 'U': factor U in the superdiagonal part of A.
90: *> If UPLO = 'L': factor L in the subdiagonal part of A.
91: *> \endverbatim
92: *>
93: *> \param[in] LDA
94: *> \verbatim
95: *> LDA is INTEGER
96: *> The leading dimension of the array A. LDA >= max(1,N).
97: *> \endverbatim
98: *>
99: *> \param[out] E
100: *> \verbatim
101: *> E is DOUBLE PRECISION array, dimension (N)
102: *> On exit, contains the superdiagonal (or subdiagonal)
103: *> elements of the symmetric block diagonal matrix D
104: *> with 1-by-1 or 2-by-2 diagonal blocks, where
105: *> If UPLO = 'U': E(i) = D(i-1,i), i=2:N, E(1) is set to 0;
106: *> If UPLO = 'L': E(i) = D(i+1,i), i=1:N-1, E(N) is set to 0.
107: *>
108: *> NOTE: For 1-by-1 diagonal block D(k), where
109: *> 1 <= k <= N, the element E(k) is set to 0 in both
110: *> UPLO = 'U' or UPLO = 'L' cases.
111: *> \endverbatim
112: *>
113: *> \param[out] IPIV
114: *> \verbatim
115: *> IPIV is INTEGER array, dimension (N)
116: *> IPIV describes the permutation matrix P in the factorization
117: *> of matrix A as follows. The absolute value of IPIV(k)
118: *> represents the index of row and column that were
119: *> interchanged with the k-th row and column. The value of UPLO
120: *> describes the order in which the interchanges were applied.
121: *> Also, the sign of IPIV represents the block structure of
122: *> the symmetric block diagonal matrix D with 1-by-1 or 2-by-2
123: *> diagonal blocks which correspond to 1 or 2 interchanges
124: *> at each factorization step. For more info see Further
125: *> Details section.
126: *>
127: *> If UPLO = 'U',
128: *> ( in factorization order, k decreases from N to 1 ):
129: *> a) A single positive entry IPIV(k) > 0 means:
130: *> D(k,k) is a 1-by-1 diagonal block.
131: *> If IPIV(k) != k, rows and columns k and IPIV(k) were
132: *> interchanged in the matrix A(1:N,1:N);
133: *> If IPIV(k) = k, no interchange occurred.
134: *>
135: *> b) A pair of consecutive negative entries
136: *> IPIV(k) < 0 and IPIV(k-1) < 0 means:
137: *> D(k-1:k,k-1:k) is a 2-by-2 diagonal block.
138: *> (NOTE: negative entries in IPIV appear ONLY in pairs).
139: *> 1) If -IPIV(k) != k, rows and columns
140: *> k and -IPIV(k) were interchanged
141: *> in the matrix A(1:N,1:N).
142: *> If -IPIV(k) = k, no interchange occurred.
143: *> 2) If -IPIV(k-1) != k-1, rows and columns
144: *> k-1 and -IPIV(k-1) were interchanged
145: *> in the matrix A(1:N,1:N).
146: *> If -IPIV(k-1) = k-1, no interchange occurred.
147: *>
148: *> c) In both cases a) and b), always ABS( IPIV(k) ) <= k.
149: *>
150: *> d) NOTE: Any entry IPIV(k) is always NONZERO on output.
151: *>
152: *> If UPLO = 'L',
153: *> ( in factorization order, k increases from 1 to N ):
154: *> a) A single positive entry IPIV(k) > 0 means:
155: *> D(k,k) is a 1-by-1 diagonal block.
156: *> If IPIV(k) != k, rows and columns k and IPIV(k) were
157: *> interchanged in the matrix A(1:N,1:N).
158: *> If IPIV(k) = k, no interchange occurred.
159: *>
160: *> b) A pair of consecutive negative entries
161: *> IPIV(k) < 0 and IPIV(k+1) < 0 means:
162: *> D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
163: *> (NOTE: negative entries in IPIV appear ONLY in pairs).
164: *> 1) If -IPIV(k) != k, rows and columns
165: *> k and -IPIV(k) were interchanged
166: *> in the matrix A(1:N,1:N).
167: *> If -IPIV(k) = k, no interchange occurred.
168: *> 2) If -IPIV(k+1) != k+1, rows and columns
169: *> k-1 and -IPIV(k-1) were interchanged
170: *> in the matrix A(1:N,1:N).
171: *> If -IPIV(k+1) = k+1, no interchange occurred.
172: *>
173: *> c) In both cases a) and b), always ABS( IPIV(k) ) >= k.
174: *>
175: *> d) NOTE: Any entry IPIV(k) is always NONZERO on output.
176: *> \endverbatim
177: *>
178: *> \param[out] WORK
179: *> \verbatim
180: *> WORK is DOUBLE PRECISION array, dimension ( MAX(1,LWORK) ).
181: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
182: *> \endverbatim
183: *>
184: *> \param[in] LWORK
185: *> \verbatim
186: *> LWORK is INTEGER
187: *> The length of WORK. LWORK >=1. For best performance
188: *> LWORK >= N*NB, where NB is the block size returned
189: *> by ILAENV.
190: *>
191: *> If LWORK = -1, then a workspace query is assumed;
192: *> the routine only calculates the optimal size of the WORK
193: *> array, returns this value as the first entry of the WORK
194: *> array, and no error message related to LWORK is issued
195: *> by XERBLA.
196: *> \endverbatim
197: *>
198: *> \param[out] INFO
199: *> \verbatim
200: *> INFO is INTEGER
201: *> = 0: successful exit
202: *>
203: *> < 0: If INFO = -k, the k-th argument had an illegal value
204: *>
205: *> > 0: If INFO = k, the matrix A is singular, because:
206: *> If UPLO = 'U': column k in the upper
207: *> triangular part of A contains all zeros.
208: *> If UPLO = 'L': column k in the lower
209: *> triangular part of A contains all zeros.
210: *>
211: *> Therefore D(k,k) is exactly zero, and superdiagonal
212: *> elements of column k of U (or subdiagonal elements of
213: *> column k of L ) are all zeros. The factorization has
214: *> been completed, but the block diagonal matrix D is
215: *> exactly singular, and division by zero will occur if
216: *> it is used to solve a system of equations.
217: *>
218: *> NOTE: INFO only stores the first occurrence of
219: *> a singularity, any subsequent occurrence of singularity
220: *> is not stored in INFO even though the factorization
221: *> always completes.
222: *> \endverbatim
223: *
224: * Authors:
225: * ========
226: *
227: *> \author Univ. of Tennessee
228: *> \author Univ. of California Berkeley
229: *> \author Univ. of Colorado Denver
230: *> \author NAG Ltd.
231: *
232: *> \date December 2016
233: *
234: *> \ingroup doubleSYcomputational
235: *
236: *> \par Further Details:
237: * =====================
238: *>
239: *> \verbatim
240: *> TODO: put correct description
241: *> \endverbatim
242: *
243: *> \par Contributors:
244: * ==================
245: *>
246: *> \verbatim
247: *>
248: *> December 2016, Igor Kozachenko,
249: *> Computer Science Division,
250: *> University of California, Berkeley
251: *>
252: *> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
253: *> School of Mathematics,
254: *> University of Manchester
255: *>
256: *> \endverbatim
257: *
258: * =====================================================================
259: SUBROUTINE DSYTRF_RK( UPLO, N, A, LDA, E, IPIV, WORK, LWORK,
260: $ INFO )
261: *
262: * -- LAPACK computational routine (version 3.7.0) --
263: * -- LAPACK is a software package provided by Univ. of Tennessee, --
264: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
265: * December 2016
266: *
267: * .. Scalar Arguments ..
268: CHARACTER UPLO
269: INTEGER INFO, LDA, LWORK, N
270: * ..
271: * .. Array Arguments ..
272: INTEGER IPIV( * )
273: DOUBLE PRECISION A( LDA, * ), E( * ), WORK( * )
274: * ..
275: *
276: * =====================================================================
277: *
278: * .. Local Scalars ..
279: LOGICAL LQUERY, UPPER
280: INTEGER I, IINFO, IP, IWS, K, KB, LDWORK, LWKOPT,
281: $ NB, NBMIN
282: * ..
283: * .. External Functions ..
284: LOGICAL LSAME
285: INTEGER ILAENV
286: EXTERNAL LSAME, ILAENV
287: * ..
288: * .. External Subroutines ..
289: EXTERNAL DLASYF_RK, DSYTF2_RK, DSWAP, XERBLA
290: * ..
291: * .. Intrinsic Functions ..
292: INTRINSIC ABS, MAX
293: * ..
294: * .. Executable Statements ..
295: *
296: * Test the input parameters.
297: *
298: INFO = 0
299: UPPER = LSAME( UPLO, 'U' )
300: LQUERY = ( LWORK.EQ.-1 )
301: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
302: INFO = -1
303: ELSE IF( N.LT.0 ) THEN
304: INFO = -2
305: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
306: INFO = -4
307: ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN
308: INFO = -8
309: END IF
310: *
311: IF( INFO.EQ.0 ) THEN
312: *
313: * Determine the block size
314: *
315: NB = ILAENV( 1, 'DSYTRF_RK', UPLO, N, -1, -1, -1 )
316: LWKOPT = N*NB
317: WORK( 1 ) = LWKOPT
318: END IF
319: *
320: IF( INFO.NE.0 ) THEN
321: CALL XERBLA( 'DSYTRF_RK', -INFO )
322: RETURN
323: ELSE IF( LQUERY ) THEN
324: RETURN
325: END IF
326: *
327: NBMIN = 2
328: LDWORK = N
329: IF( NB.GT.1 .AND. NB.LT.N ) THEN
330: IWS = LDWORK*NB
331: IF( LWORK.LT.IWS ) THEN
332: NB = MAX( LWORK / LDWORK, 1 )
333: NBMIN = MAX( 2, ILAENV( 2, 'DSYTRF_RK',
334: $ UPLO, N, -1, -1, -1 ) )
335: END IF
336: ELSE
337: IWS = 1
338: END IF
339: IF( NB.LT.NBMIN )
340: $ NB = N
341: *
342: IF( UPPER ) THEN
343: *
344: * Factorize A as U*D*U**T using the upper triangle of A
345: *
346: * K is the main loop index, decreasing from N to 1 in steps of
347: * KB, where KB is the number of columns factorized by DLASYF_RK;
348: * KB is either NB or NB-1, or K for the last block
349: *
350: K = N
351: 10 CONTINUE
352: *
353: * If K < 1, exit from loop
354: *
355: IF( K.LT.1 )
356: $ GO TO 15
357: *
358: IF( K.GT.NB ) THEN
359: *
360: * Factorize columns k-kb+1:k of A and use blocked code to
361: * update columns 1:k-kb
362: *
363: CALL DLASYF_RK( UPLO, K, NB, KB, A, LDA, E,
364: $ IPIV, WORK, LDWORK, IINFO )
365: ELSE
366: *
367: * Use unblocked code to factorize columns 1:k of A
368: *
369: CALL DSYTF2_RK( UPLO, K, A, LDA, E, IPIV, IINFO )
370: KB = K
371: END IF
372: *
373: * Set INFO on the first occurrence of a zero pivot
374: *
375: IF( INFO.EQ.0 .AND. IINFO.GT.0 )
376: $ INFO = IINFO
377: *
378: * No need to adjust IPIV
379: *
380: *
381: * Apply permutations to the leading panel 1:k-1
382: *
383: * Read IPIV from the last block factored, i.e.
384: * indices k-kb+1:k and apply row permutations to the
385: * last k+1 colunms k+1:N after that block
386: * (We can do the simple loop over IPIV with decrement -1,
387: * since the ABS value of IPIV( I ) represents the row index
388: * of the interchange with row i in both 1x1 and 2x2 pivot cases)
389: *
390: IF( K.LT.N ) THEN
391: DO I = K, ( K - KB + 1 ), -1
392: IP = ABS( IPIV( I ) )
393: IF( IP.NE.I ) THEN
394: CALL DSWAP( N-K, A( I, K+1 ), LDA,
395: $ A( IP, K+1 ), LDA )
396: END IF
397: END DO
398: END IF
399: *
400: * Decrease K and return to the start of the main loop
401: *
402: K = K - KB
403: GO TO 10
404: *
405: * This label is the exit from main loop over K decreasing
406: * from N to 1 in steps of KB
407: *
408: 15 CONTINUE
409: *
410: ELSE
411: *
412: * Factorize A as L*D*L**T using the lower triangle of A
413: *
414: * K is the main loop index, increasing from 1 to N in steps of
415: * KB, where KB is the number of columns factorized by DLASYF_RK;
416: * KB is either NB or NB-1, or N-K+1 for the last block
417: *
418: K = 1
419: 20 CONTINUE
420: *
421: * If K > N, exit from loop
422: *
423: IF( K.GT.N )
424: $ GO TO 35
425: *
426: IF( K.LE.N-NB ) THEN
427: *
428: * Factorize columns k:k+kb-1 of A and use blocked code to
429: * update columns k+kb:n
430: *
431: CALL DLASYF_RK( UPLO, N-K+1, NB, KB, A( K, K ), LDA, E( K ),
432: $ IPIV( K ), WORK, LDWORK, IINFO )
433:
434:
435: ELSE
436: *
437: * Use unblocked code to factorize columns k:n of A
438: *
439: CALL DSYTF2_RK( UPLO, N-K+1, A( K, K ), LDA, E( K ),
440: $ IPIV( K ), IINFO )
441: KB = N - K + 1
442: *
443: END IF
444: *
445: * Set INFO on the first occurrence of a zero pivot
446: *
447: IF( INFO.EQ.0 .AND. IINFO.GT.0 )
448: $ INFO = IINFO + K - 1
449: *
450: * Adjust IPIV
451: *
452: DO I = K, K + KB - 1
453: IF( IPIV( I ).GT.0 ) THEN
454: IPIV( I ) = IPIV( I ) + K - 1
455: ELSE
456: IPIV( I ) = IPIV( I ) - K + 1
457: END IF
458: END DO
459: *
460: * Apply permutations to the leading panel 1:k-1
461: *
462: * Read IPIV from the last block factored, i.e.
463: * indices k:k+kb-1 and apply row permutations to the
464: * first k-1 colunms 1:k-1 before that block
465: * (We can do the simple loop over IPIV with increment 1,
466: * since the ABS value of IPIV( I ) represents the row index
467: * of the interchange with row i in both 1x1 and 2x2 pivot cases)
468: *
469: IF( K.GT.1 ) THEN
470: DO I = K, ( K + KB - 1 ), 1
471: IP = ABS( IPIV( I ) )
472: IF( IP.NE.I ) THEN
473: CALL DSWAP( K-1, A( I, 1 ), LDA,
474: $ A( IP, 1 ), LDA )
475: END IF
476: END DO
477: END IF
478: *
479: * Increase K and return to the start of the main loop
480: *
481: K = K + KB
482: GO TO 20
483: *
484: * This label is the exit from main loop over K increasing
485: * from 1 to N in steps of KB
486: *
487: 35 CONTINUE
488: *
489: * End Lower
490: *
491: END IF
492: *
493: WORK( 1 ) = LWKOPT
494: RETURN
495: *
496: * End of DSYTRF_RK
497: *
498: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dsytrf_rk.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack