LAPACK 3.12.0
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 142 of file dlasyf_aa.f.

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*
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:162
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:110
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: