1 SUBROUTINE DPFTRI( 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 DOUBLE PRECISION A( 0: * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DPFTRI computes the inverse of a (real) symmetric positive definite
22 * matrix A using the Cholesky factorization A = U**T*U or A = L*L**T
23 * computed by DPFTRF.
24 *
25 * Arguments
26 * =========
27 *
28 * TRANSR (input) CHARACTER*1
29 * = 'N': The Normal TRANSR of RFP A is stored;
30 * = 'T': The 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) DOUBLE PRECISION array, dimension ( N*(N+1)/2 )
40 * On entry, the symmetric 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 = 'T' then RFP is
44 * the 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 * 'T'. 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 symmetric 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 Rectangular Full Packed (RFP) Format when N is
65 * even. 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 * the 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 * the transpose of the last three columns of AP lower.
84 * This covers the case N even and TRANSR = 'N'.
85 *
86 * RFP A RFP A
87 *
88 * 03 04 05 33 43 53
89 * 13 14 15 00 44 54
90 * 23 24 25 10 11 55
91 * 33 34 35 20 21 22
92 * 00 44 45 30 31 32
93 * 01 11 55 40 41 42
94 * 02 12 22 50 51 52
95 *
96 * Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
97 * transpose of RFP A above. One therefore gets:
98 *
99 *
100 * RFP A RFP A
101 *
102 * 03 13 23 33 00 01 02 33 00 10 20 30 40 50
103 * 04 14 24 34 44 11 12 43 44 11 21 31 41 51
104 * 05 15 25 35 45 55 22 53 54 55 22 32 42 52
105 *
106 *
107 * We then consider Rectangular Full Packed (RFP) Format when N is
108 * odd. We give an example where N = 5.
109 *
110 * AP is Upper AP is Lower
111 *
112 * 00 01 02 03 04 00
113 * 11 12 13 14 10 11
114 * 22 23 24 20 21 22
115 * 33 34 30 31 32 33
116 * 44 40 41 42 43 44
117 *
118 *
119 * Let TRANSR = 'N'. RFP holds AP as follows:
120 * For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
121 * three columns of AP upper. The lower triangle A(3:4,0:1) consists of
122 * the transpose of the first two columns of AP upper.
123 * For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
124 * three columns of AP lower. The upper triangle A(0:1,1:2) consists of
125 * the transpose of the last two columns of AP lower.
126 * This covers the case N odd and TRANSR = 'N'.
127 *
128 * RFP A RFP A
129 *
130 * 02 03 04 00 33 43
131 * 12 13 14 10 11 44
132 * 22 23 24 20 21 22
133 * 00 33 34 30 31 32
134 * 01 11 44 40 41 42
135 *
136 * Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
137 * transpose of RFP A above. One therefore gets:
138 *
139 * RFP A RFP A
140 *
141 * 02 12 22 00 01 00 10 20 30 40 50
142 * 03 13 23 33 11 33 11 21 31 41 51
143 * 04 14 24 34 44 43 44 22 32 42 52
144 *
145 * =====================================================================
146 *
147 * .. Parameters ..
148 DOUBLE PRECISION ONE
149 PARAMETER ( ONE = 1.0D+0 )
150 * ..
151 * .. Local Scalars ..
152 LOGICAL LOWER, NISODD, NORMALTRANSR
153 INTEGER N1, N2, K
154 * ..
155 * .. External Functions ..
156 LOGICAL LSAME
157 EXTERNAL LSAME
158 * ..
159 * .. External Subroutines ..
160 EXTERNAL XERBLA, DTFTRI, DLAUUM, DTRMM, DSYRK
161 * ..
162 * .. Intrinsic Functions ..
163 INTRINSIC MOD
164 * ..
165 * .. Executable Statements ..
166 *
167 * Test the input parameters.
168 *
169 INFO = 0
170 NORMALTRANSR = LSAME( TRANSR, 'N' )
171 LOWER = LSAME( UPLO, 'L' )
172 IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'T' ) ) THEN
173 INFO = -1
174 ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
175 INFO = -2
176 ELSE IF( N.LT.0 ) THEN
177 INFO = -3
178 END IF
179 IF( INFO.NE.0 ) THEN
180 CALL XERBLA( 'DPFTRI', -INFO )
181 RETURN
182 END IF
183 *
184 * Quick return if possible
185 *
186 IF( N.EQ.0 )
187 $ RETURN
188 *
189 * Invert the triangular Cholesky factor U or L.
190 *
191 CALL DTFTRI( TRANSR, UPLO, 'N', N, A, INFO )
192 IF( INFO.GT.0 )
193 $ RETURN
194 *
195 * If N is odd, set NISODD = .TRUE.
196 * If N is even, set K = N/2 and NISODD = .FALSE.
197 *
198 IF( MOD( N, 2 ).EQ.0 ) THEN
199 K = N / 2
200 NISODD = .FALSE.
201 ELSE
202 NISODD = .TRUE.
203 END IF
204 *
205 * Set N1 and N2 depending on LOWER
206 *
207 IF( LOWER ) THEN
208 N2 = N / 2
209 N1 = N - N2
210 ELSE
211 N1 = N / 2
212 N2 = N - N1
213 END IF
214 *
215 * Start execution of triangular matrix multiply: inv(U)*inv(U)^C or
216 * inv(L)^C*inv(L). There are eight cases.
217 *
218 IF( NISODD ) THEN
219 *
220 * N is odd
221 *
222 IF( NORMALTRANSR ) THEN
223 *
224 * N is odd and TRANSR = 'N'
225 *
226 IF( LOWER ) THEN
227 *
228 * SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:N1-1) )
229 * T1 -> a(0,0), T2 -> a(0,1), S -> a(N1,0)
230 * T1 -> a(0), T2 -> a(n), S -> a(N1)
231 *
232 CALL DLAUUM( 'L', N1, A( 0 ), N, INFO )
233 CALL DSYRK( 'L', 'T', N1, N2, ONE, A( N1 ), N, ONE,
234 $ A( 0 ), N )
235 CALL DTRMM( 'L', 'U', 'N', 'N', N2, N1, ONE, A( N ), N,
236 $ A( N1 ), N )
237 CALL DLAUUM( 'U', N2, A( N ), N, INFO )
238 *
239 ELSE
240 *
241 * SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:N2-1)
242 * T1 -> a(N1+1,0), T2 -> a(N1,0), S -> a(0,0)
243 * T1 -> a(N2), T2 -> a(N1), S -> a(0)
244 *
245 CALL DLAUUM( 'L', N1, A( N2 ), N, INFO )
246 CALL DSYRK( 'L', 'N', N1, N2, ONE, A( 0 ), N, ONE,
247 $ A( N2 ), N )
248 CALL DTRMM( 'R', 'U', 'T', 'N', N1, N2, ONE, A( N1 ), N,
249 $ A( 0 ), N )
250 CALL DLAUUM( 'U', N2, A( N1 ), N, INFO )
251 *
252 END IF
253 *
254 ELSE
255 *
256 * N is odd and TRANSR = 'T'
257 *
258 IF( LOWER ) THEN
259 *
260 * SRPA for LOWER, TRANSPOSE, and N is odd
261 * T1 -> a(0), T2 -> a(1), S -> a(0+N1*N1)
262 *
263 CALL DLAUUM( 'U', N1, A( 0 ), N1, INFO )
264 CALL DSYRK( 'U', 'N', N1, N2, ONE, A( N1*N1 ), N1, ONE,
265 $ A( 0 ), N1 )
266 CALL DTRMM( 'R', 'L', 'N', 'N', N1, N2, ONE, A( 1 ), N1,
267 $ A( N1*N1 ), N1 )
268 CALL DLAUUM( 'L', N2, A( 1 ), N1, INFO )
269 *
270 ELSE
271 *
272 * SRPA for UPPER, TRANSPOSE, and N is odd
273 * T1 -> a(0+N2*N2), T2 -> a(0+N1*N2), S -> a(0)
274 *
275 CALL DLAUUM( 'U', N1, A( N2*N2 ), N2, INFO )
276 CALL DSYRK( 'U', 'T', N1, N2, ONE, A( 0 ), N2, ONE,
277 $ A( N2*N2 ), N2 )
278 CALL DTRMM( 'L', 'L', 'T', 'N', N2, N1, ONE, A( N1*N2 ),
279 $ N2, A( 0 ), N2 )
280 CALL DLAUUM( 'L', N2, A( N1*N2 ), N2, INFO )
281 *
282 END IF
283 *
284 END IF
285 *
286 ELSE
287 *
288 * N is even
289 *
290 IF( NORMALTRANSR ) THEN
291 *
292 * N is even and TRANSR = 'N'
293 *
294 IF( LOWER ) THEN
295 *
296 * SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) )
297 * T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0)
298 * T1 -> a(1), T2 -> a(0), S -> a(k+1)
299 *
300 CALL DLAUUM( 'L', K, A( 1 ), N+1, INFO )
301 CALL DSYRK( 'L', 'T', K, K, ONE, A( K+1 ), N+1, ONE,
302 $ A( 1 ), N+1 )
303 CALL DTRMM( 'L', 'U', 'N', 'N', K, K, ONE, A( 0 ), N+1,
304 $ A( K+1 ), N+1 )
305 CALL DLAUUM( 'U', K, A( 0 ), N+1, INFO )
306 *
307 ELSE
308 *
309 * SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) )
310 * T1 -> a(k+1,0) , T2 -> a(k,0), S -> a(0,0)
311 * T1 -> a(k+1), T2 -> a(k), S -> a(0)
312 *
313 CALL DLAUUM( 'L', K, A( K+1 ), N+1, INFO )
314 CALL DSYRK( 'L', 'N', K, K, ONE, A( 0 ), N+1, ONE,
315 $ A( K+1 ), N+1 )
316 CALL DTRMM( 'R', 'U', 'T', 'N', K, K, ONE, A( K ), N+1,
317 $ A( 0 ), N+1 )
318 CALL DLAUUM( 'U', K, A( K ), N+1, INFO )
319 *
320 END IF
321 *
322 ELSE
323 *
324 * N is even and TRANSR = 'T'
325 *
326 IF( LOWER ) THEN
327 *
328 * SRPA for LOWER, TRANSPOSE, and N is even (see paper)
329 * T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1),
330 * T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k
331 *
332 CALL DLAUUM( 'U', K, A( K ), K, INFO )
333 CALL DSYRK( 'U', 'N', K, K, ONE, A( K*( K+1 ) ), K, ONE,
334 $ A( K ), K )
335 CALL DTRMM( 'R', 'L', 'N', 'N', K, K, ONE, A( 0 ), K,
336 $ A( K*( K+1 ) ), K )
337 CALL DLAUUM( 'L', K, A( 0 ), K, INFO )
338 *
339 ELSE
340 *
341 * SRPA for UPPER, TRANSPOSE, and N is even (see paper)
342 * T1 -> B(0,k+1), T2 -> B(0,k), S -> B(0,0),
343 * T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k
344 *
345 CALL DLAUUM( 'U', K, A( K*( K+1 ) ), K, INFO )
346 CALL DSYRK( 'U', 'T', K, K, ONE, A( 0 ), K, ONE,
347 $ A( K*( K+1 ) ), K )
348 CALL DTRMM( 'L', 'L', 'T', 'N', K, K, ONE, A( K*K ), K,
349 $ A( 0 ), K )
350 CALL DLAUUM( 'L', K, A( K*K ), K, INFO )
351 *
352 END IF
353 *
354 END IF
355 *
356 END IF
357 *
358 RETURN
359 *
360 * End of DPFTRI
361 *
362 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 DOUBLE PRECISION A( 0: * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DPFTRI computes the inverse of a (real) symmetric positive definite
22 * matrix A using the Cholesky factorization A = U**T*U or A = L*L**T
23 * computed by DPFTRF.
24 *
25 * Arguments
26 * =========
27 *
28 * TRANSR (input) CHARACTER*1
29 * = 'N': The Normal TRANSR of RFP A is stored;
30 * = 'T': The 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) DOUBLE PRECISION array, dimension ( N*(N+1)/2 )
40 * On entry, the symmetric 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 = 'T' then RFP is
44 * the 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 * 'T'. 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 symmetric 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 Rectangular Full Packed (RFP) Format when N is
65 * even. 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 * the 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 * the transpose of the last three columns of AP lower.
84 * This covers the case N even and TRANSR = 'N'.
85 *
86 * RFP A RFP A
87 *
88 * 03 04 05 33 43 53
89 * 13 14 15 00 44 54
90 * 23 24 25 10 11 55
91 * 33 34 35 20 21 22
92 * 00 44 45 30 31 32
93 * 01 11 55 40 41 42
94 * 02 12 22 50 51 52
95 *
96 * Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
97 * transpose of RFP A above. One therefore gets:
98 *
99 *
100 * RFP A RFP A
101 *
102 * 03 13 23 33 00 01 02 33 00 10 20 30 40 50
103 * 04 14 24 34 44 11 12 43 44 11 21 31 41 51
104 * 05 15 25 35 45 55 22 53 54 55 22 32 42 52
105 *
106 *
107 * We then consider Rectangular Full Packed (RFP) Format when N is
108 * odd. We give an example where N = 5.
109 *
110 * AP is Upper AP is Lower
111 *
112 * 00 01 02 03 04 00
113 * 11 12 13 14 10 11
114 * 22 23 24 20 21 22
115 * 33 34 30 31 32 33
116 * 44 40 41 42 43 44
117 *
118 *
119 * Let TRANSR = 'N'. RFP holds AP as follows:
120 * For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
121 * three columns of AP upper. The lower triangle A(3:4,0:1) consists of
122 * the transpose of the first two columns of AP upper.
123 * For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
124 * three columns of AP lower. The upper triangle A(0:1,1:2) consists of
125 * the transpose of the last two columns of AP lower.
126 * This covers the case N odd and TRANSR = 'N'.
127 *
128 * RFP A RFP A
129 *
130 * 02 03 04 00 33 43
131 * 12 13 14 10 11 44
132 * 22 23 24 20 21 22
133 * 00 33 34 30 31 32
134 * 01 11 44 40 41 42
135 *
136 * Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
137 * transpose of RFP A above. One therefore gets:
138 *
139 * RFP A RFP A
140 *
141 * 02 12 22 00 01 00 10 20 30 40 50
142 * 03 13 23 33 11 33 11 21 31 41 51
143 * 04 14 24 34 44 43 44 22 32 42 52
144 *
145 * =====================================================================
146 *
147 * .. Parameters ..
148 DOUBLE PRECISION ONE
149 PARAMETER ( ONE = 1.0D+0 )
150 * ..
151 * .. Local Scalars ..
152 LOGICAL LOWER, NISODD, NORMALTRANSR
153 INTEGER N1, N2, K
154 * ..
155 * .. External Functions ..
156 LOGICAL LSAME
157 EXTERNAL LSAME
158 * ..
159 * .. External Subroutines ..
160 EXTERNAL XERBLA, DTFTRI, DLAUUM, DTRMM, DSYRK
161 * ..
162 * .. Intrinsic Functions ..
163 INTRINSIC MOD
164 * ..
165 * .. Executable Statements ..
166 *
167 * Test the input parameters.
168 *
169 INFO = 0
170 NORMALTRANSR = LSAME( TRANSR, 'N' )
171 LOWER = LSAME( UPLO, 'L' )
172 IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'T' ) ) THEN
173 INFO = -1
174 ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
175 INFO = -2
176 ELSE IF( N.LT.0 ) THEN
177 INFO = -3
178 END IF
179 IF( INFO.NE.0 ) THEN
180 CALL XERBLA( 'DPFTRI', -INFO )
181 RETURN
182 END IF
183 *
184 * Quick return if possible
185 *
186 IF( N.EQ.0 )
187 $ RETURN
188 *
189 * Invert the triangular Cholesky factor U or L.
190 *
191 CALL DTFTRI( TRANSR, UPLO, 'N', N, A, INFO )
192 IF( INFO.GT.0 )
193 $ RETURN
194 *
195 * If N is odd, set NISODD = .TRUE.
196 * If N is even, set K = N/2 and NISODD = .FALSE.
197 *
198 IF( MOD( N, 2 ).EQ.0 ) THEN
199 K = N / 2
200 NISODD = .FALSE.
201 ELSE
202 NISODD = .TRUE.
203 END IF
204 *
205 * Set N1 and N2 depending on LOWER
206 *
207 IF( LOWER ) THEN
208 N2 = N / 2
209 N1 = N - N2
210 ELSE
211 N1 = N / 2
212 N2 = N - N1
213 END IF
214 *
215 * Start execution of triangular matrix multiply: inv(U)*inv(U)^C or
216 * inv(L)^C*inv(L). There are eight cases.
217 *
218 IF( NISODD ) THEN
219 *
220 * N is odd
221 *
222 IF( NORMALTRANSR ) THEN
223 *
224 * N is odd and TRANSR = 'N'
225 *
226 IF( LOWER ) THEN
227 *
228 * SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:N1-1) )
229 * T1 -> a(0,0), T2 -> a(0,1), S -> a(N1,0)
230 * T1 -> a(0), T2 -> a(n), S -> a(N1)
231 *
232 CALL DLAUUM( 'L', N1, A( 0 ), N, INFO )
233 CALL DSYRK( 'L', 'T', N1, N2, ONE, A( N1 ), N, ONE,
234 $ A( 0 ), N )
235 CALL DTRMM( 'L', 'U', 'N', 'N', N2, N1, ONE, A( N ), N,
236 $ A( N1 ), N )
237 CALL DLAUUM( 'U', N2, A( N ), N, INFO )
238 *
239 ELSE
240 *
241 * SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:N2-1)
242 * T1 -> a(N1+1,0), T2 -> a(N1,0), S -> a(0,0)
243 * T1 -> a(N2), T2 -> a(N1), S -> a(0)
244 *
245 CALL DLAUUM( 'L', N1, A( N2 ), N, INFO )
246 CALL DSYRK( 'L', 'N', N1, N2, ONE, A( 0 ), N, ONE,
247 $ A( N2 ), N )
248 CALL DTRMM( 'R', 'U', 'T', 'N', N1, N2, ONE, A( N1 ), N,
249 $ A( 0 ), N )
250 CALL DLAUUM( 'U', N2, A( N1 ), N, INFO )
251 *
252 END IF
253 *
254 ELSE
255 *
256 * N is odd and TRANSR = 'T'
257 *
258 IF( LOWER ) THEN
259 *
260 * SRPA for LOWER, TRANSPOSE, and N is odd
261 * T1 -> a(0), T2 -> a(1), S -> a(0+N1*N1)
262 *
263 CALL DLAUUM( 'U', N1, A( 0 ), N1, INFO )
264 CALL DSYRK( 'U', 'N', N1, N2, ONE, A( N1*N1 ), N1, ONE,
265 $ A( 0 ), N1 )
266 CALL DTRMM( 'R', 'L', 'N', 'N', N1, N2, ONE, A( 1 ), N1,
267 $ A( N1*N1 ), N1 )
268 CALL DLAUUM( 'L', N2, A( 1 ), N1, INFO )
269 *
270 ELSE
271 *
272 * SRPA for UPPER, TRANSPOSE, and N is odd
273 * T1 -> a(0+N2*N2), T2 -> a(0+N1*N2), S -> a(0)
274 *
275 CALL DLAUUM( 'U', N1, A( N2*N2 ), N2, INFO )
276 CALL DSYRK( 'U', 'T', N1, N2, ONE, A( 0 ), N2, ONE,
277 $ A( N2*N2 ), N2 )
278 CALL DTRMM( 'L', 'L', 'T', 'N', N2, N1, ONE, A( N1*N2 ),
279 $ N2, A( 0 ), N2 )
280 CALL DLAUUM( 'L', N2, A( N1*N2 ), N2, INFO )
281 *
282 END IF
283 *
284 END IF
285 *
286 ELSE
287 *
288 * N is even
289 *
290 IF( NORMALTRANSR ) THEN
291 *
292 * N is even and TRANSR = 'N'
293 *
294 IF( LOWER ) THEN
295 *
296 * SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) )
297 * T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0)
298 * T1 -> a(1), T2 -> a(0), S -> a(k+1)
299 *
300 CALL DLAUUM( 'L', K, A( 1 ), N+1, INFO )
301 CALL DSYRK( 'L', 'T', K, K, ONE, A( K+1 ), N+1, ONE,
302 $ A( 1 ), N+1 )
303 CALL DTRMM( 'L', 'U', 'N', 'N', K, K, ONE, A( 0 ), N+1,
304 $ A( K+1 ), N+1 )
305 CALL DLAUUM( 'U', K, A( 0 ), N+1, INFO )
306 *
307 ELSE
308 *
309 * SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) )
310 * T1 -> a(k+1,0) , T2 -> a(k,0), S -> a(0,0)
311 * T1 -> a(k+1), T2 -> a(k), S -> a(0)
312 *
313 CALL DLAUUM( 'L', K, A( K+1 ), N+1, INFO )
314 CALL DSYRK( 'L', 'N', K, K, ONE, A( 0 ), N+1, ONE,
315 $ A( K+1 ), N+1 )
316 CALL DTRMM( 'R', 'U', 'T', 'N', K, K, ONE, A( K ), N+1,
317 $ A( 0 ), N+1 )
318 CALL DLAUUM( 'U', K, A( K ), N+1, INFO )
319 *
320 END IF
321 *
322 ELSE
323 *
324 * N is even and TRANSR = 'T'
325 *
326 IF( LOWER ) THEN
327 *
328 * SRPA for LOWER, TRANSPOSE, and N is even (see paper)
329 * T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1),
330 * T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k
331 *
332 CALL DLAUUM( 'U', K, A( K ), K, INFO )
333 CALL DSYRK( 'U', 'N', K, K, ONE, A( K*( K+1 ) ), K, ONE,
334 $ A( K ), K )
335 CALL DTRMM( 'R', 'L', 'N', 'N', K, K, ONE, A( 0 ), K,
336 $ A( K*( K+1 ) ), K )
337 CALL DLAUUM( 'L', K, A( 0 ), K, INFO )
338 *
339 ELSE
340 *
341 * SRPA for UPPER, TRANSPOSE, and N is even (see paper)
342 * T1 -> B(0,k+1), T2 -> B(0,k), S -> B(0,0),
343 * T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k
344 *
345 CALL DLAUUM( 'U', K, A( K*( K+1 ) ), K, INFO )
346 CALL DSYRK( 'U', 'T', K, K, ONE, A( 0 ), K, ONE,
347 $ A( K*( K+1 ) ), K )
348 CALL DTRMM( 'L', 'L', 'T', 'N', K, K, ONE, A( K*K ), K,
349 $ A( 0 ), K )
350 CALL DLAUUM( 'L', K, A( K*K ), K, INFO )
351 *
352 END IF
353 *
354 END IF
355 *
356 END IF
357 *
358 RETURN
359 *
360 * End of DPFTRI
361 *
362 END