LAPACK  3.9.0
LAPACK: Linear Algebra PACKage

◆ 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.
Date
June 2017
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 167 of file dsytrs_3.f.

167 *
168 * -- LAPACK computational routine (version 3.7.1) --
169 * -- LAPACK is a software package provided by Univ. of Tennessee, --
170 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
171 * June 2017
172 *
173 * .. Scalar Arguments ..
174  CHARACTER UPLO
175  INTEGER INFO, LDA, LDB, N, NRHS
176 * ..
177 * .. Array Arguments ..
178  INTEGER IPIV( * )
179  DOUBLE PRECISION A( LDA, * ), B( LDB, * ), E( * )
180 * ..
181 *
182 * =====================================================================
183 *
184 * .. Parameters ..
185  DOUBLE PRECISION ONE
186  parameter( one = 1.0d+0 )
187 * ..
188 * .. Local Scalars ..
189  LOGICAL UPPER
190  INTEGER I, J, K, KP
191  DOUBLE PRECISION AK, AKM1, AKM1K, BK, BKM1, DENOM
192 * ..
193 * .. External Functions ..
194  LOGICAL LSAME
195  EXTERNAL lsame
196 * ..
197 * .. External Subroutines ..
198  EXTERNAL dscal, dswap, dtrsm, 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( 'DSYTRS_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 dswap( 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 dtrsm( '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 dscal( 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 dtrsm( '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
290  kp = abs( ipiv( k ) )
291  IF( kp.NE.k ) THEN
292  CALL dswap( 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
311  kp = abs( ipiv( k ) )
312  IF( kp.NE.k ) THEN
313  CALL dswap( 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 dtrsm( '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 dscal( 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 dtrsm('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 dswap( 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 DSYTRS_3
370 *
Here is the call graph for this function:
Here is the caller graph for this function:
dtrsm
subroutine dtrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
DTRSM
Definition: dtrsm.f:183
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
dswap
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.f:84
dscal
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:81