LAPACK 3.12.1
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 161 of file zhetrs_3.f.

163*
164* -- LAPACK computational routine --
165* -- LAPACK is a software package provided by Univ. of Tennessee, --
166* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
167*
168* .. Scalar Arguments ..
169 CHARACTER UPLO
170 INTEGER INFO, LDA, LDB, N, NRHS
171* ..
172* .. Array Arguments ..
173 INTEGER IPIV( * )
174 COMPLEX*16 A( LDA, * ), B( LDB, * ), E( * )
175* ..
176*
177* =====================================================================
178*
179* .. Parameters ..
180 COMPLEX*16 ONE
181 parameter( one = ( 1.0d+0,0.0d+0 ) )
182* ..
183* .. Local Scalars ..
184 LOGICAL UPPER
185 INTEGER I, J, K, KP
186 DOUBLE PRECISION S
187 COMPLEX*16 AK, AKM1, AKM1K, BK, BKM1, DENOM
188* ..
189* .. External Functions ..
190 LOGICAL LSAME
191 EXTERNAL lsame
192* ..
193* .. External Subroutines ..
194 EXTERNAL zdscal, zswap, ztrsm, xerbla
195* ..
196* .. Intrinsic Functions ..
197 INTRINSIC abs, dble, dconjg, max
198* ..
199* .. Executable Statements ..
200*
201 info = 0
202 upper = lsame( uplo, 'U' )
203 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
204 info = -1
205 ELSE IF( n.LT.0 ) THEN
206 info = -2
207 ELSE IF( nrhs.LT.0 ) THEN
208 info = -3
209 ELSE IF( lda.LT.max( 1, n ) ) THEN
210 info = -5
211 ELSE IF( ldb.LT.max( 1, n ) ) THEN
212 info = -9
213 END IF
214 IF( info.NE.0 ) THEN
215 CALL xerbla( 'ZHETRS_3', -info )
216 RETURN
217 END IF
218*
219* Quick return if possible
220*
221 IF( n.EQ.0 .OR. nrhs.EQ.0 )
222 $ RETURN
223*
224 IF( upper ) THEN
225*
226* Begin Upper
227*
228* Solve A*X = B, where A = U*D*U**H.
229*
230* P**T * B
231*
232* Interchange rows K and IPIV(K) of matrix B in the same order
233* that the formation order of IPIV(I) vector for Upper case.
234*
235* (We can do the simple loop over IPIV with decrement -1,
236* since the ABS value of IPIV(I) represents the row index
237* of the interchange with row i in both 1x1 and 2x2 pivot cases)
238*
239 DO k = n, 1, -1
240 kp = abs( ipiv( k ) )
241 IF( kp.NE.k ) THEN
242 CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
243 END IF
244 END DO
245*
246* Compute (U \P**T * B) -> B [ (U \P**T * B) ]
247*
248 CALL ztrsm( 'L', 'U', 'N', 'U', n, nrhs, one, a, lda, b,
249 $ 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 s = dble( one ) / dble( a( i, i ) )
257 CALL zdscal( nrhs, s, b( i, 1 ), ldb )
258 ELSE IF ( i.GT.1 ) THEN
259 akm1k = e( i )
260 akm1 = a( i-1, i-1 ) / akm1k
261 ak = a( i, i ) / dconjg( akm1k )
262 denom = akm1*ak - one
263 DO j = 1, nrhs
264 bkm1 = b( i-1, j ) / akm1k
265 bk = b( i, j ) / dconjg( akm1k )
266 b( i-1, j ) = ( ak*bkm1-bk ) / denom
267 b( i, j ) = ( akm1*bk-bkm1 ) / denom
268 END DO
269 i = i - 1
270 END IF
271 i = i - 1
272 END DO
273*
274* Compute (U**H \ B) -> B [ U**H \ (D \ (U \P**T * B) ) ]
275*
276 CALL ztrsm( 'L', 'U', 'C', 'U', n, nrhs, one, a, lda, b,
277 $ 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,
319 $ ldb )
320*
321* Compute D \ B -> B [ D \ (L \P**T * B) ]
322*
323 i = 1
324 DO WHILE ( i.LE.n )
325 IF( ipiv( i ).GT.0 ) THEN
326 s = dble( one ) / dble( a( i, i ) )
327 CALL zdscal( nrhs, s, b( i, 1 ), ldb )
328 ELSE IF( i.LT.n ) THEN
329 akm1k = e( i )
330 akm1 = a( i, i ) / dconjg( akm1k )
331 ak = a( i+1, i+1 ) / akm1k
332 denom = akm1*ak - one
333 DO j = 1, nrhs
334 bkm1 = b( i, j ) / dconjg( akm1k )
335 bk = b( i+1, j ) / akm1k
336 b( i, j ) = ( ak*bkm1-bk ) / denom
337 b( i+1, j ) = ( akm1*bk-bkm1 ) / denom
338 END DO
339 i = i + 1
340 END IF
341 i = i + 1
342 END DO
343*
344* Compute (L**H \ B) -> B [ L**H \ (D \ (L \P**T * B) ) ]
345*
346 CALL ztrsm('L', 'L', 'C', 'U', n, nrhs, one, a, lda, b,
347 $ ldb )
348*
349* P * B [ P * (L**H \ (D \ (L \P**T * B) )) ]
350*
351* Interchange rows K and IPIV(K) of matrix B in reverse order
352* from the formation order of IPIV(I) vector for Lower case.
353*
354* (We can do the simple loop over IPIV with decrement -1,
355* since the ABS value of IPIV(I) represents the row index
356* of the interchange with row i in both 1x1 and 2x2 pivot cases)
357*
358 DO k = n, 1, -1
359 kp = abs( ipiv( k ) )
360 IF( kp.NE.k ) THEN
361 CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
362 END IF
363 END DO
364*
365* END Lower
366*
367 END IF
368*
369 RETURN
370*
371* End of ZHETRS_3
372*
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: