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

SSYTRS_ROOK

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

Purpose: 
 SSYTRS_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 SSYTRF_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 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_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 SSYTRF_ROOK.
[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
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 ssytrs_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  REAL a( lda, * ), b( ldb, * )
151 * ..
152 *
153 * =====================================================================
154 *
155 * .. Parameters ..
156  REAL one
157  parameter ( one = 1.0e+0 )
158 * ..
159 * .. Local Scalars ..
160  LOGICAL upper
161  INTEGER j, k, kp
162  REAL ak, akm1, akm1k, bk, bkm1, denom
163 * ..
164 * .. External Functions ..
165  LOGICAL lsame
166  EXTERNAL lsame
167 * ..
168 * .. External Subroutines ..
169  EXTERNAL sgemv, sger, sscal, sswap, 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( 'SSYTRS_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 sswap( 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 sger( 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 sscal( 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 sswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
245 *
246  kp = -ipiv( k-1 )
247  IF( kp.NE.k-1 )
248  $ CALL sswap( 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 sger( k-2, nrhs, -one, a( 1, k ), 1, b( k, 1 ),
255  $ ldb, b( 1, 1 ), ldb )
256  CALL sger( 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 sgemv( '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 sswap( 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 sgemv( 'Transpose', k-1, nrhs, -one, b,
317  $ ldb, a( 1, k ), 1, one, b( k, 1 ), ldb )
318  CALL sgemv( '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 sswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
327 *
328  kp = -ipiv( k+1 )
329  IF( kp.NE.k+1 )
330  $ CALL sswap( 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 sswap( 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 sger( 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 sscal( 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 sswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
385 *
386  kp = -ipiv( k+1 )
387  IF( kp.NE.k+1 )
388  $ CALL sswap( 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 sger( n-k-1, nrhs, -one, a( k+2, k ), 1, b( k, 1 ),
395  $ ldb, b( k+2, 1 ), ldb )
396  CALL sger( 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 sgemv( '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 sswap( 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 sgemv( 'Transpose', n-k, nrhs, -one, b( k+1, 1 ),
457  $ ldb, a( k+1, k ), 1, one, b( k, 1 ), ldb )
458  CALL sgemv( '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 sswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
468 *
469  kp = -ipiv( k-1 )
470  IF( kp.NE.k-1 )
471  $ CALL sswap( 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 SSYTRS_ROOK
483 *
sger
subroutine sger(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
SGER
Definition: sger.f:132
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
sgemv
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
Definition: sgemv.f:158
sscal
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:55
sswap
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:53
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: