LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dlasyf_aa.f
Go to the documentation of this file.
1*> \brief \b DLASYF_AA
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DLASYF_AA + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasyf_aa.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasyf_aa.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasyf_aa.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DLASYF_AA( UPLO, J1, M, NB, A, LDA, IPIV,
22* H, LDH, WORK )
23*
24* .. Scalar Arguments ..
25* CHARACTER UPLO
26* INTEGER J1, M, NB, LDA, LDH
27* ..
28* .. Array Arguments ..
29* INTEGER IPIV( * )
30* DOUBLE PRECISION A( LDA, * ), H( LDH, * ), WORK( * )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> DLATRF_AA factorizes a panel of a real symmetric matrix A using
40*> the Aasen's algorithm. The panel consists of a set of NB rows of A
41*> when UPLO is U, or a set of NB columns when UPLO is L.
42*>
43*> In order to factorize the panel, the Aasen's algorithm requires the
44*> last row, or column, of the previous panel. The first row, or column,
45*> of A is set to be the first row, or column, of an identity matrix,
46*> which is used to factorize the first panel.
47*>
48*> The resulting J-th row of U, or J-th column of L, is stored in the
49*> (J-1)-th row, or column, of A (without the unit diagonals), while
50*> the diagonal and subdiagonal of A are overwritten by those of T.
51*>
52*> \endverbatim
53*
54* Arguments:
55* ==========
56*
57*> \param[in] UPLO
58*> \verbatim
59*> UPLO is CHARACTER*1
60*> = 'U': Upper triangle of A is stored;
61*> = 'L': Lower triangle of A is stored.
62*> \endverbatim
63*>
64*> \param[in] J1
65*> \verbatim
66*> J1 is INTEGER
67*> The location of the first row, or column, of the panel
68*> within the submatrix of A, passed to this routine, e.g.,
69*> when called by DSYTRF_AA, for the first panel, J1 is 1,
70*> while for the remaining panels, J1 is 2.
71*> \endverbatim
72*>
73*> \param[in] M
74*> \verbatim
75*> M is INTEGER
76*> The dimension of the submatrix. M >= 0.
77*> \endverbatim
78*>
79*> \param[in] NB
80*> \verbatim
81*> NB is INTEGER
82*> The dimension of the panel to be facotorized.
83*> \endverbatim
84*>
85*> \param[in,out] A
86*> \verbatim
87*> A is DOUBLE PRECISION array, dimension (LDA,M) for
88*> the first panel, while dimension (LDA,M+1) for the
89*> remaining panels.
90*>
91*> On entry, A contains the last row, or column, of
92*> the previous panel, and the trailing submatrix of A
93*> to be factorized, except for the first panel, only
94*> the panel is passed.
95*>
96*> On exit, the leading panel is factorized.
97*> \endverbatim
98*>
99*> \param[in] LDA
100*> \verbatim
101*> LDA is INTEGER
102*> The leading dimension of the array A. LDA >= max(1,M).
103*> \endverbatim
104*>
105*> \param[out] IPIV
106*> \verbatim
107*> IPIV is INTEGER array, dimension (M)
108*> Details of the row and column interchanges,
109*> the row and column k were interchanged with the row and
110*> column IPIV(k).
111*> \endverbatim
112*>
113*> \param[in,out] H
114*> \verbatim
115*> H is DOUBLE PRECISION workspace, dimension (LDH,NB).
116*>
117*> \endverbatim
118*>
119*> \param[in] LDH
120*> \verbatim
121*> LDH is INTEGER
122*> The leading dimension of the workspace H. LDH >= max(1,M).
123*> \endverbatim
124*>
125*> \param[out] WORK
126*> \verbatim
127*> WORK is DOUBLE PRECISION workspace, dimension (M).
128*> \endverbatim
129*>
130*
131* Authors:
132* ========
133*
134*> \author Univ. of Tennessee
135*> \author Univ. of California Berkeley
136*> \author Univ. of Colorado Denver
137*> \author NAG Ltd.
138*
139*> \ingroup lahef_aa
140*
141* =====================================================================
142 SUBROUTINE dlasyf_aa( UPLO, J1, M, NB, A, LDA, IPIV,
143 $ H, LDH, WORK )
144*
145* -- LAPACK computational routine --
146* -- LAPACK is a software package provided by Univ. of Tennessee, --
147* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
148*
149 IMPLICIT NONE
150*
151* .. Scalar Arguments ..
152 CHARACTER UPLO
153 INTEGER M, NB, J1, LDA, LDH
154* ..
155* .. Array Arguments ..
156 INTEGER IPIV( * )
157 DOUBLE PRECISION A( LDA, * ), H( LDH, * ), WORK( * )
158* ..
159*
160* =====================================================================
161* .. Parameters ..
162 DOUBLE PRECISION ZERO, ONE
163 parameter( zero = 0.0d+0, one = 1.0d+0 )
164*
165* .. Local Scalars ..
166 INTEGER J, K, K1, I1, I2, MJ
167 DOUBLE PRECISION PIV, ALPHA
168* ..
169* .. External Functions ..
170 LOGICAL LSAME
171 INTEGER IDAMAX, ILAENV
172 EXTERNAL lsame, ilaenv, idamax
173* ..
174* .. External Subroutines ..
175 EXTERNAL dgemv, daxpy, dcopy, dswap, dscal, dlaset,
176 $ xerbla
177* ..
178* .. Intrinsic Functions ..
179 INTRINSIC max
180* ..
181* .. Executable Statements ..
182*
183 j = 1
184*
185* K1 is the first column of the panel to be factorized
186* i.e., K1 is 2 for the first block column, and 1 for the rest of the blocks
187*
188 k1 = (2-j1)+1
189*
190 IF( lsame( uplo, 'U' ) ) THEN
191*
192* .....................................................
193* Factorize A as U**T*D*U using the upper triangle of A
194* .....................................................
195*
196 10 CONTINUE
197 IF ( j.GT.min(m, nb) )
198 $ GO TO 20
199*
200* K is the column to be factorized
201* when being called from DSYTRF_AA,
202* > for the first block column, J1 is 1, hence J1+J-1 is J,
203* > for the rest of the columns, J1 is 2, and J1+J-1 is J+1,
204*
205 k = j1+j-1
206 IF( j.EQ.m ) THEN
207*
208* Only need to compute T(J, J)
209*
210 mj = 1
211 ELSE
212 mj = m-j+1
213 END IF
214*
215* H(J:M, J) := A(J, J:M) - H(J:M, 1:(J-1)) * L(J1:(J-1), J),
216* where H(J:M, J) has been initialized to be A(J, J:M)
217*
218 IF( k.GT.2 ) THEN
219*
220* K is the column to be factorized
221* > for the first block column, K is J, skipping the first two
222* columns
223* > for the rest of the columns, K is J+1, skipping only the
224* first column
225*
226 CALL dgemv( 'No transpose', mj, j-k1,
227 $ -one, h( j, k1 ), ldh,
228 $ a( 1, j ), 1,
229 $ one, h( j, j ), 1 )
230 END IF
231*
232* Copy H(i:M, i) into WORK
233*
234 CALL dcopy( mj, h( j, j ), 1, work( 1 ), 1 )
235*
236 IF( j.GT.k1 ) THEN
237*
238* Compute WORK := WORK - L(J-1, J:M) * T(J-1,J),
239* where A(J-1, J) stores T(J-1, J) and A(J-2, J:M) stores U(J-1, J:M)
240*
241 alpha = -a( k-1, j )
242 CALL daxpy( mj, alpha, a( k-2, j ), lda, work( 1 ), 1 )
243 END IF
244*
245* Set A(J, J) = T(J, J)
246*
247 a( k, j ) = work( 1 )
248*
249 IF( j.LT.m ) THEN
250*
251* Compute WORK(2:M) = T(J, J) L(J, (J+1):M)
252* where A(J, J) stores T(J, J) and A(J-1, (J+1):M) stores U(J, (J+1):M)
253*
254 IF( k.GT.1 ) THEN
255 alpha = -a( k, j )
256 CALL daxpy( m-j, alpha, a( k-1, j+1 ), lda,
257 $ work( 2 ), 1 )
258 ENDIF
259*
260* Find max(|WORK(2:M)|)
261*
262 i2 = idamax( m-j, work( 2 ), 1 ) + 1
263 piv = work( i2 )
264*
265* Apply symmetric pivot
266*
267 IF( (i2.NE.2) .AND. (piv.NE.0) ) THEN
268*
269* Swap WORK(I1) and WORK(I2)
270*
271 i1 = 2
272 work( i2 ) = work( i1 )
273 work( i1 ) = piv
274*
275* Swap A(I1, I1+1:M) with A(I1+1:M, I2)
276*
277 i1 = i1+j-1
278 i2 = i2+j-1
279 CALL dswap( i2-i1-1, a( j1+i1-1, i1+1 ), lda,
280 $ a( j1+i1, i2 ), 1 )
281*
282* Swap A(I1, I2+1:M) with A(I2, I2+1:M)
283*
284 IF( i2.LT.m )
285 $ CALL dswap( m-i2, a( j1+i1-1, i2+1 ), lda,
286 $ a( j1+i2-1, i2+1 ), lda )
287*
288* Swap A(I1, I1) with A(I2,I2)
289*
290 piv = a( i1+j1-1, i1 )
291 a( j1+i1-1, i1 ) = a( j1+i2-1, i2 )
292 a( j1+i2-1, i2 ) = piv
293*
294* Swap H(I1, 1:J1) with H(I2, 1:J1)
295*
296 CALL dswap( i1-1, h( i1, 1 ), ldh, h( i2, 1 ), ldh )
297 ipiv( i1 ) = i2
298*
299 IF( i1.GT.(k1-1) ) THEN
300*
301* Swap L(1:I1-1, I1) with L(1:I1-1, I2),
302* skipping the first column
303*
304 CALL dswap( i1-k1+1, a( 1, i1 ), 1,
305 $ a( 1, i2 ), 1 )
306 END IF
307 ELSE
308 ipiv( j+1 ) = j+1
309 ENDIF
310*
311* Set A(J, J+1) = T(J, J+1)
312*
313 a( k, j+1 ) = work( 2 )
314*
315 IF( j.LT.nb ) THEN
316*
317* Copy A(J+1:M, J+1) into H(J:M, J),
318*
319 CALL dcopy( m-j, a( k+1, j+1 ), lda,
320 $ h( j+1, j+1 ), 1 )
321 END IF
322*
323* Compute L(J+2, J+1) = WORK( 3:M ) / T(J, J+1),
324* where A(J, J+1) = T(J, J+1) and A(J+2:M, J) = L(J+2:M, J+1)
325*
326 IF( j.LT.(m-1) ) THEN
327 IF( a( k, j+1 ).NE.zero ) THEN
328 alpha = one / a( k, j+1 )
329 CALL dcopy( m-j-1, work( 3 ), 1, a( k, j+2 ), lda )
330 CALL dscal( m-j-1, alpha, a( k, j+2 ), lda )
331 ELSE
332 CALL dlaset( 'Full', 1, m-j-1, zero, zero,
333 $ a( k, j+2 ), lda)
334 END IF
335 END IF
336 END IF
337 j = j + 1
338 GO TO 10
339 20 CONTINUE
340*
341 ELSE
342*
343* .....................................................
344* Factorize A as L*D*L**T using the lower triangle of A
345* .....................................................
346*
347 30 CONTINUE
348 IF( j.GT.min( m, nb ) )
349 $ GO TO 40
350*
351* K is the column to be factorized
352* when being called from DSYTRF_AA,
353* > for the first block column, J1 is 1, hence J1+J-1 is J,
354* > for the rest of the columns, J1 is 2, and J1+J-1 is J+1,
355*
356 k = j1+j-1
357 IF( j.EQ.m ) THEN
358*
359* Only need to compute T(J, J)
360*
361 mj = 1
362 ELSE
363 mj = m-j+1
364 END IF
365*
366* H(J:M, J) := A(J:M, J) - H(J:M, 1:(J-1)) * L(J, J1:(J-1))^T,
367* where H(J:M, J) has been initialized to be A(J:M, J)
368*
369 IF( k.GT.2 ) THEN
370*
371* K is the column to be factorized
372* > for the first block column, K is J, skipping the first two
373* columns
374* > for the rest of the columns, K is J+1, skipping only the
375* first column
376*
377 CALL dgemv( 'No transpose', mj, j-k1,
378 $ -one, h( j, k1 ), ldh,
379 $ a( j, 1 ), lda,
380 $ one, h( j, j ), 1 )
381 END IF
382*
383* Copy H(J:M, J) into WORK
384*
385 CALL dcopy( mj, h( j, j ), 1, work( 1 ), 1 )
386*
387 IF( j.GT.k1 ) THEN
388*
389* Compute WORK := WORK - L(J:M, J-1) * T(J-1,J),
390* where A(J-1, J) = T(J-1, J) and A(J, J-2) = L(J, J-1)
391*
392 alpha = -a( j, k-1 )
393 CALL daxpy( mj, alpha, a( j, k-2 ), 1, work( 1 ), 1 )
394 END IF
395*
396* Set A(J, J) = T(J, J)
397*
398 a( j, k ) = work( 1 )
399*
400 IF( j.LT.m ) THEN
401*
402* Compute WORK(2:M) = T(J, J) L((J+1):M, J)
403* where A(J, J) = T(J, J) and A((J+1):M, J-1) = L((J+1):M, J)
404*
405 IF( k.GT.1 ) THEN
406 alpha = -a( j, k )
407 CALL daxpy( m-j, alpha, a( j+1, k-1 ), 1,
408 $ work( 2 ), 1 )
409 ENDIF
410*
411* Find max(|WORK(2:M)|)
412*
413 i2 = idamax( m-j, work( 2 ), 1 ) + 1
414 piv = work( i2 )
415*
416* Apply symmetric pivot
417*
418 IF( (i2.NE.2) .AND. (piv.NE.0) ) THEN
419*
420* Swap WORK(I1) and WORK(I2)
421*
422 i1 = 2
423 work( i2 ) = work( i1 )
424 work( i1 ) = piv
425*
426* Swap A(I1+1:M, I1) with A(I2, I1+1:M)
427*
428 i1 = i1+j-1
429 i2 = i2+j-1
430 CALL dswap( i2-i1-1, a( i1+1, j1+i1-1 ), 1,
431 $ a( i2, j1+i1 ), lda )
432*
433* Swap A(I2+1:M, I1) with A(I2+1:M, I2)
434*
435 IF( i2.LT.m )
436 $ CALL dswap( m-i2, a( i2+1, j1+i1-1 ), 1,
437 $ a( i2+1, j1+i2-1 ), 1 )
438*
439* Swap A(I1, I1) with A(I2, I2)
440*
441 piv = a( i1, j1+i1-1 )
442 a( i1, j1+i1-1 ) = a( i2, j1+i2-1 )
443 a( i2, j1+i2-1 ) = piv
444*
445* Swap H(I1, I1:J1) with H(I2, I2:J1)
446*
447 CALL dswap( i1-1, h( i1, 1 ), ldh, h( i2, 1 ), ldh )
448 ipiv( i1 ) = i2
449*
450 IF( i1.GT.(k1-1) ) THEN
451*
452* Swap L(1:I1-1, I1) with L(1:I1-1, I2),
453* skipping the first column
454*
455 CALL dswap( i1-k1+1, a( i1, 1 ), lda,
456 $ a( i2, 1 ), lda )
457 END IF
458 ELSE
459 ipiv( j+1 ) = j+1
460 ENDIF
461*
462* Set A(J+1, J) = T(J+1, J)
463*
464 a( j+1, k ) = work( 2 )
465*
466 IF( j.LT.nb ) THEN
467*
468* Copy A(J+1:M, J+1) into H(J+1:M, J),
469*
470 CALL dcopy( m-j, a( j+1, k+1 ), 1,
471 $ h( j+1, j+1 ), 1 )
472 END IF
473*
474* Compute L(J+2, J+1) = WORK( 3:M ) / T(J, J+1),
475* where A(J, J+1) = T(J, J+1) and A(J+2:M, J) = L(J+2:M, J+1)
476*
477 IF( j.LT.(m-1) ) THEN
478 IF( a( j+1, k ).NE.zero ) THEN
479 alpha = one / a( j+1, k )
480 CALL dcopy( m-j-1, work( 3 ), 1, a( j+2, k ), 1 )
481 CALL dscal( m-j-1, alpha, a( j+2, k ), 1 )
482 ELSE
483 CALL dlaset( 'Full', m-j-1, 1, zero, zero,
484 $ a( j+2, k ), lda )
485 END IF
486 END IF
487 END IF
488 j = j + 1
489 GO TO 30
490 40 CONTINUE
491 END IF
492 RETURN
493*
494* End of DLASYF_AA
495*
496 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
DGEMV
Definition dgemv.f:158
subroutine dlasyf_aa(uplo, j1, m, nb, a, lda, ipiv, h, ldh, work)
DLASYF_AA
Definition dlasyf_aa.f:144
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:110
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82