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

CSYTRS_ROOK

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

Purpose: 
 CSYTRS_ROOK solves a system of linear equations A*X = B with
 a complex symmetric matrix A using the factorization A = U*D*U**T or
 A = L*D*L**T computed by CSYTRF_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 COMPLEX array, dimension (LDA,N)
          The block diagonal matrix D and the multipliers used to
          obtain the factor U or L as computed by CSYTRF_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 CSYTRF_ROOK.
[in,out]B
          B is COMPLEX 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
November 2015
Contributors: 
   November 2015, 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 csytrs_rook.f.

138 *
139 * -- LAPACK computational routine (version 3.6.0) --
140 * -- LAPACK is a software package provided by Univ. of Tennessee, --
141 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
142 * November 2015
143 *
144 * .. Scalar Arguments ..
145  CHARACTER uplo
146  INTEGER info, lda, ldb, n, nrhs
147 * ..
148 * .. Array Arguments ..
149  INTEGER ipiv( * )
150  COMPLEX a( lda, * ), b( ldb, * )
151 * ..
152 *
153 * =====================================================================
154 *
155 * .. Parameters ..
156  COMPLEX cone
157  parameter ( cone = ( 1.0e+0, 0.0e+0 ) )
158 * ..
159 * .. Local Scalars ..
160  LOGICAL upper
161  INTEGER j, k, kp
162  COMPLEX ak, akm1, akm1k, bk, bkm1, denom
163 * ..
164 * .. External Functions ..
165  LOGICAL lsame
166  EXTERNAL lsame
167 * ..
168 * .. External Subroutines ..
169  EXTERNAL cgemv, cgeru, cscal, cswap, 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( 'CSYTRS_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 cswap( 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 cgeru( k-1, nrhs, -cone, 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 cscal( nrhs, cone / 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 cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
245 *
246  kp = -ipiv( k-1 )
247  IF( kp.NE.k-1 )
248  $ CALL cswap( 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 cgeru( k-2, nrhs,-cone, a( 1, k ), 1, b( k, 1 ),
255  $ ldb, b( 1, 1 ), ldb )
256  CALL cgeru( k-2, nrhs,-cone, 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 - cone
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 cgemv( 'Transpose', k-1, nrhs, -cone, b,
300  $ ldb, a( 1, k ), 1, cone, b( k, 1 ), ldb )
301 *
302 * Interchange rows K and IPIV(K).
303 *
304  kp = ipiv( k )
305  IF( kp.NE.k )
306  $ CALL cswap( 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 cgemv( 'Transpose', k-1, nrhs, -cone, b,
317  $ ldb, a( 1, k ), 1, cone, b( k, 1 ), ldb )
318  CALL cgemv( 'Transpose', k-1, nrhs, -cone, b,
319  $ ldb, a( 1, k+1 ), 1, cone, 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 cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
327 *
328  kp = -ipiv( k+1 )
329  IF( kp.NE.k+1 )
330  $ CALL cswap( 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 cswap( 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 cgeru( n-k, nrhs, -cone, 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 cscal( nrhs, cone / 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 cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
385 *
386  kp = -ipiv( k+1 )
387  IF( kp.NE.k+1 )
388  $ CALL cswap( 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 cgeru( n-k-1, nrhs,-cone, a( k+2, k ), 1, b( k, 1 ),
395  $ ldb, b( k+2, 1 ), ldb )
396  CALL cgeru( n-k-1, nrhs,-cone, 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 - cone
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 cgemv( 'Transpose', n-k, nrhs, -cone, b( k+1, 1 ),
440  $ ldb, a( k+1, k ), 1, cone, b( k, 1 ), ldb )
441 *
442 * Interchange rows K and IPIV(K).
443 *
444  kp = ipiv( k )
445  IF( kp.NE.k )
446  $ CALL cswap( 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 cgemv( 'Transpose', n-k, nrhs, -cone, b( k+1, 1 ),
457  $ ldb, a( k+1, k ), 1, cone, b( k, 1 ), ldb )
458  CALL cgemv( 'Transpose', n-k, nrhs, -cone, b( k+1, 1 ),
459  $ ldb, a( k+1, k-1 ), 1, cone, 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 cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
468 *
469  kp = -ipiv( k-1 )
470  IF( kp.NE.k-1 )
471  $ CALL cswap( 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 CSYTRS_ROOK
483 *
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
cscal
subroutine cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:54
cgemv
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
Definition: cgemv.f:160
cswap
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:52
cgeru
subroutine cgeru(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
CGERU
Definition: cgeru.f:132
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: