Annotation of rpl/lapack/lapack/dsytrs2.f, revision 1.13
1.4 bertrand 1: *> \brief \b DSYTRS2
2: *
3: * =========== DOCUMENTATION ===========
4: *
1.12 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.4 bertrand 7: *
8: *> \htmlonly
1.12 bertrand 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">
1.4 bertrand 15: *> [TXT]</a>
1.12 bertrand 16: *> \endhtmlonly
1.4 bertrand 17: *
18: * Definition:
19: * ===========
20: *
1.12 bertrand 21: * SUBROUTINE DSYTRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB,
1.4 bertrand 22: * WORK, INFO )
1.12 bertrand 23: *
1.4 bertrand 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: * ..
1.12 bertrand 32: *
1.4 bertrand 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: *>
1.9 bertrand 69: *> \param[in,out] A
1.4 bertrand 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.
1.9 bertrand 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.
1.4 bertrand 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
1.10 bertrand 109: *> WORK is DOUBLE PRECISION array, dimension (N)
1.4 bertrand 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: *
1.12 bertrand 122: *> \author Univ. of Tennessee
123: *> \author Univ. of California Berkeley
124: *> \author Univ. of Colorado Denver
125: *> \author NAG Ltd.
1.4 bertrand 126: *
1.10 bertrand 127: *> \date June 2016
1.4 bertrand 128: *
129: *> \ingroup doubleSYcomputational
130: *
131: * =====================================================================
1.12 bertrand 132: SUBROUTINE DSYTRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB,
1.1 bertrand 133: $ WORK, INFO )
134: *
1.12 bertrand 135: * -- LAPACK computational routine (version 3.7.0) --
1.1 bertrand 136: * -- LAPACK is a software package provided by Univ. of Tennessee, --
137: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.10 bertrand 138: * June 2016
1.1 bertrand 139: *
140: * .. Scalar Arguments ..
141: CHARACTER UPLO
142: INTEGER INFO, LDA, LDB, N, NRHS
143: * ..
144: * .. Array Arguments ..
145: INTEGER IPIV( * )
146: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * )
147: * ..
148: *
149: * =====================================================================
150: *
151: * .. Parameters ..
152: DOUBLE PRECISION ONE
153: PARAMETER ( ONE = 1.0D+0 )
154: * ..
155: * .. Local Scalars ..
156: LOGICAL UPPER
157: INTEGER I, IINFO, J, K, KP
158: DOUBLE PRECISION AK, AKM1, AKM1K, BK, BKM1, DENOM
159: * ..
160: * .. External Functions ..
161: LOGICAL LSAME
162: EXTERNAL LSAME
163: * ..
164: * .. External Subroutines ..
165: EXTERNAL DSCAL, DSYCONV, DSWAP, DTRSM, XERBLA
166: * ..
167: * .. Intrinsic Functions ..
168: INTRINSIC MAX
169: * ..
170: * .. Executable Statements ..
171: *
172: INFO = 0
173: UPPER = LSAME( UPLO, 'U' )
174: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
175: INFO = -1
176: ELSE IF( N.LT.0 ) THEN
177: INFO = -2
178: ELSE IF( NRHS.LT.0 ) THEN
179: INFO = -3
180: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
181: INFO = -5
182: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
183: INFO = -8
184: END IF
185: IF( INFO.NE.0 ) THEN
186: CALL XERBLA( 'DSYTRS2', -INFO )
187: RETURN
188: END IF
189: *
190: * Quick return if possible
191: *
192: IF( N.EQ.0 .OR. NRHS.EQ.0 )
193: $ RETURN
194: *
195: * Convert A
196: *
197: CALL DSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO )
198: *
199: IF( UPPER ) THEN
200: *
1.3 bertrand 201: * Solve A*X = B, where A = U*D*U**T.
1.1 bertrand 202: *
1.12 bertrand 203: * P**T * B
1.1 bertrand 204: K=N
205: DO WHILE ( K .GE. 1 )
206: IF( IPIV( K ).GT.0 ) THEN
207: * 1 x 1 diagonal block
208: * Interchange rows K and IPIV(K).
209: KP = IPIV( K )
210: IF( KP.NE.K )
211: $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
212: K=K-1
213: ELSE
214: * 2 x 2 diagonal block
215: * Interchange rows K-1 and -IPIV(K).
216: KP = -IPIV( K )
217: IF( KP.EQ.-IPIV( K-1 ) )
218: $ CALL DSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
219: K=K-2
220: END IF
221: END DO
222: *
1.3 bertrand 223: * Compute (U \P**T * B) -> B [ (U \P**T * B) ]
1.1 bertrand 224: *
1.3 bertrand 225: CALL DTRSM('L','U','N','U',N,NRHS,ONE,A,LDA,B,LDB)
1.1 bertrand 226: *
1.3 bertrand 227: * Compute D \ B -> B [ D \ (U \P**T * B) ]
1.12 bertrand 228: *
1.1 bertrand 229: I=N
230: DO WHILE ( I .GE. 1 )
231: IF( IPIV(I) .GT. 0 ) THEN
1.3 bertrand 232: CALL DSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), LDB )
1.1 bertrand 233: ELSEIF ( I .GT. 1) THEN
234: IF ( IPIV(I-1) .EQ. IPIV(I) ) THEN
235: AKM1K = WORK(I)
236: AKM1 = A( I-1, I-1 ) / AKM1K
237: AK = A( I, I ) / AKM1K
238: DENOM = AKM1*AK - ONE
239: DO 15 J = 1, NRHS
240: BKM1 = B( I-1, J ) / AKM1K
241: BK = B( I, J ) / AKM1K
242: B( I-1, J ) = ( AK*BKM1-BK ) / DENOM
243: B( I, J ) = ( AKM1*BK-BKM1 ) / DENOM
244: 15 CONTINUE
245: I = I - 1
246: ENDIF
247: ENDIF
248: I = I - 1
249: END DO
250: *
1.3 bertrand 251: * Compute (U**T \ B) -> B [ U**T \ (D \ (U \P**T * B) ) ]
1.1 bertrand 252: *
1.3 bertrand 253: CALL DTRSM('L','U','T','U',N,NRHS,ONE,A,LDA,B,LDB)
1.1 bertrand 254: *
1.3 bertrand 255: * P * B [ P * (U**T \ (D \ (U \P**T * B) )) ]
1.1 bertrand 256: *
257: K=1
258: DO WHILE ( K .LE. N )
259: IF( IPIV( K ).GT.0 ) THEN
260: * 1 x 1 diagonal block
261: * Interchange rows K and IPIV(K).
262: KP = IPIV( K )
263: IF( KP.NE.K )
264: $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
265: K=K+1
266: ELSE
267: * 2 x 2 diagonal block
268: * Interchange rows K-1 and -IPIV(K).
269: KP = -IPIV( K )
270: IF( K .LT. N .AND. KP.EQ.-IPIV( K+1 ) )
271: $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
272: K=K+2
273: ENDIF
274: END DO
275: *
276: ELSE
277: *
1.3 bertrand 278: * Solve A*X = B, where A = L*D*L**T.
1.1 bertrand 279: *
1.12 bertrand 280: * P**T * B
1.1 bertrand 281: K=1
282: DO WHILE ( K .LE. N )
283: IF( IPIV( K ).GT.0 ) THEN
284: * 1 x 1 diagonal block
285: * Interchange rows K and IPIV(K).
286: KP = IPIV( K )
287: IF( KP.NE.K )
288: $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
289: K=K+1
290: ELSE
291: * 2 x 2 diagonal block
292: * Interchange rows K and -IPIV(K+1).
293: KP = -IPIV( K+1 )
294: IF( KP.EQ.-IPIV( K ) )
295: $ CALL DSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
296: K=K+2
297: ENDIF
298: END DO
299: *
1.3 bertrand 300: * Compute (L \P**T * B) -> B [ (L \P**T * B) ]
1.1 bertrand 301: *
1.3 bertrand 302: CALL DTRSM('L','L','N','U',N,NRHS,ONE,A,LDA,B,LDB)
1.1 bertrand 303: *
1.3 bertrand 304: * Compute D \ B -> B [ D \ (L \P**T * B) ]
1.12 bertrand 305: *
1.1 bertrand 306: I=1
307: DO WHILE ( I .LE. N )
308: IF( IPIV(I) .GT. 0 ) THEN
1.3 bertrand 309: CALL DSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), LDB )
1.1 bertrand 310: ELSE
311: AKM1K = WORK(I)
312: AKM1 = A( I, I ) / AKM1K
313: AK = A( I+1, I+1 ) / AKM1K
314: DENOM = AKM1*AK - ONE
315: DO 25 J = 1, NRHS
316: BKM1 = B( I, J ) / AKM1K
317: BK = B( I+1, J ) / AKM1K
318: B( I, J ) = ( AK*BKM1-BK ) / DENOM
319: B( I+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
320: 25 CONTINUE
321: I = I + 1
322: ENDIF
323: I = I + 1
324: END DO
325: *
1.3 bertrand 326: * Compute (L**T \ B) -> B [ L**T \ (D \ (L \P**T * B) ) ]
1.12 bertrand 327: *
1.3 bertrand 328: CALL DTRSM('L','L','T','U',N,NRHS,ONE,A,LDA,B,LDB)
1.1 bertrand 329: *
1.3 bertrand 330: * P * B [ P * (L**T \ (D \ (L \P**T * B) )) ]
1.1 bertrand 331: *
332: K=N
333: DO WHILE ( K .GE. 1 )
334: IF( IPIV( K ).GT.0 ) THEN
335: * 1 x 1 diagonal block
336: * Interchange rows K and IPIV(K).
337: KP = IPIV( K )
338: IF( KP.NE.K )
339: $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
340: K=K-1
341: ELSE
342: * 2 x 2 diagonal block
343: * Interchange rows K-1 and -IPIV(K).
344: KP = -IPIV( K )
345: IF( K.GT.1 .AND. KP.EQ.-IPIV( K-1 ) )
346: $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
347: K=K-2
348: ENDIF
349: END DO
350: *
351: END IF
352: *
353: * Revert A
354: *
355: CALL DSYCONV( UPLO, 'R', N, A, LDA, IPIV, WORK, IINFO )
356: *
357: RETURN
358: *
359: * End of DSYTRS2
360: *
361: 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