LAPACK  3.9.0
LAPACK: Linear Algebra PACKage

◆ zlasyf_aa()

subroutine zlasyf_aa ( character  UPLO,
integer  J1,
integer  M,
integer  NB,
complex*16, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
complex*16, dimension( ldh, * )  H,
integer  LDH,
complex*16, dimension( * )  WORK 
)

ZLASYF_AA

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

Purpose:
 DLATRF_AA factorizes a panel of a complex 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 ZSYTRF_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 COMPLEX*16 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 COMPLEX*16 workspace, dimension (LDH,NB).
[in]LDH
          LDH is INTEGER
          The leading dimension of the workspace H. LDH >= max(1,M).
[out]WORK
          WORK is COMPLEX*16 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 zlasyf_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  COMPLEX*16 A( LDA, * ), H( LDH, * ), WORK( * )
161 * ..
162 *
163 * =====================================================================
164 * .. Parameters ..
165  COMPLEX*16 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  COMPLEX*16 PIV, ALPHA
171 * ..
172 * .. External Functions ..
173  LOGICAL LSAME
174  INTEGER IZAMAX, ILAENV
175  EXTERNAL lsame, ilaenv, izamax
176 * ..
177 * .. External Subroutines ..
178  EXTERNAL zgemv, zaxpy, zscal, zcopy, zswap, zlaset,
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 ZSYTRF_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 zgemv( '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 zcopy( 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 zaxpy( 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 zaxpy( 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 = izamax( 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 zswap( 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 zswap( 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 zswap( 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 zswap( 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 zcopy( 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 zcopy( m-j-1, work( 3 ), 1, a( k, j+2 ), lda )
333  CALL zscal( m-j-1, alpha, a( k, j+2 ), lda )
334  ELSE
335  CALL zlaset( '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 ZSYTRF_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 zgemv( '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 zcopy( 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 zaxpy( 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 zaxpy( 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 = izamax( 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 zswap( 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 zswap( 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 zswap( 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 zswap( 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 zcopy( 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 zcopy( m-j-1, work( 3 ), 1, a( j+2, k ), 1 )
484  CALL zscal( m-j-1, alpha, a( j+2, k ), 1 )
485  ELSE
486  CALL zlaset( '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 ZLASYF_AA
498 *
Here is the call graph for this function:
Here is the caller graph for this function:
zaxpy
subroutine zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY
Definition: zaxpy.f:90
zgemv
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
Definition: zgemv.f:160
zcopy
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:83
izamax
integer function izamax(N, ZX, INCX)
IZAMAX
Definition: izamax.f:73
zlaset
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: zlaset.f:108
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
zswap
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:83
ilaenv
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: tstiee.f:83
zscal
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:80