LAPACK  3.9.0
LAPACK: Linear Algebra PACKage

◆ clasyf_aa()

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

CLASYF_AA

Download CLASYF_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 CSYTRF_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 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 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 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 clasyf_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 A( LDA, * ), H( LDH, * ), WORK( * )
161 * ..
162 *
163 * =====================================================================
164 * .. Parameters ..
165  COMPLEX ZERO, ONE
166  parameter( zero = 0.0e+0, one = 1.0e+0 )
167 *
168 * .. Local Scalars ..
169  INTEGER J, K, K1, I1, I2, MJ
170  COMPLEX PIV, ALPHA
171 * ..
172 * .. External Functions ..
173  LOGICAL LSAME
174  INTEGER ICAMAX, ILAENV
175  EXTERNAL lsame, ilaenv, icamax
176 * ..
177 * .. External Subroutines ..
178  EXTERNAL caxpy, cgemv, cscal, ccopy, cswap, claset,
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 CSYTRF_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 cgemv( '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 ccopy( 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 caxpy( 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 caxpy( 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 = icamax( 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 cswap( 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 cswap( 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 cswap( 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 cswap( 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 ccopy( 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 ccopy( m-j-1, work( 3 ), 1, a( k, j+2 ), lda )
333  CALL cscal( m-j-1, alpha, a( k, j+2 ), lda )
334  ELSE
335  CALL claset( '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 CSYTRF_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 cgemv( '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 ccopy( 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 caxpy( 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 caxpy( 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 = icamax( 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 cswap( 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 cswap( 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 cswap( 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 cswap( 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 ccopy( 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 ccopy( m-j-1, work( 3 ), 1, a( j+2, k ), 1 )
484  CALL cscal( m-j-1, alpha, a( j+2, k ), 1 )
485  ELSE
486  CALL claset( '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 CLASYF_AA
498 *
Here is the call graph for this function:
Here is the caller graph for this function:
cgemv
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
Definition: cgemv.f:160
cscal
subroutine cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:80
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
cswap
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:83
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
claset
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: claset.f:108
ccopy
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:83
ilaenv
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: tstiee.f:83
caxpy
subroutine caxpy(N, CA, CX, INCX, CY, INCY)
CAXPY
Definition: caxpy.f:90
icamax
integer function icamax(N, CX, INCX)
ICAMAX
Definition: icamax.f:73