LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ csytrs_rook()

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.
Contributors:
   December 2016, 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 csytrs_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 COMPLEX A( LDA, * ), B( LDB, * )
148* ..
149*
150* =====================================================================
151*
152* .. Parameters ..
153 COMPLEX CONE
154 parameter( cone = ( 1.0e+0, 0.0e+0 ) )
155* ..
156* .. Local Scalars ..
157 LOGICAL UPPER
158 INTEGER J, K, KP
159 COMPLEX AK, AKM1, AKM1K, BK, BKM1, DENOM
160* ..
161* .. External Functions ..
162 LOGICAL LSAME
163 EXTERNAL lsame
164* ..
165* .. External Subroutines ..
166 EXTERNAL cgemv, cgeru, cscal, cswap, 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( 'CSYTRS_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 cswap( 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 cgeru( k-1, nrhs, -cone, 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 cscal( nrhs, cone / 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 cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
242*
243 kp = -ipiv( k-1 )
244 IF( kp.NE.k-1 )
245 $ CALL cswap( 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 cgeru( k-2, nrhs,-cone, a( 1, k ), 1, b( k, 1 ),
252 $ ldb, b( 1, 1 ), ldb )
253 CALL cgeru( k-2, nrhs,-cone, 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 - cone
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 cgemv( 'Transpose', k-1, nrhs, -cone, b,
297 $ ldb, a( 1, k ), 1, cone, b( k, 1 ), ldb )
298*
299* Interchange rows K and IPIV(K).
300*
301 kp = ipiv( k )
302 IF( kp.NE.k )
303 $ CALL cswap( 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 cgemv( 'Transpose', k-1, nrhs, -cone, b,
314 $ ldb, a( 1, k ), 1, cone, b( k, 1 ), ldb )
315 CALL cgemv( 'Transpose', k-1, nrhs, -cone, b,
316 $ ldb, a( 1, k+1 ), 1, cone, 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 cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
324*
325 kp = -ipiv( k+1 )
326 IF( kp.NE.k+1 )
327 $ CALL cswap( 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 cswap( 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 cgeru( n-k, nrhs, -cone, 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 cscal( nrhs, cone / 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 cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
382*
383 kp = -ipiv( k+1 )
384 IF( kp.NE.k+1 )
385 $ CALL cswap( 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 cgeru( n-k-1, nrhs,-cone, a( k+2, k ), 1, b( k, 1 ),
392 $ ldb, b( k+2, 1 ), ldb )
393 CALL cgeru( n-k-1, nrhs,-cone, 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 - cone
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 cgemv( 'Transpose', n-k, nrhs, -cone, b( k+1, 1 ),
437 $ ldb, a( k+1, k ), 1, cone, b( k, 1 ), ldb )
438*
439* Interchange rows K and IPIV(K).
440*
441 kp = ipiv( k )
442 IF( kp.NE.k )
443 $ CALL cswap( 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 cgemv( 'Transpose', n-k, nrhs, -cone, b( k+1, 1 ),
454 $ ldb, a( k+1, k ), 1, cone, b( k, 1 ), ldb )
455 CALL cgemv( 'Transpose', n-k, nrhs, -cone, b( k+1, 1 ),
456 $ ldb, a( k+1, k-1 ), 1, cone, 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 cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
465*
466 kp = -ipiv( k-1 )
467 IF( kp.NE.k-1 )
468 $ CALL cswap( 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 CSYTRS_ROOK
480*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine cgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
CGEMV
Definition cgemv.f:160
subroutine cgeru(m, n, alpha, x, incx, y, incy, a, lda)
CGERU
Definition cgeru.f:130
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine cscal(n, ca, cx, incx)
CSCAL
Definition cscal.f:78
subroutine cswap(n, cx, incx, cy, incy)
CSWAP
Definition cswap.f:81
Here is the call graph for this function:
Here is the caller graph for this function: