LAPACK  3.9.0
LAPACK: Linear Algebra PACKage

◆ slasyf_aa()

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

SLASYF_AA

Download SLASYF_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 SSYTRF_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 REAL 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 REAL workspace, dimension (LDH,NB).
[in]LDH
          LDH is INTEGER
          The leading dimension of the workspace H. LDH >= max(1,M).
[out]WORK
          WORK is REAL 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 slasyf_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  REAL A( LDA, * ), H( LDH, * ), WORK( * )
161 * ..
162 *
163 * =====================================================================
164 * .. Parameters ..
165  REAL 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  REAL PIV, ALPHA
171 * ..
172 * .. External Functions ..
173  LOGICAL LSAME
174  INTEGER ISAMAX, ILAENV
175  EXTERNAL lsame, ilaenv, isamax
176 * ..
177 * .. External Subroutines ..
178  EXTERNAL saxpy, sgemv, sscal, scopy, sswap, slaset,
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 SSYTRF_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 sgemv( '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 scopy( 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 saxpy( 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 saxpy( 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 = isamax( 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 sswap( 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 sswap( 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 sswap( 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 sswap( 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 scopy( 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 scopy( m-j-1, work( 3 ), 1, a( k, j+2 ), lda )
333  CALL sscal( m-j-1, alpha, a( k, j+2 ), lda )
334  ELSE
335  CALL slaset( '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 SSYTRF_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 sgemv( '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 scopy( 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 saxpy( 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 saxpy( 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 = isamax( 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 sswap( 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 sswap( 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 sswap( 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 sswap( 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 scopy( 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 scopy( m-j-1, work( 3 ), 1, a( j+2, k ), 1 )
484  CALL sscal( m-j-1, alpha, a( j+2, k ), 1 )
485  ELSE
486  CALL slaset( '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 SLASYF_AA
498 *
Here is the call graph for this function:
Here is the caller graph for this function:
sswap
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:84
isamax
integer function isamax(N, SX, INCX)
ISAMAX
Definition: isamax.f:73
sgemv
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
Definition: sgemv.f:158
scopy
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:84
sscal
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:81
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
slaset
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: slaset.f:112
ilaenv
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: tstiee.f:83
saxpy
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
Definition: saxpy.f:91