LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dlasyf_aa()

subroutine dlasyf_aa ( character uplo,
integer j1,
integer m,
integer nb,
double precision, dimension( lda, * ) a,
integer lda,
integer, dimension( * ) ipiv,
double precision, dimension( ldh, * ) h,
integer ldh,
double precision, dimension( * ) work )

DLASYF_AA

Download DLASYF_AA + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> DLATRF_AA factorizes a panel of a real symmetric matrix A using
!> the Aasen's algorithm. The panel consists of a set of NB rows of A
!> when UPLO is U, or a set of NB columns when UPLO is L.
!>
!> In order to factorize the panel, the Aasen's algorithm requires the
!> last row, or column, of the previous panel. The first row, or column,
!> of A is set to be the first row, or column, of an identity matrix,
!> which is used to factorize the first panel.
!>
!> The resulting J-th row of U, or J-th column of L, is stored in the
!> (J-1)-th row, or column, of A (without the unit diagonals), while
!> the diagonal and subdiagonal of A are overwritten by those of T.
!>
!>
Parameters
[in]UPLO
!>          UPLO is CHARACTER*1
!>          = 'U':  Upper triangle of A is stored;
!>          = 'L':  Lower triangle of A is stored.
!>
[in]J1
!>          J1 is INTEGER
!>          The location of the first row, or column, of the panel
!>          within the submatrix of A, passed to this routine, e.g.,
!>          when called by DSYTRF_AA, for the first panel, J1 is 1,
!>          while for the remaining panels, J1 is 2.
!>
[in]M
!>          M is INTEGER
!>          The dimension of the submatrix. M >= 0.
!>
[in]NB
!>          NB is INTEGER
!>          The dimension of the panel to be facotorized.
!>
[in,out]A
!>          A is DOUBLE PRECISION array, dimension (LDA,M) for
!>          the first panel, while dimension (LDA,M+1) for the
!>          remaining panels.
!>
!>          On entry, A contains the last row, or column, of
!>          the previous panel, and the trailing submatrix of A
!>          to be factorized, except for the first panel, only
!>          the panel is passed.
!>
!>          On exit, the leading panel is factorized.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,M).
!>
[out]IPIV
!>          IPIV is INTEGER array, dimension (M)
!>          Details of the row and column interchanges,
!>          the row and column k were interchanged with the row and
!>          column IPIV(k).
!>
[in,out]H
!>          H is DOUBLE PRECISION workspace, dimension (LDH,NB).
!>
!>
[in]LDH
!>          LDH is INTEGER
!>          The leading dimension of the workspace H. LDH >= max(1,M).
!>
[out]WORK
!>          WORK is DOUBLE PRECISION workspace, dimension (M).
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 140 of file dlasyf_aa.f.

142*
143* -- LAPACK computational routine --
144* -- LAPACK is a software package provided by Univ. of Tennessee, --
145* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
146*
147 IMPLICIT NONE
148*
149* .. Scalar Arguments ..
150 CHARACTER UPLO
151 INTEGER M, NB, J1, LDA, LDH
152* ..
153* .. Array Arguments ..
154 INTEGER IPIV( * )
155 DOUBLE PRECISION A( LDA, * ), H( LDH, * ), WORK( * )
156* ..
157*
158* =====================================================================
159* .. Parameters ..
160 DOUBLE PRECISION ZERO, ONE
161 parameter( zero = 0.0d+0, one = 1.0d+0 )
162*
163* .. Local Scalars ..
164 INTEGER J, K, K1, I1, I2, MJ
165 DOUBLE PRECISION PIV, ALPHA
166* ..
167* .. External Functions ..
168 LOGICAL LSAME
169 INTEGER IDAMAX, ILAENV
170 EXTERNAL lsame, ilaenv, idamax
171* ..
172* .. External Subroutines ..
173 EXTERNAL dgemv, daxpy, dcopy, dswap, dscal,
174 $ dlaset,
175 $ xerbla
176* ..
177* .. Intrinsic Functions ..
178 INTRINSIC max
179* ..
180* .. Executable Statements ..
181*
182 j = 1
183*
184* K1 is the first column of the panel to be factorized
185* i.e., K1 is 2 for the first block column, and 1 for the rest of the blocks
186*
187 k1 = (2-j1)+1
188*
189 IF( lsame( uplo, 'U' ) ) THEN
190*
191* .....................................................
192* Factorize A as U**T*D*U using the upper triangle of A
193* .....................................................
194*
195 10 CONTINUE
196 IF ( j.GT.min(m, nb) )
197 $ GO TO 20
198*
199* K is the column to be factorized
200* when being called from DSYTRF_AA,
201* > for the first block column, J1 is 1, hence J1+J-1 is J,
202* > for the rest of the columns, J1 is 2, and J1+J-1 is J+1,
203*
204 k = j1+j-1
205 IF( j.EQ.m ) THEN
206*
207* Only need to compute T(J, J)
208*
209 mj = 1
210 ELSE
211 mj = m-j+1
212 END IF
213*
214* H(J:M, J) := A(J, J:M) - H(J:M, 1:(J-1)) * L(J1:(J-1), J),
215* where H(J:M, J) has been initialized to be A(J, J:M)
216*
217 IF( k.GT.2 ) THEN
218*
219* K is the column to be factorized
220* > for the first block column, K is J, skipping the first two
221* columns
222* > for the rest of the columns, K is J+1, skipping only the
223* first column
224*
225 CALL dgemv( 'No transpose', mj, j-k1,
226 $ -one, h( j, k1 ), ldh,
227 $ a( 1, j ), 1,
228 $ one, h( j, j ), 1 )
229 END IF
230*
231* Copy H(i:M, i) into WORK
232*
233 CALL dcopy( mj, h( j, j ), 1, work( 1 ), 1 )
234*
235 IF( j.GT.k1 ) THEN
236*
237* Compute WORK := WORK - L(J-1, J:M) * T(J-1,J),
238* where A(J-1, J) stores T(J-1, J) and A(J-2, J:M) stores U(J-1, J:M)
239*
240 alpha = -a( k-1, j )
241 CALL daxpy( mj, alpha, a( k-2, j ), lda, work( 1 ), 1 )
242 END IF
243*
244* Set A(J, J) = T(J, J)
245*
246 a( k, j ) = work( 1 )
247*
248 IF( j.LT.m ) THEN
249*
250* Compute WORK(2:M) = T(J, J) L(J, (J+1):M)
251* where A(J, J) stores T(J, J) and A(J-1, (J+1):M) stores U(J, (J+1):M)
252*
253 IF( k.GT.1 ) THEN
254 alpha = -a( k, j )
255 CALL daxpy( m-j, alpha, a( k-1, j+1 ), lda,
256 $ work( 2 ), 1 )
257 ENDIF
258*
259* Find max(|WORK(2:M)|)
260*
261 i2 = idamax( m-j, work( 2 ), 1 ) + 1
262 piv = work( i2 )
263*
264* Apply symmetric pivot
265*
266 IF( (i2.NE.2) .AND. (piv.NE.0) ) THEN
267*
268* Swap WORK(I1) and WORK(I2)
269*
270 i1 = 2
271 work( i2 ) = work( i1 )
272 work( i1 ) = piv
273*
274* Swap A(I1, I1+1:M) with A(I1+1:M, I2)
275*
276 i1 = i1+j-1
277 i2 = i2+j-1
278 CALL dswap( i2-i1-1, a( j1+i1-1, i1+1 ), lda,
279 $ a( j1+i1, i2 ), 1 )
280*
281* Swap A(I1, I2+1:M) with A(I2, I2+1:M)
282*
283 IF( i2.LT.m )
284 $ CALL dswap( m-i2, a( j1+i1-1, i2+1 ), lda,
285 $ a( j1+i2-1, i2+1 ), lda )
286*
287* Swap A(I1, I1) with A(I2,I2)
288*
289 piv = a( i1+j1-1, i1 )
290 a( j1+i1-1, i1 ) = a( j1+i2-1, i2 )
291 a( j1+i2-1, i2 ) = piv
292*
293* Swap H(I1, 1:J1) with H(I2, 1:J1)
294*
295 CALL dswap( i1-1, h( i1, 1 ), ldh, h( i2, 1 ), ldh )
296 ipiv( i1 ) = i2
297*
298 IF( i1.GT.(k1-1) ) THEN
299*
300* Swap L(1:I1-1, I1) with L(1:I1-1, I2),
301* skipping the first column
302*
303 CALL dswap( i1-k1+1, a( 1, i1 ), 1,
304 $ a( 1, i2 ), 1 )
305 END IF
306 ELSE
307 ipiv( j+1 ) = j+1
308 ENDIF
309*
310* Set A(J, J+1) = T(J, J+1)
311*
312 a( k, j+1 ) = work( 2 )
313*
314 IF( j.LT.nb ) THEN
315*
316* Copy A(J+1:M, J+1) into H(J:M, J),
317*
318 CALL dcopy( m-j, a( k+1, j+1 ), lda,
319 $ h( j+1, j+1 ), 1 )
320 END IF
321*
322* Compute L(J+2, J+1) = WORK( 3:M ) / T(J, J+1),
323* where A(J, J+1) = T(J, J+1) and A(J+2:M, J) = L(J+2:M, J+1)
324*
325 IF( j.LT.(m-1) ) THEN
326 IF( a( k, j+1 ).NE.zero ) THEN
327 alpha = one / a( k, j+1 )
328 CALL dcopy( m-j-1, work( 3 ), 1, a( k, j+2 ), lda )
329 CALL dscal( m-j-1, alpha, a( k, j+2 ), lda )
330 ELSE
331 CALL dlaset( 'Full', 1, m-j-1, zero, zero,
332 $ a( k, j+2 ), lda)
333 END IF
334 END IF
335 END IF
336 j = j + 1
337 GO TO 10
338 20 CONTINUE
339*
340 ELSE
341*
342* .....................................................
343* Factorize A as L*D*L**T using the lower triangle of A
344* .....................................................
345*
346 30 CONTINUE
347 IF( j.GT.min( m, nb ) )
348 $ GO TO 40
349*
350* K is the column to be factorized
351* when being called from DSYTRF_AA,
352* > for the first block column, J1 is 1, hence J1+J-1 is J,
353* > for the rest of the columns, J1 is 2, and J1+J-1 is J+1,
354*
355 k = j1+j-1
356 IF( j.EQ.m ) THEN
357*
358* Only need to compute T(J, J)
359*
360 mj = 1
361 ELSE
362 mj = m-j+1
363 END IF
364*
365* H(J:M, J) := A(J:M, J) - H(J:M, 1:(J-1)) * L(J, J1:(J-1))^T,
366* where H(J:M, J) has been initialized to be A(J:M, J)
367*
368 IF( k.GT.2 ) THEN
369*
370* K is the column to be factorized
371* > for the first block column, K is J, skipping the first two
372* columns
373* > for the rest of the columns, K is J+1, skipping only the
374* first column
375*
376 CALL dgemv( 'No transpose', mj, j-k1,
377 $ -one, h( j, k1 ), ldh,
378 $ a( j, 1 ), lda,
379 $ one, h( j, j ), 1 )
380 END IF
381*
382* Copy H(J:M, J) into WORK
383*
384 CALL dcopy( mj, h( j, j ), 1, work( 1 ), 1 )
385*
386 IF( j.GT.k1 ) THEN
387*
388* Compute WORK := WORK - L(J:M, J-1) * T(J-1,J),
389* where A(J-1, J) = T(J-1, J) and A(J, J-2) = L(J, J-1)
390*
391 alpha = -a( j, k-1 )
392 CALL daxpy( mj, alpha, a( j, k-2 ), 1, work( 1 ), 1 )
393 END IF
394*
395* Set A(J, J) = T(J, J)
396*
397 a( j, k ) = work( 1 )
398*
399 IF( j.LT.m ) THEN
400*
401* Compute WORK(2:M) = T(J, J) L((J+1):M, J)
402* where A(J, J) = T(J, J) and A((J+1):M, J-1) = L((J+1):M, J)
403*
404 IF( k.GT.1 ) THEN
405 alpha = -a( j, k )
406 CALL daxpy( m-j, alpha, a( j+1, k-1 ), 1,
407 $ work( 2 ), 1 )
408 ENDIF
409*
410* Find max(|WORK(2:M)|)
411*
412 i2 = idamax( m-j, work( 2 ), 1 ) + 1
413 piv = work( i2 )
414*
415* Apply symmetric pivot
416*
417 IF( (i2.NE.2) .AND. (piv.NE.0) ) THEN
418*
419* Swap WORK(I1) and WORK(I2)
420*
421 i1 = 2
422 work( i2 ) = work( i1 )
423 work( i1 ) = piv
424*
425* Swap A(I1+1:M, I1) with A(I2, I1+1:M)
426*
427 i1 = i1+j-1
428 i2 = i2+j-1
429 CALL dswap( i2-i1-1, a( i1+1, j1+i1-1 ), 1,
430 $ a( i2, j1+i1 ), lda )
431*
432* Swap A(I2+1:M, I1) with A(I2+1:M, I2)
433*
434 IF( i2.LT.m )
435 $ CALL dswap( m-i2, a( i2+1, j1+i1-1 ), 1,
436 $ a( i2+1, j1+i2-1 ), 1 )
437*
438* Swap A(I1, I1) with A(I2, I2)
439*
440 piv = a( i1, j1+i1-1 )
441 a( i1, j1+i1-1 ) = a( i2, j1+i2-1 )
442 a( i2, j1+i2-1 ) = piv
443*
444* Swap H(I1, I1:J1) with H(I2, I2:J1)
445*
446 CALL dswap( i1-1, h( i1, 1 ), ldh, h( i2, 1 ), ldh )
447 ipiv( i1 ) = i2
448*
449 IF( i1.GT.(k1-1) ) THEN
450*
451* Swap L(1:I1-1, I1) with L(1:I1-1, I2),
452* skipping the first column
453*
454 CALL dswap( i1-k1+1, a( i1, 1 ), lda,
455 $ a( i2, 1 ), lda )
456 END IF
457 ELSE
458 ipiv( j+1 ) = j+1
459 ENDIF
460*
461* Set A(J+1, J) = T(J+1, J)
462*
463 a( j+1, k ) = work( 2 )
464*
465 IF( j.LT.nb ) THEN
466*
467* Copy A(J+1:M, J+1) into H(J+1:M, J),
468*
469 CALL dcopy( m-j, a( j+1, k+1 ), 1,
470 $ h( j+1, j+1 ), 1 )
471 END IF
472*
473* Compute L(J+2, J+1) = WORK( 3:M ) / T(J, J+1),
474* where A(J, J+1) = T(J, J+1) and A(J+2:M, J) = L(J+2:M, J+1)
475*
476 IF( j.LT.(m-1) ) THEN
477 IF( a( j+1, k ).NE.zero ) THEN
478 alpha = one / a( j+1, k )
479 CALL dcopy( m-j-1, work( 3 ), 1, a( j+2, k ), 1 )
480 CALL dscal( m-j-1, alpha, a( j+2, k ), 1 )
481 ELSE
482 CALL dlaset( 'Full', m-j-1, 1, zero, zero,
483 $ a( j+2, k ), lda )
484 END IF
485 END IF
486 END IF
487 j = j + 1
488 GO TO 30
489 40 CONTINUE
490 END IF
491 RETURN
492*
493* End of DLASYF_AA
494*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
daxpy
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
dgemv
subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
DGEMV
Definition dgemv.f:158
idamax
integer function idamax(n, dx, incx)
IDAMAX
Definition idamax.f:71
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:160
dlaset
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:108
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
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
Here is the call graph for this function:
Here is the caller graph for this function: