Annotation of rpl/lapack/lapack/zsytrs_rook.f, revision 1.6
1.1 bertrand 1: *> \brief \b ZSYTRS_ROOK
2: *
3: * =========== DOCUMENTATION ===========
4: *
1.5 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.1 bertrand 7: *
8: *> \htmlonly
1.5 bertrand 9: *> Download ZSYTRS_ROOK + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zsytrs_rook.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zsytrs_rook.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zsytrs_rook.f">
1.1 bertrand 15: *> [TXT]</a>
1.5 bertrand 16: *> \endhtmlonly
1.1 bertrand 17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZSYTRS_ROOK( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO )
1.5 bertrand 22: *
1.1 bertrand 23: * .. Scalar Arguments ..
24: * CHARACTER UPLO
25: * INTEGER INFO, LDA, LDB, N, NRHS
26: * ..
27: * .. Array Arguments ..
28: * INTEGER IPIV( * )
29: * COMPLEX*16 A( LDA, * ), B( LDB, * )
30: * ..
1.5 bertrand 31: *
1.1 bertrand 32: *
33: *> \par Purpose:
34: * =============
35: *>
36: *> \verbatim
37: *>
38: *> ZSYTRS_ROOK solves a system of linear equations A*X = B with
39: *> a complex symmetric matrix A using the factorization A = U*D*U**T or
40: *> A = L*D*L**T computed by ZSYTRF_ROOK.
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] A
69: *> \verbatim
70: *> A is COMPLEX*16 array, dimension (LDA,N)
71: *> The block diagonal matrix D and the multipliers used to
72: *> obtain the factor U or L as computed by ZSYTRF_ROOK.
73: *> \endverbatim
74: *>
75: *> \param[in] LDA
76: *> \verbatim
77: *> LDA is INTEGER
78: *> The leading dimension of the array A. LDA >= max(1,N).
79: *> \endverbatim
80: *>
81: *> \param[in] IPIV
82: *> \verbatim
83: *> IPIV is INTEGER array, dimension (N)
84: *> Details of the interchanges and the block structure of D
85: *> as determined by ZSYTRF_ROOK.
86: *> \endverbatim
87: *>
88: *> \param[in,out] B
89: *> \verbatim
90: *> B is COMPLEX*16 array, dimension (LDB,NRHS)
91: *> On entry, the right hand side matrix B.
92: *> On exit, the solution matrix X.
93: *> \endverbatim
94: *>
95: *> \param[in] LDB
96: *> \verbatim
97: *> LDB is INTEGER
98: *> The leading dimension of the array B. LDB >= max(1,N).
99: *> \endverbatim
100: *>
101: *> \param[out] INFO
102: *> \verbatim
103: *> INFO is INTEGER
104: *> = 0: successful exit
105: *> < 0: if INFO = -i, the i-th argument had an illegal value
106: *> \endverbatim
107: *
108: * Authors:
109: * ========
110: *
1.5 bertrand 111: *> \author Univ. of Tennessee
112: *> \author Univ. of California Berkeley
113: *> \author Univ. of Colorado Denver
114: *> \author NAG Ltd.
1.1 bertrand 115: *
1.5 bertrand 116: *> \date December 2016
1.1 bertrand 117: *
118: *> \ingroup complex16SYcomputational
119: *
120: *> \par Contributors:
121: * ==================
122: *>
123: *> \verbatim
124: *>
1.5 bertrand 125: *> December 2016, Igor Kozachenko,
1.1 bertrand 126: *> Computer Science Division,
127: *> University of California, Berkeley
128: *>
129: *> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
130: *> School of Mathematics,
131: *> University of Manchester
132: *>
133: *> \endverbatim
134: *
135: * =====================================================================
136: SUBROUTINE ZSYTRS_ROOK( UPLO, N, NRHS, A, LDA, IPIV, B, LDB,
137: $ INFO )
138: *
1.5 bertrand 139: * -- LAPACK computational routine (version 3.7.0) --
1.1 bertrand 140: * -- LAPACK is a software package provided by Univ. of Tennessee, --
141: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.5 bertrand 142: * December 2016
1.1 bertrand 143: *
144: * .. Scalar Arguments ..
145: CHARACTER UPLO
146: INTEGER INFO, LDA, LDB, N, NRHS
147: * ..
148: * .. Array Arguments ..
149: INTEGER IPIV( * )
150: COMPLEX*16 A( LDA, * ), B( LDB, * )
151: * ..
152: *
153: * =====================================================================
154: *
155: * .. Parameters ..
156: COMPLEX*16 CONE
157: PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) )
158: * ..
159: * .. Local Scalars ..
160: LOGICAL UPPER
161: INTEGER J, K, KP
162: COMPLEX*16 AK, AKM1, AKM1K, BK, BKM1, DENOM
163: * ..
164: * .. External Functions ..
165: LOGICAL LSAME
166: EXTERNAL LSAME
167: * ..
168: * .. External Subroutines ..
169: EXTERNAL ZGEMV, ZGERU, ZSCAL, ZSWAP, XERBLA
170: * ..
171: * .. Intrinsic Functions ..
172: INTRINSIC MAX
173: * ..
174: * .. Executable Statements ..
175: *
176: INFO = 0
177: UPPER = LSAME( UPLO, 'U' )
178: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
179: INFO = -1
180: ELSE IF( N.LT.0 ) THEN
181: INFO = -2
182: ELSE IF( NRHS.LT.0 ) THEN
183: INFO = -3
184: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
185: INFO = -5
186: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
187: INFO = -8
188: END IF
189: IF( INFO.NE.0 ) THEN
190: CALL XERBLA( 'ZSYTRS_ROOK', -INFO )
191: RETURN
192: END IF
193: *
194: * Quick return if possible
195: *
196: IF( N.EQ.0 .OR. NRHS.EQ.0 )
197: $ RETURN
198: *
199: IF( UPPER ) THEN
200: *
201: * Solve A*X = B, where A = U*D*U**T.
202: *
203: * First solve U*D*X = B, overwriting B with X.
204: *
205: * K is the main loop index, decreasing from N to 1 in steps of
206: * 1 or 2, depending on the size of the diagonal blocks.
207: *
208: K = N
209: 10 CONTINUE
210: *
211: * If K < 1, exit from loop.
212: *
213: IF( K.LT.1 )
214: $ GO TO 30
215: *
216: IF( IPIV( K ).GT.0 ) THEN
217: *
218: * 1 x 1 diagonal block
219: *
220: * Interchange rows K and IPIV(K).
221: *
222: KP = IPIV( K )
223: IF( KP.NE.K )
224: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
225: *
226: * Multiply by inv(U(K)), where U(K) is the transformation
227: * stored in column K of A.
228: *
229: CALL ZGERU( K-1, NRHS, -CONE, A( 1, K ), 1, B( K, 1 ), LDB,
230: $ B( 1, 1 ), LDB )
231: *
232: * Multiply by the inverse of the diagonal block.
233: *
234: CALL ZSCAL( NRHS, CONE / A( K, K ), B( K, 1 ), LDB )
235: K = K - 1
236: ELSE
237: *
238: * 2 x 2 diagonal block
239: *
240: * Interchange rows K and -IPIV(K) THEN K-1 and -IPIV(K-1)
241: *
242: KP = -IPIV( K )
243: IF( KP.NE.K )
244: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
245: *
246: KP = -IPIV( K-1 )
247: IF( KP.NE.K-1 )
248: $ CALL ZSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
249: *
250: * Multiply by inv(U(K)), where U(K) is the transformation
251: * stored in columns K-1 and K of A.
252: *
253: IF( K.GT.2 ) THEN
254: CALL ZGERU( K-2, NRHS,-CONE, A( 1, K ), 1, B( K, 1 ),
255: $ LDB, B( 1, 1 ), LDB )
256: CALL ZGERU( K-2, NRHS,-CONE, A( 1, K-1 ), 1, B( K-1, 1 ),
257: $ LDB, B( 1, 1 ), LDB )
258: END IF
259: *
260: * Multiply by the inverse of the diagonal block.
261: *
262: AKM1K = A( K-1, K )
263: AKM1 = A( K-1, K-1 ) / AKM1K
264: AK = A( K, K ) / AKM1K
265: DENOM = AKM1*AK - CONE
266: DO 20 J = 1, NRHS
267: BKM1 = B( K-1, J ) / AKM1K
268: BK = B( K, J ) / AKM1K
269: B( K-1, J ) = ( AK*BKM1-BK ) / DENOM
270: B( K, J ) = ( AKM1*BK-BKM1 ) / DENOM
271: 20 CONTINUE
272: K = K - 2
273: END IF
274: *
275: GO TO 10
276: 30 CONTINUE
277: *
278: * Next solve U**T *X = B, overwriting B with X.
279: *
280: * K is the main loop index, increasing from 1 to N in steps of
281: * 1 or 2, depending on the size of the diagonal blocks.
282: *
283: K = 1
284: 40 CONTINUE
285: *
286: * If K > N, exit from loop.
287: *
288: IF( K.GT.N )
289: $ GO TO 50
290: *
291: IF( IPIV( K ).GT.0 ) THEN
292: *
293: * 1 x 1 diagonal block
294: *
295: * Multiply by inv(U**T(K)), where U(K) is the transformation
296: * stored in column K of A.
297: *
298: IF( K.GT.1 )
299: $ CALL ZGEMV( 'Transpose', K-1, NRHS, -CONE, B,
300: $ LDB, A( 1, K ), 1, CONE, B( K, 1 ), LDB )
301: *
302: * Interchange rows K and IPIV(K).
303: *
304: KP = IPIV( K )
305: IF( KP.NE.K )
306: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
307: K = K + 1
308: ELSE
309: *
310: * 2 x 2 diagonal block
311: *
312: * Multiply by inv(U**T(K+1)), where U(K+1) is the transformation
313: * stored in columns K and K+1 of A.
314: *
315: IF( K.GT.1 ) THEN
316: CALL ZGEMV( 'Transpose', K-1, NRHS, -CONE, B,
317: $ LDB, A( 1, K ), 1, CONE, B( K, 1 ), LDB )
318: CALL ZGEMV( 'Transpose', K-1, NRHS, -CONE, B,
319: $ LDB, A( 1, K+1 ), 1, CONE, B( K+1, 1 ), LDB )
320: END IF
321: *
322: * Interchange rows K and -IPIV(K) THEN K+1 and -IPIV(K+1).
323: *
324: KP = -IPIV( K )
325: IF( KP.NE.K )
326: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
327: *
328: KP = -IPIV( K+1 )
329: IF( KP.NE.K+1 )
330: $ CALL ZSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
331: *
332: K = K + 2
333: END IF
334: *
335: GO TO 40
336: 50 CONTINUE
337: *
338: ELSE
339: *
340: * Solve A*X = B, where A = L*D*L**T.
341: *
342: * First solve L*D*X = B, overwriting B with X.
343: *
344: * K is the main loop index, increasing from 1 to N in steps of
345: * 1 or 2, depending on the size of the diagonal blocks.
346: *
347: K = 1
348: 60 CONTINUE
349: *
350: * If K > N, exit from loop.
351: *
352: IF( K.GT.N )
353: $ GO TO 80
354: *
355: IF( IPIV( K ).GT.0 ) THEN
356: *
357: * 1 x 1 diagonal block
358: *
359: * Interchange rows K and IPIV(K).
360: *
361: KP = IPIV( K )
362: IF( KP.NE.K )
363: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
364: *
365: * Multiply by inv(L(K)), where L(K) is the transformation
366: * stored in column K of A.
367: *
368: IF( K.LT.N )
369: $ CALL ZGERU( N-K, NRHS, -CONE, A( K+1, K ), 1, B( K, 1 ),
370: $ LDB, B( K+1, 1 ), LDB )
371: *
372: * Multiply by the inverse of the diagonal block.
373: *
374: CALL ZSCAL( NRHS, CONE / A( K, K ), B( K, 1 ), LDB )
375: K = K + 1
376: ELSE
377: *
378: * 2 x 2 diagonal block
379: *
380: * Interchange rows K and -IPIV(K) THEN K+1 and -IPIV(K+1)
381: *
382: KP = -IPIV( K )
383: IF( KP.NE.K )
384: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
385: *
386: KP = -IPIV( K+1 )
387: IF( KP.NE.K+1 )
388: $ CALL ZSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
389: *
390: * Multiply by inv(L(K)), where L(K) is the transformation
391: * stored in columns K and K+1 of A.
392: *
393: IF( K.LT.N-1 ) THEN
394: CALL ZGERU( N-K-1, NRHS,-CONE, A( K+2, K ), 1, B( K, 1 ),
395: $ LDB, B( K+2, 1 ), LDB )
396: CALL ZGERU( N-K-1, NRHS,-CONE, A( K+2, K+1 ), 1,
397: $ B( K+1, 1 ), LDB, B( K+2, 1 ), LDB )
398: END IF
399: *
400: * Multiply by the inverse of the diagonal block.
401: *
402: AKM1K = A( K+1, K )
403: AKM1 = A( K, K ) / AKM1K
404: AK = A( K+1, K+1 ) / AKM1K
405: DENOM = AKM1*AK - CONE
406: DO 70 J = 1, NRHS
407: BKM1 = B( K, J ) / AKM1K
408: BK = B( K+1, J ) / AKM1K
409: B( K, J ) = ( AK*BKM1-BK ) / DENOM
410: B( K+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
411: 70 CONTINUE
412: K = K + 2
413: END IF
414: *
415: GO TO 60
416: 80 CONTINUE
417: *
418: * Next solve L**T *X = B, overwriting B with X.
419: *
420: * K is the main loop index, decreasing from N to 1 in steps of
421: * 1 or 2, depending on the size of the diagonal blocks.
422: *
423: K = N
424: 90 CONTINUE
425: *
426: * If K < 1, exit from loop.
427: *
428: IF( K.LT.1 )
429: $ GO TO 100
430: *
431: IF( IPIV( K ).GT.0 ) THEN
432: *
433: * 1 x 1 diagonal block
434: *
435: * Multiply by inv(L**T(K)), where L(K) is the transformation
436: * stored in column K of A.
437: *
438: IF( K.LT.N )
439: $ CALL ZGEMV( 'Transpose', N-K, NRHS, -CONE, B( K+1, 1 ),
440: $ LDB, A( K+1, K ), 1, CONE, B( K, 1 ), LDB )
441: *
442: * Interchange rows K and IPIV(K).
443: *
444: KP = IPIV( K )
445: IF( KP.NE.K )
446: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
447: K = K - 1
448: ELSE
449: *
450: * 2 x 2 diagonal block
451: *
452: * Multiply by inv(L**T(K-1)), where L(K-1) is the transformation
453: * stored in columns K-1 and K of A.
454: *
455: IF( K.LT.N ) THEN
456: CALL ZGEMV( 'Transpose', N-K, NRHS, -CONE, B( K+1, 1 ),
457: $ LDB, A( K+1, K ), 1, CONE, B( K, 1 ), LDB )
458: CALL ZGEMV( 'Transpose', N-K, NRHS, -CONE, B( K+1, 1 ),
459: $ LDB, A( K+1, K-1 ), 1, CONE, B( K-1, 1 ),
460: $ LDB )
461: END IF
462: *
463: * Interchange rows K and -IPIV(K) THEN K-1 and -IPIV(K-1)
464: *
465: KP = -IPIV( K )
466: IF( KP.NE.K )
467: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
468: *
469: KP = -IPIV( K-1 )
470: IF( KP.NE.K-1 )
471: $ CALL ZSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
472: *
473: K = K - 2
474: END IF
475: *
476: GO TO 90
477: 100 CONTINUE
478: END IF
479: *
480: RETURN
481: *
482: * End of ZSYTRS_ROOK
483: *
484: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zsytrs_rook.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack