LAPACK  3.7.0
LAPACK: Linear Algebra PACKage
subroutine ssytrs_3 ( character  UPLO,
integer  N,
integer  NRHS,
real, dimension( lda, * )  A,
integer  LDA,
real, dimension( * )  E,
integer, dimension( * )  IPIV,
real, dimension( ldb, * )  B,
integer  LDB,
integer  INFO 
)

SSYTRS_3

Download SSYTRS_3 + dependencies [TGZ] [ZIP] [TXT]

Purpose: 
 SSYTRS_3 solves a system of linear equations A * X = B with a real
 symmetric matrix A using the factorization computed
 by SSYTRF_RK or SSYTRF_BK:

    A = P*U*D*(U**T)*(P**T) or A = P*L*D*(L**T)*(P**T),

 where U (or L) is unit upper (or lower) triangular matrix,
 U**T (or L**T) is the transpose of U (or L), P is a permutation
 matrix, P**T is the transpose of P, and D is symmetric and block
 diagonal with 1-by-1 and 2-by-2 diagonal blocks.

 This algorithm is using Level 3 BLAS.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          Specifies whether the details of the factorization are
          stored as an upper or lower triangular matrix:
          = 'U':  Upper triangular, form is A = P*U*D*(U**T)*(P**T);
          = 'L':  Lower triangular, form is A = P*L*D*(L**T)*(P**T).
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in]NRHS
          NRHS is INTEGER
          The number of right hand sides, i.e., the number of columns
          of the matrix B.  NRHS >= 0.
[in]A
          A is REAL array, dimension (LDA,N)
          Diagonal of the block diagonal matrix D and factors U or L
          as computed by SSYTRF_RK and SSYTRF_BK:
            a) ONLY diagonal elements of the symmetric block diagonal
               matrix D on the diagonal of A, i.e. D(k,k) = A(k,k);
               (superdiagonal (or subdiagonal) elements of D
                should be provided on entry in array E), and
            b) If UPLO = 'U': factor U in the superdiagonal part of A.
               If UPLO = 'L': factor L in the subdiagonal part of A.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[in]E
          E is REAL array, dimension (N)
          On entry, contains the superdiagonal (or subdiagonal)
          elements of the symmetric block diagonal matrix D
          with 1-by-1 or 2-by-2 diagonal blocks, where
          If UPLO = 'U': E(i) = D(i-1,i),i=2:N, E(1) not refernced;
          If UPLO = 'L': E(i) = D(i+1,i),i=1:N-1, E(N) not referenced.

          NOTE: For 1-by-1 diagonal block D(k), where
          1 <= k <= N, the element E(k) is not referenced in both
          UPLO = 'U' or UPLO = 'L' cases.
[in]IPIV
          IPIV is INTEGER array, dimension (N)
          Details of the interchanges and the block structure of D
          as determined by SSYTRF_RK or SSYTRF_BK.
[in,out]B
          B is REAL array, dimension (LDB,NRHS)
          On entry, the right hand side matrix B.
          On exit, the solution matrix X.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B.  LDB >= max(1,N).
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
December 2016
Contributors: 
  December 2016,  Igor Kozachenko,
                  Computer Science Division,
                  University of California, Berkeley

  September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
                  School of Mathematics,
                  University of Manchester

Definition at line 167 of file ssytrs_3.f.

167 *
168 * -- LAPACK computational routine (version 3.7.0) --
169 * -- LAPACK is a software package provided by Univ. of Tennessee, --
170 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
171 * December 2016
172 *
173 * .. Scalar Arguments ..
174  CHARACTER uplo
175  INTEGER info, lda, ldb, n, nrhs
176 * ..
177 * .. Array Arguments ..
178  INTEGER ipiv( * )
179  REAL a( lda, * ), b( ldb, * ), e( * )
180 * ..
181 *
182 * =====================================================================
183 *
184 * .. Parameters ..
185  REAL one
186  parameter ( one = 1.0e+0 )
187 * ..
188 * .. Local Scalars ..
189  LOGICAL upper
190  INTEGER i, j, k, kp
191  REAL ak, akm1, akm1k, bk, bkm1, denom
192 * ..
193 * .. External Functions ..
194  LOGICAL lsame
195  EXTERNAL lsame
196 * ..
197 * .. External Subroutines ..
198  EXTERNAL sscal, sswap, strsm, xerbla
199 * ..
200 * .. Intrinsic Functions ..
201  INTRINSIC abs, max
202 * ..
203 * .. Executable Statements ..
204 *
205  info = 0
206  upper = lsame( uplo, 'U' )
207  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
208  info = -1
209  ELSE IF( n.LT.0 ) THEN
210  info = -2
211  ELSE IF( nrhs.LT.0 ) THEN
212  info = -3
213  ELSE IF( lda.LT.max( 1, n ) ) THEN
214  info = -5
215  ELSE IF( ldb.LT.max( 1, n ) ) THEN
216  info = -9
217  END IF
218  IF( info.NE.0 ) THEN
219  CALL xerbla( 'SSYTRS_3', -info )
220  RETURN
221  END IF
222 *
223 * Quick return if possible
224 *
225  IF( n.EQ.0 .OR. nrhs.EQ.0 )
226  $ RETURN
227 *
228  IF( upper ) THEN
229 *
230 * Begin Upper
231 *
232 * Solve A*X = B, where A = U*D*U**T.
233 *
234 * P**T * B
235 *
236 * Interchange rows K and IPIV(K) of matrix B in the same order
237 * that the formation order of IPIV(I) vector for Upper case.
238 *
239 * (We can do the simple loop over IPIV with decrement -1,
240 * since the ABS value of IPIV(I) represents the row index
241 * of the interchange with row i in both 1x1 and 2x2 pivot cases)
242 *
243  DO k = n, 1, -1
244  kp = abs( ipiv( k ) )
245  IF( kp.NE.k ) THEN
246  CALL sswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
247  END IF
248  END DO
249 *
250 * Compute (U \P**T * B) -> B [ (U \P**T * B) ]
251 *
252  CALL strsm( 'L', 'U', 'N', 'U', n, nrhs, one, a, lda, b, ldb )
253 *
254 * Compute D \ B -> B [ D \ (U \P**T * B) ]
255 *
256  i = n
257  DO WHILE ( i.GE.1 )
258  IF( ipiv( i ).GT.0 ) THEN
259  CALL sscal( nrhs, one / a( i, i ), b( i, 1 ), ldb )
260  ELSE IF ( i.GT.1 ) THEN
261  akm1k = e( i )
262  akm1 = a( i-1, i-1 ) / akm1k
263  ak = a( i, i ) / akm1k
264  denom = akm1*ak - one
265  DO j = 1, nrhs
266  bkm1 = b( i-1, j ) / akm1k
267  bk = b( i, j ) / akm1k
268  b( i-1, j ) = ( ak*bkm1-bk ) / denom
269  b( i, j ) = ( akm1*bk-bkm1 ) / denom
270  END DO
271  i = i - 1
272  END IF
273  i = i - 1
274  END DO
275 *
276 * Compute (U**T \ B) -> B [ U**T \ (D \ (U \P**T * B) ) ]
277 *
278  CALL strsm( 'L', 'U', 'T', 'U', n, nrhs, one, a, lda, b, ldb )
279 *
280 * P * B [ P * (U**T \ (D \ (U \P**T * B) )) ]
281 *
282 * Interchange rows K and IPIV(K) of matrix B in reverse order
283 * from the formation order of IPIV(I) vector for Upper case.
284 *
285 * (We can do the simple loop over IPIV with increment 1,
286 * since the ABS value of IPIV(I) represents the row index
287 * of the interchange with row i in both 1x1 and 2x2 pivot cases)
288 *
289  DO k = 1, n, 1
290  kp = abs( ipiv( k ) )
291  IF( kp.NE.k ) THEN
292  CALL sswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
293  END IF
294  END DO
295 *
296  ELSE
297 *
298 * Begin Lower
299 *
300 * Solve A*X = B, where A = L*D*L**T.
301 *
302 * P**T * B
303 * Interchange rows K and IPIV(K) of matrix B in the same order
304 * that the formation order of IPIV(I) vector for Lower case.
305 *
306 * (We can do the simple loop over IPIV with increment 1,
307 * since the ABS value of IPIV(I) represents the row index
308 * of the interchange with row i in both 1x1 and 2x2 pivot cases)
309 *
310  DO k = 1, n, 1
311  kp = abs( ipiv( k ) )
312  IF( kp.NE.k ) THEN
313  CALL sswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
314  END IF
315  END DO
316 *
317 * Compute (L \P**T * B) -> B [ (L \P**T * B) ]
318 *
319  CALL strsm( 'L', 'L', 'N', 'U', n, nrhs, one, a, lda, b, ldb )
320 *
321 * Compute D \ B -> B [ D \ (L \P**T * B) ]
322 *
323  i = 1
324  DO WHILE ( i.LE.n )
325  IF( ipiv( i ).GT.0 ) THEN
326  CALL sscal( nrhs, one / a( i, i ), b( i, 1 ), ldb )
327  ELSE IF( i.LT.n ) THEN
328  akm1k = e( i )
329  akm1 = a( i, i ) / akm1k
330  ak = a( i+1, i+1 ) / akm1k
331  denom = akm1*ak - one
332  DO j = 1, nrhs
333  bkm1 = b( i, j ) / akm1k
334  bk = b( i+1, j ) / akm1k
335  b( i, j ) = ( ak*bkm1-bk ) / denom
336  b( i+1, j ) = ( akm1*bk-bkm1 ) / denom
337  END DO
338  i = i + 1
339  END IF
340  i = i + 1
341  END DO
342 *
343 * Compute (L**T \ B) -> B [ L**T \ (D \ (L \P**T * B) ) ]
344 *
345  CALL strsm('L', 'L', 'T', 'U', n, nrhs, one, a, lda, b, ldb )
346 *
347 * P * B [ P * (L**T \ (D \ (L \P**T * B) )) ]
348 *
349 * Interchange rows K and IPIV(K) of matrix B in reverse order
350 * from the formation order of IPIV(I) vector for Lower case.
351 *
352 * (We can do the simple loop over IPIV with decrement -1,
353 * since the ABS value of IPIV(I) represents the row index
354 * of the interchange with row i in both 1x1 and 2x2 pivot cases)
355 *
356  DO k = n, 1, -1
357  kp = abs( ipiv( k ) )
358  IF( kp.NE.k ) THEN
359  CALL sswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
360  END IF
361  END DO
362 *
363 * END Lower
364 *
365  END IF
366 *
367  RETURN
368 *
369 * End of SSYTRS_3
370 *
strsm
subroutine strsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
STRSM
Definition: strsm.f:183
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
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

Here is the call graph for this function:

Here is the caller graph for this function: