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