LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ zhetrs_3()

subroutine zhetrs_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 
)

ZHETRS_3

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

Purpose:
 ZHETRS_3 solves a system of linear equations A * X = B with a complex
 Hermitian matrix A using the factorization computed
 by ZHETRF_RK or ZHETRF_BK:

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

 where U (or L) is unit upper (or lower) triangular matrix,
 U**H (or L**H) is the conjugate of U (or L), P is a permutation
 matrix, P**T is the transpose of P, and D is Hermitian 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**H)*(P**T);
          = 'L':  Lower triangular, form is A = P*L*D*(L**H)*(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 ZHETRF_RK and ZHETRF_BK:
            a) ONLY diagonal elements of the Hermitian 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 Hermitian 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 ZHETRF_RK or ZHETRF_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 zhetrs_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  DOUBLE PRECISION S
189  COMPLEX*16 AK, AKM1, AKM1K, BK, BKM1, DENOM
190 * ..
191 * .. External Functions ..
192  LOGICAL LSAME
193  EXTERNAL lsame
194 * ..
195 * .. External Subroutines ..
196  EXTERNAL zdscal, zswap, ztrsm, xerbla
197 * ..
198 * .. Intrinsic Functions ..
199  INTRINSIC abs, dble, dconjg, max
200 * ..
201 * .. Executable Statements ..
202 *
203  info = 0
204  upper = lsame( uplo, 'U' )
205  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
206  info = -1
207  ELSE IF( n.LT.0 ) THEN
208  info = -2
209  ELSE IF( nrhs.LT.0 ) THEN
210  info = -3
211  ELSE IF( lda.LT.max( 1, n ) ) THEN
212  info = -5
213  ELSE IF( ldb.LT.max( 1, n ) ) THEN
214  info = -9
215  END IF
216  IF( info.NE.0 ) THEN
217  CALL xerbla( 'ZHETRS_3', -info )
218  RETURN
219  END IF
220 *
221 * Quick return if possible
222 *
223  IF( n.EQ.0 .OR. nrhs.EQ.0 )
224  $ RETURN
225 *
226  IF( upper ) THEN
227 *
228 * Begin Upper
229 *
230 * Solve A*X = B, where A = U*D*U**H.
231 *
232 * P**T * B
233 *
234 * Interchange rows K and IPIV(K) of matrix B in the same order
235 * that the formation order of IPIV(I) vector for Upper case.
236 *
237 * (We can do the simple loop over IPIV with decrement -1,
238 * since the ABS value of IPIV(I) represents the row index
239 * of the interchange with row i in both 1x1 and 2x2 pivot cases)
240 *
241  DO k = n, 1, -1
242  kp = abs( ipiv( k ) )
243  IF( kp.NE.k ) THEN
244  CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
245  END IF
246  END DO
247 *
248 * Compute (U \P**T * B) -> B [ (U \P**T * B) ]
249 *
250  CALL ztrsm( 'L', 'U', 'N', 'U', n, nrhs, one, a, lda, b, ldb )
251 *
252 * Compute D \ B -> B [ D \ (U \P**T * B) ]
253 *
254  i = n
255  DO WHILE ( i.GE.1 )
256  IF( ipiv( i ).GT.0 ) THEN
257  s = dble( one ) / dble( a( i, i ) )
258  CALL zdscal( nrhs, s, b( i, 1 ), ldb )
259  ELSE IF ( i.GT.1 ) THEN
260  akm1k = e( i )
261  akm1 = a( i-1, i-1 ) / akm1k
262  ak = a( i, i ) / dconjg( akm1k )
263  denom = akm1*ak - one
264  DO j = 1, nrhs
265  bkm1 = b( i-1, j ) / akm1k
266  bk = b( i, j ) / dconjg( akm1k )
267  b( i-1, j ) = ( ak*bkm1-bk ) / denom
268  b( i, j ) = ( akm1*bk-bkm1 ) / denom
269  END DO
270  i = i - 1
271  END IF
272  i = i - 1
273  END DO
274 *
275 * Compute (U**H \ B) -> B [ U**H \ (D \ (U \P**T * B) ) ]
276 *
277  CALL ztrsm( 'L', 'U', 'C', 'U', n, nrhs, one, a, lda, b, ldb )
278 *
279 * P * B [ P * (U**H \ (D \ (U \P**T * B) )) ]
280 *
281 * Interchange rows K and IPIV(K) of matrix B in reverse order
282 * from the formation order of IPIV(I) vector for Upper case.
283 *
284 * (We can do the simple loop over IPIV with increment 1,
285 * since the ABS value of IPIV(I) represents the row index
286 * of the interchange with row i in both 1x1 and 2x2 pivot cases)
287 *
288  DO k = 1, n, 1
289  kp = abs( ipiv( k ) )
290  IF( kp.NE.k ) THEN
291  CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
292  END IF
293  END DO
294 *
295  ELSE
296 *
297 * Begin Lower
298 *
299 * Solve A*X = B, where A = L*D*L**H.
300 *
301 * P**T * B
302 * Interchange rows K and IPIV(K) of matrix B in the same order
303 * that the formation order of IPIV(I) vector for Lower case.
304 *
305 * (We can do the simple loop over IPIV with increment 1,
306 * since the ABS value of IPIV(I) represents the row index
307 * of the interchange with row i in both 1x1 and 2x2 pivot cases)
308 *
309  DO k = 1, n, 1
310  kp = abs( ipiv( k ) )
311  IF( kp.NE.k ) THEN
312  CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
313  END IF
314  END DO
315 *
316 * Compute (L \P**T * B) -> B [ (L \P**T * B) ]
317 *
318  CALL ztrsm( 'L', 'L', 'N', 'U', n, nrhs, one, a, lda, b, ldb )
319 *
320 * Compute D \ B -> B [ D \ (L \P**T * B) ]
321 *
322  i = 1
323  DO WHILE ( i.LE.n )
324  IF( ipiv( i ).GT.0 ) THEN
325  s = dble( one ) / dble( a( i, i ) )
326  CALL zdscal( nrhs, s, b( i, 1 ), ldb )
327  ELSE IF( i.LT.n ) THEN
328  akm1k = e( i )
329  akm1 = a( i, i ) / dconjg( akm1k )
330  ak = a( i+1, i+1 ) / akm1k
331  denom = akm1*ak - one
332  DO j = 1, nrhs
333  bkm1 = b( i, j ) / dconjg( akm1k )
334  bk = b( i+1, j ) / akm1k
335  b( i, j ) = ( ak*bkm1-bk ) / denom
336  b( i+1, j ) = ( akm1*bk-bkm1 ) / denom
337  END DO
338  i = i + 1
339  END IF
340  i = i + 1
341  END DO
342 *
343 * Compute (L**H \ B) -> B [ L**H \ (D \ (L \P**T * B) ) ]
344 *
345  CALL ztrsm('L', 'L', 'C', 'U', n, nrhs, one, a, lda, b, ldb )
346 *
347 * P * B [ P * (L**H \ (D \ (L \P**T * B) )) ]
348 *
349 * Interchange rows K and IPIV(K) of matrix B in reverse order
350 * from the formation order of IPIV(I) vector for Lower case.
351 *
352 * (We can do the simple loop over IPIV with decrement -1,
353 * since the ABS value of IPIV(I) represents the row index
354 * of the interchange with row i in both 1x1 and 2x2 pivot cases)
355 *
356  DO k = n, 1, -1
357  kp = abs( ipiv( k ) )
358  IF( kp.NE.k ) THEN
359  CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
360  END IF
361  END DO
362 *
363 * END Lower
364 *
365  END IF
366 *
367  RETURN
368 *
369 * End of ZHETRS_3
370 *
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 zdscal(N, DA, ZX, INCX)
ZDSCAL
Definition: zdscal.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: