LAPACK  3.8.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  CALL zswap( m-i2, a( j1+i1-1, i2+1 ), lda,
288  $ a( j1+i2-1, i2+1 ), lda )
289 *
290 * Swap A(I1, I1) with A(I2,I2)
291 *
292  piv = a( i1+j1-1, i1 )
293  a( j1+i1-1, i1 ) = a( j1+i2-1, i2 )
294  a( j1+i2-1, i2 ) = piv
295 *
296 * Swap H(I1, 1:J1) with H(I2, 1:J1)
297 *
298  CALL zswap( i1-1, h( i1, 1 ), ldh, h( i2, 1 ), ldh )
299  ipiv( i1 ) = i2
300 *
301  IF( i1.GT.(k1-1) ) THEN
302 *
303 * Swap L(1:I1-1, I1) with L(1:I1-1, I2),
304 * skipping the first column
305 *
306  CALL zswap( i1-k1+1, a( 1, i1 ), 1,
307  $ a( 1, i2 ), 1 )
308  END IF
309  ELSE
310  ipiv( j+1 ) = j+1
311  ENDIF
312 *
313 * Set A(J, J+1) = T(J, J+1)
314 *
315  a( k, j+1 ) = work( 2 )
316 *
317  IF( j.LT.nb ) THEN
318 *
319 * Copy A(J+1:M, J+1) into H(J:M, J),
320 *
321  CALL zcopy( m-j, a( k+1, j+1 ), lda,
322  $ h( j+1, j+1 ), 1 )
323  END IF
324 *
325 * Compute L(J+2, J+1) = WORK( 3:M ) / T(J, J+1),
326 * where A(J, J+1) = T(J, J+1) and A(J+2:M, J) = L(J+2:M, J+1)
327 *
328  IF( a( k, j+1 ).NE.zero ) THEN
329  alpha = one / a( k, j+1 )
330  CALL zcopy( m-j-1, work( 3 ), 1, a( k, j+2 ), lda )
331  CALL zscal( m-j-1, alpha, a( k, j+2 ), lda )
332  ELSE
333  CALL zlaset( 'Full', 1, m-j-1, zero, zero,
334  $ a( k, j+2 ), lda)
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 ZSYTRF_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 zgemv( '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 zcopy( 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 zaxpy( 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 zaxpy( 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 = izamax( 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 zswap( 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  CALL zswap( 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 zswap( 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 zswap( 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 zcopy( 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( a( j+1, k ).NE.zero ) THEN
477  alpha = one / a( j+1, k )
478  CALL zcopy( m-j-1, work( 3 ), 1, a( j+2, k ), 1 )
479  CALL zscal( m-j-1, alpha, a( j+2, k ), 1 )
480  ELSE
481  CALL zlaset( 'Full', m-j-1, 1, zero, zero,
482  $ a( j+2, k ), lda )
483  END IF
484  END IF
485  j = j + 1
486  GO TO 30
487  40 CONTINUE
488  END IF
489  RETURN
490 *
491 * End of ZLASYF_AA
492 *
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:83
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
Definition: zgemv.f:160
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:83
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: tstiee.f:83
integer function izamax(N, ZX, INCX)
IZAMAX
Definition: izamax.f:73
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
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
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY
Definition: zaxpy.f:90
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:80
Here is the call graph for this function:
Here is the caller graph for this function: