LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dsptrs.f
Go to the documentation of this file.
1*> \brief \b DSPTRS
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DSPTRS + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsptrs.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsptrs.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsptrs.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DSPTRS( UPLO, N, NRHS, AP, IPIV, B, LDB, INFO )
22*
23* .. Scalar Arguments ..
24* CHARACTER UPLO
25* INTEGER INFO, LDB, N, NRHS
26* ..
27* .. Array Arguments ..
28* INTEGER IPIV( * )
29* DOUBLE PRECISION AP( * ), B( LDB, * )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> DSPTRS solves a system of linear equations A*X = B with a real
39*> symmetric matrix A stored in packed format using the factorization
40*> A = U*D*U**T or A = L*D*L**T computed by DSPTRF.
41*> \endverbatim
42*
43* Arguments:
44* ==========
45*
46*> \param[in] UPLO
47*> \verbatim
48*> UPLO is CHARACTER*1
49*> Specifies whether the details of the factorization are stored
50*> as an upper or lower triangular matrix.
51*> = 'U': Upper triangular, form is A = U*D*U**T;
52*> = 'L': Lower triangular, form is A = L*D*L**T.
53*> \endverbatim
54*>
55*> \param[in] N
56*> \verbatim
57*> N is INTEGER
58*> The order of the matrix A. N >= 0.
59*> \endverbatim
60*>
61*> \param[in] NRHS
62*> \verbatim
63*> NRHS is INTEGER
64*> The number of right hand sides, i.e., the number of columns
65*> of the matrix B. NRHS >= 0.
66*> \endverbatim
67*>
68*> \param[in] AP
69*> \verbatim
70*> AP is DOUBLE PRECISION array, dimension (N*(N+1)/2)
71*> The block diagonal matrix D and the multipliers used to
72*> obtain the factor U or L as computed by DSPTRF, stored as a
73*> packed triangular matrix.
74*> \endverbatim
75*>
76*> \param[in] IPIV
77*> \verbatim
78*> IPIV is INTEGER array, dimension (N)
79*> Details of the interchanges and the block structure of D
80*> as determined by DSPTRF.
81*> \endverbatim
82*>
83*> \param[in,out] B
84*> \verbatim
85*> B is DOUBLE PRECISION array, dimension (LDB,NRHS)
86*> On entry, the right hand side matrix B.
87*> On exit, the solution matrix X.
88*> \endverbatim
89*>
90*> \param[in] LDB
91*> \verbatim
92*> LDB is INTEGER
93*> The leading dimension of the array B. LDB >= max(1,N).
94*> \endverbatim
95*>
96*> \param[out] INFO
97*> \verbatim
98*> INFO is INTEGER
99*> = 0: successful exit
100*> < 0: if INFO = -i, the i-th argument had an illegal value
101*> \endverbatim
102*
103* Authors:
104* ========
105*
106*> \author Univ. of Tennessee
107*> \author Univ. of California Berkeley
108*> \author Univ. of Colorado Denver
109*> \author NAG Ltd.
110*
111*> \ingroup hptrs
112*
113* =====================================================================
114 SUBROUTINE dsptrs( UPLO, N, NRHS, AP, IPIV, B, LDB, INFO )
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 DOUBLE PRECISION AP( * ), B( LDB, * )
127* ..
128*
129* =====================================================================
130*
131* .. Parameters ..
132 DOUBLE PRECISION ONE
133 parameter( one = 1.0d+0 )
134* ..
135* .. Local Scalars ..
136 LOGICAL UPPER
137 INTEGER J, K, KC, KP
138 DOUBLE PRECISION AK, AKM1, AKM1K, BK, BKM1, DENOM
139* ..
140* .. External Functions ..
141 LOGICAL LSAME
142 EXTERNAL lsame
143* ..
144* .. External Subroutines ..
145 EXTERNAL dgemv, dger, dscal, dswap, xerbla
146* ..
147* .. Intrinsic Functions ..
148 INTRINSIC max
149* ..
150* .. Executable Statements ..
151*
152 info = 0
153 upper = lsame( uplo, 'U' )
154 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
155 info = -1
156 ELSE IF( n.LT.0 ) THEN
157 info = -2
158 ELSE IF( nrhs.LT.0 ) THEN
159 info = -3
160 ELSE IF( ldb.LT.max( 1, n ) ) THEN
161 info = -7
162 END IF
163 IF( info.NE.0 ) THEN
164 CALL xerbla( 'DSPTRS', -info )
165 RETURN
166 END IF
167*
168* Quick return if possible
169*
170 IF( n.EQ.0 .OR. nrhs.EQ.0 )
171 $ RETURN
172*
173 IF( upper ) THEN
174*
175* Solve A*X = B, where A = U*D*U**T.
176*
177* First solve U*D*X = B, overwriting B with X.
178*
179* K is the main loop index, decreasing from N to 1 in steps of
180* 1 or 2, depending on the size of the diagonal blocks.
181*
182 k = n
183 kc = n*( n+1 ) / 2 + 1
184 10 CONTINUE
185*
186* If K < 1, exit from loop.
187*
188 IF( k.LT.1 )
189 $ GO TO 30
190*
191 kc = kc - k
192 IF( ipiv( k ).GT.0 ) THEN
193*
194* 1 x 1 diagonal block
195*
196* Interchange rows K and IPIV(K).
197*
198 kp = ipiv( k )
199 IF( kp.NE.k )
200 $ CALL dswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
201*
202* Multiply by inv(U(K)), where U(K) is the transformation
203* stored in column K of A.
204*
205 CALL dger( k-1, nrhs, -one, ap( kc ), 1, b( k, 1 ), ldb,
206 $ b( 1, 1 ), ldb )
207*
208* Multiply by the inverse of the diagonal block.
209*
210 CALL dscal( nrhs, one / ap( kc+k-1 ), b( k, 1 ), ldb )
211 k = k - 1
212 ELSE
213*
214* 2 x 2 diagonal block
215*
216* Interchange rows K-1 and -IPIV(K).
217*
218 kp = -ipiv( k )
219 IF( kp.NE.k-1 )
220 $ CALL dswap( nrhs, b( k-1, 1 ), ldb, b( kp, 1 ), ldb )
221*
222* Multiply by inv(U(K)), where U(K) is the transformation
223* stored in columns K-1 and K of A.
224*
225 CALL dger( k-2, nrhs, -one, ap( kc ), 1, b( k, 1 ), ldb,
226 $ b( 1, 1 ), ldb )
227 CALL dger( k-2, nrhs, -one, ap( kc-( k-1 ) ), 1,
228 $ b( k-1, 1 ), ldb, b( 1, 1 ), ldb )
229*
230* Multiply by the inverse of the diagonal block.
231*
232 akm1k = ap( kc+k-2 )
233 akm1 = ap( kc-1 ) / akm1k
234 ak = ap( kc+k-1 ) / akm1k
235 denom = akm1*ak - one
236 DO 20 j = 1, nrhs
237 bkm1 = b( k-1, j ) / akm1k
238 bk = b( k, j ) / akm1k
239 b( k-1, j ) = ( ak*bkm1-bk ) / denom
240 b( k, j ) = ( akm1*bk-bkm1 ) / denom
241 20 CONTINUE
242 kc = kc - k + 1
243 k = k - 2
244 END IF
245*
246 GO TO 10
247 30 CONTINUE
248*
249* Next solve U**T*X = B, overwriting B with X.
250*
251* K is the main loop index, increasing from 1 to N in steps of
252* 1 or 2, depending on the size of the diagonal blocks.
253*
254 k = 1
255 kc = 1
256 40 CONTINUE
257*
258* If K > N, exit from loop.
259*
260 IF( k.GT.n )
261 $ GO TO 50
262*
263 IF( ipiv( k ).GT.0 ) THEN
264*
265* 1 x 1 diagonal block
266*
267* Multiply by inv(U**T(K)), where U(K) is the transformation
268* stored in column K of A.
269*
270 CALL dgemv( 'Transpose', k-1, nrhs, -one, b, ldb, ap( kc ),
271 $ 1, one, b( k, 1 ), ldb )
272*
273* Interchange rows K and IPIV(K).
274*
275 kp = ipiv( k )
276 IF( kp.NE.k )
277 $ CALL dswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
278 kc = kc + k
279 k = k + 1
280 ELSE
281*
282* 2 x 2 diagonal block
283*
284* Multiply by inv(U**T(K+1)), where U(K+1) is the transformation
285* stored in columns K and K+1 of A.
286*
287 CALL dgemv( 'Transpose', k-1, nrhs, -one, b, ldb, 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, 1 ),
359 $ ldb, b( k+2, 1 ), ldb )
360 CALL dger( n-k-1, nrhs, -one, ap( kc+n-k+2 ), 1,
361 $ b( k+1, 1 ), ldb, b( k+2, 1 ), ldb )
362 END IF
363*
364* Multiply by the inverse of the diagonal block.
365*
366 akm1k = ap( kc+1 )
367 akm1 = ap( kc ) / akm1k
368 ak = ap( kc+n-k+1 ) / akm1k
369 denom = akm1*ak - one
370 DO 70 j = 1, nrhs
371 bkm1 = b( k, j ) / akm1k
372 bk = b( k+1, j ) / akm1k
373 b( k, j ) = ( ak*bkm1-bk ) / denom
374 b( k+1, j ) = ( akm1*bk-bkm1 ) / denom
375 70 CONTINUE
376 kc = kc + 2*( n-k ) + 1
377 k = k + 2
378 END IF
379*
380 GO TO 60
381 80 CONTINUE
382*
383* Next solve L**T*X = B, overwriting B with X.
384*
385* K is the main loop index, decreasing from N to 1 in steps of
386* 1 or 2, depending on the size of the diagonal blocks.
387*
388 k = n
389 kc = n*( n+1 ) / 2 + 1
390 90 CONTINUE
391*
392* If K < 1, exit from loop.
393*
394 IF( k.LT.1 )
395 $ GO TO 100
396*
397 kc = kc - ( n-k+1 )
398 IF( ipiv( k ).GT.0 ) THEN
399*
400* 1 x 1 diagonal block
401*
402* Multiply by inv(L**T(K)), where L(K) is the transformation
403* stored in column K of A.
404*
405 IF( k.LT.n )
406 $ CALL dgemv( 'Transpose', n-k, nrhs, -one, b( k+1, 1 ),
407 $ ldb, ap( kc+1 ), 1, one, b( k, 1 ), ldb )
408*
409* Interchange rows K and IPIV(K).
410*
411 kp = ipiv( k )
412 IF( kp.NE.k )
413 $ CALL dswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
414 k = k - 1
415 ELSE
416*
417* 2 x 2 diagonal block
418*
419* Multiply by inv(L**T(K-1)), where L(K-1) is the transformation
420* stored in columns K-1 and K of A.
421*
422 IF( k.LT.n ) THEN
423 CALL dgemv( 'Transpose', n-k, nrhs, -one, b( k+1, 1 ),
424 $ ldb, ap( kc+1 ), 1, one, b( k, 1 ), ldb )
425 CALL dgemv( 'Transpose', n-k, nrhs, -one, b( k+1, 1 ),
426 $ ldb, ap( kc-( n-k ) ), 1, one, b( k-1, 1 ),
427 $ ldb )
428 END IF
429*
430* Interchange rows K and -IPIV(K).
431*
432 kp = -ipiv( k )
433 IF( kp.NE.k )
434 $ CALL dswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
435 kc = kc - ( n-k+2 )
436 k = k - 2
437 END IF
438*
439 GO TO 90
440 100 CONTINUE
441 END IF
442*
443 RETURN
444*
445* End of DSPTRS
446*
447 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
DGEMV
Definition dgemv.f:158
subroutine dger(m, n, alpha, x, incx, y, incy, a, lda)
DGER
Definition dger.f:130
subroutine dsptrs(uplo, n, nrhs, ap, ipiv, b, ldb, info)
DSPTRS
Definition dsptrs.f:115
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82