LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
zlasyf.f
Go to the documentation of this file.
1*> \brief \b ZLASYF computes a partial factorization of a complex symmetric matrix using the Bunch-Kaufman diagonal pivoting method.
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download ZLASYF + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlasyf.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlasyf.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlasyf.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE ZLASYF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO )
22*
23* .. Scalar Arguments ..
24* CHARACTER UPLO
25* INTEGER INFO, KB, LDA, LDW, N, NB
26* ..
27* .. Array Arguments ..
28* INTEGER IPIV( * )
29* COMPLEX*16 A( LDA, * ), W( LDW, * )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> ZLASYF computes a partial factorization of a complex symmetric matrix
39*> A using the Bunch-Kaufman diagonal pivoting method. The partial
40*> factorization has the form:
41*>
42*> A = ( I U12 ) ( A11 0 ) ( I 0 ) if UPLO = 'U', or:
43*> ( 0 U22 ) ( 0 D ) ( U12**T U22**T )
44*>
45*> A = ( L11 0 ) ( D 0 ) ( L11**T L21**T ) if UPLO = 'L'
46*> ( L21 I ) ( 0 A22 ) ( 0 I )
47*>
48*> where the order of D is at most NB. The actual order is returned in
49*> the argument KB, and is either NB or NB-1, or N if N <= NB.
50*> Note that U**T denotes the transpose of U.
51*>
52*> ZLASYF is an auxiliary routine called by ZSYTRF. It uses blocked code
53*> (calling Level 3 BLAS) to update the submatrix A11 (if UPLO = 'U') or
54*> A22 (if UPLO = 'L').
55*> \endverbatim
56*
57* Arguments:
58* ==========
59*
60*> \param[in] UPLO
61*> \verbatim
62*> UPLO is CHARACTER*1
63*> Specifies whether the upper or lower triangular part of the
64*> symmetric matrix A is stored:
65*> = 'U': Upper triangular
66*> = 'L': Lower triangular
67*> \endverbatim
68*>
69*> \param[in] N
70*> \verbatim
71*> N is INTEGER
72*> The order of the matrix A. N >= 0.
73*> \endverbatim
74*>
75*> \param[in] NB
76*> \verbatim
77*> NB is INTEGER
78*> The maximum number of columns of the matrix A that should be
79*> factored. NB should be at least 2 to allow for 2-by-2 pivot
80*> blocks.
81*> \endverbatim
82*>
83*> \param[out] KB
84*> \verbatim
85*> KB is INTEGER
86*> The number of columns of A that were actually factored.
87*> KB is either NB-1 or NB, or N if N <= NB.
88*> \endverbatim
89*>
90*> \param[in,out] A
91*> \verbatim
92*> A is COMPLEX*16 array, dimension (LDA,N)
93*> On entry, the symmetric matrix A. If UPLO = 'U', the leading
94*> n-by-n upper triangular part of A contains the upper
95*> triangular part of the matrix A, and the strictly lower
96*> triangular part of A is not referenced. If UPLO = 'L', the
97*> leading n-by-n lower triangular part of A contains the lower
98*> triangular part of the matrix A, and the strictly upper
99*> triangular part of A is not referenced.
100*> On exit, A contains details of the partial factorization.
101*> \endverbatim
102*>
103*> \param[in] LDA
104*> \verbatim
105*> LDA is INTEGER
106*> The leading dimension of the array A. LDA >= max(1,N).
107*> \endverbatim
108*>
109*> \param[out] IPIV
110*> \verbatim
111*> IPIV is INTEGER array, dimension (N)
112*> Details of the interchanges and the block structure of D.
113*>
114*> If UPLO = 'U':
115*> Only the last KB elements of IPIV are set.
116*>
117*> If IPIV(k) > 0, then rows and columns k and IPIV(k) were
118*> interchanged and D(k,k) is a 1-by-1 diagonal block.
119*>
120*> If IPIV(k) = IPIV(k-1) < 0, then rows and columns
121*> k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
122*> is a 2-by-2 diagonal block.
123*>
124*> If UPLO = 'L':
125*> Only the first KB elements of IPIV are set.
126*>
127*> If IPIV(k) > 0, then rows and columns k and IPIV(k) were
128*> interchanged and D(k,k) is a 1-by-1 diagonal block.
129*>
130*> If IPIV(k) = IPIV(k+1) < 0, then rows and columns
131*> k+1 and -IPIV(k) were interchanged and D(k:k+1,k:k+1)
132*> is a 2-by-2 diagonal block.
133*> \endverbatim
134*>
135*> \param[out] W
136*> \verbatim
137*> W is COMPLEX*16 array, dimension (LDW,NB)
138*> \endverbatim
139*>
140*> \param[in] LDW
141*> \verbatim
142*> LDW is INTEGER
143*> The leading dimension of the array W. LDW >= max(1,N).
144*> \endverbatim
145*>
146*> \param[out] INFO
147*> \verbatim
148*> INFO is INTEGER
149*> = 0: successful exit
150*> > 0: if INFO = k, D(k,k) is exactly zero. The factorization
151*> has been completed, but the block diagonal matrix D is
152*> exactly singular.
153*> \endverbatim
154*
155* Authors:
156* ========
157*
158*> \author Univ. of Tennessee
159*> \author Univ. of California Berkeley
160*> \author Univ. of Colorado Denver
161*> \author NAG Ltd.
162*
163*> \ingroup lahef
164*
165*> \par Contributors:
166* ==================
167*>
168*> \verbatim
169*>
170*> November 2013, Igor Kozachenko,
171*> Computer Science Division,
172*> University of California, Berkeley
173*> \endverbatim
174*
175* =====================================================================
176 SUBROUTINE zlasyf( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO )
177*
178* -- LAPACK computational routine --
179* -- LAPACK is a software package provided by Univ. of Tennessee, --
180* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
181*
182* .. Scalar Arguments ..
183 CHARACTER UPLO
184 INTEGER INFO, KB, LDA, LDW, N, NB
185* ..
186* .. Array Arguments ..
187 INTEGER IPIV( * )
188 COMPLEX*16 A( LDA, * ), W( LDW, * )
189* ..
190*
191* =====================================================================
192*
193* .. Parameters ..
194 DOUBLE PRECISION ZERO, ONE
195 parameter( zero = 0.0d+0, one = 1.0d+0 )
196 DOUBLE PRECISION EIGHT, SEVTEN
197 parameter( eight = 8.0d+0, sevten = 17.0d+0 )
198 COMPLEX*16 CONE
199 parameter( cone = ( 1.0d+0, 0.0d+0 ) )
200* ..
201* .. Local Scalars ..
202 INTEGER IMAX, J, JB, JJ, JMAX, JP, K, KK, KKW, KP,
203 $ KSTEP, KW
204 DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, ROWMAX
205 COMPLEX*16 D11, D21, D22, R1, T, Z
206* ..
207* .. External Functions ..
208 LOGICAL LSAME
209 INTEGER IZAMAX
210 EXTERNAL lsame, izamax
211* ..
212* .. External Subroutines ..
213 EXTERNAL zcopy, zgemm, zgemv, zscal, zswap
214* ..
215* .. Intrinsic Functions ..
216 INTRINSIC abs, dble, dimag, max, min, sqrt
217* ..
218* .. Statement Functions ..
219 DOUBLE PRECISION CABS1
220* ..
221* .. Statement Function definitions ..
222 cabs1( z ) = abs( dble( z ) ) + abs( dimag( z ) )
223* ..
224* .. Executable Statements ..
225*
226 info = 0
227*
228* Initialize ALPHA for use in choosing pivot block size.
229*
230 alpha = ( one+sqrt( sevten ) ) / eight
231*
232 IF( lsame( uplo, 'U' ) ) THEN
233*
234* Factorize the trailing columns of A using the upper triangle
235* of A and working backwards, and compute the matrix W = U12*D
236* for use in updating A11
237*
238* K is the main loop index, decreasing from N in steps of 1 or 2
239*
240* KW is the column of W which corresponds to column K of A
241*
242 k = n
243 10 CONTINUE
244 kw = nb + k - n
245*
246* Exit from loop
247*
248 IF( ( k.LE.n-nb+1 .AND. nb.LT.n ) .OR. k.LT.1 )
249 $ GO TO 30
250*
251* Copy column K of A to column KW of W and update it
252*
253 CALL zcopy( k, a( 1, k ), 1, w( 1, kw ), 1 )
254 IF( k.LT.n )
255 $ CALL zgemv( 'No transpose', k, n-k, -cone, a( 1, k+1 ), lda,
256 $ w( k, kw+1 ), ldw, cone, w( 1, kw ), 1 )
257*
258 kstep = 1
259*
260* Determine rows and columns to be interchanged and whether
261* a 1-by-1 or 2-by-2 pivot block will be used
262*
263 absakk = cabs1( w( k, kw ) )
264*
265* IMAX is the row-index of the largest off-diagonal element in
266
267*
268 IF( k.GT.1 ) THEN
269 imax = izamax( k-1, w( 1, kw ), 1 )
270 colmax = cabs1( w( imax, kw ) )
271 ELSE
272 colmax = zero
273 END IF
274*
275 IF( max( absakk, colmax ).EQ.zero ) THEN
276*
277* Column K is zero or underflow: set INFO and continue
278*
279 IF( info.EQ.0 )
280 $ info = k
281 kp = k
282 ELSE
283 IF( absakk.GE.alpha*colmax ) THEN
284*
285* no interchange, use 1-by-1 pivot block
286*
287 kp = k
288 ELSE
289*
290* Copy column IMAX to column KW-1 of W and update it
291*
292 CALL zcopy( imax, a( 1, imax ), 1, w( 1, kw-1 ), 1 )
293 CALL zcopy( k-imax, a( imax, imax+1 ), lda,
294 $ w( imax+1, kw-1 ), 1 )
295 IF( k.LT.n )
296 $ CALL zgemv( 'No transpose', k, n-k, -cone,
297 $ a( 1, k+1 ), lda, w( imax, kw+1 ), ldw,
298 $ cone, w( 1, kw-1 ), 1 )
299*
300* JMAX is the column-index of the largest off-diagonal
301* element in row IMAX, and ROWMAX is its absolute value
302*
303 jmax = imax + izamax( k-imax, w( imax+1, kw-1 ), 1 )
304 rowmax = cabs1( w( jmax, kw-1 ) )
305 IF( imax.GT.1 ) THEN
306 jmax = izamax( imax-1, w( 1, kw-1 ), 1 )
307 rowmax = max( rowmax, cabs1( w( jmax, kw-1 ) ) )
308 END IF
309*
310 IF( absakk.GE.alpha*colmax*( colmax / rowmax ) ) THEN
311*
312* no interchange, use 1-by-1 pivot block
313*
314 kp = k
315 ELSE IF( cabs1( w( imax, kw-1 ) ).GE.alpha*rowmax ) THEN
316*
317* interchange rows and columns K and IMAX, use 1-by-1
318* pivot block
319*
320 kp = imax
321*
322* copy column KW-1 of W to column KW of W
323*
324 CALL zcopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
325 ELSE
326*
327* interchange rows and columns K-1 and IMAX, use 2-by-2
328* pivot block
329*
330 kp = imax
331 kstep = 2
332 END IF
333 END IF
334*
335* ============================================================
336*
337* KK is the column of A where pivoting step stopped
338*
339 kk = k - kstep + 1
340*
341* KKW is the column of W which corresponds to column KK of A
342*
343 kkw = nb + kk - n
344*
345* Interchange rows and columns KP and KK.
346* Updated column KP is already stored in column KKW of W.
347*
348 IF( kp.NE.kk ) THEN
349*
350* Copy non-updated column KK to column KP of submatrix A
351* at step K. No need to copy element into column K
352* (or K and K-1 for 2-by-2 pivot) of A, since these columns
353* will be later overwritten.
354*
355 a( kp, kp ) = a( kk, kk )
356 CALL zcopy( kk-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
357 $ lda )
358 IF( kp.GT.1 )
359 $ CALL zcopy( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
360*
361* Interchange rows KK and KP in last K+1 to N columns of A
362* (columns K (or K and K-1 for 2-by-2 pivot) of A will be
363* later overwritten). Interchange rows KK and KP
364* in last KKW to NB columns of W.
365*
366 IF( k.LT.n )
367 $ CALL zswap( n-k, a( kk, k+1 ), lda, a( kp, k+1 ),
368 $ lda )
369 CALL zswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),
370 $ ldw )
371 END IF
372*
373 IF( kstep.EQ.1 ) THEN
374*
375* 1-by-1 pivot block D(k): column kw of W now holds
376*
377* W(kw) = U(k)*D(k),
378*
379* where U(k) is the k-th column of U
380*
381* Store subdiag. elements of column U(k)
382* and 1-by-1 block D(k) in column k of A.
383* NOTE: Diagonal element U(k,k) is a UNIT element
384* and not stored.
385* A(k,k) := D(k,k) = W(k,kw)
386* A(1:k-1,k) := U(1:k-1,k) = W(1:k-1,kw)/D(k,k)
387*
388 CALL zcopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
389 r1 = cone / a( k, k )
390 CALL zscal( k-1, r1, a( 1, k ), 1 )
391*
392 ELSE
393*
394* 2-by-2 pivot block D(k): columns kw and kw-1 of W now hold
395*
396* ( W(kw-1) W(kw) ) = ( U(k-1) U(k) )*D(k)
397*
398* where U(k) and U(k-1) are the k-th and (k-1)-th columns
399* of U
400*
401* Store U(1:k-2,k-1) and U(1:k-2,k) and 2-by-2
402* block D(k-1:k,k-1:k) in columns k-1 and k of A.
403* NOTE: 2-by-2 diagonal block U(k-1:k,k-1:k) is a UNIT
404* block and not stored.
405* A(k-1:k,k-1:k) := D(k-1:k,k-1:k) = W(k-1:k,kw-1:kw)
406* A(1:k-2,k-1:k) := U(1:k-2,k:k-1:k) =
407* = W(1:k-2,kw-1:kw) * ( D(k-1:k,k-1:k)**(-1) )
408*
409 IF( k.GT.2 ) THEN
410*
411* Compose the columns of the inverse of 2-by-2 pivot
412* block D in the following way to reduce the number
413* of FLOPS when we myltiply panel ( W(kw-1) W(kw) ) by
414* this inverse
415*
416* D**(-1) = ( d11 d21 )**(-1) =
417* ( d21 d22 )
418*
419* = 1/(d11*d22-d21**2) * ( ( d22 ) (-d21 ) ) =
420* ( (-d21 ) ( d11 ) )
421*
422* = 1/d21 * 1/((d11/d21)*(d22/d21)-1) *
423*
424* * ( ( d22/d21 ) ( -1 ) ) =
425* ( ( -1 ) ( d11/d21 ) )
426*
427* = 1/d21 * 1/(D22*D11-1) * ( ( D11 ) ( -1 ) ) =
428* ( ( -1 ) ( D22 ) )
429*
430* = 1/d21 * T * ( ( D11 ) ( -1 ) )
431* ( ( -1 ) ( D22 ) )
432*
433* = D21 * ( ( D11 ) ( -1 ) )
434* ( ( -1 ) ( D22 ) )
435*
436 d21 = w( k-1, kw )
437 d11 = w( k, kw ) / d21
438 d22 = w( k-1, kw-1 ) / d21
439 t = cone / ( d11*d22-cone )
440 d21 = t / d21
441*
442* Update elements in columns A(k-1) and A(k) as
443* dot products of rows of ( W(kw-1) W(kw) ) and columns
444* of D**(-1)
445*
446 DO 20 j = 1, k - 2
447 a( j, k-1 ) = d21*( d11*w( j, kw-1 )-w( j, kw ) )
448 a( j, k ) = d21*( d22*w( j, kw )-w( j, kw-1 ) )
449 20 CONTINUE
450 END IF
451*
452* Copy D(k) to A
453*
454 a( k-1, k-1 ) = w( k-1, kw-1 )
455 a( k-1, k ) = w( k-1, kw )
456 a( k, k ) = w( k, kw )
457*
458 END IF
459*
460 END IF
461*
462* Store details of the interchanges in IPIV
463*
464 IF( kstep.EQ.1 ) THEN
465 ipiv( k ) = kp
466 ELSE
467 ipiv( k ) = -kp
468 ipiv( k-1 ) = -kp
469 END IF
470*
471* Decrease K and return to the start of the main loop
472*
473 k = k - kstep
474 GO TO 10
475*
476 30 CONTINUE
477*
478* Update the upper triangle of A11 (= A(1:k,1:k)) as
479*
480* A11 := A11 - U12*D*U12**T = A11 - U12*W**T
481*
482* computing blocks of NB columns at a time
483*
484 DO 50 j = ( ( k-1 ) / nb )*nb + 1, 1, -nb
485 jb = min( nb, k-j+1 )
486*
487* Update the upper triangle of the diagonal block
488*
489 DO 40 jj = j, j + jb - 1
490 CALL zgemv( 'No transpose', jj-j+1, n-k, -cone,
491 $ a( j, k+1 ), lda, w( jj, kw+1 ), ldw, cone,
492 $ a( j, jj ), 1 )
493 40 CONTINUE
494*
495* Update the rectangular superdiagonal block
496*
497 CALL zgemm( 'No transpose', 'Transpose', j-1, jb, n-k,
498 $ -cone, a( 1, k+1 ), lda, w( j, kw+1 ), ldw,
499 $ cone, a( 1, j ), lda )
500 50 CONTINUE
501*
502* Put U12 in standard form by partially undoing the interchanges
503* in columns k+1:n looping backwards from k+1 to n
504*
505 j = k + 1
506 60 CONTINUE
507*
508* Undo the interchanges (if any) of rows JJ and JP at each
509* step J
510*
511* (Here, J is a diagonal index)
512 jj = j
513 jp = ipiv( j )
514 IF( jp.LT.0 ) THEN
515 jp = -jp
516* (Here, J is a diagonal index)
517 j = j + 1
518 END IF
519* (NOTE: Here, J is used to determine row length. Length N-J+1
520* of the rows to swap back doesn't include diagonal element)
521 j = j + 1
522 IF( jp.NE.jj .AND. j.LE.n )
523 $ CALL zswap( n-j+1, a( jp, j ), lda, a( jj, j ), lda )
524 IF( j.LT.n )
525 $ GO TO 60
526*
527* Set KB to the number of columns factorized
528*
529 kb = n - k
530*
531 ELSE
532*
533* Factorize the leading columns of A using the lower triangle
534* of A and working forwards, and compute the matrix W = L21*D
535* for use in updating A22
536*
537* K is the main loop index, increasing from 1 in steps of 1 or 2
538*
539 k = 1
540 70 CONTINUE
541*
542* Exit from loop
543*
544 IF( ( k.GE.nb .AND. nb.LT.n ) .OR. k.GT.n )
545 $ GO TO 90
546*
547* Copy column K of A to column K of W and update it
548*
549 CALL zcopy( n-k+1, a( k, k ), 1, w( k, k ), 1 )
550 CALL zgemv( 'No transpose', n-k+1, k-1, -cone, a( k, 1 ), lda,
551 $ w( k, 1 ), ldw, cone, w( k, k ), 1 )
552*
553 kstep = 1
554*
555* Determine rows and columns to be interchanged and whether
556* a 1-by-1 or 2-by-2 pivot block will be used
557*
558 absakk = cabs1( w( k, k ) )
559*
560* IMAX is the row-index of the largest off-diagonal element in
561
562*
563 IF( k.LT.n ) THEN
564 imax = k + izamax( n-k, w( k+1, k ), 1 )
565 colmax = cabs1( w( imax, k ) )
566 ELSE
567 colmax = zero
568 END IF
569*
570 IF( max( absakk, colmax ).EQ.zero ) THEN
571*
572* Column K is zero or underflow: set INFO and continue
573*
574 IF( info.EQ.0 )
575 $ info = k
576 kp = k
577 ELSE
578 IF( absakk.GE.alpha*colmax ) THEN
579*
580* no interchange, use 1-by-1 pivot block
581*
582 kp = k
583 ELSE
584*
585* Copy column IMAX to column K+1 of W and update it
586*
587 CALL zcopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1 )
588 CALL zcopy( n-imax+1, a( imax, imax ), 1, w( imax, k+1 ),
589 $ 1 )
590 CALL zgemv( 'No transpose', n-k+1, k-1, -cone, a( k, 1 ),
591 $ lda, w( imax, 1 ), ldw, cone, w( k, k+1 ),
592 $ 1 )
593*
594* JMAX is the column-index of the largest off-diagonal
595* element in row IMAX, and ROWMAX is its absolute value
596*
597 jmax = k - 1 + izamax( imax-k, w( k, k+1 ), 1 )
598 rowmax = cabs1( w( jmax, k+1 ) )
599 IF( imax.LT.n ) THEN
600 jmax = imax + izamax( n-imax, w( imax+1, k+1 ), 1 )
601 rowmax = max( rowmax, cabs1( w( jmax, k+1 ) ) )
602 END IF
603*
604 IF( absakk.GE.alpha*colmax*( colmax / rowmax ) ) THEN
605*
606* no interchange, use 1-by-1 pivot block
607*
608 kp = k
609 ELSE IF( cabs1( w( imax, k+1 ) ).GE.alpha*rowmax ) THEN
610*
611* interchange rows and columns K and IMAX, use 1-by-1
612* pivot block
613*
614 kp = imax
615*
616* copy column K+1 of W to column K of W
617*
618 CALL zcopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
619 ELSE
620*
621* interchange rows and columns K+1 and IMAX, use 2-by-2
622* pivot block
623*
624 kp = imax
625 kstep = 2
626 END IF
627 END IF
628*
629* ============================================================
630*
631* KK is the column of A where pivoting step stopped
632*
633 kk = k + kstep - 1
634*
635* Interchange rows and columns KP and KK.
636* Updated column KP is already stored in column KK of W.
637*
638 IF( kp.NE.kk ) THEN
639*
640* Copy non-updated column KK to column KP of submatrix A
641* at step K. No need to copy element into column K
642* (or K and K+1 for 2-by-2 pivot) of A, since these columns
643* will be later overwritten.
644*
645 a( kp, kp ) = a( kk, kk )
646 CALL zcopy( kp-kk-1, a( kk+1, kk ), 1, a( kp, kk+1 ),
647 $ lda )
648 IF( kp.LT.n )
649 $ CALL zcopy( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ), 1 )
650*
651* Interchange rows KK and KP in first K-1 columns of A
652* (columns K (or K and K+1 for 2-by-2 pivot) of A will be
653* later overwritten). Interchange rows KK and KP
654* in first KK columns of W.
655*
656 IF( k.GT.1 )
657 $ CALL zswap( k-1, a( kk, 1 ), lda, a( kp, 1 ), lda )
658 CALL zswap( kk, w( kk, 1 ), ldw, w( kp, 1 ), ldw )
659 END IF
660*
661 IF( kstep.EQ.1 ) THEN
662*
663* 1-by-1 pivot block D(k): column k of W now holds
664*
665* W(k) = L(k)*D(k),
666*
667* where L(k) is the k-th column of L
668*
669* Store subdiag. elements of column L(k)
670* and 1-by-1 block D(k) in column k of A.
671* (NOTE: Diagonal element L(k,k) is a UNIT element
672* and not stored)
673* A(k,k) := D(k,k) = W(k,k)
674* A(k+1:N,k) := L(k+1:N,k) = W(k+1:N,k)/D(k,k)
675*
676 CALL zcopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
677 IF( k.LT.n ) THEN
678 r1 = cone / a( k, k )
679 CALL zscal( n-k, r1, a( k+1, k ), 1 )
680 END IF
681*
682 ELSE
683*
684* 2-by-2 pivot block D(k): columns k and k+1 of W now hold
685*
686* ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
687*
688* where L(k) and L(k+1) are the k-th and (k+1)-th columns
689* of L
690*
691* Store L(k+2:N,k) and L(k+2:N,k+1) and 2-by-2
692* block D(k:k+1,k:k+1) in columns k and k+1 of A.
693* (NOTE: 2-by-2 diagonal block L(k:k+1,k:k+1) is a UNIT
694* block and not stored)
695* A(k:k+1,k:k+1) := D(k:k+1,k:k+1) = W(k:k+1,k:k+1)
696* A(k+2:N,k:k+1) := L(k+2:N,k:k+1) =
697* = W(k+2:N,k:k+1) * ( D(k:k+1,k:k+1)**(-1) )
698*
699 IF( k.LT.n-1 ) THEN
700*
701* Compose the columns of the inverse of 2-by-2 pivot
702* block D in the following way to reduce the number
703* of FLOPS when we myltiply panel ( W(k) W(k+1) ) by
704* this inverse
705*
706* D**(-1) = ( d11 d21 )**(-1) =
707* ( d21 d22 )
708*
709* = 1/(d11*d22-d21**2) * ( ( d22 ) (-d21 ) ) =
710* ( (-d21 ) ( d11 ) )
711*
712* = 1/d21 * 1/((d11/d21)*(d22/d21)-1) *
713*
714* * ( ( d22/d21 ) ( -1 ) ) =
715* ( ( -1 ) ( d11/d21 ) )
716*
717* = 1/d21 * 1/(D22*D11-1) * ( ( D11 ) ( -1 ) ) =
718* ( ( -1 ) ( D22 ) )
719*
720* = 1/d21 * T * ( ( D11 ) ( -1 ) )
721* ( ( -1 ) ( D22 ) )
722*
723* = D21 * ( ( D11 ) ( -1 ) )
724* ( ( -1 ) ( D22 ) )
725*
726 d21 = w( k+1, k )
727 d11 = w( k+1, k+1 ) / d21
728 d22 = w( k, k ) / d21
729 t = cone / ( d11*d22-cone )
730 d21 = t / d21
731*
732* Update elements in columns A(k) and A(k+1) as
733* dot products of rows of ( W(k) W(k+1) ) and columns
734* of D**(-1)
735*
736 DO 80 j = k + 2, n
737 a( j, k ) = d21*( d11*w( j, k )-w( j, k+1 ) )
738 a( j, k+1 ) = d21*( d22*w( j, k+1 )-w( j, k ) )
739 80 CONTINUE
740 END IF
741*
742* Copy D(k) to A
743*
744 a( k, k ) = w( k, k )
745 a( k+1, k ) = w( k+1, k )
746 a( k+1, k+1 ) = w( k+1, k+1 )
747*
748 END IF
749*
750 END IF
751*
752* Store details of the interchanges in IPIV
753*
754 IF( kstep.EQ.1 ) THEN
755 ipiv( k ) = kp
756 ELSE
757 ipiv( k ) = -kp
758 ipiv( k+1 ) = -kp
759 END IF
760*
761* Increase K and return to the start of the main loop
762*
763 k = k + kstep
764 GO TO 70
765*
766 90 CONTINUE
767*
768* Update the lower triangle of A22 (= A(k:n,k:n)) as
769*
770* A22 := A22 - L21*D*L21**T = A22 - L21*W**T
771*
772* computing blocks of NB columns at a time
773*
774 DO 110 j = k, n, nb
775 jb = min( nb, n-j+1 )
776*
777* Update the lower triangle of the diagonal block
778*
779 DO 100 jj = j, j + jb - 1
780 CALL zgemv( 'No transpose', j+jb-jj, k-1, -cone,
781 $ a( jj, 1 ), lda, w( jj, 1 ), ldw, cone,
782 $ a( jj, jj ), 1 )
783 100 CONTINUE
784*
785* Update the rectangular subdiagonal block
786*
787 IF( j+jb.LE.n )
788 $ CALL zgemm( 'No transpose', 'Transpose', n-j-jb+1, jb,
789 $ k-1, -cone, a( j+jb, 1 ), lda, w( j, 1 ),
790 $ ldw, cone, a( j+jb, j ), lda )
791 110 CONTINUE
792*
793* Put L21 in standard form by partially undoing the interchanges
794* of rows in columns 1:k-1 looping backwards from k-1 to 1
795*
796 j = k - 1
797 120 CONTINUE
798*
799* Undo the interchanges (if any) of rows JJ and JP at each
800* step J
801*
802* (Here, J is a diagonal index)
803 jj = j
804 jp = ipiv( j )
805 IF( jp.LT.0 ) THEN
806 jp = -jp
807* (Here, J is a diagonal index)
808 j = j - 1
809 END IF
810* (NOTE: Here, J is used to determine row length. Length J
811* of the rows to swap back doesn't include diagonal element)
812 j = j - 1
813 IF( jp.NE.jj .AND. j.GE.1 )
814 $ CALL zswap( j, a( jp, 1 ), lda, a( jj, 1 ), lda )
815 IF( j.GT.1 )
816 $ GO TO 120
817*
818* Set KB to the number of columns factorized
819*
820 kb = k - 1
821*
822 END IF
823 RETURN
824*
825* End of ZLASYF
826*
827 END
zcopy
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
Definition zcopy.f:81
zgemm
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
Definition zgemm.f:188
zgemv
subroutine zgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
ZGEMV
Definition zgemv.f:160
zlasyf
subroutine zlasyf(uplo, n, nb, kb, a, lda, ipiv, w, ldw, info)
ZLASYF computes a partial factorization of a complex symmetric matrix using the Bunch-Kaufman diagona...
Definition zlasyf.f:177
zscal
subroutine zscal(n, za, zx, incx)
ZSCAL
Definition zscal.f:78
zswap
subroutine zswap(n, zx, incx, zy, incy)
ZSWAP
Definition zswap.f:81