LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
clasyf.f
Go to the documentation of this file.
1*> \brief \b CLASYF 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 CLASYF + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clasyf.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clasyf.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clasyf.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CLASYF( 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 A( LDA, * ), W( LDW, * )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> CLASYF 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*> CLASYF is an auxiliary routine called by CSYTRF. 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 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 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 complexSYcomputational
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 clasyf( 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 A( LDA, * ), W( LDW, * )
189* ..
190*
191* =====================================================================
192*
193* .. Parameters ..
194 REAL ZERO, ONE
195 parameter( zero = 0.0e+0, one = 1.0e+0 )
196 REAL EIGHT, SEVTEN
197 parameter( eight = 8.0e+0, sevten = 17.0e+0 )
198 COMPLEX CONE
199 parameter( cone = ( 1.0e+0, 0.0e+0 ) )
200* ..
201* .. Local Scalars ..
202 INTEGER IMAX, J, JB, JJ, JMAX, JP, K, KK, KKW, KP,
203 $ KSTEP, KW
204 REAL ABSAKK, ALPHA, COLMAX, ROWMAX
205 COMPLEX D11, D21, D22, R1, T, Z
206* ..
207* .. External Functions ..
208 LOGICAL LSAME
209 INTEGER ICAMAX
210 EXTERNAL lsame, icamax
211* ..
212* .. External Subroutines ..
213 EXTERNAL ccopy, cgemm, cgemv, cscal, cswap
214* ..
215* .. Intrinsic Functions ..
216 INTRINSIC abs, aimag, max, min, real, sqrt
217* ..
218* .. Statement Functions ..
219 REAL CABS1
220* ..
221* .. Statement Function definitions ..
222 cabs1( z ) = abs( real( z ) ) + abs( aimag( 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 ccopy( k, a( 1, k ), 1, w( 1, kw ), 1 )
254 IF( k.LT.n )
255 $ CALL cgemv( '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* column K, and COLMAX is its absolute value.
267* Determine both COLMAX and IMAX.
268*
269 IF( k.GT.1 ) THEN
270 imax = icamax( k-1, w( 1, kw ), 1 )
271 colmax = cabs1( w( imax, kw ) )
272 ELSE
273 colmax = zero
274 END IF
275*
276 IF( max( absakk, colmax ).EQ.zero ) THEN
277*
278* Column K is zero or underflow: set INFO and continue
279*
280 IF( info.EQ.0 )
281 $ info = k
282 kp = k
283 ELSE
284 IF( absakk.GE.alpha*colmax ) THEN
285*
286* no interchange, use 1-by-1 pivot block
287*
288 kp = k
289 ELSE
290*
291* Copy column IMAX to column KW-1 of W and update it
292*
293 CALL ccopy( imax, a( 1, imax ), 1, w( 1, kw-1 ), 1 )
294 CALL ccopy( k-imax, a( imax, imax+1 ), lda,
295 $ w( imax+1, kw-1 ), 1 )
296 IF( k.LT.n )
297 $ CALL cgemv( 'No transpose', k, n-k, -cone,
298 $ a( 1, k+1 ), lda, w( imax, kw+1 ), ldw,
299 $ cone, w( 1, kw-1 ), 1 )
300*
301* JMAX is the column-index of the largest off-diagonal
302* element in row IMAX, and ROWMAX is its absolute value
303*
304 jmax = imax + icamax( k-imax, w( imax+1, kw-1 ), 1 )
305 rowmax = cabs1( w( jmax, kw-1 ) )
306 IF( imax.GT.1 ) THEN
307 jmax = icamax( imax-1, w( 1, kw-1 ), 1 )
308 rowmax = max( rowmax, cabs1( w( jmax, kw-1 ) ) )
309 END IF
310*
311 IF( absakk.GE.alpha*colmax*( colmax / rowmax ) ) THEN
312*
313* no interchange, use 1-by-1 pivot block
314*
315 kp = k
316 ELSE IF( cabs1( w( imax, kw-1 ) ).GE.alpha*rowmax ) THEN
317*
318* interchange rows and columns K and IMAX, use 1-by-1
319* pivot block
320*
321 kp = imax
322*
323* copy column KW-1 of W to column KW of W
324*
325 CALL ccopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
326 ELSE
327*
328* interchange rows and columns K-1 and IMAX, use 2-by-2
329* pivot block
330*
331 kp = imax
332 kstep = 2
333 END IF
334 END IF
335*
336* ============================================================
337*
338* KK is the column of A where pivoting step stopped
339*
340 kk = k - kstep + 1
341*
342* KKW is the column of W which corresponds to column KK of A
343*
344 kkw = nb + kk - n
345*
346* Interchange rows and columns KP and KK.
347* Updated column KP is already stored in column KKW of W.
348*
349 IF( kp.NE.kk ) THEN
350*
351* Copy non-updated column KK to column KP of submatrix A
352* at step K. No need to copy element into column K
353* (or K and K-1 for 2-by-2 pivot) of A, since these columns
354* will be later overwritten.
355*
356 a( kp, kp ) = a( kk, kk )
357 CALL ccopy( kk-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
358 $ lda )
359 IF( kp.GT.1 )
360 $ CALL ccopy( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
361*
362* Interchange rows KK and KP in last K+1 to N columns of A
363* (columns K (or K and K-1 for 2-by-2 pivot) of A will be
364* later overwritten). Interchange rows KK and KP
365* in last KKW to NB columns of W.
366*
367 IF( k.LT.n )
368 $ CALL cswap( n-k, a( kk, k+1 ), lda, a( kp, k+1 ),
369 $ lda )
370 CALL cswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),
371 $ ldw )
372 END IF
373*
374 IF( kstep.EQ.1 ) THEN
375*
376* 1-by-1 pivot block D(k): column kw of W now holds
377*
378* W(kw) = U(k)*D(k),
379*
380* where U(k) is the k-th column of U
381*
382* Store subdiag. elements of column U(k)
383* and 1-by-1 block D(k) in column k of A.
384* NOTE: Diagonal element U(k,k) is a UNIT element
385* and not stored.
386* A(k,k) := D(k,k) = W(k,kw)
387* A(1:k-1,k) := U(1:k-1,k) = W(1:k-1,kw)/D(k,k)
388*
389 CALL ccopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
390 r1 = cone / a( k, k )
391 CALL cscal( k-1, r1, a( 1, k ), 1 )
392*
393 ELSE
394*
395* 2-by-2 pivot block D(k): columns kw and kw-1 of W now hold
396*
397* ( W(kw-1) W(kw) ) = ( U(k-1) U(k) )*D(k)
398*
399* where U(k) and U(k-1) are the k-th and (k-1)-th columns
400* of U
401*
402* Store U(1:k-2,k-1) and U(1:k-2,k) and 2-by-2
403* block D(k-1:k,k-1:k) in columns k-1 and k of A.
404* NOTE: 2-by-2 diagonal block U(k-1:k,k-1:k) is a UNIT
405* block and not stored.
406* A(k-1:k,k-1:k) := D(k-1:k,k-1:k) = W(k-1:k,kw-1:kw)
407* A(1:k-2,k-1:k) := U(1:k-2,k:k-1:k) =
408* = W(1:k-2,kw-1:kw) * ( D(k-1:k,k-1:k)**(-1) )
409*
410 IF( k.GT.2 ) THEN
411*
412* Compose the columns of the inverse of 2-by-2 pivot
413* block D in the following way to reduce the number
414* of FLOPS when we myltiply panel ( W(kw-1) W(kw) ) by
415* this inverse
416*
417* D**(-1) = ( d11 d21 )**(-1) =
418* ( d21 d22 )
419*
420* = 1/(d11*d22-d21**2) * ( ( d22 ) (-d21 ) ) =
421* ( (-d21 ) ( d11 ) )
422*
423* = 1/d21 * 1/((d11/d21)*(d22/d21)-1) *
424*
425* * ( ( d22/d21 ) ( -1 ) ) =
426* ( ( -1 ) ( d11/d21 ) )
427*
428* = 1/d21 * 1/(D22*D11-1) * ( ( D11 ) ( -1 ) ) =
429* ( ( -1 ) ( D22 ) )
430*
431* = 1/d21 * T * ( ( D11 ) ( -1 ) )
432* ( ( -1 ) ( D22 ) )
433*
434* = D21 * ( ( D11 ) ( -1 ) )
435* ( ( -1 ) ( D22 ) )
436*
437 d21 = w( k-1, kw )
438 d11 = w( k, kw ) / d21
439 d22 = w( k-1, kw-1 ) / d21
440 t = cone / ( d11*d22-cone )
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 d21 = t / d21
447 DO 20 j = 1, k - 2
448 a( j, k-1 ) = d21*( d11*w( j, kw-1 )-w( j, kw ) )
449 a( j, k ) = d21*( d22*w( j, kw )-w( j, kw-1 ) )
450 20 CONTINUE
451 END IF
452*
453* Copy D(k) to A
454*
455 a( k-1, k-1 ) = w( k-1, kw-1 )
456 a( k-1, k ) = w( k-1, kw )
457 a( k, k ) = w( k, kw )
458*
459 END IF
460*
461 END IF
462*
463* Store details of the interchanges in IPIV
464*
465 IF( kstep.EQ.1 ) THEN
466 ipiv( k ) = kp
467 ELSE
468 ipiv( k ) = -kp
469 ipiv( k-1 ) = -kp
470 END IF
471*
472* Decrease K and return to the start of the main loop
473*
474 k = k - kstep
475 GO TO 10
476*
477 30 CONTINUE
478*
479* Update the upper triangle of A11 (= A(1:k,1:k)) as
480*
481* A11 := A11 - U12*D*U12**T = A11 - U12*W**T
482*
483* computing blocks of NB columns at a time
484*
485 DO 50 j = ( ( k-1 ) / nb )*nb + 1, 1, -nb
486 jb = min( nb, k-j+1 )
487*
488* Update the upper triangle of the diagonal block
489*
490 DO 40 jj = j, j + jb - 1
491 CALL cgemv( 'No transpose', jj-j+1, n-k, -cone,
492 $ a( j, k+1 ), lda, w( jj, kw+1 ), ldw, cone,
493 $ a( j, jj ), 1 )
494 40 CONTINUE
495*
496* Update the rectangular superdiagonal block
497*
498 CALL cgemm( 'No transpose', 'Transpose', j-1, jb, n-k,
499 $ -cone, a( 1, k+1 ), lda, w( j, kw+1 ), ldw,
500 $ cone, a( 1, j ), lda )
501 50 CONTINUE
502*
503* Put U12 in standard form by partially undoing the interchanges
504* in columns k+1:n looping backwards from k+1 to n
505*
506 j = k + 1
507 60 CONTINUE
508*
509* Undo the interchanges (if any) of rows JJ and JP at each
510* step J
511*
512* (Here, J is a diagonal index)
513 jj = j
514 jp = ipiv( j )
515 IF( jp.LT.0 ) THEN
516 jp = -jp
517* (Here, J is a diagonal index)
518 j = j + 1
519 END IF
520* (NOTE: Here, J is used to determine row length. Length N-J+1
521* of the rows to swap back doesn't include diagonal element)
522 j = j + 1
523 IF( jp.NE.jj .AND. j.LE.n )
524 $ CALL cswap( n-j+1, a( jp, j ), lda, a( jj, j ), lda )
525 IF( j.LT.n )
526 $ GO TO 60
527*
528* Set KB to the number of columns factorized
529*
530 kb = n - k
531*
532 ELSE
533*
534* Factorize the leading columns of A using the lower triangle
535* of A and working forwards, and compute the matrix W = L21*D
536* for use in updating A22
537*
538* K is the main loop index, increasing from 1 in steps of 1 or 2
539*
540 k = 1
541 70 CONTINUE
542*
543* Exit from loop
544*
545 IF( ( k.GE.nb .AND. nb.LT.n ) .OR. k.GT.n )
546 $ GO TO 90
547*
548* Copy column K of A to column K of W and update it
549*
550 CALL ccopy( n-k+1, a( k, k ), 1, w( k, k ), 1 )
551 CALL cgemv( 'No transpose', n-k+1, k-1, -cone, a( k, 1 ), lda,
552 $ w( k, 1 ), ldw, cone, w( k, k ), 1 )
553*
554 kstep = 1
555*
556* Determine rows and columns to be interchanged and whether
557* a 1-by-1 or 2-by-2 pivot block will be used
558*
559 absakk = cabs1( w( k, k ) )
560*
561* IMAX is the row-index of the largest off-diagonal element in
562* column K, and COLMAX is its absolute value.
563* Determine both COLMAX and IMAX.
564*
565 IF( k.LT.n ) THEN
566 imax = k + icamax( n-k, w( k+1, k ), 1 )
567 colmax = cabs1( w( imax, k ) )
568 ELSE
569 colmax = zero
570 END IF
571*
572 IF( max( absakk, colmax ).EQ.zero ) THEN
573*
574* Column K is zero or underflow: set INFO and continue
575*
576 IF( info.EQ.0 )
577 $ info = k
578 kp = k
579 ELSE
580 IF( absakk.GE.alpha*colmax ) THEN
581*
582* no interchange, use 1-by-1 pivot block
583*
584 kp = k
585 ELSE
586*
587* Copy column IMAX to column K+1 of W and update it
588*
589 CALL ccopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1 )
590 CALL ccopy( n-imax+1, a( imax, imax ), 1, w( imax, k+1 ),
591 $ 1 )
592 CALL cgemv( 'No transpose', n-k+1, k-1, -cone, a( k, 1 ),
593 $ lda, w( imax, 1 ), ldw, cone, w( k, k+1 ),
594 $ 1 )
595*
596* JMAX is the column-index of the largest off-diagonal
597* element in row IMAX, and ROWMAX is its absolute value
598*
599 jmax = k - 1 + icamax( imax-k, w( k, k+1 ), 1 )
600 rowmax = cabs1( w( jmax, k+1 ) )
601 IF( imax.LT.n ) THEN
602 jmax = imax + icamax( n-imax, w( imax+1, k+1 ), 1 )
603 rowmax = max( rowmax, cabs1( w( jmax, k+1 ) ) )
604 END IF
605*
606 IF( absakk.GE.alpha*colmax*( colmax / rowmax ) ) THEN
607*
608* no interchange, use 1-by-1 pivot block
609*
610 kp = k
611 ELSE IF( cabs1( w( imax, k+1 ) ).GE.alpha*rowmax ) THEN
612*
613* interchange rows and columns K and IMAX, use 1-by-1
614* pivot block
615*
616 kp = imax
617*
618* copy column K+1 of W to column K of W
619*
620 CALL ccopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
621 ELSE
622*
623* interchange rows and columns K+1 and IMAX, use 2-by-2
624* pivot block
625*
626 kp = imax
627 kstep = 2
628 END IF
629 END IF
630*
631* ============================================================
632*
633* KK is the column of A where pivoting step stopped
634*
635 kk = k + kstep - 1
636*
637* Interchange rows and columns KP and KK.
638* Updated column KP is already stored in column KK of W.
639*
640 IF( kp.NE.kk ) THEN
641*
642* Copy non-updated column KK to column KP of submatrix A
643* at step K. No need to copy element into column K
644* (or K and K+1 for 2-by-2 pivot) of A, since these columns
645* will be later overwritten.
646*
647 a( kp, kp ) = a( kk, kk )
648 CALL ccopy( kp-kk-1, a( kk+1, kk ), 1, a( kp, kk+1 ),
649 $ lda )
650 IF( kp.LT.n )
651 $ CALL ccopy( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ), 1 )
652*
653* Interchange rows KK and KP in first K-1 columns of A
654* (columns K (or K and K+1 for 2-by-2 pivot) of A will be
655* later overwritten). Interchange rows KK and KP
656* in first KK columns of W.
657*
658 IF( k.GT.1 )
659 $ CALL cswap( k-1, a( kk, 1 ), lda, a( kp, 1 ), lda )
660 CALL cswap( kk, w( kk, 1 ), ldw, w( kp, 1 ), ldw )
661 END IF
662*
663 IF( kstep.EQ.1 ) THEN
664*
665* 1-by-1 pivot block D(k): column k of W now holds
666*
667* W(k) = L(k)*D(k),
668*
669* where L(k) is the k-th column of L
670*
671* Store subdiag. elements of column L(k)
672* and 1-by-1 block D(k) in column k of A.
673* (NOTE: Diagonal element L(k,k) is a UNIT element
674* and not stored)
675* A(k,k) := D(k,k) = W(k,k)
676* A(k+1:N,k) := L(k+1:N,k) = W(k+1:N,k)/D(k,k)
677*
678 CALL ccopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
679 IF( k.LT.n ) THEN
680 r1 = cone / a( k, k )
681 CALL cscal( n-k, r1, a( k+1, k ), 1 )
682 END IF
683*
684 ELSE
685*
686* 2-by-2 pivot block D(k): columns k and k+1 of W now hold
687*
688* ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
689*
690* where L(k) and L(k+1) are the k-th and (k+1)-th columns
691* of L
692*
693* Store L(k+2:N,k) and L(k+2:N,k+1) and 2-by-2
694* block D(k:k+1,k:k+1) in columns k and k+1 of A.
695* (NOTE: 2-by-2 diagonal block L(k:k+1,k:k+1) is a UNIT
696* block and not stored)
697* A(k:k+1,k:k+1) := D(k:k+1,k:k+1) = W(k:k+1,k:k+1)
698* A(k+2:N,k:k+1) := L(k+2:N,k:k+1) =
699* = W(k+2:N,k:k+1) * ( D(k:k+1,k:k+1)**(-1) )
700*
701 IF( k.LT.n-1 ) THEN
702*
703* Compose the columns of the inverse of 2-by-2 pivot
704* block D in the following way to reduce the number
705* of FLOPS when we myltiply panel ( W(k) W(k+1) ) by
706* this inverse
707*
708* D**(-1) = ( d11 d21 )**(-1) =
709* ( d21 d22 )
710*
711* = 1/(d11*d22-d21**2) * ( ( d22 ) (-d21 ) ) =
712* ( (-d21 ) ( d11 ) )
713*
714* = 1/d21 * 1/((d11/d21)*(d22/d21)-1) *
715*
716* * ( ( d22/d21 ) ( -1 ) ) =
717* ( ( -1 ) ( d11/d21 ) )
718*
719* = 1/d21 * 1/(D22*D11-1) * ( ( D11 ) ( -1 ) ) =
720* ( ( -1 ) ( D22 ) )
721*
722* = 1/d21 * T * ( ( D11 ) ( -1 ) )
723* ( ( -1 ) ( D22 ) )
724*
725* = D21 * ( ( D11 ) ( -1 ) )
726* ( ( -1 ) ( D22 ) )
727*
728 d21 = w( k+1, k )
729 d11 = w( k+1, k+1 ) / d21
730 d22 = w( k, k ) / d21
731 t = cone / ( d11*d22-cone )
732 d21 = t / d21
733*
734* Update elements in columns A(k) and A(k+1) as
735* dot products of rows of ( W(k) W(k+1) ) and columns
736* of D**(-1)
737*
738 DO 80 j = k + 2, n
739 a( j, k ) = d21*( d11*w( j, k )-w( j, k+1 ) )
740 a( j, k+1 ) = d21*( d22*w( j, k+1 )-w( j, k ) )
741 80 CONTINUE
742 END IF
743*
744* Copy D(k) to A
745*
746 a( k, k ) = w( k, k )
747 a( k+1, k ) = w( k+1, k )
748 a( k+1, k+1 ) = w( k+1, k+1 )
749*
750 END IF
751*
752 END IF
753*
754* Store details of the interchanges in IPIV
755*
756 IF( kstep.EQ.1 ) THEN
757 ipiv( k ) = kp
758 ELSE
759 ipiv( k ) = -kp
760 ipiv( k+1 ) = -kp
761 END IF
762*
763* Increase K and return to the start of the main loop
764*
765 k = k + kstep
766 GO TO 70
767*
768 90 CONTINUE
769*
770* Update the lower triangle of A22 (= A(k:n,k:n)) as
771*
772* A22 := A22 - L21*D*L21**T = A22 - L21*W**T
773*
774* computing blocks of NB columns at a time
775*
776 DO 110 j = k, n, nb
777 jb = min( nb, n-j+1 )
778*
779* Update the lower triangle of the diagonal block
780*
781 DO 100 jj = j, j + jb - 1
782 CALL cgemv( 'No transpose', j+jb-jj, k-1, -cone,
783 $ a( jj, 1 ), lda, w( jj, 1 ), ldw, cone,
784 $ a( jj, jj ), 1 )
785 100 CONTINUE
786*
787* Update the rectangular subdiagonal block
788*
789 IF( j+jb.LE.n )
790 $ CALL cgemm( 'No transpose', 'Transpose', n-j-jb+1, jb,
791 $ k-1, -cone, a( j+jb, 1 ), lda, w( j, 1 ),
792 $ ldw, cone, a( j+jb, j ), lda )
793 110 CONTINUE
794*
795* Put L21 in standard form by partially undoing the interchanges
796* of rows in columns 1:k-1 looping backwards from k-1 to 1
797*
798 j = k - 1
799 120 CONTINUE
800*
801* Undo the interchanges (if any) of rows JJ and JP at each
802* step J
803*
804* (Here, J is a diagonal index)
805 jj = j
806 jp = ipiv( j )
807 IF( jp.LT.0 ) THEN
808 jp = -jp
809* (Here, J is a diagonal index)
810 j = j - 1
811 END IF
812* (NOTE: Here, J is used to determine row length. Length J
813* of the rows to swap back doesn't include diagonal element)
814 j = j - 1
815 IF( jp.NE.jj .AND. j.GE.1 )
816 $ CALL cswap( j, a( jp, 1 ), lda, a( jj, 1 ), lda )
817 IF( j.GT.1 )
818 $ GO TO 120
819*
820* Set KB to the number of columns factorized
821*
822 kb = k - 1
823*
824 END IF
825 RETURN
826*
827* End of CLASYF
828*
829 END
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:81
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:81
subroutine cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:78
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
Definition: cgemv.f:158
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:187
subroutine clasyf(UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO)
CLASYF computes a partial factorization of a complex symmetric matrix using the Bunch-Kaufman diagona...
Definition: clasyf.f:177