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