Annotation of rpl/lapack/lapack/dsytrs2.f, revision 1.7
1.4 bertrand 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] 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: *> \endverbatim
75: *>
76: *> \param[in] LDA
77: *> \verbatim
78: *> LDA is INTEGER
79: *> The leading dimension of the array A. LDA >= max(1,N).
80: *> \endverbatim
81: *>
82: *> \param[in] IPIV
83: *> \verbatim
84: *> IPIV is INTEGER array, dimension (N)
85: *> Details of the interchanges and the block structure of D
86: *> as determined by DSYTRF.
87: *> \endverbatim
88: *>
89: *> \param[in,out] B
90: *> \verbatim
91: *> B is DOUBLE PRECISION array, dimension (LDB,NRHS)
92: *> On entry, the right hand side matrix B.
93: *> On exit, the solution matrix X.
94: *> \endverbatim
95: *>
96: *> \param[in] LDB
97: *> \verbatim
98: *> LDB is INTEGER
99: *> The leading dimension of the array B. LDB >= max(1,N).
100: *> \endverbatim
101: *>
102: *> \param[out] WORK
103: *> \verbatim
104: *> WORK is REAL array, dimension (N)
105: *> \endverbatim
106: *>
107: *> \param[out] INFO
108: *> \verbatim
109: *> INFO is INTEGER
110: *> = 0: successful exit
111: *> < 0: if INFO = -i, the i-th argument had an illegal value
112: *> \endverbatim
113: *
114: * Authors:
115: * ========
116: *
117: *> \author Univ. of Tennessee
118: *> \author Univ. of California Berkeley
119: *> \author Univ. of Colorado Denver
120: *> \author NAG Ltd.
121: *
122: *> \date November 2011
123: *
124: *> \ingroup doubleSYcomputational
125: *
126: * =====================================================================
1.1 bertrand 127: SUBROUTINE DSYTRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB,
128: $ WORK, INFO )
129: *
1.4 bertrand 130: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 131: * -- LAPACK is a software package provided by Univ. of Tennessee, --
132: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.4 bertrand 133: * November 2011
1.1 bertrand 134: *
135: * .. Scalar Arguments ..
136: CHARACTER UPLO
137: INTEGER INFO, LDA, LDB, N, NRHS
138: * ..
139: * .. Array Arguments ..
140: INTEGER IPIV( * )
141: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * )
142: * ..
143: *
144: * =====================================================================
145: *
146: * .. Parameters ..
147: DOUBLE PRECISION ONE
148: PARAMETER ( ONE = 1.0D+0 )
149: * ..
150: * .. Local Scalars ..
151: LOGICAL UPPER
152: INTEGER I, IINFO, J, K, KP
153: DOUBLE PRECISION AK, AKM1, AKM1K, BK, BKM1, DENOM
154: * ..
155: * .. External Functions ..
156: LOGICAL LSAME
157: EXTERNAL LSAME
158: * ..
159: * .. External Subroutines ..
160: EXTERNAL DSCAL, DSYCONV, DSWAP, DTRSM, XERBLA
161: * ..
162: * .. Intrinsic Functions ..
163: INTRINSIC MAX
164: * ..
165: * .. Executable Statements ..
166: *
167: INFO = 0
168: UPPER = LSAME( UPLO, 'U' )
169: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
170: INFO = -1
171: ELSE IF( N.LT.0 ) THEN
172: INFO = -2
173: ELSE IF( NRHS.LT.0 ) THEN
174: INFO = -3
175: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
176: INFO = -5
177: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
178: INFO = -8
179: END IF
180: IF( INFO.NE.0 ) THEN
181: CALL XERBLA( 'DSYTRS2', -INFO )
182: RETURN
183: END IF
184: *
185: * Quick return if possible
186: *
187: IF( N.EQ.0 .OR. NRHS.EQ.0 )
188: $ RETURN
189: *
190: * Convert A
191: *
192: CALL DSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO )
193: *
194: IF( UPPER ) THEN
195: *
1.3 bertrand 196: * Solve A*X = B, where A = U*D*U**T.
1.1 bertrand 197: *
1.3 bertrand 198: * P**T * B
1.1 bertrand 199: K=N
200: DO WHILE ( K .GE. 1 )
201: IF( IPIV( K ).GT.0 ) THEN
202: * 1 x 1 diagonal block
203: * Interchange rows K and IPIV(K).
204: KP = IPIV( K )
205: IF( KP.NE.K )
206: $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
207: K=K-1
208: ELSE
209: * 2 x 2 diagonal block
210: * Interchange rows K-1 and -IPIV(K).
211: KP = -IPIV( K )
212: IF( KP.EQ.-IPIV( K-1 ) )
213: $ CALL DSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
214: K=K-2
215: END IF
216: END DO
217: *
1.3 bertrand 218: * Compute (U \P**T * B) -> B [ (U \P**T * B) ]
1.1 bertrand 219: *
1.3 bertrand 220: CALL DTRSM('L','U','N','U',N,NRHS,ONE,A,LDA,B,LDB)
1.1 bertrand 221: *
1.3 bertrand 222: * Compute D \ B -> B [ D \ (U \P**T * B) ]
1.1 bertrand 223: *
224: I=N
225: DO WHILE ( I .GE. 1 )
226: IF( IPIV(I) .GT. 0 ) THEN
1.3 bertrand 227: CALL DSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), LDB )
1.1 bertrand 228: ELSEIF ( I .GT. 1) THEN
229: IF ( IPIV(I-1) .EQ. IPIV(I) ) THEN
230: AKM1K = WORK(I)
231: AKM1 = A( I-1, I-1 ) / AKM1K
232: AK = A( I, I ) / AKM1K
233: DENOM = AKM1*AK - ONE
234: DO 15 J = 1, NRHS
235: BKM1 = B( I-1, J ) / AKM1K
236: BK = B( I, J ) / AKM1K
237: B( I-1, J ) = ( AK*BKM1-BK ) / DENOM
238: B( I, J ) = ( AKM1*BK-BKM1 ) / DENOM
239: 15 CONTINUE
240: I = I - 1
241: ENDIF
242: ENDIF
243: I = I - 1
244: END DO
245: *
1.3 bertrand 246: * Compute (U**T \ B) -> B [ U**T \ (D \ (U \P**T * B) ) ]
1.1 bertrand 247: *
1.3 bertrand 248: CALL DTRSM('L','U','T','U',N,NRHS,ONE,A,LDA,B,LDB)
1.1 bertrand 249: *
1.3 bertrand 250: * P * B [ P * (U**T \ (D \ (U \P**T * B) )) ]
1.1 bertrand 251: *
252: K=1
253: DO WHILE ( K .LE. N )
254: IF( IPIV( K ).GT.0 ) THEN
255: * 1 x 1 diagonal block
256: * Interchange rows K and IPIV(K).
257: KP = IPIV( K )
258: IF( KP.NE.K )
259: $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
260: K=K+1
261: ELSE
262: * 2 x 2 diagonal block
263: * Interchange rows K-1 and -IPIV(K).
264: KP = -IPIV( K )
265: IF( K .LT. N .AND. KP.EQ.-IPIV( K+1 ) )
266: $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
267: K=K+2
268: ENDIF
269: END DO
270: *
271: ELSE
272: *
1.3 bertrand 273: * Solve A*X = B, where A = L*D*L**T.
1.1 bertrand 274: *
1.3 bertrand 275: * P**T * B
1.1 bertrand 276: K=1
277: DO WHILE ( K .LE. N )
278: IF( IPIV( K ).GT.0 ) THEN
279: * 1 x 1 diagonal block
280: * Interchange rows K and IPIV(K).
281: KP = IPIV( K )
282: IF( KP.NE.K )
283: $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
284: K=K+1
285: ELSE
286: * 2 x 2 diagonal block
287: * Interchange rows K and -IPIV(K+1).
288: KP = -IPIV( K+1 )
289: IF( KP.EQ.-IPIV( K ) )
290: $ CALL DSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
291: K=K+2
292: ENDIF
293: END DO
294: *
1.3 bertrand 295: * Compute (L \P**T * B) -> B [ (L \P**T * B) ]
1.1 bertrand 296: *
1.3 bertrand 297: CALL DTRSM('L','L','N','U',N,NRHS,ONE,A,LDA,B,LDB)
1.1 bertrand 298: *
1.3 bertrand 299: * Compute D \ B -> B [ D \ (L \P**T * B) ]
1.1 bertrand 300: *
301: I=1
302: DO WHILE ( I .LE. N )
303: IF( IPIV(I) .GT. 0 ) THEN
1.3 bertrand 304: CALL DSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), LDB )
1.1 bertrand 305: ELSE
306: AKM1K = WORK(I)
307: AKM1 = A( I, I ) / AKM1K
308: AK = A( I+1, I+1 ) / AKM1K
309: DENOM = AKM1*AK - ONE
310: DO 25 J = 1, NRHS
311: BKM1 = B( I, J ) / AKM1K
312: BK = B( I+1, J ) / AKM1K
313: B( I, J ) = ( AK*BKM1-BK ) / DENOM
314: B( I+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
315: 25 CONTINUE
316: I = I + 1
317: ENDIF
318: I = I + 1
319: END DO
320: *
1.3 bertrand 321: * Compute (L**T \ B) -> B [ L**T \ (D \ (L \P**T * B) ) ]
1.1 bertrand 322: *
1.3 bertrand 323: CALL DTRSM('L','L','T','U',N,NRHS,ONE,A,LDA,B,LDB)
1.1 bertrand 324: *
1.3 bertrand 325: * P * B [ P * (L**T \ (D \ (L \P**T * B) )) ]
1.1 bertrand 326: *
327: K=N
328: DO WHILE ( K .GE. 1 )
329: IF( IPIV( K ).GT.0 ) THEN
330: * 1 x 1 diagonal block
331: * Interchange rows K and IPIV(K).
332: KP = IPIV( K )
333: IF( KP.NE.K )
334: $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
335: K=K-1
336: ELSE
337: * 2 x 2 diagonal block
338: * Interchange rows K-1 and -IPIV(K).
339: KP = -IPIV( K )
340: IF( K.GT.1 .AND. KP.EQ.-IPIV( K-1 ) )
341: $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
342: K=K-2
343: ENDIF
344: END DO
345: *
346: END IF
347: *
348: * Revert A
349: *
350: CALL DSYCONV( UPLO, 'R', N, A, LDA, IPIV, WORK, IINFO )
351: *
352: RETURN
353: *
354: * End of DSYTRS2
355: *
356: 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