LAPACK  3.9.0
LAPACK: Linear Algebra PACKage

◆ 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.
Date
November 2017

Definition at line 146 of file dlasyf_aa.f.

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