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

◆ dsytrs_3()

subroutine dsytrs_3 ( character  uplo,
integer  n,
integer  nrhs,
double precision, dimension( lda, * )  a,
integer  lda,
double precision, dimension( * )  e,
integer, dimension( * )  ipiv,
double precision, dimension( ldb, * )  b,
integer  ldb,
integer  info 
)

DSYTRS_3

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

Purpose:
 DSYTRS_3 solves a system of linear equations A * X = B with a real
 symmetric matrix A using the factorization computed
 by DSYTRF_RK or DSYTRF_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 DOUBLE PRECISION array, dimension (LDA,N)
          Diagonal of the block diagonal matrix D and factors U or L
          as computed by DSYTRF_RK and DSYTRF_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 DOUBLE PRECISION 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 DSYTRF_RK or DSYTRF_BK.
[in,out]B
          B is DOUBLE PRECISION 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.
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 dsytrs_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 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), E( * )
177* ..
178*
179* =====================================================================
180*
181* .. Parameters ..
182 DOUBLE PRECISION ONE
183 parameter( one = 1.0d+0 )
184* ..
185* .. Local Scalars ..
186 LOGICAL UPPER
187 INTEGER I, J, K, KP
188 DOUBLE PRECISION AK, AKM1, AKM1K, BK, BKM1, DENOM
189* ..
190* .. External Functions ..
191 LOGICAL LSAME
192 EXTERNAL lsame
193* ..
194* .. External Subroutines ..
195 EXTERNAL dscal, dswap, dtrsm, 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( 'DSYTRS_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 dswap( 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 dtrsm( '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 dscal( 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 dtrsm( '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
287 kp = abs( ipiv( k ) )
288 IF( kp.NE.k ) THEN
289 CALL dswap( 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
308 kp = abs( ipiv( k ) )
309 IF( kp.NE.k ) THEN
310 CALL dswap( 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 dtrsm( '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 dscal( 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 dtrsm('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 dswap( 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 DSYTRS_3
367*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
dswap
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82
dtrsm
subroutine dtrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
DTRSM
Definition dtrsm.f:181
Here is the call graph for this function:
Here is the caller graph for this function: