LAPACK  3.7.0
LAPACK: Linear Algebra PACKage
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,
integer  INFO 
)

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,N).
[out]IPIV
          IPIV is INTEGER array, dimension (N)
          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).
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 0:  if INFO = i, D(i,i) is exactly zero.  The factorization
                has been completed, but the block diagonal matrix D is
                exactly singular, and division by zero will occur if it
                is used to solve a system of equations.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
December 2016

Definition at line 156 of file slasyf_aa.f.

156 *
157 * -- LAPACK computational routine (version 3.7.0) --
158 * -- LAPACK is a software package provided by Univ. of Tennessee, --
159 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
160 * December 2016
161 *
162  IMPLICIT NONE
163 *
164 * .. Scalar Arguments ..
165  CHARACTER uplo
166  INTEGER m, nb, j1, lda, ldh, info
167 * ..
168 * .. Array Arguments ..
169  INTEGER ipiv( * )
170  REAL a( lda, * ), h( ldh, * ), work( * )
171 * ..
172 *
173 * =====================================================================
174 * .. Parameters ..
175  REAL zero, one
176  parameter ( zero = 0.0e+0, one = 1.0e+0 )
177 *
178 * .. Local Scalars ..
179  INTEGER j, k, k1, i1, i2
180  REAL piv, alpha
181 * ..
182 * .. External Functions ..
183  LOGICAL lsame
184  INTEGER isamax, ilaenv
185  EXTERNAL lsame, ilaenv, isamax
186 * ..
187 * .. External Subroutines ..
188  EXTERNAL xerbla
189 * ..
190 * .. Intrinsic Functions ..
191  INTRINSIC max
192 * ..
193 * .. Executable Statements ..
194 *
195  info = 0
196  j = 1
197 *
198 * K1 is the first column of the panel to be factorized
199 * i.e., K1 is 2 for the first block column, and 1 for the rest of the blocks
200 *
201  k1 = (2-j1)+1
202 *
203  IF( lsame( uplo, 'U' ) ) THEN
204 *
205 * .....................................................
206 * Factorize A as U**T*D*U using the upper triangle of A
207 * .....................................................
208 *
209  10 CONTINUE
210  IF ( j.GT.min(m, nb) )
211  $ GO TO 20
212 *
213 * K is the column to be factorized
214 * when being called from SSYTRF_AA,
215 * > for the first block column, J1 is 1, hence J1+J-1 is J,
216 * > for the rest of the columns, J1 is 2, and J1+J-1 is J+1,
217 *
218  k = j1+j-1
219 *
220 * H(J:N, J) := A(J, J:N) - H(J:N, 1:(J-1)) * L(J1:(J-1), J),
221 * where H(J:N, J) has been initialized to be A(J, J:N)
222 *
223  IF( k.GT.2 ) THEN
224 *
225 * K is the column to be factorized
226 * > for the first block column, K is J, skipping the first two
227 * columns
228 * > for the rest of the columns, K is J+1, skipping only the
229 * first column
230 *
231  CALL sgemv( 'No transpose', m-j+1, j-k1,
232  $ -one, h( j, k1 ), ldh,
233  $ a( 1, j ), 1,
234  $ one, h( j, j ), 1 )
235  END IF
236 *
237 * Copy H(i:n, i) into WORK
238 *
239  CALL scopy( m-j+1, h( j, j ), 1, work( 1 ), 1 )
240 *
241  IF( j.GT.k1 ) THEN
242 *
243 * Compute WORK := WORK - L(J-1, J:N) * T(J-1,J),
244 * where A(J-1, J) stores T(J-1, J) and A(J-2, J:N) stores U(J-1, J:N)
245 *
246  alpha = -a( k-1, j )
247  CALL saxpy( m-j+1, alpha, a( k-2, j ), lda, work( 1 ), 1 )
248  END IF
249 *
250 * Set A(J, J) = T(J, J)
251 *
252  a( k, j ) = work( 1 )
253 *
254  IF( j.LT.m ) THEN
255 *
256 * Compute WORK(2:N) = T(J, J) L(J, (J+1):N)
257 * where A(J, J) stores T(J, J) and A(J-1, (J+1):N) stores U(J, (J+1):N)
258 *
259  IF( k.GT.1 ) THEN
260  alpha = -a( k, j )
261  CALL saxpy( m-j, alpha, a( k-1, j+1 ), lda,
262  $ work( 2 ), 1 )
263  ENDIF
264 *
265 * Find max(|WORK(2:n)|)
266 *
267  i2 = isamax( m-j, work( 2 ), 1 ) + 1
268  piv = work( i2 )
269 *
270 * Apply symmetric pivot
271 *
272  IF( (i2.NE.2) .AND. (piv.NE.0) ) THEN
273 *
274 * Swap WORK(I1) and WORK(I2)
275 *
276  i1 = 2
277  work( i2 ) = work( i1 )
278  work( i1 ) = piv
279 *
280 * Swap A(I1, I1+1:N) with A(I1+1:N, I2)
281 *
282  i1 = i1+j-1
283  i2 = i2+j-1
284  CALL sswap( i2-i1-1, a( j1+i1-1, i1+1 ), lda,
285  $ a( j1+i1, i2 ), 1 )
286 *
287 * Swap A(I1, I2+1:N) with A(I2, I2+1:N)
288 *
289  CALL sswap( m-i2, a( j1+i1-1, i2+1 ), lda,
290  $ a( j1+i2-1, i2+1 ), lda )
291 *
292 * Swap A(I1, I1) with A(I2,I2)
293 *
294  piv = a( i1+j1-1, i1 )
295  a( j1+i1-1, i1 ) = a( j1+i2-1, i2 )
296  a( j1+i2-1, i2 ) = piv
297 *
298 * Swap H(I1, 1:J1) with H(I2, 1:J1)
299 *
300  CALL sswap( i1-1, h( i1, 1 ), ldh, h( i2, 1 ), ldh )
301  ipiv( i1 ) = i2
302 *
303  IF( i1.GT.(k1-1) ) THEN
304 *
305 * Swap L(1:I1-1, I1) with L(1:I1-1, I2),
306 * skipping the first column
307 *
308  CALL sswap( i1-k1+1, a( 1, i1 ), 1,
309  $ a( 1, i2 ), 1 )
310  END IF
311  ELSE
312  ipiv( j+1 ) = j+1
313  ENDIF
314 *
315 * Set A(J, J+1) = T(J, J+1)
316 *
317  a( k, j+1 ) = work( 2 )
318  IF( (a( k, j ).EQ.zero ) .AND.
319  $ ( (j.EQ.m) .OR. (a( k, j+1 ).EQ.zero))) THEN
320  IF(info .EQ. 0) THEN
321  info = j
322  ENDIF
323  END IF
324 *
325  IF( j.LT.nb ) THEN
326 *
327 * Copy A(J+1:N, J+1) into H(J:N, J),
328 *
329  CALL scopy( m-j, a( k+1, j+1 ), lda,
330  $ h( j+1, j+1 ), 1 )
331  END IF
332 *
333 * Compute L(J+2, J+1) = WORK( 3:N ) / T(J, J+1),
334 * where A(J, J+1) = T(J, J+1) and A(J+2:N, J) = L(J+2:N, J+1)
335 *
336  IF( a( k, j+1 ).NE.zero ) THEN
337  alpha = one / a( k, j+1 )
338  CALL scopy( m-j-1, work( 3 ), 1, a( k, j+2 ), lda )
339  CALL sscal( m-j-1, alpha, a( k, j+2 ), lda )
340  ELSE
341  CALL slaset( 'Full', 1, m-j-1, zero, zero,
342  $ a( k, j+2 ), lda)
343  END IF
344  ELSE
345  IF( (a( k, j ).EQ.zero) .AND. (info.EQ.0) ) THEN
346  info = j
347  END IF
348  END IF
349  j = j + 1
350  GO TO 10
351  20 CONTINUE
352 *
353  ELSE
354 *
355 * .....................................................
356 * Factorize A as L*D*L**T using the lower triangle of A
357 * .....................................................
358 *
359  30 CONTINUE
360  IF( j.GT.min( m, nb ) )
361  $ GO TO 40
362 *
363 * K is the column to be factorized
364 * when being called from SSYTRF_AA,
365 * > for the first block column, J1 is 1, hence J1+J-1 is J,
366 * > for the rest of the columns, J1 is 2, and J1+J-1 is J+1,
367 *
368  k = j1+j-1
369 *
370 * H(J:N, J) := A(J:N, J) - H(J:N, 1:(J-1)) * L(J, J1:(J-1))^T,
371 * where H(J:N, J) has been initialized to be A(J:N, J)
372 *
373  IF( k.GT.2 ) THEN
374 *
375 * K is the column to be factorized
376 * > for the first block column, K is J, skipping the first two
377 * columns
378 * > for the rest of the columns, K is J+1, skipping only the
379 * first column
380 *
381  CALL sgemv( 'No transpose', m-j+1, j-k1,
382  $ -one, h( j, k1 ), ldh,
383  $ a( j, 1 ), lda,
384  $ one, h( j, j ), 1 )
385  END IF
386 *
387 * Copy H(J:N, J) into WORK
388 *
389  CALL scopy( m-j+1, h( j, j ), 1, work( 1 ), 1 )
390 *
391  IF( j.GT.k1 ) THEN
392 *
393 * Compute WORK := WORK - L(J:N, J-1) * T(J-1,J),
394 * where A(J-1, J) = T(J-1, J) and A(J, J-2) = L(J, J-1)
395 *
396  alpha = -a( j, k-1 )
397  CALL saxpy( m-j+1, alpha, a( j, k-2 ), 1, work( 1 ), 1 )
398  END IF
399 *
400 * Set A(J, J) = T(J, J)
401 *
402  a( j, k ) = work( 1 )
403 *
404  IF( j.LT.m ) THEN
405 *
406 * Compute WORK(2:N) = T(J, J) L((J+1):N, J)
407 * where A(J, J) = T(J, J) and A((J+1):N, J-1) = L((J+1):N, J)
408 *
409  IF( k.GT.1 ) THEN
410  alpha = -a( j, k )
411  CALL saxpy( m-j, alpha, a( j+1, k-1 ), 1,
412  $ work( 2 ), 1 )
413  ENDIF
414 *
415 * Find max(|WORK(2:n)|)
416 *
417  i2 = isamax( m-j, work( 2 ), 1 ) + 1
418  piv = work( i2 )
419 *
420 * Apply symmetric pivot
421 *
422  IF( (i2.NE.2) .AND. (piv.NE.0) ) THEN
423 *
424 * Swap WORK(I1) and WORK(I2)
425 *
426  i1 = 2
427  work( i2 ) = work( i1 )
428  work( i1 ) = piv
429 *
430 * Swap A(I1+1:N, I1) with A(I2, I1+1:N)
431 *
432  i1 = i1+j-1
433  i2 = i2+j-1
434  CALL sswap( i2-i1-1, a( i1+1, j1+i1-1 ), 1,
435  $ a( i2, j1+i1 ), lda )
436 *
437 * Swap A(I2+1:N, I1) with A(I2+1:N, I2)
438 *
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  IF( (a( j, k ).EQ.zero) .AND.
469  $ ( (j.EQ.m) .OR. (a( j+1, k ).EQ.zero)) ) THEN
470  IF (info .EQ. 0)
471  $ info = j
472  END IF
473 *
474  IF( j.LT.nb ) THEN
475 *
476 * Copy A(J+1:N, J+1) into H(J+1:N, J),
477 *
478  CALL scopy( m-j, a( j+1, k+1 ), 1,
479  $ h( j+1, j+1 ), 1 )
480  END IF
481 *
482 * Compute L(J+2, J+1) = WORK( 3:N ) / T(J, J+1),
483 * where A(J, J+1) = T(J, J+1) and A(J+2:N, J) = L(J+2:N, J+1)
484 *
485  IF( a( j+1, k ).NE.zero ) THEN
486  alpha = one / a( j+1, k )
487  CALL scopy( m-j-1, work( 3 ), 1, a( j+2, k ), 1 )
488  CALL sscal( m-j-1, alpha, a( j+2, k ), 1 )
489  ELSE
490  CALL slaset( 'Full', m-j-1, 1, zero, zero,
491  $ a( j+2, k ), lda )
492  END IF
493  ELSE
494  IF( (a( j, k ).EQ.zero) .AND. (info.EQ.0) ) THEN
495  info = j
496  END IF
497  END IF
498  j = j + 1
499  GO TO 30
500  40 CONTINUE
501  END IF
502  RETURN
503 *
504 * End of SLASYF_AA
505 *
ilaenv
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: tstiee.f:83
isamax
integer function isamax(N, SX, INCX)
ISAMAX
Definition: isamax.f:53
sgemv
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
Definition: sgemv.f:158
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
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
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
saxpy
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
Definition: saxpy.f:54
sscal
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:55
sswap
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:53
scopy
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:53

Here is the call graph for this function:

Here is the caller graph for this function: