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

◆ zhptrs()

subroutine zhptrs ( character  uplo,
integer  n,
integer  nrhs,
complex*16, dimension( * )  ap,
integer, dimension( * )  ipiv,
complex*16, dimension( ldb, * )  b,
integer  ldb,
integer  info 
)

ZHPTRS

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

Purpose:
 ZHPTRS solves a system of linear equations A*X = B with a complex
 Hermitian matrix A stored in packed format using the factorization
 A = U*D*U**H or A = L*D*L**H computed by ZHPTRF.
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 = U*D*U**H;
          = 'L':  Lower triangular, form is A = L*D*L**H.
[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]AP
          AP is COMPLEX*16 array, dimension (N*(N+1)/2)
          The block diagonal matrix D and the multipliers used to
          obtain the factor U or L as computed by ZHPTRF, stored as a
          packed triangular matrix.
[in]IPIV
          IPIV is INTEGER array, dimension (N)
          Details of the interchanges and the block structure of D
          as determined by ZHPTRF.
[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.

Definition at line 114 of file zhptrs.f.

115*
116* -- LAPACK computational routine --
117* -- LAPACK is a software package provided by Univ. of Tennessee, --
118* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
119*
120* .. Scalar Arguments ..
121 CHARACTER UPLO
122 INTEGER INFO, LDB, N, NRHS
123* ..
124* .. Array Arguments ..
125 INTEGER IPIV( * )
126 COMPLEX*16 AP( * ), B( LDB, * )
127* ..
128*
129* =====================================================================
130*
131* .. Parameters ..
132 COMPLEX*16 ONE
133 parameter( one = ( 1.0d+0, 0.0d+0 ) )
134* ..
135* .. Local Scalars ..
136 LOGICAL UPPER
137 INTEGER J, K, KC, KP
138 DOUBLE PRECISION S
139 COMPLEX*16 AK, AKM1, AKM1K, BK, BKM1, DENOM
140* ..
141* .. External Functions ..
142 LOGICAL LSAME
143 EXTERNAL lsame
144* ..
145* .. External Subroutines ..
146 EXTERNAL xerbla, zdscal, zgemv, zgeru, zlacgv, zswap
147* ..
148* .. Intrinsic Functions ..
149 INTRINSIC dble, dconjg, max
150* ..
151* .. Executable Statements ..
152*
153 info = 0
154 upper = lsame( uplo, 'U' )
155 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
156 info = -1
157 ELSE IF( n.LT.0 ) THEN
158 info = -2
159 ELSE IF( nrhs.LT.0 ) THEN
160 info = -3
161 ELSE IF( ldb.LT.max( 1, n ) ) THEN
162 info = -7
163 END IF
164 IF( info.NE.0 ) THEN
165 CALL xerbla( 'ZHPTRS', -info )
166 RETURN
167 END IF
168*
169* Quick return if possible
170*
171 IF( n.EQ.0 .OR. nrhs.EQ.0 )
172 $ RETURN
173*
174 IF( upper ) THEN
175*
176* Solve A*X = B, where A = U*D*U**H.
177*
178* First solve U*D*X = B, overwriting B with X.
179*
180* K is the main loop index, decreasing from N to 1 in steps of
181* 1 or 2, depending on the size of the diagonal blocks.
182*
183 k = n
184 kc = n*( n+1 ) / 2 + 1
185 10 CONTINUE
186*
187* If K < 1, exit from loop.
188*
189 IF( k.LT.1 )
190 $ GO TO 30
191*
192 kc = kc - k
193 IF( ipiv( k ).GT.0 ) THEN
194*
195* 1 x 1 diagonal block
196*
197* Interchange rows K and IPIV(K).
198*
199 kp = ipiv( k )
200 IF( kp.NE.k )
201 $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
202*
203* Multiply by inv(U(K)), where U(K) is the transformation
204* stored in column K of A.
205*
206 CALL zgeru( k-1, nrhs, -one, ap( kc ), 1, b( k, 1 ), ldb,
207 $ b( 1, 1 ), ldb )
208*
209* Multiply by the inverse of the diagonal block.
210*
211 s = dble( one ) / dble( ap( kc+k-1 ) )
212 CALL zdscal( nrhs, s, b( k, 1 ), ldb )
213 k = k - 1
214 ELSE
215*
216* 2 x 2 diagonal block
217*
218* Interchange rows K-1 and -IPIV(K).
219*
220 kp = -ipiv( k )
221 IF( kp.NE.k-1 )
222 $ CALL zswap( nrhs, b( k-1, 1 ), ldb, b( kp, 1 ), ldb )
223*
224* Multiply by inv(U(K)), where U(K) is the transformation
225* stored in columns K-1 and K of A.
226*
227 CALL zgeru( k-2, nrhs, -one, ap( kc ), 1, b( k, 1 ), ldb,
228 $ b( 1, 1 ), ldb )
229 CALL zgeru( k-2, nrhs, -one, ap( kc-( k-1 ) ), 1,
230 $ b( k-1, 1 ), ldb, b( 1, 1 ), ldb )
231*
232* Multiply by the inverse of the diagonal block.
233*
234 akm1k = ap( kc+k-2 )
235 akm1 = ap( kc-1 ) / akm1k
236 ak = ap( kc+k-1 ) / dconjg( akm1k )
237 denom = akm1*ak - one
238 DO 20 j = 1, nrhs
239 bkm1 = b( k-1, j ) / akm1k
240 bk = b( k, j ) / dconjg( akm1k )
241 b( k-1, j ) = ( ak*bkm1-bk ) / denom
242 b( k, j ) = ( akm1*bk-bkm1 ) / denom
243 20 CONTINUE
244 kc = kc - k + 1
245 k = k - 2
246 END IF
247*
248 GO TO 10
249 30 CONTINUE
250*
251* Next solve U**H *X = B, overwriting B with X.
252*
253* K is the main loop index, increasing from 1 to N in steps of
254* 1 or 2, depending on the size of the diagonal blocks.
255*
256 k = 1
257 kc = 1
258 40 CONTINUE
259*
260* If K > N, exit from loop.
261*
262 IF( k.GT.n )
263 $ GO TO 50
264*
265 IF( ipiv( k ).GT.0 ) THEN
266*
267* 1 x 1 diagonal block
268*
269* Multiply by inv(U**H(K)), where U(K) is the transformation
270* stored in column K of A.
271*
272 IF( k.GT.1 ) THEN
273 CALL zlacgv( nrhs, b( k, 1 ), ldb )
274 CALL zgemv( 'Conjugate transpose', k-1, nrhs, -one, b,
275 $ ldb, ap( kc ), 1, one, b( k, 1 ), ldb )
276 CALL zlacgv( nrhs, b( k, 1 ), ldb )
277 END IF
278*
279* Interchange rows K and IPIV(K).
280*
281 kp = ipiv( k )
282 IF( kp.NE.k )
283 $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
284 kc = kc + k
285 k = k + 1
286 ELSE
287*
288* 2 x 2 diagonal block
289*
290* Multiply by inv(U**H(K+1)), where U(K+1) is the transformation
291* stored in columns K and K+1 of A.
292*
293 IF( k.GT.1 ) THEN
294 CALL zlacgv( nrhs, b( k, 1 ), ldb )
295 CALL zgemv( 'Conjugate transpose', k-1, nrhs, -one, b,
296 $ ldb, ap( kc ), 1, one, b( k, 1 ), ldb )
297 CALL zlacgv( nrhs, b( k, 1 ), ldb )
298*
299 CALL zlacgv( nrhs, b( k+1, 1 ), ldb )
300 CALL zgemv( 'Conjugate transpose', k-1, nrhs, -one, b,
301 $ ldb, ap( kc+k ), 1, one, b( k+1, 1 ), ldb )
302 CALL zlacgv( nrhs, b( k+1, 1 ), ldb )
303 END IF
304*
305* Interchange rows K and -IPIV(K).
306*
307 kp = -ipiv( k )
308 IF( kp.NE.k )
309 $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
310 kc = kc + 2*k + 1
311 k = k + 2
312 END IF
313*
314 GO TO 40
315 50 CONTINUE
316*
317 ELSE
318*
319* Solve A*X = B, where A = L*D*L**H.
320*
321* First solve L*D*X = B, overwriting B with X.
322*
323* K is the main loop index, increasing from 1 to N in steps of
324* 1 or 2, depending on the size of the diagonal blocks.
325*
326 k = 1
327 kc = 1
328 60 CONTINUE
329*
330* If K > N, exit from loop.
331*
332 IF( k.GT.n )
333 $ GO TO 80
334*
335 IF( ipiv( k ).GT.0 ) THEN
336*
337* 1 x 1 diagonal block
338*
339* Interchange rows K and IPIV(K).
340*
341 kp = ipiv( k )
342 IF( kp.NE.k )
343 $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
344*
345* Multiply by inv(L(K)), where L(K) is the transformation
346* stored in column K of A.
347*
348 IF( k.LT.n )
349 $ CALL zgeru( n-k, nrhs, -one, ap( kc+1 ), 1, b( k, 1 ),
350 $ ldb, b( k+1, 1 ), ldb )
351*
352* Multiply by the inverse of the diagonal block.
353*
354 s = dble( one ) / dble( ap( kc ) )
355 CALL zdscal( nrhs, s, b( k, 1 ), ldb )
356 kc = kc + n - k + 1
357 k = k + 1
358 ELSE
359*
360* 2 x 2 diagonal block
361*
362* Interchange rows K+1 and -IPIV(K).
363*
364 kp = -ipiv( k )
365 IF( kp.NE.k+1 )
366 $ CALL zswap( nrhs, b( k+1, 1 ), ldb, b( kp, 1 ), ldb )
367*
368* Multiply by inv(L(K)), where L(K) is the transformation
369* stored in columns K and K+1 of A.
370*
371 IF( k.LT.n-1 ) THEN
372 CALL zgeru( n-k-1, nrhs, -one, ap( kc+2 ), 1, b( k, 1 ),
373 $ ldb, b( k+2, 1 ), ldb )
374 CALL zgeru( n-k-1, nrhs, -one, ap( kc+n-k+2 ), 1,
375 $ b( k+1, 1 ), ldb, b( k+2, 1 ), ldb )
376 END IF
377*
378* Multiply by the inverse of the diagonal block.
379*
380 akm1k = ap( kc+1 )
381 akm1 = ap( kc ) / dconjg( akm1k )
382 ak = ap( kc+n-k+1 ) / akm1k
383 denom = akm1*ak - one
384 DO 70 j = 1, nrhs
385 bkm1 = b( k, j ) / dconjg( akm1k )
386 bk = b( k+1, j ) / akm1k
387 b( k, j ) = ( ak*bkm1-bk ) / denom
388 b( k+1, j ) = ( akm1*bk-bkm1 ) / denom
389 70 CONTINUE
390 kc = kc + 2*( n-k ) + 1
391 k = k + 2
392 END IF
393*
394 GO TO 60
395 80 CONTINUE
396*
397* Next solve L**H *X = B, overwriting B with X.
398*
399* K is the main loop index, decreasing from N to 1 in steps of
400* 1 or 2, depending on the size of the diagonal blocks.
401*
402 k = n
403 kc = n*( n+1 ) / 2 + 1
404 90 CONTINUE
405*
406* If K < 1, exit from loop.
407*
408 IF( k.LT.1 )
409 $ GO TO 100
410*
411 kc = kc - ( n-k+1 )
412 IF( ipiv( k ).GT.0 ) THEN
413*
414* 1 x 1 diagonal block
415*
416* Multiply by inv(L**H(K)), where L(K) is the transformation
417* stored in column K of A.
418*
419 IF( k.LT.n ) THEN
420 CALL zlacgv( nrhs, b( k, 1 ), ldb )
421 CALL zgemv( 'Conjugate transpose', n-k, nrhs, -one,
422 $ b( k+1, 1 ), ldb, ap( kc+1 ), 1, one,
423 $ b( k, 1 ), ldb )
424 CALL zlacgv( nrhs, b( k, 1 ), ldb )
425 END IF
426*
427* Interchange rows K and IPIV(K).
428*
429 kp = ipiv( k )
430 IF( kp.NE.k )
431 $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
432 k = k - 1
433 ELSE
434*
435* 2 x 2 diagonal block
436*
437* Multiply by inv(L**H(K-1)), where L(K-1) is the transformation
438* stored in columns K-1 and K of A.
439*
440 IF( k.LT.n ) THEN
441 CALL zlacgv( nrhs, b( k, 1 ), ldb )
442 CALL zgemv( 'Conjugate transpose', n-k, nrhs, -one,
443 $ b( k+1, 1 ), ldb, ap( kc+1 ), 1, one,
444 $ b( k, 1 ), ldb )
445 CALL zlacgv( nrhs, b( k, 1 ), ldb )
446*
447 CALL zlacgv( nrhs, b( k-1, 1 ), ldb )
448 CALL zgemv( 'Conjugate transpose', n-k, nrhs, -one,
449 $ b( k+1, 1 ), ldb, ap( kc-( n-k ) ), 1, one,
450 $ b( k-1, 1 ), ldb )
451 CALL zlacgv( nrhs, b( k-1, 1 ), ldb )
452 END IF
453*
454* Interchange rows K and -IPIV(K).
455*
456 kp = -ipiv( k )
457 IF( kp.NE.k )
458 $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
459 kc = kc - ( n-k+2 )
460 k = k - 2
461 END IF
462*
463 GO TO 90
464 100 CONTINUE
465 END IF
466*
467 RETURN
468*
469* End of ZHPTRS
470*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
zgemv
subroutine zgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
ZGEMV
Definition zgemv.f:160
zgeru
subroutine zgeru(m, n, alpha, x, incx, y, incy, a, lda)
ZGERU
Definition zgeru.f:130
zlacgv
subroutine zlacgv(n, x, incx)
ZLACGV conjugates a complex vector.
Definition zlacgv.f:74
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
Here is the call graph for this function:
Here is the caller graph for this function: