1 SUBROUTINE ZPFTRI( TRANSR, UPLO, N, A, INFO )
2 *
3 * -- LAPACK routine (version 3.3.1) --
4 *
5 * -- Contributed by Fred Gustavson of the IBM Watson Research Center --
6 * -- April 2011 --
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, N
14 * .. Array Arguments ..
15 COMPLEX*16 A( 0: * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * ZPFTRI computes the inverse of a complex Hermitian positive definite
22 * matrix A using the Cholesky factorization A = U**H*U or A = L*L**H
23 * computed by ZPFTRF.
24 *
25 * Arguments
26 * =========
27 *
28 * TRANSR (input) CHARACTER*1
29 * = 'N': The Normal TRANSR of RFP A is stored;
30 * = 'C': The Conjugate-transpose TRANSR of RFP A is stored.
31 *
32 * UPLO (input) CHARACTER*1
33 * = 'U': Upper triangle of A is stored;
34 * = 'L': Lower triangle of A is stored.
35 *
36 * N (input) INTEGER
37 * The order of the matrix A. N >= 0.
38 *
39 * A (input/output) COMPLEX*16 array, dimension ( N*(N+1)/2 );
40 * On entry, the Hermitian matrix A in RFP format. RFP format is
41 * described by TRANSR, UPLO, and N as follows: If TRANSR = 'N'
42 * then RFP A is (0:N,0:k-1) when N is even; k=N/2. RFP A is
43 * (0:N-1,0:k) when N is odd; k=N/2. IF TRANSR = 'C' then RFP is
44 * the Conjugate-transpose of RFP A as defined when
45 * TRANSR = 'N'. The contents of RFP A are defined by UPLO as
46 * follows: If UPLO = 'U' the RFP A contains the nt elements of
47 * upper packed A. If UPLO = 'L' the RFP A contains the elements
48 * of lower packed A. The LDA of RFP A is (N+1)/2 when TRANSR =
49 * 'C'. When TRANSR is 'N' the LDA is N+1 when N is even and N
50 * is odd. See the Note below for more details.
51 *
52 * On exit, the Hermitian inverse of the original matrix, in the
53 * same storage format.
54 *
55 * INFO (output) INTEGER
56 * = 0: successful exit
57 * < 0: if INFO = -i, the i-th argument had an illegal value
58 * > 0: if INFO = i, the (i,i) element of the factor U or L is
59 * zero, and the inverse could not be computed.
60 *
61 * Further Details
62 * ===============
63 *
64 * We first consider Standard Packed Format when N is even.
65 * We give an example where N = 6.
66 *
67 * AP is Upper AP is Lower
68 *
69 * 00 01 02 03 04 05 00
70 * 11 12 13 14 15 10 11
71 * 22 23 24 25 20 21 22
72 * 33 34 35 30 31 32 33
73 * 44 45 40 41 42 43 44
74 * 55 50 51 52 53 54 55
75 *
76 *
77 * Let TRANSR = 'N'. RFP holds AP as follows:
78 * For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
79 * three columns of AP upper. The lower triangle A(4:6,0:2) consists of
80 * conjugate-transpose of the first three columns of AP upper.
81 * For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
82 * three columns of AP lower. The upper triangle A(0:2,0:2) consists of
83 * conjugate-transpose of the last three columns of AP lower.
84 * To denote conjugate we place -- above the element. This covers the
85 * case N even and TRANSR = 'N'.
86 *
87 * RFP A RFP A
88 *
89 * -- -- --
90 * 03 04 05 33 43 53
91 * -- --
92 * 13 14 15 00 44 54
93 * --
94 * 23 24 25 10 11 55
95 *
96 * 33 34 35 20 21 22
97 * --
98 * 00 44 45 30 31 32
99 * -- --
100 * 01 11 55 40 41 42
101 * -- -- --
102 * 02 12 22 50 51 52
103 *
104 * Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
105 * transpose of RFP A above. One therefore gets:
106 *
107 *
108 * RFP A RFP A
109 *
110 * -- -- -- -- -- -- -- -- -- --
111 * 03 13 23 33 00 01 02 33 00 10 20 30 40 50
112 * -- -- -- -- -- -- -- -- -- --
113 * 04 14 24 34 44 11 12 43 44 11 21 31 41 51
114 * -- -- -- -- -- -- -- -- -- --
115 * 05 15 25 35 45 55 22 53 54 55 22 32 42 52
116 *
117 *
118 * We next consider Standard Packed Format when N is odd.
119 * We give an example where N = 5.
120 *
121 * AP is Upper AP is Lower
122 *
123 * 00 01 02 03 04 00
124 * 11 12 13 14 10 11
125 * 22 23 24 20 21 22
126 * 33 34 30 31 32 33
127 * 44 40 41 42 43 44
128 *
129 *
130 * Let TRANSR = 'N'. RFP holds AP as follows:
131 * For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
132 * three columns of AP upper. The lower triangle A(3:4,0:1) consists of
133 * conjugate-transpose of the first two columns of AP upper.
134 * For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
135 * three columns of AP lower. The upper triangle A(0:1,1:2) consists of
136 * conjugate-transpose of the last two columns of AP lower.
137 * To denote conjugate we place -- above the element. This covers the
138 * case N odd and TRANSR = 'N'.
139 *
140 * RFP A RFP A
141 *
142 * -- --
143 * 02 03 04 00 33 43
144 * --
145 * 12 13 14 10 11 44
146 *
147 * 22 23 24 20 21 22
148 * --
149 * 00 33 34 30 31 32
150 * -- --
151 * 01 11 44 40 41 42
152 *
153 * Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
154 * transpose of RFP A above. One therefore gets:
155 *
156 *
157 * RFP A RFP A
158 *
159 * -- -- -- -- -- -- -- -- --
160 * 02 12 22 00 01 00 10 20 30 40 50
161 * -- -- -- -- -- -- -- -- --
162 * 03 13 23 33 11 33 11 21 31 41 51
163 * -- -- -- -- -- -- -- -- --
164 * 04 14 24 34 44 43 44 22 32 42 52
165 *
166 * =====================================================================
167 *
168 * .. Parameters ..
169 DOUBLE PRECISION ONE
170 COMPLEX*16 CONE
171 PARAMETER ( ONE = 1.D0, CONE = ( 1.D0, 0.D0 ) )
172 * ..
173 * .. Local Scalars ..
174 LOGICAL LOWER, NISODD, NORMALTRANSR
175 INTEGER N1, N2, K
176 * ..
177 * .. External Functions ..
178 LOGICAL LSAME
179 EXTERNAL LSAME
180 * ..
181 * .. External Subroutines ..
182 EXTERNAL XERBLA, ZTFTRI, ZLAUUM, ZTRMM, ZHERK
183 * ..
184 * .. Intrinsic Functions ..
185 INTRINSIC MOD
186 * ..
187 * .. Executable Statements ..
188 *
189 * Test the input parameters.
190 *
191 INFO = 0
192 NORMALTRANSR = LSAME( TRANSR, 'N' )
193 LOWER = LSAME( UPLO, 'L' )
194 IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'C' ) ) THEN
195 INFO = -1
196 ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
197 INFO = -2
198 ELSE IF( N.LT.0 ) THEN
199 INFO = -3
200 END IF
201 IF( INFO.NE.0 ) THEN
202 CALL XERBLA( 'ZPFTRI', -INFO )
203 RETURN
204 END IF
205 *
206 * Quick return if possible
207 *
208 IF( N.EQ.0 )
209 $ RETURN
210 *
211 * Invert the triangular Cholesky factor U or L.
212 *
213 CALL ZTFTRI( TRANSR, UPLO, 'N', N, A, INFO )
214 IF( INFO.GT.0 )
215 $ RETURN
216 *
217 * If N is odd, set NISODD = .TRUE.
218 * If N is even, set K = N/2 and NISODD = .FALSE.
219 *
220 IF( MOD( N, 2 ).EQ.0 ) THEN
221 K = N / 2
222 NISODD = .FALSE.
223 ELSE
224 NISODD = .TRUE.
225 END IF
226 *
227 * Set N1 and N2 depending on LOWER
228 *
229 IF( LOWER ) THEN
230 N2 = N / 2
231 N1 = N - N2
232 ELSE
233 N1 = N / 2
234 N2 = N - N1
235 END IF
236 *
237 * Start execution of triangular matrix multiply: inv(U)*inv(U)^C or
238 * inv(L)^C*inv(L). There are eight cases.
239 *
240 IF( NISODD ) THEN
241 *
242 * N is odd
243 *
244 IF( NORMALTRANSR ) THEN
245 *
246 * N is odd and TRANSR = 'N'
247 *
248 IF( LOWER ) THEN
249 *
250 * SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:N1-1) )
251 * T1 -> a(0,0), T2 -> a(0,1), S -> a(N1,0)
252 * T1 -> a(0), T2 -> a(n), S -> a(N1)
253 *
254 CALL ZLAUUM( 'L', N1, A( 0 ), N, INFO )
255 CALL ZHERK( 'L', 'C', N1, N2, ONE, A( N1 ), N, ONE,
256 $ A( 0 ), N )
257 CALL ZTRMM( 'L', 'U', 'N', 'N', N2, N1, CONE, A( N ), N,
258 $ A( N1 ), N )
259 CALL ZLAUUM( 'U', N2, A( N ), N, INFO )
260 *
261 ELSE
262 *
263 * SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:N2-1)
264 * T1 -> a(N1+1,0), T2 -> a(N1,0), S -> a(0,0)
265 * T1 -> a(N2), T2 -> a(N1), S -> a(0)
266 *
267 CALL ZLAUUM( 'L', N1, A( N2 ), N, INFO )
268 CALL ZHERK( 'L', 'N', N1, N2, ONE, A( 0 ), N, ONE,
269 $ A( N2 ), N )
270 CALL ZTRMM( 'R', 'U', 'C', 'N', N1, N2, CONE, A( N1 ), N,
271 $ A( 0 ), N )
272 CALL ZLAUUM( 'U', N2, A( N1 ), N, INFO )
273 *
274 END IF
275 *
276 ELSE
277 *
278 * N is odd and TRANSR = 'C'
279 *
280 IF( LOWER ) THEN
281 *
282 * SRPA for LOWER, TRANSPOSE, and N is odd
283 * T1 -> a(0), T2 -> a(1), S -> a(0+N1*N1)
284 *
285 CALL ZLAUUM( 'U', N1, A( 0 ), N1, INFO )
286 CALL ZHERK( 'U', 'N', N1, N2, ONE, A( N1*N1 ), N1, ONE,
287 $ A( 0 ), N1 )
288 CALL ZTRMM( 'R', 'L', 'N', 'N', N1, N2, CONE, A( 1 ), N1,
289 $ A( N1*N1 ), N1 )
290 CALL ZLAUUM( 'L', N2, A( 1 ), N1, INFO )
291 *
292 ELSE
293 *
294 * SRPA for UPPER, TRANSPOSE, and N is odd
295 * T1 -> a(0+N2*N2), T2 -> a(0+N1*N2), S -> a(0)
296 *
297 CALL ZLAUUM( 'U', N1, A( N2*N2 ), N2, INFO )
298 CALL ZHERK( 'U', 'C', N1, N2, ONE, A( 0 ), N2, ONE,
299 $ A( N2*N2 ), N2 )
300 CALL ZTRMM( 'L', 'L', 'C', 'N', N2, N1, CONE, A( N1*N2 ),
301 $ N2, A( 0 ), N2 )
302 CALL ZLAUUM( 'L', N2, A( N1*N2 ), N2, INFO )
303 *
304 END IF
305 *
306 END IF
307 *
308 ELSE
309 *
310 * N is even
311 *
312 IF( NORMALTRANSR ) THEN
313 *
314 * N is even and TRANSR = 'N'
315 *
316 IF( LOWER ) THEN
317 *
318 * SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) )
319 * T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0)
320 * T1 -> a(1), T2 -> a(0), S -> a(k+1)
321 *
322 CALL ZLAUUM( 'L', K, A( 1 ), N+1, INFO )
323 CALL ZHERK( 'L', 'C', K, K, ONE, A( K+1 ), N+1, ONE,
324 $ A( 1 ), N+1 )
325 CALL ZTRMM( 'L', 'U', 'N', 'N', K, K, CONE, A( 0 ), N+1,
326 $ A( K+1 ), N+1 )
327 CALL ZLAUUM( 'U', K, A( 0 ), N+1, INFO )
328 *
329 ELSE
330 *
331 * SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) )
332 * T1 -> a(k+1,0) , T2 -> a(k,0), S -> a(0,0)
333 * T1 -> a(k+1), T2 -> a(k), S -> a(0)
334 *
335 CALL ZLAUUM( 'L', K, A( K+1 ), N+1, INFO )
336 CALL ZHERK( 'L', 'N', K, K, ONE, A( 0 ), N+1, ONE,
337 $ A( K+1 ), N+1 )
338 CALL ZTRMM( 'R', 'U', 'C', 'N', K, K, CONE, A( K ), N+1,
339 $ A( 0 ), N+1 )
340 CALL ZLAUUM( 'U', K, A( K ), N+1, INFO )
341 *
342 END IF
343 *
344 ELSE
345 *
346 * N is even and TRANSR = 'C'
347 *
348 IF( LOWER ) THEN
349 *
350 * SRPA for LOWER, TRANSPOSE, and N is even (see paper)
351 * T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1),
352 * T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k
353 *
354 CALL ZLAUUM( 'U', K, A( K ), K, INFO )
355 CALL ZHERK( 'U', 'N', K, K, ONE, A( K*( K+1 ) ), K, ONE,
356 $ A( K ), K )
357 CALL ZTRMM( 'R', 'L', 'N', 'N', K, K, CONE, A( 0 ), K,
358 $ A( K*( K+1 ) ), K )
359 CALL ZLAUUM( 'L', K, A( 0 ), K, INFO )
360 *
361 ELSE
362 *
363 * SRPA for UPPER, TRANSPOSE, and N is even (see paper)
364 * T1 -> B(0,k+1), T2 -> B(0,k), S -> B(0,0),
365 * T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k
366 *
367 CALL ZLAUUM( 'U', K, A( K*( K+1 ) ), K, INFO )
368 CALL ZHERK( 'U', 'C', K, K, ONE, A( 0 ), K, ONE,
369 $ A( K*( K+1 ) ), K )
370 CALL ZTRMM( 'L', 'L', 'C', 'N', K, K, CONE, A( K*K ), K,
371 $ A( 0 ), K )
372 CALL ZLAUUM( 'L', K, A( K*K ), K, INFO )
373 *
374 END IF
375 *
376 END IF
377 *
378 END IF
379 *
380 RETURN
381 *
382 * End of ZPFTRI
383 *
384 END
2 *
3 * -- LAPACK routine (version 3.3.1) --
4 *
5 * -- Contributed by Fred Gustavson of the IBM Watson Research Center --
6 * -- April 2011 --
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, N
14 * .. Array Arguments ..
15 COMPLEX*16 A( 0: * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * ZPFTRI computes the inverse of a complex Hermitian positive definite
22 * matrix A using the Cholesky factorization A = U**H*U or A = L*L**H
23 * computed by ZPFTRF.
24 *
25 * Arguments
26 * =========
27 *
28 * TRANSR (input) CHARACTER*1
29 * = 'N': The Normal TRANSR of RFP A is stored;
30 * = 'C': The Conjugate-transpose TRANSR of RFP A is stored.
31 *
32 * UPLO (input) CHARACTER*1
33 * = 'U': Upper triangle of A is stored;
34 * = 'L': Lower triangle of A is stored.
35 *
36 * N (input) INTEGER
37 * The order of the matrix A. N >= 0.
38 *
39 * A (input/output) COMPLEX*16 array, dimension ( N*(N+1)/2 );
40 * On entry, the Hermitian matrix A in RFP format. RFP format is
41 * described by TRANSR, UPLO, and N as follows: If TRANSR = 'N'
42 * then RFP A is (0:N,0:k-1) when N is even; k=N/2. RFP A is
43 * (0:N-1,0:k) when N is odd; k=N/2. IF TRANSR = 'C' then RFP is
44 * the Conjugate-transpose of RFP A as defined when
45 * TRANSR = 'N'. The contents of RFP A are defined by UPLO as
46 * follows: If UPLO = 'U' the RFP A contains the nt elements of
47 * upper packed A. If UPLO = 'L' the RFP A contains the elements
48 * of lower packed A. The LDA of RFP A is (N+1)/2 when TRANSR =
49 * 'C'. When TRANSR is 'N' the LDA is N+1 when N is even and N
50 * is odd. See the Note below for more details.
51 *
52 * On exit, the Hermitian inverse of the original matrix, in the
53 * same storage format.
54 *
55 * INFO (output) INTEGER
56 * = 0: successful exit
57 * < 0: if INFO = -i, the i-th argument had an illegal value
58 * > 0: if INFO = i, the (i,i) element of the factor U or L is
59 * zero, and the inverse could not be computed.
60 *
61 * Further Details
62 * ===============
63 *
64 * We first consider Standard Packed Format when N is even.
65 * We give an example where N = 6.
66 *
67 * AP is Upper AP is Lower
68 *
69 * 00 01 02 03 04 05 00
70 * 11 12 13 14 15 10 11
71 * 22 23 24 25 20 21 22
72 * 33 34 35 30 31 32 33
73 * 44 45 40 41 42 43 44
74 * 55 50 51 52 53 54 55
75 *
76 *
77 * Let TRANSR = 'N'. RFP holds AP as follows:
78 * For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
79 * three columns of AP upper. The lower triangle A(4:6,0:2) consists of
80 * conjugate-transpose of the first three columns of AP upper.
81 * For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
82 * three columns of AP lower. The upper triangle A(0:2,0:2) consists of
83 * conjugate-transpose of the last three columns of AP lower.
84 * To denote conjugate we place -- above the element. This covers the
85 * case N even and TRANSR = 'N'.
86 *
87 * RFP A RFP A
88 *
89 * -- -- --
90 * 03 04 05 33 43 53
91 * -- --
92 * 13 14 15 00 44 54
93 * --
94 * 23 24 25 10 11 55
95 *
96 * 33 34 35 20 21 22
97 * --
98 * 00 44 45 30 31 32
99 * -- --
100 * 01 11 55 40 41 42
101 * -- -- --
102 * 02 12 22 50 51 52
103 *
104 * Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
105 * transpose of RFP A above. One therefore gets:
106 *
107 *
108 * RFP A RFP A
109 *
110 * -- -- -- -- -- -- -- -- -- --
111 * 03 13 23 33 00 01 02 33 00 10 20 30 40 50
112 * -- -- -- -- -- -- -- -- -- --
113 * 04 14 24 34 44 11 12 43 44 11 21 31 41 51
114 * -- -- -- -- -- -- -- -- -- --
115 * 05 15 25 35 45 55 22 53 54 55 22 32 42 52
116 *
117 *
118 * We next consider Standard Packed Format when N is odd.
119 * We give an example where N = 5.
120 *
121 * AP is Upper AP is Lower
122 *
123 * 00 01 02 03 04 00
124 * 11 12 13 14 10 11
125 * 22 23 24 20 21 22
126 * 33 34 30 31 32 33
127 * 44 40 41 42 43 44
128 *
129 *
130 * Let TRANSR = 'N'. RFP holds AP as follows:
131 * For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
132 * three columns of AP upper. The lower triangle A(3:4,0:1) consists of
133 * conjugate-transpose of the first two columns of AP upper.
134 * For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
135 * three columns of AP lower. The upper triangle A(0:1,1:2) consists of
136 * conjugate-transpose of the last two columns of AP lower.
137 * To denote conjugate we place -- above the element. This covers the
138 * case N odd and TRANSR = 'N'.
139 *
140 * RFP A RFP A
141 *
142 * -- --
143 * 02 03 04 00 33 43
144 * --
145 * 12 13 14 10 11 44
146 *
147 * 22 23 24 20 21 22
148 * --
149 * 00 33 34 30 31 32
150 * -- --
151 * 01 11 44 40 41 42
152 *
153 * Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
154 * transpose of RFP A above. One therefore gets:
155 *
156 *
157 * RFP A RFP A
158 *
159 * -- -- -- -- -- -- -- -- --
160 * 02 12 22 00 01 00 10 20 30 40 50
161 * -- -- -- -- -- -- -- -- --
162 * 03 13 23 33 11 33 11 21 31 41 51
163 * -- -- -- -- -- -- -- -- --
164 * 04 14 24 34 44 43 44 22 32 42 52
165 *
166 * =====================================================================
167 *
168 * .. Parameters ..
169 DOUBLE PRECISION ONE
170 COMPLEX*16 CONE
171 PARAMETER ( ONE = 1.D0, CONE = ( 1.D0, 0.D0 ) )
172 * ..
173 * .. Local Scalars ..
174 LOGICAL LOWER, NISODD, NORMALTRANSR
175 INTEGER N1, N2, K
176 * ..
177 * .. External Functions ..
178 LOGICAL LSAME
179 EXTERNAL LSAME
180 * ..
181 * .. External Subroutines ..
182 EXTERNAL XERBLA, ZTFTRI, ZLAUUM, ZTRMM, ZHERK
183 * ..
184 * .. Intrinsic Functions ..
185 INTRINSIC MOD
186 * ..
187 * .. Executable Statements ..
188 *
189 * Test the input parameters.
190 *
191 INFO = 0
192 NORMALTRANSR = LSAME( TRANSR, 'N' )
193 LOWER = LSAME( UPLO, 'L' )
194 IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'C' ) ) THEN
195 INFO = -1
196 ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
197 INFO = -2
198 ELSE IF( N.LT.0 ) THEN
199 INFO = -3
200 END IF
201 IF( INFO.NE.0 ) THEN
202 CALL XERBLA( 'ZPFTRI', -INFO )
203 RETURN
204 END IF
205 *
206 * Quick return if possible
207 *
208 IF( N.EQ.0 )
209 $ RETURN
210 *
211 * Invert the triangular Cholesky factor U or L.
212 *
213 CALL ZTFTRI( TRANSR, UPLO, 'N', N, A, INFO )
214 IF( INFO.GT.0 )
215 $ RETURN
216 *
217 * If N is odd, set NISODD = .TRUE.
218 * If N is even, set K = N/2 and NISODD = .FALSE.
219 *
220 IF( MOD( N, 2 ).EQ.0 ) THEN
221 K = N / 2
222 NISODD = .FALSE.
223 ELSE
224 NISODD = .TRUE.
225 END IF
226 *
227 * Set N1 and N2 depending on LOWER
228 *
229 IF( LOWER ) THEN
230 N2 = N / 2
231 N1 = N - N2
232 ELSE
233 N1 = N / 2
234 N2 = N - N1
235 END IF
236 *
237 * Start execution of triangular matrix multiply: inv(U)*inv(U)^C or
238 * inv(L)^C*inv(L). There are eight cases.
239 *
240 IF( NISODD ) THEN
241 *
242 * N is odd
243 *
244 IF( NORMALTRANSR ) THEN
245 *
246 * N is odd and TRANSR = 'N'
247 *
248 IF( LOWER ) THEN
249 *
250 * SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:N1-1) )
251 * T1 -> a(0,0), T2 -> a(0,1), S -> a(N1,0)
252 * T1 -> a(0), T2 -> a(n), S -> a(N1)
253 *
254 CALL ZLAUUM( 'L', N1, A( 0 ), N, INFO )
255 CALL ZHERK( 'L', 'C', N1, N2, ONE, A( N1 ), N, ONE,
256 $ A( 0 ), N )
257 CALL ZTRMM( 'L', 'U', 'N', 'N', N2, N1, CONE, A( N ), N,
258 $ A( N1 ), N )
259 CALL ZLAUUM( 'U', N2, A( N ), N, INFO )
260 *
261 ELSE
262 *
263 * SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:N2-1)
264 * T1 -> a(N1+1,0), T2 -> a(N1,0), S -> a(0,0)
265 * T1 -> a(N2), T2 -> a(N1), S -> a(0)
266 *
267 CALL ZLAUUM( 'L', N1, A( N2 ), N, INFO )
268 CALL ZHERK( 'L', 'N', N1, N2, ONE, A( 0 ), N, ONE,
269 $ A( N2 ), N )
270 CALL ZTRMM( 'R', 'U', 'C', 'N', N1, N2, CONE, A( N1 ), N,
271 $ A( 0 ), N )
272 CALL ZLAUUM( 'U', N2, A( N1 ), N, INFO )
273 *
274 END IF
275 *
276 ELSE
277 *
278 * N is odd and TRANSR = 'C'
279 *
280 IF( LOWER ) THEN
281 *
282 * SRPA for LOWER, TRANSPOSE, and N is odd
283 * T1 -> a(0), T2 -> a(1), S -> a(0+N1*N1)
284 *
285 CALL ZLAUUM( 'U', N1, A( 0 ), N1, INFO )
286 CALL ZHERK( 'U', 'N', N1, N2, ONE, A( N1*N1 ), N1, ONE,
287 $ A( 0 ), N1 )
288 CALL ZTRMM( 'R', 'L', 'N', 'N', N1, N2, CONE, A( 1 ), N1,
289 $ A( N1*N1 ), N1 )
290 CALL ZLAUUM( 'L', N2, A( 1 ), N1, INFO )
291 *
292 ELSE
293 *
294 * SRPA for UPPER, TRANSPOSE, and N is odd
295 * T1 -> a(0+N2*N2), T2 -> a(0+N1*N2), S -> a(0)
296 *
297 CALL ZLAUUM( 'U', N1, A( N2*N2 ), N2, INFO )
298 CALL ZHERK( 'U', 'C', N1, N2, ONE, A( 0 ), N2, ONE,
299 $ A( N2*N2 ), N2 )
300 CALL ZTRMM( 'L', 'L', 'C', 'N', N2, N1, CONE, A( N1*N2 ),
301 $ N2, A( 0 ), N2 )
302 CALL ZLAUUM( 'L', N2, A( N1*N2 ), N2, INFO )
303 *
304 END IF
305 *
306 END IF
307 *
308 ELSE
309 *
310 * N is even
311 *
312 IF( NORMALTRANSR ) THEN
313 *
314 * N is even and TRANSR = 'N'
315 *
316 IF( LOWER ) THEN
317 *
318 * SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) )
319 * T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0)
320 * T1 -> a(1), T2 -> a(0), S -> a(k+1)
321 *
322 CALL ZLAUUM( 'L', K, A( 1 ), N+1, INFO )
323 CALL ZHERK( 'L', 'C', K, K, ONE, A( K+1 ), N+1, ONE,
324 $ A( 1 ), N+1 )
325 CALL ZTRMM( 'L', 'U', 'N', 'N', K, K, CONE, A( 0 ), N+1,
326 $ A( K+1 ), N+1 )
327 CALL ZLAUUM( 'U', K, A( 0 ), N+1, INFO )
328 *
329 ELSE
330 *
331 * SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) )
332 * T1 -> a(k+1,0) , T2 -> a(k,0), S -> a(0,0)
333 * T1 -> a(k+1), T2 -> a(k), S -> a(0)
334 *
335 CALL ZLAUUM( 'L', K, A( K+1 ), N+1, INFO )
336 CALL ZHERK( 'L', 'N', K, K, ONE, A( 0 ), N+1, ONE,
337 $ A( K+1 ), N+1 )
338 CALL ZTRMM( 'R', 'U', 'C', 'N', K, K, CONE, A( K ), N+1,
339 $ A( 0 ), N+1 )
340 CALL ZLAUUM( 'U', K, A( K ), N+1, INFO )
341 *
342 END IF
343 *
344 ELSE
345 *
346 * N is even and TRANSR = 'C'
347 *
348 IF( LOWER ) THEN
349 *
350 * SRPA for LOWER, TRANSPOSE, and N is even (see paper)
351 * T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1),
352 * T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k
353 *
354 CALL ZLAUUM( 'U', K, A( K ), K, INFO )
355 CALL ZHERK( 'U', 'N', K, K, ONE, A( K*( K+1 ) ), K, ONE,
356 $ A( K ), K )
357 CALL ZTRMM( 'R', 'L', 'N', 'N', K, K, CONE, A( 0 ), K,
358 $ A( K*( K+1 ) ), K )
359 CALL ZLAUUM( 'L', K, A( 0 ), K, INFO )
360 *
361 ELSE
362 *
363 * SRPA for UPPER, TRANSPOSE, and N is even (see paper)
364 * T1 -> B(0,k+1), T2 -> B(0,k), S -> B(0,0),
365 * T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k
366 *
367 CALL ZLAUUM( 'U', K, A( K*( K+1 ) ), K, INFO )
368 CALL ZHERK( 'U', 'C', K, K, ONE, A( 0 ), K, ONE,
369 $ A( K*( K+1 ) ), K )
370 CALL ZTRMM( 'L', 'L', 'C', 'N', K, K, CONE, A( K*K ), K,
371 $ A( 0 ), K )
372 CALL ZLAUUM( 'L', K, A( K*K ), K, INFO )
373 *
374 END IF
375 *
376 END IF
377 *
378 END IF
379 *
380 RETURN
381 *
382 * End of ZPFTRI
383 *
384 END