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: *> \ingroup doubleSYcomputational
233: *
234: *> \par Further Details:
235: * =====================
236: *>
237: *> \verbatim
238: *> TODO: put correct description
239: *> \endverbatim
240: *
241: *> \par Contributors:
242: * ==================
243: *>
244: *> \verbatim
245: *>
246: *> December 2016, Igor Kozachenko,
247: *> Computer Science Division,
248: *> University of California, Berkeley
249: *>
250: *> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
251: *> School of Mathematics,
252: *> University of Manchester
253: *>
254: *> \endverbatim
255: *
256: * =====================================================================
257: SUBROUTINE DSYTRF_RK( UPLO, N, A, LDA, E, IPIV, WORK, LWORK,
258: $ INFO )
259: *
260: * -- LAPACK computational routine --
261: * -- LAPACK is a software package provided by Univ. of Tennessee, --
262: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
263: *
264: * .. Scalar Arguments ..
265: CHARACTER UPLO
266: INTEGER INFO, LDA, LWORK, N
267: * ..
268: * .. Array Arguments ..
269: INTEGER IPIV( * )
270: DOUBLE PRECISION A( LDA, * ), E( * ), WORK( * )
271: * ..
272: *
273: * =====================================================================
274: *
275: * .. Local Scalars ..
276: LOGICAL LQUERY, UPPER
277: INTEGER I, IINFO, IP, IWS, K, KB, LDWORK, LWKOPT,
278: $ NB, NBMIN
279: * ..
280: * .. External Functions ..
281: LOGICAL LSAME
282: INTEGER ILAENV
283: EXTERNAL LSAME, ILAENV
284: * ..
285: * .. External Subroutines ..
286: EXTERNAL DLASYF_RK, DSYTF2_RK, DSWAP, XERBLA
287: * ..
288: * .. Intrinsic Functions ..
289: INTRINSIC ABS, MAX
290: * ..
291: * .. Executable Statements ..
292: *
293: * Test the input parameters.
294: *
295: INFO = 0
296: UPPER = LSAME( UPLO, 'U' )
297: LQUERY = ( LWORK.EQ.-1 )
298: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
299: INFO = -1
300: ELSE IF( N.LT.0 ) THEN
301: INFO = -2
302: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
303: INFO = -4
304: ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN
305: INFO = -8
306: END IF
307: *
308: IF( INFO.EQ.0 ) THEN
309: *
310: * Determine the block size
311: *
312: NB = ILAENV( 1, 'DSYTRF_RK', UPLO, N, -1, -1, -1 )
313: LWKOPT = N*NB
314: WORK( 1 ) = LWKOPT
315: END IF
316: *
317: IF( INFO.NE.0 ) THEN
318: CALL XERBLA( 'DSYTRF_RK', -INFO )
319: RETURN
320: ELSE IF( LQUERY ) THEN
321: RETURN
322: END IF
323: *
324: NBMIN = 2
325: LDWORK = N
326: IF( NB.GT.1 .AND. NB.LT.N ) THEN
327: IWS = LDWORK*NB
328: IF( LWORK.LT.IWS ) THEN
329: NB = MAX( LWORK / LDWORK, 1 )
330: NBMIN = MAX( 2, ILAENV( 2, 'DSYTRF_RK',
331: $ UPLO, N, -1, -1, -1 ) )
332: END IF
333: ELSE
334: IWS = 1
335: END IF
336: IF( NB.LT.NBMIN )
337: $ NB = N
338: *
339: IF( UPPER ) THEN
340: *
341: * Factorize A as U*D*U**T using the upper triangle of A
342: *
343: * K is the main loop index, decreasing from N to 1 in steps of
344: * KB, where KB is the number of columns factorized by DLASYF_RK;
345: * KB is either NB or NB-1, or K for the last block
346: *
347: K = N
348: 10 CONTINUE
349: *
350: * If K < 1, exit from loop
351: *
352: IF( K.LT.1 )
353: $ GO TO 15
354: *
355: IF( K.GT.NB ) THEN
356: *
357: * Factorize columns k-kb+1:k of A and use blocked code to
358: * update columns 1:k-kb
359: *
360: CALL DLASYF_RK( UPLO, K, NB, KB, A, LDA, E,
361: $ IPIV, WORK, LDWORK, IINFO )
362: ELSE
363: *
364: * Use unblocked code to factorize columns 1:k of A
365: *
366: CALL DSYTF2_RK( UPLO, K, A, LDA, E, IPIV, IINFO )
367: KB = K
368: END IF
369: *
370: * Set INFO on the first occurrence of a zero pivot
371: *
372: IF( INFO.EQ.0 .AND. IINFO.GT.0 )
373: $ INFO = IINFO
374: *
375: * No need to adjust IPIV
376: *
377: *
378: * Apply permutations to the leading panel 1:k-1
379: *
380: * Read IPIV from the last block factored, i.e.
381: * indices k-kb+1:k and apply row permutations to the
382: * last k+1 colunms k+1:N after that block
383: * (We can do the simple loop over IPIV with decrement -1,
384: * since the ABS value of IPIV( I ) represents the row index
385: * of the interchange with row i in both 1x1 and 2x2 pivot cases)
386: *
387: IF( K.LT.N ) THEN
388: DO I = K, ( K - KB + 1 ), -1
389: IP = ABS( IPIV( I ) )
390: IF( IP.NE.I ) THEN
391: CALL DSWAP( N-K, A( I, K+1 ), LDA,
392: $ A( IP, K+1 ), LDA )
393: END IF
394: END DO
395: END IF
396: *
397: * Decrease K and return to the start of the main loop
398: *
399: K = K - KB
400: GO TO 10
401: *
402: * This label is the exit from main loop over K decreasing
403: * from N to 1 in steps of KB
404: *
405: 15 CONTINUE
406: *
407: ELSE
408: *
409: * Factorize A as L*D*L**T using the lower triangle of A
410: *
411: * K is the main loop index, increasing from 1 to N in steps of
412: * KB, where KB is the number of columns factorized by DLASYF_RK;
413: * KB is either NB or NB-1, or N-K+1 for the last block
414: *
415: K = 1
416: 20 CONTINUE
417: *
418: * If K > N, exit from loop
419: *
420: IF( K.GT.N )
421: $ GO TO 35
422: *
423: IF( K.LE.N-NB ) THEN
424: *
425: * Factorize columns k:k+kb-1 of A and use blocked code to
426: * update columns k+kb:n
427: *
428: CALL DLASYF_RK( UPLO, N-K+1, NB, KB, A( K, K ), LDA, E( K ),
429: $ IPIV( K ), WORK, LDWORK, IINFO )
430:
431:
432: ELSE
433: *
434: * Use unblocked code to factorize columns k:n of A
435: *
436: CALL DSYTF2_RK( UPLO, N-K+1, A( K, K ), LDA, E( K ),
437: $ IPIV( K ), IINFO )
438: KB = N - K + 1
439: *
440: END IF
441: *
442: * Set INFO on the first occurrence of a zero pivot
443: *
444: IF( INFO.EQ.0 .AND. IINFO.GT.0 )
445: $ INFO = IINFO + K - 1
446: *
447: * Adjust IPIV
448: *
449: DO I = K, K + KB - 1
450: IF( IPIV( I ).GT.0 ) THEN
451: IPIV( I ) = IPIV( I ) + K - 1
452: ELSE
453: IPIV( I ) = IPIV( I ) - K + 1
454: END IF
455: END DO
456: *
457: * Apply permutations to the leading panel 1:k-1
458: *
459: * Read IPIV from the last block factored, i.e.
460: * indices k:k+kb-1 and apply row permutations to the
461: * first k-1 colunms 1:k-1 before that block
462: * (We can do the simple loop over IPIV with increment 1,
463: * since the ABS value of IPIV( I ) represents the row index
464: * of the interchange with row i in both 1x1 and 2x2 pivot cases)
465: *
466: IF( K.GT.1 ) THEN
467: DO I = K, ( K + KB - 1 ), 1
468: IP = ABS( IPIV( I ) )
469: IF( IP.NE.I ) THEN
470: CALL DSWAP( K-1, A( I, 1 ), LDA,
471: $ A( IP, 1 ), LDA )
472: END IF
473: END DO
474: END IF
475: *
476: * Increase K and return to the start of the main loop
477: *
478: K = K + KB
479: GO TO 20
480: *
481: * This label is the exit from main loop over K increasing
482: * from 1 to N in steps of KB
483: *
484: 35 CONTINUE
485: *
486: * End Lower
487: *
488: END IF
489: *
490: WORK( 1 ) = LWKOPT
491: RETURN
492: *
493: * End of DSYTRF_RK
494: *
495: 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