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

◆ csytrs_3()

subroutine csytrs_3 ( character  uplo,
integer  n,
integer  nrhs,
complex, dimension( lda, * )  a,
integer  lda,
complex, dimension( * )  e,
integer, dimension( * )  ipiv,
complex, dimension( ldb, * )  b,
integer  ldb,
integer  info 
)

CSYTRS_3

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

Purpose:
 CSYTRS_3 solves a system of linear equations A * X = B with a complex
 symmetric matrix A using the factorization computed
 by CSYTRF_RK or CSYTRF_BK:

    A = P*U*D*(U**T)*(P**T) or A = P*L*D*(L**T)*(P**T),

 where U (or L) is unit upper (or lower) triangular matrix,
 U**T (or L**T) is the transpose of U (or L), P is a permutation
 matrix, P**T is the transpose of P, and D is symmetric and block
 diagonal with 1-by-1 and 2-by-2 diagonal blocks.

 This algorithm is using Level 3 BLAS.
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 = P*U*D*(U**T)*(P**T);
          = 'L':  Lower triangular, form is A = P*L*D*(L**T)*(P**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)
          Diagonal of the block diagonal matrix D and factors U or L
          as computed by CSYTRF_RK and CSYTRF_BK:
            a) ONLY diagonal elements of the symmetric block diagonal
               matrix D on the diagonal of A, i.e. D(k,k) = A(k,k);
               (superdiagonal (or subdiagonal) elements of D
                should be provided on entry in array E), and
            b) If UPLO = 'U': factor U in the superdiagonal part of A.
               If UPLO = 'L': factor L in the subdiagonal part of A.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[in]E
          E is COMPLEX array, dimension (N)
          On entry, contains the superdiagonal (or subdiagonal)
          elements of the symmetric block diagonal matrix D
          with 1-by-1 or 2-by-2 diagonal blocks, where
          If UPLO = 'U': E(i) = D(i-1,i),i=2:N, E(1) not referenced;
          If UPLO = 'L': E(i) = D(i+1,i),i=1:N-1, E(N) not referenced.

          NOTE: For 1-by-1 diagonal block D(k), where
          1 <= k <= N, the element E(k) is not referenced in both
          UPLO = 'U' or UPLO = 'L' cases.
[in]IPIV
          IPIV is INTEGER array, dimension (N)
          Details of the interchanges and the block structure of D
          as determined by CSYTRF_RK or CSYTRF_BK.
[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:
  June 2017,  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 163 of file csytrs_3.f.

165*
166* -- LAPACK computational routine --
167* -- LAPACK is a software package provided by Univ. of Tennessee, --
168* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
169*
170* .. Scalar Arguments ..
171 CHARACTER UPLO
172 INTEGER INFO, LDA, LDB, N, NRHS
173* ..
174* .. Array Arguments ..
175 INTEGER IPIV( * )
176 COMPLEX A( LDA, * ), B( LDB, * ), E( * )
177* ..
178*
179* =====================================================================
180*
181* .. Parameters ..
182 COMPLEX ONE
183 parameter( one = ( 1.0e+0,0.0e+0 ) )
184* ..
185* .. Local Scalars ..
186 LOGICAL UPPER
187 INTEGER I, J, K, KP
188 COMPLEX AK, AKM1, AKM1K, BK, BKM1, DENOM
189* ..
190* .. External Functions ..
191 LOGICAL LSAME
192 EXTERNAL lsame
193* ..
194* .. External Subroutines ..
195 EXTERNAL cscal, cswap, ctrsm, xerbla
196* ..
197* .. Intrinsic Functions ..
198 INTRINSIC abs, max
199* ..
200* .. Executable Statements ..
201*
202 info = 0
203 upper = lsame( uplo, 'U' )
204 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
205 info = -1
206 ELSE IF( n.LT.0 ) THEN
207 info = -2
208 ELSE IF( nrhs.LT.0 ) THEN
209 info = -3
210 ELSE IF( lda.LT.max( 1, n ) ) THEN
211 info = -5
212 ELSE IF( ldb.LT.max( 1, n ) ) THEN
213 info = -9
214 END IF
215 IF( info.NE.0 ) THEN
216 CALL xerbla( 'CSYTRS_3', -info )
217 RETURN
218 END IF
219*
220* Quick return if possible
221*
222 IF( n.EQ.0 .OR. nrhs.EQ.0 )
223 $ RETURN
224*
225 IF( upper ) THEN
226*
227* Begin Upper
228*
229* Solve A*X = B, where A = U*D*U**T.
230*
231* P**T * B
232*
233* Interchange rows K and IPIV(K) of matrix B in the same order
234* that the formation order of IPIV(I) vector for Upper case.
235*
236* (We can do the simple loop over IPIV with decrement -1,
237* since the ABS value of IPIV(I) represents the row index
238* of the interchange with row i in both 1x1 and 2x2 pivot cases)
239*
240 DO k = n, 1, -1
241 kp = abs( ipiv( k ) )
242 IF( kp.NE.k ) THEN
243 CALL cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
244 END IF
245 END DO
246*
247* Compute (U \P**T * B) -> B [ (U \P**T * B) ]
248*
249 CALL ctrsm( 'L', 'U', 'N', 'U', n, nrhs, one, a, lda, b, ldb )
250*
251* Compute D \ B -> B [ D \ (U \P**T * B) ]
252*
253 i = n
254 DO WHILE ( i.GE.1 )
255 IF( ipiv( i ).GT.0 ) THEN
256 CALL cscal( nrhs, one / a( i, i ), b( i, 1 ), ldb )
257 ELSE IF ( i.GT.1 ) THEN
258 akm1k = e( i )
259 akm1 = a( i-1, i-1 ) / akm1k
260 ak = a( i, i ) / akm1k
261 denom = akm1*ak - one
262 DO j = 1, nrhs
263 bkm1 = b( i-1, j ) / akm1k
264 bk = b( i, j ) / akm1k
265 b( i-1, j ) = ( ak*bkm1-bk ) / denom
266 b( i, j ) = ( akm1*bk-bkm1 ) / denom
267 END DO
268 i = i - 1
269 END IF
270 i = i - 1
271 END DO
272*
273* Compute (U**T \ B) -> B [ U**T \ (D \ (U \P**T * B) ) ]
274*
275 CALL ctrsm( 'L', 'U', 'T', 'U', n, nrhs, one, a, lda, b, ldb )
276*
277* P * B [ P * (U**T \ (D \ (U \P**T * B) )) ]
278*
279* Interchange rows K and IPIV(K) of matrix B in reverse order
280* from the formation order of IPIV(I) vector for Upper case.
281*
282* (We can do the simple loop over IPIV with increment 1,
283* since the ABS value of IPIV( I ) represents the row index
284* of the interchange with row i in both 1x1 and 2x2 pivot cases)
285*
286 DO k = 1, n, 1
287 kp = abs( ipiv( k ) )
288 IF( kp.NE.k ) THEN
289 CALL cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
290 END IF
291 END DO
292*
293 ELSE
294*
295* Begin Lower
296*
297* Solve A*X = B, where A = L*D*L**T.
298*
299* P**T * B
300* Interchange rows K and IPIV(K) of matrix B in the same order
301* that the formation order of IPIV(I) vector for Lower case.
302*
303* (We can do the simple loop over IPIV with increment 1,
304* since the ABS value of IPIV(I) represents the row index
305* of the interchange with row i in both 1x1 and 2x2 pivot cases)
306*
307 DO k = 1, n, 1
308 kp = abs( ipiv( k ) )
309 IF( kp.NE.k ) THEN
310 CALL cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
311 END IF
312 END DO
313*
314* Compute (L \P**T * B) -> B [ (L \P**T * B) ]
315*
316 CALL ctrsm( 'L', 'L', 'N', 'U', n, nrhs, one, a, lda, b, ldb )
317*
318* Compute D \ B -> B [ D \ (L \P**T * B) ]
319*
320 i = 1
321 DO WHILE ( i.LE.n )
322 IF( ipiv( i ).GT.0 ) THEN
323 CALL cscal( nrhs, one / a( i, i ), b( i, 1 ), ldb )
324 ELSE IF( i.LT.n ) THEN
325 akm1k = e( i )
326 akm1 = a( i, i ) / akm1k
327 ak = a( i+1, i+1 ) / akm1k
328 denom = akm1*ak - one
329 DO j = 1, nrhs
330 bkm1 = b( i, j ) / akm1k
331 bk = b( i+1, j ) / akm1k
332 b( i, j ) = ( ak*bkm1-bk ) / denom
333 b( i+1, j ) = ( akm1*bk-bkm1 ) / denom
334 END DO
335 i = i + 1
336 END IF
337 i = i + 1
338 END DO
339*
340* Compute (L**T \ B) -> B [ L**T \ (D \ (L \P**T * B) ) ]
341*
342 CALL ctrsm('L', 'L', 'T', 'U', n, nrhs, one, a, lda, b, ldb )
343*
344* P * B [ P * (L**T \ (D \ (L \P**T * B) )) ]
345*
346* Interchange rows K and IPIV(K) of matrix B in reverse order
347* from the formation order of IPIV(I) vector for Lower case.
348*
349* (We can do the simple loop over IPIV with decrement -1,
350* since the ABS value of IPIV(I) represents the row index
351* of the interchange with row i in both 1x1 and 2x2 pivot cases)
352*
353 DO k = n, 1, -1
354 kp = abs( ipiv( k ) )
355 IF( kp.NE.k ) THEN
356 CALL cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
357 END IF
358 END DO
359*
360* END Lower
361*
362 END IF
363*
364 RETURN
365*
366* End of CSYTRS_3
367*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
cscal
subroutine cscal(n, ca, cx, incx)
CSCAL
Definition cscal.f:78
cswap
subroutine cswap(n, cx, incx, cy, incy)
CSWAP
Definition cswap.f:81
ctrsm
subroutine ctrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
CTRSM
Definition ctrsm.f:180
Here is the call graph for this function:
Here is the caller graph for this function: