 LAPACK  3.10.0 LAPACK: Linear Algebra PACKage

◆ ssytrs_3()

 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

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 referenced; 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
Contributors:
June 2017,  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 163 of file ssytrs_3.f.

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