LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ ssytrs2()

subroutine ssytrs2 ( character uplo,
integer n,
integer nrhs,
real, dimension( lda, * ) a,
integer lda,
integer, dimension( * ) ipiv,
real, dimension( ldb, * ) b,
integer ldb,
real, dimension( * ) work,
integer info )

SSYTRS2

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

Purpose:
!>
!> SSYTRS2 solves a system of linear equations A*X = B with a real
!> symmetric matrix A using the factorization A = U*D*U**T or
!> A = L*D*L**T computed by SSYTRF and converted by SSYCONV.
!>
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 = U*D*U**T;
!>          = 'L':  Lower triangular, form is A = L*D*L**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,out]A
!>          A is REAL array, dimension (LDA,N)
!>          The block diagonal matrix D and the multipliers used to
!>          obtain the factor U or L as computed by SSYTRF.
!>          Note that A is input / output. This might be counter-intuitive,
!>          and one may think that A is input only. A is input / output. This
!>          is because, at the start of the subroutine, we permute A in a
!>           form and then we permute A back to its original form at
!>          the end.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[in]IPIV
!>          IPIV is INTEGER array, dimension (N)
!>          Details of the interchanges and the block structure of D
!>          as determined by SSYTRF.
!>
[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]WORK
!>          WORK is REAL array, dimension (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.

Definition at line 128 of file ssytrs2.f.

130*
131* -- LAPACK computational routine --
132* -- LAPACK is a software package provided by Univ. of Tennessee, --
133* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
134*
135* .. Scalar Arguments ..
136 CHARACTER UPLO
137 INTEGER INFO, LDA, LDB, N, NRHS
138* ..
139* .. Array Arguments ..
140 INTEGER IPIV( * )
141 REAL A( LDA, * ), B( LDB, * ), WORK( * )
142* ..
143*
144* =====================================================================
145*
146* .. Parameters ..
147 REAL ONE
148 parameter( one = 1.0e+0 )
149* ..
150* .. Local Scalars ..
151 LOGICAL UPPER
152 INTEGER I, IINFO, J, K, KP
153 REAL AK, AKM1, AKM1K, BK, BKM1, DENOM
154* ..
155* .. External Functions ..
156 LOGICAL LSAME
157 EXTERNAL lsame
158* ..
159* .. External Subroutines ..
160 EXTERNAL sscal, ssyconv, sswap, strsm,
161 $ xerbla
162* ..
163* .. Intrinsic Functions ..
164 INTRINSIC max
165* ..
166* .. Executable Statements ..
167*
168 info = 0
169 upper = lsame( uplo, 'U' )
170 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
171 info = -1
172 ELSE IF( n.LT.0 ) THEN
173 info = -2
174 ELSE IF( nrhs.LT.0 ) THEN
175 info = -3
176 ELSE IF( lda.LT.max( 1, n ) ) THEN
177 info = -5
178 ELSE IF( ldb.LT.max( 1, n ) ) THEN
179 info = -8
180 END IF
181 IF( info.NE.0 ) THEN
182 CALL xerbla( 'SSYTRS2', -info )
183 RETURN
184 END IF
185*
186* Quick return if possible
187*
188 IF( n.EQ.0 .OR. nrhs.EQ.0 )
189 $ RETURN
190*
191* Convert A
192*
193 CALL ssyconv( uplo, 'C', n, a, lda, ipiv, work, iinfo )
194*
195 IF( upper ) THEN
196*
197* Solve A*X = B, where A = U*D*U**T.
198*
199* P**T * B
200 k=n
201 DO WHILE ( k .GE. 1 )
202 IF( ipiv( k ).GT.0 ) THEN
203* 1 x 1 diagonal block
204* Interchange rows K and IPIV(K).
205 kp = ipiv( k )
206 IF( kp.NE.k )
207 $ CALL sswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
208 k=k-1
209 ELSE
210* 2 x 2 diagonal block
211* Interchange rows K-1 and -IPIV(K).
212 kp = -ipiv( k )
213 IF( kp.EQ.-ipiv( k-1 ) )
214 $ CALL sswap( nrhs, b( k-1, 1 ), ldb, b( kp, 1 ), ldb )
215 k=k-2
216 END IF
217 END DO
218*
219* Compute (U \P**T * B) -> B [ (U \P**T * B) ]
220*
221 CALL strsm('L','U','N','U',n,nrhs,one,a,lda,b,ldb)
222*
223* Compute D \ B -> B [ D \ (U \P**T * B) ]
224*
225 i=n
226 DO WHILE ( i .GE. 1 )
227 IF( ipiv(i) .GT. 0 ) THEN
228 CALL sscal( nrhs, one / a( i, i ), b( i, 1 ), ldb )
229 ELSEIF ( i .GT. 1) THEN
230 IF ( ipiv(i-1) .EQ. ipiv(i) ) THEN
231 akm1k = work(i)
232 akm1 = a( i-1, i-1 ) / akm1k
233 ak = a( i, i ) / akm1k
234 denom = akm1*ak - one
235 DO 15 j = 1, nrhs
236 bkm1 = b( i-1, j ) / akm1k
237 bk = b( i, j ) / akm1k
238 b( i-1, j ) = ( ak*bkm1-bk ) / denom
239 b( i, j ) = ( akm1*bk-bkm1 ) / denom
240 15 CONTINUE
241 i = i - 1
242 ENDIF
243 ENDIF
244 i = i - 1
245 END DO
246*
247* Compute (U**T \ B) -> B [ U**T \ (D \ (U \P**T * B) ) ]
248*
249 CALL strsm('L','U','T','U',n,nrhs,one,a,lda,b,ldb)
250*
251* P * B [ P * (U**T \ (D \ (U \P**T * B) )) ]
252*
253 k=1
254 DO WHILE ( k .LE. n )
255 IF( ipiv( k ).GT.0 ) THEN
256* 1 x 1 diagonal block
257* Interchange rows K and IPIV(K).
258 kp = ipiv( k )
259 IF( kp.NE.k )
260 $ CALL sswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
261 k=k+1
262 ELSE
263* 2 x 2 diagonal block
264* Interchange rows K-1 and -IPIV(K).
265 kp = -ipiv( k )
266 IF( k .LT. n .AND. kp.EQ.-ipiv( k+1 ) )
267 $ CALL sswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
268 k=k+2
269 ENDIF
270 END DO
271*
272 ELSE
273*
274* Solve A*X = B, where A = L*D*L**T.
275*
276* P**T * B
277 k=1
278 DO WHILE ( k .LE. n )
279 IF( ipiv( k ).GT.0 ) THEN
280* 1 x 1 diagonal block
281* Interchange rows K and IPIV(K).
282 kp = ipiv( k )
283 IF( kp.NE.k )
284 $ CALL sswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
285 k=k+1
286 ELSE
287* 2 x 2 diagonal block
288* Interchange rows K and -IPIV(K+1).
289 kp = -ipiv( k+1 )
290 IF( kp.EQ.-ipiv( k ) )
291 $ CALL sswap( nrhs, b( k+1, 1 ), ldb, b( kp, 1 ), ldb )
292 k=k+2
293 ENDIF
294 END DO
295*
296* Compute (L \P**T * B) -> B [ (L \P**T * B) ]
297*
298 CALL strsm('L','L','N','U',n,nrhs,one,a,lda,b,ldb)
299*
300* Compute D \ B -> B [ D \ (L \P**T * B) ]
301*
302 i=1
303 DO WHILE ( i .LE. n )
304 IF( ipiv(i) .GT. 0 ) THEN
305 CALL sscal( nrhs, one / a( i, i ), b( i, 1 ), ldb )
306 ELSE
307 akm1k = work(i)
308 akm1 = a( i, i ) / akm1k
309 ak = a( i+1, i+1 ) / akm1k
310 denom = akm1*ak - one
311 DO 25 j = 1, nrhs
312 bkm1 = b( i, j ) / akm1k
313 bk = b( i+1, j ) / akm1k
314 b( i, j ) = ( ak*bkm1-bk ) / denom
315 b( i+1, j ) = ( akm1*bk-bkm1 ) / denom
316 25 CONTINUE
317 i = i + 1
318 ENDIF
319 i = i + 1
320 END DO
321*
322* Compute (L**T \ B) -> B [ L**T \ (D \ (L \P**T * B) ) ]
323*
324 CALL strsm('L','L','T','U',n,nrhs,one,a,lda,b,ldb)
325*
326* P * B [ P * (L**T \ (D \ (L \P**T * B) )) ]
327*
328 k=n
329 DO WHILE ( k .GE. 1 )
330 IF( ipiv( k ).GT.0 ) THEN
331* 1 x 1 diagonal block
332* Interchange rows K and IPIV(K).
333 kp = ipiv( k )
334 IF( kp.NE.k )
335 $ CALL sswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
336 k=k-1
337 ELSE
338* 2 x 2 diagonal block
339* Interchange rows K-1 and -IPIV(K).
340 kp = -ipiv( k )
341 IF( k.GT.1 .AND. kp.EQ.-ipiv( k-1 ) )
342 $ CALL sswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
343 k=k-2
344 ENDIF
345 END DO
346*
347 END IF
348*
349* Revert A
350*
351 CALL ssyconv( uplo, 'R', n, a, lda, ipiv, work, iinfo )
352*
353 RETURN
354*
355* End of SSYTRS2
356*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
sscal
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
sswap
subroutine sswap(n, sx, incx, sy, incy)
SSWAP
Definition sswap.f:82
ssyconv
subroutine ssyconv(uplo, way, n, a, lda, ipiv, e, info)
SSYCONV
Definition ssyconv.f:112
strsm
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: