LAPACK  3.7.0
LAPACK: Linear Algebra PACKage
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,
integer  INFO 
)

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,N).
[out]IPIV
          IPIV is INTEGER array, dimension (N)
          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).
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 0:  if INFO = i, D(i,i) is exactly zero.  The factorization
                has been completed, but the block diagonal matrix D is
                exactly singular, and division by zero will occur if it
                is used to solve a system of equations.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
December 2016

Definition at line 156 of file dlasyf_aa.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: