1: *> \brief \b ZSYTRF_AA
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download ZSYTRF_AA + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zsytrf_aa.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zsytrf_aa.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zsytrf_aa.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZSYTRF_AA( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO )
22: *
23: * .. Scalar Arguments ..
24: * CHARACTER UPLO
25: * INTEGER N, LDA, LWORK, INFO
26: * ..
27: * .. Array Arguments ..
28: * INTEGER IPIV( * )
29: * COMPLEX*16 A( LDA, * ), WORK( * )
30: * ..
31: *
32: *> \par Purpose:
33: * =============
34: *>
35: *> \verbatim
36: *>
37: *> ZSYTRF_AA computes the factorization of a complex symmetric matrix A
38: *> using the Aasen's algorithm. The form of the factorization is
39: *>
40: *> A = U*T*U**T or A = L*T*L**T
41: *>
42: *> where U (or L) is a product of permutation and unit upper (lower)
43: *> triangular matrices, and T is a complex symmetric tridiagonal matrix.
44: *>
45: *> This is the blocked version of the algorithm, calling Level 3 BLAS.
46: *> \endverbatim
47: *
48: * Arguments:
49: * ==========
50: *
51: *> \param[in] UPLO
52: *> \verbatim
53: *> UPLO is CHARACTER*1
54: *> = 'U': Upper triangle of A is stored;
55: *> = 'L': Lower triangle of A is stored.
56: *> \endverbatim
57: *>
58: *> \param[in] N
59: *> \verbatim
60: *> N is INTEGER
61: *> The order of the matrix A. N >= 0.
62: *> \endverbatim
63: *>
64: *> \param[in,out] A
65: *> \verbatim
66: *> A is COMPLEX*16 array, dimension (LDA,N)
67: *> On entry, the symmetric matrix A. If UPLO = 'U', the leading
68: *> N-by-N upper triangular part of A contains the upper
69: *> triangular part of the matrix A, and the strictly lower
70: *> triangular part of A is not referenced. If UPLO = 'L', the
71: *> leading N-by-N lower triangular part of A contains the lower
72: *> triangular part of the matrix A, and the strictly upper
73: *> triangular part of A is not referenced.
74: *>
75: *> On exit, the tridiagonal matrix is stored in the diagonals
76: *> and the subdiagonals of A just below (or above) the diagonals,
77: *> and L is stored below (or above) the subdiaonals, when UPLO
78: *> is 'L' (or 'U').
79: *> \endverbatim
80: *>
81: *> \param[in] LDA
82: *> \verbatim
83: *> LDA is INTEGER
84: *> The leading dimension of the array A. LDA >= max(1,N).
85: *> \endverbatim
86: *>
87: *> \param[out] IPIV
88: *> \verbatim
89: *> IPIV is INTEGER array, dimension (N)
90: *> On exit, it contains the details of the interchanges, i.e.,
91: *> the row and column k of A were interchanged with the
92: *> row and column IPIV(k).
93: *> \endverbatim
94: *>
95: *> \param[out] WORK
96: *> \verbatim
97: *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
98: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
99: *> \endverbatim
100: *>
101: *> \param[in] LWORK
102: *> \verbatim
103: *> LWORK is INTEGER
104: *> The length of WORK. LWORK >=MAX(1,2*N). For optimum performance
105: *> LWORK >= N*(1+NB), where NB is the optimal blocksize.
106: *>
107: *> If LWORK = -1, then a workspace query is assumed; the routine
108: *> only calculates the optimal size of the WORK array, returns
109: *> this value as the first entry of the WORK array, and no error
110: *> message related to LWORK is issued by XERBLA.
111: *> \endverbatim
112: *>
113: *> \param[out] INFO
114: *> \verbatim
115: *> INFO is INTEGER
116: *> = 0: successful exit
117: *> < 0: if INFO = -i, the i-th argument had an illegal value.
118: *> \endverbatim
119: *
120: * Authors:
121: * ========
122: *
123: *> \author Univ. of Tennessee
124: *> \author Univ. of California Berkeley
125: *> \author Univ. of Colorado Denver
126: *> \author NAG Ltd.
127: *
128: *> \date November 2017
129: *
130: *> \ingroup complex16SYcomputational
131: *
132: * =====================================================================
133: SUBROUTINE ZSYTRF_AA( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO)
134: *
135: * -- LAPACK computational routine (version 3.8.0) --
136: * -- LAPACK is a software package provided by Univ. of Tennessee, --
137: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
138: * November 2017
139: *
140: IMPLICIT NONE
141: *
142: * .. Scalar Arguments ..
143: CHARACTER UPLO
144: INTEGER N, LDA, LWORK, INFO
145: * ..
146: * .. Array Arguments ..
147: INTEGER IPIV( * )
148: COMPLEX*16 A( LDA, * ), WORK( * )
149: * ..
150: *
151: * =====================================================================
152: * .. Parameters ..
153: COMPLEX*16 ZERO, ONE
154: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
155: *
156: * .. Local Scalars ..
157: LOGICAL LQUERY, UPPER
158: INTEGER J, LWKOPT
159: INTEGER NB, MJ, NJ, K1, K2, J1, J2, J3, JB
160: COMPLEX*16 ALPHA
161: * ..
162: * .. External Functions ..
163: LOGICAL LSAME
164: INTEGER ILAENV
165: EXTERNAL LSAME, ILAENV
166: * ..
167: * .. External Subroutines ..
168: EXTERNAL ZLASYF_AA, ZGEMM, ZGEMV, ZSCAL, ZCOPY,
169: $ ZSWAP, XERBLA
170: * ..
171: * .. Intrinsic Functions ..
172: INTRINSIC MAX
173: * ..
174: * .. Executable Statements ..
175: *
176: * Determine the block size
177: *
178: NB = ILAENV( 1, 'ZSYTRF_AA', UPLO, N, -1, -1, -1 )
179: *
180: * Test the input parameters.
181: *
182: INFO = 0
183: UPPER = LSAME( UPLO, 'U' )
184: LQUERY = ( LWORK.EQ.-1 )
185: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
186: INFO = -1
187: ELSE IF( N.LT.0 ) THEN
188: INFO = -2
189: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
190: INFO = -4
191: ELSE IF( LWORK.LT.MAX( 1, 2*N ) .AND. .NOT.LQUERY ) THEN
192: INFO = -7
193: END IF
194: *
195: IF( INFO.EQ.0 ) THEN
196: LWKOPT = (NB+1)*N
197: WORK( 1 ) = LWKOPT
198: END IF
199: *
200: IF( INFO.NE.0 ) THEN
201: CALL XERBLA( 'ZSYTRF_AA', -INFO )
202: RETURN
203: ELSE IF( LQUERY ) THEN
204: RETURN
205: END IF
206: *
207: * Quick return
208: *
209: IF ( N.EQ.0 ) THEN
210: RETURN
211: ENDIF
212: IPIV( 1 ) = 1
213: IF ( N.EQ.1 ) THEN
214: RETURN
215: END IF
216: *
217: * Adjust block size based on the workspace size
218: *
219: IF( LWORK.LT.((1+NB)*N) ) THEN
220: NB = ( LWORK-N ) / N
221: END IF
222: *
223: IF( UPPER ) THEN
224: *
225: * .....................................................
226: * Factorize A as L*D*L**T using the upper triangle of A
227: * .....................................................
228: *
229: * Copy first row A(1, 1:N) into H(1:n) (stored in WORK(1:N))
230: *
231: CALL ZCOPY( N, A( 1, 1 ), LDA, WORK( 1 ), 1 )
232: *
233: * J is the main loop index, increasing from 1 to N in steps of
234: * JB, where JB is the number of columns factorized by ZLASYF;
235: * JB is either NB, or N-J+1 for the last block
236: *
237: J = 0
238: 10 CONTINUE
239: IF( J.GE.N )
240: $ GO TO 20
241: *
242: * each step of the main loop
243: * J is the last column of the previous panel
244: * J1 is the first column of the current panel
245: * K1 identifies if the previous column of the panel has been
246: * explicitly stored, e.g., K1=1 for the first panel, and
247: * K1=0 for the rest
248: *
249: J1 = J + 1
250: JB = MIN( N-J1+1, NB )
251: K1 = MAX(1, J)-J
252: *
253: * Panel factorization
254: *
255: CALL ZLASYF_AA( UPLO, 2-K1, N-J, JB,
256: $ A( MAX(1, J), J+1 ), LDA,
257: $ IPIV( J+1 ), WORK, N, WORK( N*NB+1 ) )
258: *
259: * Ajust IPIV and apply it back (J-th step picks (J+1)-th pivot)
260: *
261: DO J2 = J+2, MIN(N, J+JB+1)
262: IPIV( J2 ) = IPIV( J2 ) + J
263: IF( (J2.NE.IPIV(J2)) .AND. ((J1-K1).GT.2) ) THEN
264: CALL ZSWAP( J1-K1-2, A( 1, J2 ), 1,
265: $ A( 1, IPIV(J2) ), 1 )
266: END IF
267: END DO
268: J = J + JB
269: *
270: * Trailing submatrix update, where
271: * the row A(J1-1, J2-1:N) stores U(J1, J2+1:N) and
272: * WORK stores the current block of the auxiriarly matrix H
273: *
274: IF( J.LT.N ) THEN
275: *
276: * If first panel and JB=1 (NB=1), then nothing to do
277: *
278: IF( J1.GT.1 .OR. JB.GT.1 ) THEN
279: *
280: * Merge rank-1 update with BLAS-3 update
281: *
282: ALPHA = A( J, J+1 )
283: A( J, J+1 ) = ONE
284: CALL ZCOPY( N-J, A( J-1, J+1 ), LDA,
285: $ WORK( (J+1-J1+1)+JB*N ), 1 )
286: CALL ZSCAL( N-J, ALPHA, WORK( (J+1-J1+1)+JB*N ), 1 )
287: *
288: * K1 identifies if the previous column of the panel has been
289: * explicitly stored, e.g., K1=1 and K2= 0 for the first panel,
290: * while K1=0 and K2=1 for the rest
291: *
292: IF( J1.GT.1 ) THEN
293: *
294: * Not first panel
295: *
296: K2 = 1
297: ELSE
298: *
299: * First panel
300: *
301: K2 = 0
302: *
303: * First update skips the first column
304: *
305: JB = JB - 1
306: END IF
307: *
308: DO J2 = J+1, N, NB
309: NJ = MIN( NB, N-J2+1 )
310: *
311: * Update (J2, J2) diagonal block with ZGEMV
312: *
313: J3 = J2
314: DO MJ = NJ-1, 1, -1
315: CALL ZGEMV( 'No transpose', MJ, JB+1,
316: $ -ONE, WORK( J3-J1+1+K1*N ), N,
317: $ A( J1-K2, J3 ), 1,
318: $ ONE, A( J3, J3 ), LDA )
319: J3 = J3 + 1
320: END DO
321: *
322: * Update off-diagonal block of J2-th block row with ZGEMM
323: *
324: CALL ZGEMM( 'Transpose', 'Transpose',
325: $ NJ, N-J3+1, JB+1,
326: $ -ONE, A( J1-K2, J2 ), LDA,
327: $ WORK( J3-J1+1+K1*N ), N,
328: $ ONE, A( J2, J3 ), LDA )
329: END DO
330: *
331: * Recover T( J, J+1 )
332: *
333: A( J, J+1 ) = ALPHA
334: END IF
335: *
336: * WORK(J+1, 1) stores H(J+1, 1)
337: *
338: CALL ZCOPY( N-J, A( J+1, J+1 ), LDA, WORK( 1 ), 1 )
339: END IF
340: GO TO 10
341: ELSE
342: *
343: * .....................................................
344: * Factorize A as L*D*L**T using the lower triangle of A
345: * .....................................................
346: *
347: * copy first column A(1:N, 1) into H(1:N, 1)
348: * (stored in WORK(1:N))
349: *
350: CALL ZCOPY( N, A( 1, 1 ), 1, WORK( 1 ), 1 )
351: *
352: * J is the main loop index, increasing from 1 to N in steps of
353: * JB, where JB is the number of columns factorized by ZLASYF;
354: * JB is either NB, or N-J+1 for the last block
355: *
356: J = 0
357: 11 CONTINUE
358: IF( J.GE.N )
359: $ GO TO 20
360: *
361: * each step of the main loop
362: * J is the last column of the previous panel
363: * J1 is the first column of the current panel
364: * K1 identifies if the previous column of the panel has been
365: * explicitly stored, e.g., K1=1 for the first panel, and
366: * K1=0 for the rest
367: *
368: J1 = J+1
369: JB = MIN( N-J1+1, NB )
370: K1 = MAX(1, J)-J
371: *
372: * Panel factorization
373: *
374: CALL ZLASYF_AA( UPLO, 2-K1, N-J, JB,
375: $ A( J+1, MAX(1, J) ), LDA,
376: $ IPIV( J+1 ), WORK, N, WORK( N*NB+1 ) )
377: *
378: * Ajust IPIV and apply it back (J-th step picks (J+1)-th pivot)
379: *
380: DO J2 = J+2, MIN(N, J+JB+1)
381: IPIV( J2 ) = IPIV( J2 ) + J
382: IF( (J2.NE.IPIV(J2)) .AND. ((J1-K1).GT.2) ) THEN
383: CALL ZSWAP( J1-K1-2, A( J2, 1 ), LDA,
384: $ A( IPIV(J2), 1 ), LDA )
385: END IF
386: END DO
387: J = J + JB
388: *
389: * Trailing submatrix update, where
390: * A(J2+1, J1-1) stores L(J2+1, J1) and
391: * WORK(J2+1, 1) stores H(J2+1, 1)
392: *
393: IF( J.LT.N ) THEN
394: *
395: * if first panel and JB=1 (NB=1), then nothing to do
396: *
397: IF( J1.GT.1 .OR. JB.GT.1 ) THEN
398: *
399: * Merge rank-1 update with BLAS-3 update
400: *
401: ALPHA = A( J+1, J )
402: A( J+1, J ) = ONE
403: CALL ZCOPY( N-J, A( J+1, J-1 ), 1,
404: $ WORK( (J+1-J1+1)+JB*N ), 1 )
405: CALL ZSCAL( N-J, ALPHA, WORK( (J+1-J1+1)+JB*N ), 1 )
406: *
407: * K1 identifies if the previous column of the panel has been
408: * explicitly stored, e.g., K1=1 and K2= 0 for the first panel,
409: * while K1=0 and K2=1 for the rest
410: *
411: IF( J1.GT.1 ) THEN
412: *
413: * Not first panel
414: *
415: K2 = 1
416: ELSE
417: *
418: * First panel
419: *
420: K2 = 0
421: *
422: * First update skips the first column
423: *
424: JB = JB - 1
425: END IF
426: *
427: DO J2 = J+1, N, NB
428: NJ = MIN( NB, N-J2+1 )
429: *
430: * Update (J2, J2) diagonal block with ZGEMV
431: *
432: J3 = J2
433: DO MJ = NJ-1, 1, -1
434: CALL ZGEMV( 'No transpose', MJ, JB+1,
435: $ -ONE, WORK( J3-J1+1+K1*N ), N,
436: $ A( J3, J1-K2 ), LDA,
437: $ ONE, A( J3, J3 ), 1 )
438: J3 = J3 + 1
439: END DO
440: *
441: * Update off-diagonal block in J2-th block column with ZGEMM
442: *
443: CALL ZGEMM( 'No transpose', 'Transpose',
444: $ N-J3+1, NJ, JB+1,
445: $ -ONE, WORK( J3-J1+1+K1*N ), N,
446: $ A( J2, J1-K2 ), LDA,
447: $ ONE, A( J3, J2 ), LDA )
448: END DO
449: *
450: * Recover T( J+1, J )
451: *
452: A( J+1, J ) = ALPHA
453: END IF
454: *
455: * WORK(J+1, 1) stores H(J+1, 1)
456: *
457: CALL ZCOPY( N-J, A( J+1, J+1 ), 1, WORK( 1 ), 1 )
458: END IF
459: GO TO 11
460: END IF
461: *
462: 20 CONTINUE
463: RETURN
464: *
465: * End of ZSYTRF_AA
466: *
467: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zsytrf_aa.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack