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 REAL 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*> \date November 2015
128*
129*> \ingroup doubleSYcomputational
130*
131*  =====================================================================
132      SUBROUTINE DSYTRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB,
133     $                    WORK, INFO )
134*
135*  -- LAPACK computational routine (version 3.6.0) --
136*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
137*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
138*     November 2015
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*
201*        Solve A*X = B, where A = U*D*U**T.
202*
203*       P**T * B
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*
223*  Compute (U \P**T * B) -> B    [ (U \P**T * B) ]
224*
225        CALL DTRSM('L','U','N','U',N,NRHS,ONE,A,LDA,B,LDB)
226*
227*  Compute D \ B -> B   [ D \ (U \P**T * B) ]
228*
229         I=N
230         DO WHILE ( I .GE. 1 )
231            IF( IPIV(I) .GT. 0 ) THEN
232              CALL DSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), LDB )
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*
251*      Compute (U**T \ B) -> B   [ U**T \ (D \ (U \P**T * B) ) ]
252*
253         CALL DTRSM('L','U','T','U',N,NRHS,ONE,A,LDA,B,LDB)
254*
255*       P * B  [ P * (U**T \ (D \ (U \P**T * B) )) ]
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*
278*        Solve A*X = B, where A = L*D*L**T.
279*
280*       P**T * B
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*
300*  Compute (L \P**T * B) -> B    [ (L \P**T * B) ]
301*
302        CALL DTRSM('L','L','N','U',N,NRHS,ONE,A,LDA,B,LDB)
303*
304*  Compute D \ B -> B   [ D \ (L \P**T * B) ]
305*
306         I=1
307         DO WHILE ( I .LE. N )
308            IF( IPIV(I) .GT. 0 ) THEN
309              CALL DSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), LDB )
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*
326*  Compute (L**T \ B) -> B   [ L**T \ (D \ (L \P**T * B) ) ]
327*
328        CALL DTRSM('L','L','T','U',N,NRHS,ONE,A,LDA,B,LDB)
329*
330*       P * B  [ P * (L**T \ (D \ (L \P**T * B) )) ]
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
362