LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ zsytrs_3()

subroutine zsytrs_3 ( character  UPLO,
integer  N,
integer  NRHS,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( * )  E,
integer, dimension( * )  IPIV,
complex*16, dimension( ldb, * )  B,
integer  LDB,
integer  INFO 
)

ZSYTRS_3

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

Purpose:
 ZSYTRS_3 solves a system of linear equations A * X = B with a complex
 symmetric matrix A using the factorization computed
 by ZSYTRF_RK or ZSYTRF_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*16 array, dimension (LDA,N)
          Diagonal of the block diagonal matrix D and factors U or L
          as computed by ZSYTRF_RK and ZSYTRF_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*16 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 ZSYTRF_RK or ZSYTRF_BK.
[in,out]B
          B is COMPLEX*16 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 zsytrs_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*16 A( LDA, * ), B( LDB, * ), E( * )
177 * ..
178 *
179 * =====================================================================
180 *
181 * .. Parameters ..
182  COMPLEX*16 ONE
183  parameter( one = ( 1.0d+0,0.0d+0 ) )
184 * ..
185 * .. Local Scalars ..
186  LOGICAL UPPER
187  INTEGER I, J, K, KP
188  COMPLEX*16 AK, AKM1, AKM1K, BK, BKM1, DENOM
189 * ..
190 * .. External Functions ..
191  LOGICAL LSAME
192  EXTERNAL lsame
193 * ..
194 * .. External Subroutines ..
195  EXTERNAL zscal, zswap, ztrsm, 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( 'ZSYTRS_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 zswap( 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 ztrsm( '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 zscal( 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 ztrsm( '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 zswap( 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 zswap( 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 ztrsm( '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 zscal( 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 ztrsm('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 zswap( 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 ZSYTRS_3
367 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:81
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:78
subroutine ztrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
ZTRSM
Definition: ztrsm.f:180
Here is the call graph for this function:
Here is the caller graph for this function: