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

◆ dsptrs()

subroutine dsptrs ( character uplo,
integer n,
integer nrhs,
double precision, dimension( * ) ap,
integer, dimension( * ) ipiv,
double precision, dimension( ldb, * ) b,
integer ldb,
integer info )

DSPTRS

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

Purpose:
!>
!> DSPTRS solves a system of linear equations A*X = B with a real
!> symmetric matrix A stored in packed format using the factorization
!> A = U*D*U**T or A = L*D*L**T computed by DSPTRF.
!>
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**T;
!>          = 'L':  Lower triangular, form is A = L*D*L**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]AP
!>          AP is DOUBLE PRECISION 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 DSPTRF, 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 DSPTRF.
!>
[in,out]B
!>          B is DOUBLE PRECISION 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 112 of file dsptrs.f.

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