LAPACK  3.10.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.

Definition at line 142 of file zlasyf_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  COMPLEX*16 A( LDA, * ), H( LDH, * ), WORK( * )
158 * ..
159 *
160 * =====================================================================
161 * .. Parameters ..
162  COMPLEX*16 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  COMPLEX*16 PIV, ALPHA
168 * ..
169 * .. External Functions ..
170  LOGICAL LSAME
171  INTEGER IZAMAX, ILAENV
172  EXTERNAL lsame, ilaenv, izamax
173 * ..
174 * .. External Subroutines ..
175  EXTERNAL zgemv, zaxpy, zscal, zcopy, zswap, zlaset,
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 ZSYTRF_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 zgemv( '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 zcopy( 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 zaxpy( 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 zaxpy( 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 = izamax( 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 zswap( 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 zswap( 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 zswap( 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 zswap( 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 zcopy( 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 zcopy( m-j-1, work( 3 ), 1, a( k, j+2 ), lda )
330  CALL zscal( m-j-1, alpha, a( k, j+2 ), lda )
331  ELSE
332  CALL zlaset( '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 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  IF( i2.LT.m )
436  $ CALL zswap( 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 zswap( 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 zswap( 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 zcopy( 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 zcopy( m-j-1, work( 3 ), 1, a( j+2, k ), 1 )
481  CALL zscal( m-j-1, alpha, a( j+2, k ), 1 )
482  ELSE
483  CALL zlaset( '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 ZLASYF_AA
495 *
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: ilaenv.f:162
integer function izamax(N, ZX, INCX)
IZAMAX
Definition: izamax.f:71
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:81
subroutine zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY
Definition: zaxpy.f:88
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:78
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:81
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
Definition: zgemv.f:158
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:106
Here is the call graph for this function:
Here is the caller graph for this function: