1: *> \brief \b DSYTRS2
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DSYTRS2 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsytrs2.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsytrs2.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsytrs2.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DSYTRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB,
22: * WORK, INFO )
23: *
24: * .. Scalar Arguments ..
25: * CHARACTER UPLO
26: * INTEGER INFO, LDA, LDB, N, NRHS
27: * ..
28: * .. Array Arguments ..
29: * INTEGER IPIV( * )
30: * DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * )
31: * ..
32: *
33: *
34: *> \par Purpose:
35: * =============
36: *>
37: *> \verbatim
38: *>
39: *> DSYTRS2 solves a system of linear equations A*X = B with a real
40: *> symmetric matrix A using the factorization A = U*D*U**T or
41: *> A = L*D*L**T computed by DSYTRF and converted by DSYCONV.
42: *> \endverbatim
43: *
44: * Arguments:
45: * ==========
46: *
47: *> \param[in] UPLO
48: *> \verbatim
49: *> UPLO is CHARACTER*1
50: *> Specifies whether the details of the factorization are stored
51: *> as an upper or lower triangular matrix.
52: *> = 'U': Upper triangular, form is A = U*D*U**T;
53: *> = 'L': Lower triangular, form is A = L*D*L**T.
54: *> \endverbatim
55: *>
56: *> \param[in] N
57: *> \verbatim
58: *> N is INTEGER
59: *> The order of the matrix A. N >= 0.
60: *> \endverbatim
61: *>
62: *> \param[in] NRHS
63: *> \verbatim
64: *> NRHS is INTEGER
65: *> The number of right hand sides, i.e., the number of columns
66: *> of the matrix B. NRHS >= 0.
67: *> \endverbatim
68: *>
69: *> \param[in,out] A
70: *> \verbatim
71: *> A is DOUBLE PRECISION array, dimension (LDA,N)
72: *> The block diagonal matrix D and the multipliers used to
73: *> obtain the factor U or L as computed by DSYTRF.
74: *> Note that A is input / output. This might be counter-intuitive,
75: *> and one may think that A is input only. A is input / output. This
76: *> is because, at the start of the subroutine, we permute A in a
77: *> "better" form and then we permute A back to its original form at
78: *> the end.
79: *> \endverbatim
80: *>
81: *> \param[in] LDA
82: *> \verbatim
83: *> LDA is INTEGER
84: *> The leading dimension of the array A. LDA >= max(1,N).
85: *> \endverbatim
86: *>
87: *> \param[in] IPIV
88: *> \verbatim
89: *> IPIV is INTEGER array, dimension (N)
90: *> Details of the interchanges and the block structure of D
91: *> as determined by DSYTRF.
92: *> \endverbatim
93: *>
94: *> \param[in,out] B
95: *> \verbatim
96: *> B is DOUBLE PRECISION array, dimension (LDB,NRHS)
97: *> On entry, the right hand side matrix B.
98: *> On exit, the solution matrix X.
99: *> \endverbatim
100: *>
101: *> \param[in] LDB
102: *> \verbatim
103: *> LDB is INTEGER
104: *> The leading dimension of the array B. LDB >= max(1,N).
105: *> \endverbatim
106: *>
107: *> \param[out] WORK
108: *> \verbatim
109: *> WORK is DOUBLE PRECISION array, dimension (N)
110: *> \endverbatim
111: *>
112: *> \param[out] INFO
113: *> \verbatim
114: *> INFO is INTEGER
115: *> = 0: successful exit
116: *> < 0: if INFO = -i, the i-th argument had an illegal value
117: *> \endverbatim
118: *
119: * Authors:
120: * ========
121: *
122: *> \author Univ. of Tennessee
123: *> \author Univ. of California Berkeley
124: *> \author Univ. of Colorado Denver
125: *> \author NAG Ltd.
126: *
127: *> \ingroup doubleSYcomputational
128: *
129: * =====================================================================
130: SUBROUTINE DSYTRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB,
131: $ WORK, INFO )
132: *
133: * -- LAPACK computational routine --
134: * -- LAPACK is a software package provided by Univ. of Tennessee, --
135: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
136: *
137: * .. Scalar Arguments ..
138: CHARACTER UPLO
139: INTEGER INFO, LDA, LDB, N, NRHS
140: * ..
141: * .. Array Arguments ..
142: INTEGER IPIV( * )
143: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * )
144: * ..
145: *
146: * =====================================================================
147: *
148: * .. Parameters ..
149: DOUBLE PRECISION ONE
150: PARAMETER ( ONE = 1.0D+0 )
151: * ..
152: * .. Local Scalars ..
153: LOGICAL UPPER
154: INTEGER I, IINFO, J, K, KP
155: DOUBLE PRECISION AK, AKM1, AKM1K, BK, BKM1, DENOM
156: * ..
157: * .. External Functions ..
158: LOGICAL LSAME
159: EXTERNAL LSAME
160: * ..
161: * .. External Subroutines ..
162: EXTERNAL DSCAL, DSYCONV, DSWAP, DTRSM, XERBLA
163: * ..
164: * .. Intrinsic Functions ..
165: INTRINSIC MAX
166: * ..
167: * .. Executable Statements ..
168: *
169: INFO = 0
170: UPPER = LSAME( UPLO, 'U' )
171: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
172: INFO = -1
173: ELSE IF( N.LT.0 ) THEN
174: INFO = -2
175: ELSE IF( NRHS.LT.0 ) THEN
176: INFO = -3
177: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
178: INFO = -5
179: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
180: INFO = -8
181: END IF
182: IF( INFO.NE.0 ) THEN
183: CALL XERBLA( 'DSYTRS2', -INFO )
184: RETURN
185: END IF
186: *
187: * Quick return if possible
188: *
189: IF( N.EQ.0 .OR. NRHS.EQ.0 )
190: $ RETURN
191: *
192: * Convert A
193: *
194: CALL DSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO )
195: *
196: IF( UPPER ) THEN
197: *
198: * Solve A*X = B, where A = U*D*U**T.
199: *
200: * P**T * B
201: K=N
202: DO WHILE ( K .GE. 1 )
203: IF( IPIV( K ).GT.0 ) THEN
204: * 1 x 1 diagonal block
205: * Interchange rows K and IPIV(K).
206: KP = IPIV( K )
207: IF( KP.NE.K )
208: $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
209: K=K-1
210: ELSE
211: * 2 x 2 diagonal block
212: * Interchange rows K-1 and -IPIV(K).
213: KP = -IPIV( K )
214: IF( KP.EQ.-IPIV( K-1 ) )
215: $ CALL DSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
216: K=K-2
217: END IF
218: END DO
219: *
220: * Compute (U \P**T * B) -> B [ (U \P**T * B) ]
221: *
222: CALL DTRSM('L','U','N','U',N,NRHS,ONE,A,LDA,B,LDB)
223: *
224: * Compute D \ B -> B [ D \ (U \P**T * B) ]
225: *
226: I=N
227: DO WHILE ( I .GE. 1 )
228: IF( IPIV(I) .GT. 0 ) THEN
229: CALL DSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), LDB )
230: ELSEIF ( I .GT. 1) THEN
231: IF ( IPIV(I-1) .EQ. IPIV(I) ) THEN
232: AKM1K = WORK(I)
233: AKM1 = A( I-1, I-1 ) / AKM1K
234: AK = A( I, I ) / AKM1K
235: DENOM = AKM1*AK - ONE
236: DO 15 J = 1, NRHS
237: BKM1 = B( I-1, J ) / AKM1K
238: BK = B( I, J ) / AKM1K
239: B( I-1, J ) = ( AK*BKM1-BK ) / DENOM
240: B( I, J ) = ( AKM1*BK-BKM1 ) / DENOM
241: 15 CONTINUE
242: I = I - 1
243: ENDIF
244: ENDIF
245: I = I - 1
246: END DO
247: *
248: * Compute (U**T \ B) -> B [ U**T \ (D \ (U \P**T * B) ) ]
249: *
250: CALL DTRSM('L','U','T','U',N,NRHS,ONE,A,LDA,B,LDB)
251: *
252: * P * B [ P * (U**T \ (D \ (U \P**T * B) )) ]
253: *
254: K=1
255: DO WHILE ( K .LE. N )
256: IF( IPIV( K ).GT.0 ) THEN
257: * 1 x 1 diagonal block
258: * Interchange rows K and IPIV(K).
259: KP = IPIV( K )
260: IF( KP.NE.K )
261: $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
262: K=K+1
263: ELSE
264: * 2 x 2 diagonal block
265: * Interchange rows K-1 and -IPIV(K).
266: KP = -IPIV( K )
267: IF( K .LT. N .AND. KP.EQ.-IPIV( K+1 ) )
268: $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
269: K=K+2
270: ENDIF
271: END DO
272: *
273: ELSE
274: *
275: * Solve A*X = B, where A = L*D*L**T.
276: *
277: * P**T * B
278: K=1
279: DO WHILE ( K .LE. N )
280: IF( IPIV( K ).GT.0 ) THEN
281: * 1 x 1 diagonal block
282: * Interchange rows K and IPIV(K).
283: KP = IPIV( K )
284: IF( KP.NE.K )
285: $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
286: K=K+1
287: ELSE
288: * 2 x 2 diagonal block
289: * Interchange rows K and -IPIV(K+1).
290: KP = -IPIV( K+1 )
291: IF( KP.EQ.-IPIV( K ) )
292: $ CALL DSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
293: K=K+2
294: ENDIF
295: END DO
296: *
297: * Compute (L \P**T * B) -> B [ (L \P**T * B) ]
298: *
299: CALL DTRSM('L','L','N','U',N,NRHS,ONE,A,LDA,B,LDB)
300: *
301: * Compute D \ B -> B [ D \ (L \P**T * B) ]
302: *
303: I=1
304: DO WHILE ( I .LE. N )
305: IF( IPIV(I) .GT. 0 ) THEN
306: CALL DSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), LDB )
307: ELSE
308: AKM1K = WORK(I)
309: AKM1 = A( I, I ) / AKM1K
310: AK = A( I+1, I+1 ) / AKM1K
311: DENOM = AKM1*AK - ONE
312: DO 25 J = 1, NRHS
313: BKM1 = B( I, J ) / AKM1K
314: BK = B( I+1, J ) / AKM1K
315: B( I, J ) = ( AK*BKM1-BK ) / DENOM
316: B( I+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
317: 25 CONTINUE
318: I = I + 1
319: ENDIF
320: I = I + 1
321: END DO
322: *
323: * Compute (L**T \ B) -> B [ L**T \ (D \ (L \P**T * B) ) ]
324: *
325: CALL DTRSM('L','L','T','U',N,NRHS,ONE,A,LDA,B,LDB)
326: *
327: * P * B [ P * (L**T \ (D \ (L \P**T * B) )) ]
328: *
329: K=N
330: DO WHILE ( K .GE. 1 )
331: IF( IPIV( K ).GT.0 ) THEN
332: * 1 x 1 diagonal block
333: * Interchange rows K and IPIV(K).
334: KP = IPIV( K )
335: IF( KP.NE.K )
336: $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
337: K=K-1
338: ELSE
339: * 2 x 2 diagonal block
340: * Interchange rows K-1 and -IPIV(K).
341: KP = -IPIV( K )
342: IF( K.GT.1 .AND. KP.EQ.-IPIV( K-1 ) )
343: $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
344: K=K-2
345: ENDIF
346: END DO
347: *
348: END IF
349: *
350: * Revert A
351: *
352: CALL DSYCONV( UPLO, 'R', N, A, LDA, IPIV, WORK, IINFO )
353: *
354: RETURN
355: *
356: * End of DSYTRS2
357: *
358: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dsytrs2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack