LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine dsytrs_rook ( character  UPLO,
integer  N,
integer  NRHS,
double precision, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
double precision, dimension( ldb, * )  B,
integer  LDB,
integer  INFO 
)

DSYTRS_ROOK

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

Purpose: 
 DSYTRS_ROOK 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 DSYTRF_ROOK.
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]A
          A is DOUBLE PRECISION array, dimension (LDA,N)
          The block diagonal matrix D and the multipliers used to
          obtain the factor U or L as computed by DSYTRF_ROOK.
[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 DSYTRF_ROOK.
[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
April 2012
Contributors: 
   April 2012, 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 138 of file dsytrs_rook.f.

138 *
139 * -- LAPACK computational routine (version 3.4.1) --
140 * -- LAPACK is a software package provided by Univ. of Tennessee, --
141 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
142 * April 2012
143 *
144 * .. Scalar Arguments ..
145  CHARACTER uplo
146  INTEGER info, lda, ldb, n, nrhs
147 * ..
148 * .. Array Arguments ..
149  INTEGER ipiv( * )
150  DOUBLE PRECISION a( lda, * ), b( ldb, * )
151 * ..
152 *
153 * =====================================================================
154 *
155 * .. Parameters ..
156  DOUBLE PRECISION one
157  parameter ( one = 1.0d+0 )
158 * ..
159 * .. Local Scalars ..
160  LOGICAL upper
161  INTEGER j, k, kp
162  DOUBLE PRECISION ak, akm1, akm1k, bk, bkm1, denom
163 * ..
164 * .. External Functions ..
165  LOGICAL lsame
166  EXTERNAL lsame
167 * ..
168 * .. External Subroutines ..
169  EXTERNAL dgemv, dger, dscal, dswap, xerbla
170 * ..
171 * .. Intrinsic Functions ..
172  INTRINSIC max
173 * ..
174 * .. Executable Statements ..
175 *
176  info = 0
177  upper = lsame( uplo, 'U' )
178  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
179  info = -1
180  ELSE IF( n.LT.0 ) THEN
181  info = -2
182  ELSE IF( nrhs.LT.0 ) THEN
183  info = -3
184  ELSE IF( lda.LT.max( 1, n ) ) THEN
185  info = -5
186  ELSE IF( ldb.LT.max( 1, n ) ) THEN
187  info = -8
188  END IF
189  IF( info.NE.0 ) THEN
190  CALL xerbla( 'DSYTRS_ROOK', -info )
191  RETURN
192  END IF
193 *
194 * Quick return if possible
195 *
196  IF( n.EQ.0 .OR. nrhs.EQ.0 )
197  $ RETURN
198 *
199  IF( upper ) THEN
200 *
201 * Solve A*X = B, where A = U*D*U**T.
202 *
203 * First solve U*D*X = B, overwriting B with X.
204 *
205 * K is the main loop index, decreasing from N to 1 in steps of
206 * 1 or 2, depending on the size of the diagonal blocks.
207 *
208  k = n
209  10 CONTINUE
210 *
211 * If K < 1, exit from loop.
212 *
213  IF( k.LT.1 )
214  $ GO TO 30
215 *
216  IF( ipiv( k ).GT.0 ) THEN
217 *
218 * 1 x 1 diagonal block
219 *
220 * Interchange rows K and IPIV(K).
221 *
222  kp = ipiv( k )
223  IF( kp.NE.k )
224  $ CALL dswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
225 *
226 * Multiply by inv(U(K)), where U(K) is the transformation
227 * stored in column K of A.
228 *
229  CALL dger( k-1, nrhs, -one, a( 1, k ), 1, b( k, 1 ), ldb,
230  $ b( 1, 1 ), ldb )
231 *
232 * Multiply by the inverse of the diagonal block.
233 *
234  CALL dscal( nrhs, one / a( k, k ), b( k, 1 ), ldb )
235  k = k - 1
236  ELSE
237 *
238 * 2 x 2 diagonal block
239 *
240 * Interchange rows K and -IPIV(K) THEN K-1 and -IPIV(K-1)
241 *
242  kp = -ipiv( k )
243  IF( kp.NE.k )
244  $ CALL dswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
245 *
246  kp = -ipiv( k-1 )
247  IF( kp.NE.k-1 )
248  $ CALL dswap( nrhs, b( k-1, 1 ), ldb, b( kp, 1 ), ldb )
249 *
250 * Multiply by inv(U(K)), where U(K) is the transformation
251 * stored in columns K-1 and K of A.
252 *
253  IF( k.GT.2 ) THEN
254  CALL dger( k-2, nrhs, -one, a( 1, k ), 1, b( k, 1 ),
255  $ ldb, b( 1, 1 ), ldb )
256  CALL dger( k-2, nrhs, -one, a( 1, k-1 ), 1, b( k-1, 1 ),
257  $ ldb, b( 1, 1 ), ldb )
258  END IF
259 *
260 * Multiply by the inverse of the diagonal block.
261 *
262  akm1k = a( k-1, k )
263  akm1 = a( k-1, k-1 ) / akm1k
264  ak = a( k, k ) / akm1k
265  denom = akm1*ak - one
266  DO 20 j = 1, nrhs
267  bkm1 = b( k-1, j ) / akm1k
268  bk = b( k, j ) / akm1k
269  b( k-1, j ) = ( ak*bkm1-bk ) / denom
270  b( k, j ) = ( akm1*bk-bkm1 ) / denom
271  20 CONTINUE
272  k = k - 2
273  END IF
274 *
275  GO TO 10
276  30 CONTINUE
277 *
278 * Next solve U**T *X = B, overwriting B with X.
279 *
280 * K is the main loop index, increasing from 1 to N in steps of
281 * 1 or 2, depending on the size of the diagonal blocks.
282 *
283  k = 1
284  40 CONTINUE
285 *
286 * If K > N, exit from loop.
287 *
288  IF( k.GT.n )
289  $ GO TO 50
290 *
291  IF( ipiv( k ).GT.0 ) THEN
292 *
293 * 1 x 1 diagonal block
294 *
295 * Multiply by inv(U**T(K)), where U(K) is the transformation
296 * stored in column K of A.
297 *
298  IF( k.GT.1 )
299  $ CALL dgemv( 'Transpose', k-1, nrhs, -one, b,
300  $ ldb, a( 1, k ), 1, one, b( k, 1 ), ldb )
301 *
302 * Interchange rows K and IPIV(K).
303 *
304  kp = ipiv( k )
305  IF( kp.NE.k )
306  $ CALL dswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
307  k = k + 1
308  ELSE
309 *
310 * 2 x 2 diagonal block
311 *
312 * Multiply by inv(U**T(K+1)), where U(K+1) is the transformation
313 * stored in columns K and K+1 of A.
314 *
315  IF( k.GT.1 ) THEN
316  CALL dgemv( 'Transpose', k-1, nrhs, -one, b,
317  $ ldb, a( 1, k ), 1, one, b( k, 1 ), ldb )
318  CALL dgemv( 'Transpose', k-1, nrhs, -one, b,
319  $ ldb, a( 1, k+1 ), 1, one, b( k+1, 1 ), ldb )
320  END IF
321 *
322 * Interchange rows K and -IPIV(K) THEN K+1 and -IPIV(K+1).
323 *
324  kp = -ipiv( k )
325  IF( kp.NE.k )
326  $ CALL dswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
327 *
328  kp = -ipiv( k+1 )
329  IF( kp.NE.k+1 )
330  $ CALL dswap( nrhs, b( k+1, 1 ), ldb, b( kp, 1 ), ldb )
331 *
332  k = k + 2
333  END IF
334 *
335  GO TO 40
336  50 CONTINUE
337 *
338  ELSE
339 *
340 * Solve A*X = B, where A = L*D*L**T.
341 *
342 * First solve L*D*X = B, overwriting B with X.
343 *
344 * K is the main loop index, increasing from 1 to N in steps of
345 * 1 or 2, depending on the size of the diagonal blocks.
346 *
347  k = 1
348  60 CONTINUE
349 *
350 * If K > N, exit from loop.
351 *
352  IF( k.GT.n )
353  $ GO TO 80
354 *
355  IF( ipiv( k ).GT.0 ) THEN
356 *
357 * 1 x 1 diagonal block
358 *
359 * Interchange rows K and IPIV(K).
360 *
361  kp = ipiv( k )
362  IF( kp.NE.k )
363  $ CALL dswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
364 *
365 * Multiply by inv(L(K)), where L(K) is the transformation
366 * stored in column K of A.
367 *
368  IF( k.LT.n )
369  $ CALL dger( n-k, nrhs, -one, a( k+1, k ), 1, b( k, 1 ),
370  $ ldb, b( k+1, 1 ), ldb )
371 *
372 * Multiply by the inverse of the diagonal block.
373 *
374  CALL dscal( nrhs, one / a( k, k ), b( k, 1 ), ldb )
375  k = k + 1
376  ELSE
377 *
378 * 2 x 2 diagonal block
379 *
380 * Interchange rows K and -IPIV(K) THEN K+1 and -IPIV(K+1)
381 *
382  kp = -ipiv( k )
383  IF( kp.NE.k )
384  $ CALL dswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
385 *
386  kp = -ipiv( k+1 )
387  IF( kp.NE.k+1 )
388  $ CALL dswap( nrhs, b( k+1, 1 ), ldb, b( kp, 1 ), ldb )
389 *
390 * Multiply by inv(L(K)), where L(K) is the transformation
391 * stored in columns K and K+1 of A.
392 *
393  IF( k.LT.n-1 ) THEN
394  CALL dger( n-k-1, nrhs, -one, a( k+2, k ), 1, b( k, 1 ),
395  $ ldb, b( k+2, 1 ), ldb )
396  CALL dger( n-k-1, nrhs, -one, a( k+2, k+1 ), 1,
397  $ b( k+1, 1 ), ldb, b( k+2, 1 ), ldb )
398  END IF
399 *
400 * Multiply by the inverse of the diagonal block.
401 *
402  akm1k = a( k+1, k )
403  akm1 = a( k, k ) / akm1k
404  ak = a( k+1, k+1 ) / akm1k
405  denom = akm1*ak - one
406  DO 70 j = 1, nrhs
407  bkm1 = b( k, j ) / akm1k
408  bk = b( k+1, j ) / akm1k
409  b( k, j ) = ( ak*bkm1-bk ) / denom
410  b( k+1, j ) = ( akm1*bk-bkm1 ) / denom
411  70 CONTINUE
412  k = k + 2
413  END IF
414 *
415  GO TO 60
416  80 CONTINUE
417 *
418 * Next solve L**T *X = B, overwriting B with X.
419 *
420 * K is the main loop index, decreasing from N to 1 in steps of
421 * 1 or 2, depending on the size of the diagonal blocks.
422 *
423  k = n
424  90 CONTINUE
425 *
426 * If K < 1, exit from loop.
427 *
428  IF( k.LT.1 )
429  $ GO TO 100
430 *
431  IF( ipiv( k ).GT.0 ) THEN
432 *
433 * 1 x 1 diagonal block
434 *
435 * Multiply by inv(L**T(K)), where L(K) is the transformation
436 * stored in column K of A.
437 *
438  IF( k.LT.n )
439  $ CALL dgemv( 'Transpose', n-k, nrhs, -one, b( k+1, 1 ),
440  $ ldb, a( k+1, k ), 1, one, b( k, 1 ), ldb )
441 *
442 * Interchange rows K and IPIV(K).
443 *
444  kp = ipiv( k )
445  IF( kp.NE.k )
446  $ CALL dswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
447  k = k - 1
448  ELSE
449 *
450 * 2 x 2 diagonal block
451 *
452 * Multiply by inv(L**T(K-1)), where L(K-1) is the transformation
453 * stored in columns K-1 and K of A.
454 *
455  IF( k.LT.n ) THEN
456  CALL dgemv( 'Transpose', n-k, nrhs, -one, b( k+1, 1 ),
457  $ ldb, a( k+1, k ), 1, one, b( k, 1 ), ldb )
458  CALL dgemv( 'Transpose', n-k, nrhs, -one, b( k+1, 1 ),
459  $ ldb, a( k+1, k-1 ), 1, one, b( k-1, 1 ),
460  $ ldb )
461  END IF
462 *
463 * Interchange rows K and -IPIV(K) THEN K-1 and -IPIV(K-1)
464 *
465  kp = -ipiv( k )
466  IF( kp.NE.k )
467  $ CALL dswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
468 *
469  kp = -ipiv( k-1 )
470  IF( kp.NE.k-1 )
471  $ CALL dswap( nrhs, b( k-1, 1 ), ldb, b( kp, 1 ), ldb )
472 *
473  k = k - 2
474  END IF
475 *
476  GO TO 90
477  100 CONTINUE
478  END IF
479 *
480  RETURN
481 *
482 * End of DSYTRS_ROOK
483 *
dgemv
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
Definition: dgemv.f:158
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
dswap
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.f:53
dger
subroutine dger(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
DGER
Definition: dger.f:132
dscal
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:55
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the call graph for this function:

Here is the caller graph for this function: