1: SUBROUTINE ZPFTRS( TRANSR, UPLO, N, NRHS, A, B, LDB, INFO )
2: *
3: * -- LAPACK routine (version 3.3.0) --
4: *
5: * -- Contributed by Fred Gustavson of the IBM Watson Research Center --
6: * November 2010
7: *
8: * -- LAPACK is a software package provided by Univ. of Tennessee, --
9: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
10: *
11: * .. Scalar Arguments ..
12: CHARACTER TRANSR, UPLO
13: INTEGER INFO, LDB, N, NRHS
14: * ..
15: * .. Array Arguments ..
16: COMPLEX*16 A( 0: * ), B( LDB, * )
17: * ..
18: *
19: * Purpose
20: * =======
21: *
22: * ZPFTRS solves a system of linear equations A*X = B with a Hermitian
23: * positive definite matrix A using the Cholesky factorization
24: * A = U**H*U or A = L*L**H computed by ZPFTRF.
25: *
26: * Arguments
27: * =========
28: *
29: * TRANSR (input) CHARACTER*1
30: * = 'N': The Normal TRANSR of RFP A is stored;
31: * = 'C': The Conjugate-transpose TRANSR of RFP A is stored.
32: *
33: * UPLO (input) CHARACTER*1
34: * = 'U': Upper triangle of RFP A is stored;
35: * = 'L': Lower triangle of RFP A is stored.
36: *
37: * N (input) INTEGER
38: * The order of the matrix A. N >= 0.
39: *
40: * NRHS (input) INTEGER
41: * The number of right hand sides, i.e., the number of columns
42: * of the matrix B. NRHS >= 0.
43: *
44: * A (input) COMPLEX*16 array, dimension ( N*(N+1)/2 );
45: * The triangular factor U or L from the Cholesky factorization
46: * of RFP A = U**H*U or RFP A = L*L**H, as computed by ZPFTRF.
47: * See note below for more details about RFP A.
48: *
49: * B (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
50: * On entry, the right hand side matrix B.
51: * On exit, the solution matrix X.
52: *
53: * LDB (input) INTEGER
54: * The leading dimension of the array B. LDB >= max(1,N).
55: *
56: * INFO (output) INTEGER
57: * = 0: successful exit
58: * < 0: if INFO = -i, the i-th argument had an illegal value
59: *
60: * Further Details
61: * ===============
62: *
63: * We first consider Standard Packed Format when N is even.
64: * We give an example where N = 6.
65: *
66: * AP is Upper AP is Lower
67: *
68: * 00 01 02 03 04 05 00
69: * 11 12 13 14 15 10 11
70: * 22 23 24 25 20 21 22
71: * 33 34 35 30 31 32 33
72: * 44 45 40 41 42 43 44
73: * 55 50 51 52 53 54 55
74: *
75: *
76: * Let TRANSR = 'N'. RFP holds AP as follows:
77: * For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
78: * three columns of AP upper. The lower triangle A(4:6,0:2) consists of
79: * conjugate-transpose of the first three columns of AP upper.
80: * For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
81: * three columns of AP lower. The upper triangle A(0:2,0:2) consists of
82: * conjugate-transpose of the last three columns of AP lower.
83: * To denote conjugate we place -- above the element. This covers the
84: * case N even and TRANSR = 'N'.
85: *
86: * RFP A RFP A
87: *
88: * -- -- --
89: * 03 04 05 33 43 53
90: * -- --
91: * 13 14 15 00 44 54
92: * --
93: * 23 24 25 10 11 55
94: *
95: * 33 34 35 20 21 22
96: * --
97: * 00 44 45 30 31 32
98: * -- --
99: * 01 11 55 40 41 42
100: * -- -- --
101: * 02 12 22 50 51 52
102: *
103: * Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
104: * transpose of RFP A above. One therefore gets:
105: *
106: *
107: * RFP A RFP A
108: *
109: * -- -- -- -- -- -- -- -- -- --
110: * 03 13 23 33 00 01 02 33 00 10 20 30 40 50
111: * -- -- -- -- -- -- -- -- -- --
112: * 04 14 24 34 44 11 12 43 44 11 21 31 41 51
113: * -- -- -- -- -- -- -- -- -- --
114: * 05 15 25 35 45 55 22 53 54 55 22 32 42 52
115: *
116: *
117: * We next consider Standard Packed Format when N is odd.
118: * We give an example where N = 5.
119: *
120: * AP is Upper AP is Lower
121: *
122: * 00 01 02 03 04 00
123: * 11 12 13 14 10 11
124: * 22 23 24 20 21 22
125: * 33 34 30 31 32 33
126: * 44 40 41 42 43 44
127: *
128: *
129: * Let TRANSR = 'N'. RFP holds AP as follows:
130: * For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
131: * three columns of AP upper. The lower triangle A(3:4,0:1) consists of
132: * conjugate-transpose of the first two columns of AP upper.
133: * For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
134: * three columns of AP lower. The upper triangle A(0:1,1:2) consists of
135: * conjugate-transpose of the last two columns of AP lower.
136: * To denote conjugate we place -- above the element. This covers the
137: * case N odd and TRANSR = 'N'.
138: *
139: * RFP A RFP A
140: *
141: * -- --
142: * 02 03 04 00 33 43
143: * --
144: * 12 13 14 10 11 44
145: *
146: * 22 23 24 20 21 22
147: * --
148: * 00 33 34 30 31 32
149: * -- --
150: * 01 11 44 40 41 42
151: *
152: * Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
153: * transpose of RFP A above. One therefore gets:
154: *
155: *
156: * RFP A RFP A
157: *
158: * -- -- -- -- -- -- -- -- --
159: * 02 12 22 00 01 00 10 20 30 40 50
160: * -- -- -- -- -- -- -- -- --
161: * 03 13 23 33 11 33 11 21 31 41 51
162: * -- -- -- -- -- -- -- -- --
163: * 04 14 24 34 44 43 44 22 32 42 52
164: *
165: * =====================================================================
166: *
167: * .. Parameters ..
168: COMPLEX*16 CONE
169: PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) )
170: * ..
171: * .. Local Scalars ..
172: LOGICAL LOWER, NORMALTRANSR
173: * ..
174: * .. External Functions ..
175: LOGICAL LSAME
176: EXTERNAL LSAME
177: * ..
178: * .. External Subroutines ..
179: EXTERNAL XERBLA, ZTFSM
180: * ..
181: * .. Intrinsic Functions ..
182: INTRINSIC MAX
183: * ..
184: * .. Executable Statements ..
185: *
186: * Test the input parameters.
187: *
188: INFO = 0
189: NORMALTRANSR = LSAME( TRANSR, 'N' )
190: LOWER = LSAME( UPLO, 'L' )
191: IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'C' ) ) THEN
192: INFO = -1
193: ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
194: INFO = -2
195: ELSE IF( N.LT.0 ) THEN
196: INFO = -3
197: ELSE IF( NRHS.LT.0 ) THEN
198: INFO = -4
199: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
200: INFO = -7
201: END IF
202: IF( INFO.NE.0 ) THEN
203: CALL XERBLA( 'ZPFTRS', -INFO )
204: RETURN
205: END IF
206: *
207: * Quick return if possible
208: *
209: IF( N.EQ.0 .OR. NRHS.EQ.0 )
210: + RETURN
211: *
212: * start execution: there are two triangular solves
213: *
214: IF( LOWER ) THEN
215: CALL ZTFSM( TRANSR, 'L', UPLO, 'N', 'N', N, NRHS, CONE, A, B,
216: + LDB )
217: CALL ZTFSM( TRANSR, 'L', UPLO, 'C', 'N', N, NRHS, CONE, A, B,
218: + LDB )
219: ELSE
220: CALL ZTFSM( TRANSR, 'L', UPLO, 'C', 'N', N, NRHS, CONE, A, B,
221: + LDB )
222: CALL ZTFSM( TRANSR, 'L', UPLO, 'N', 'N', N, NRHS, CONE, A, B,
223: + LDB )
224: END IF
225: *
226: RETURN
227: *
228: * End of ZPFTRS
229: *
230: END
CVSweb interface <joel.bertrand@systella.fr>