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

◆ 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*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
zdscal
subroutine zdscal(n, da, zx, incx)
ZDSCAL
Definition zdscal.f:78
zswap
subroutine zswap(n, zx, incx, zy, incy)
ZSWAP
Definition zswap.f:81
ztrsm
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: