LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ ssytrs_rook()

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.
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 134 of file ssytrs_rook.f.

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