001:       SUBROUTINE ZPFTRI( TRANSR, UPLO, N, A, INFO )
002: *
003: *  -- LAPACK routine (version 3.2)                                    --
004: *
005: *  -- Contributed by Fred Gustavson of the IBM Watson Research Center --
006: *  -- November 2008                                                   --
007: *
008: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
009: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
010: *
011: *     .. Scalar Arguments ..
012:       CHARACTER          TRANSR, UPLO
013:       INTEGER            INFO, N
014: *     .. Array Arguments ..
015:       COMPLEX*16         A( 0: * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  ZPFTRI computes the inverse of a complex Hermitian positive definite
022: *  matrix A using the Cholesky factorization A = U**H*U or A = L*L**H
023: *  computed by ZPFTRF.
024: *
025: *  Arguments
026: *  =========
027: *
028: *  TRANSR    (input) CHARACTER
029: *          = 'N':  The Normal TRANSR of RFP A is stored;
030: *          = 'C':  The Conjugate-transpose TRANSR of RFP A is stored.
031: *
032: *  UPLO    (input) CHARACTER
033: *          = 'U':  Upper triangle of A is stored;
034: *          = 'L':  Lower triangle of A is stored.
035: *
036: *  N       (input) INTEGER
037: *          The order of the matrix A.  N >= 0.
038: *
039: *  A       (input/output) COMPLEX*16 array, dimension ( N*(N+1)/2 );
040: *          On entry, the Hermitian matrix A in RFP format. RFP format is
041: *          described by TRANSR, UPLO, and N as follows: If TRANSR = 'N'
042: *          then RFP A is (0:N,0:k-1) when N is even; k=N/2. RFP A is
043: *          (0:N-1,0:k) when N is odd; k=N/2. IF TRANSR = 'C' then RFP is
044: *          the Conjugate-transpose of RFP A as defined when
045: *          TRANSR = 'N'. The contents of RFP A are defined by UPLO as
046: *          follows: If UPLO = 'U' the RFP A contains the nt elements of
047: *          upper packed A. If UPLO = 'L' the RFP A contains the elements
048: *          of lower packed A. The LDA of RFP A is (N+1)/2 when TRANSR =
049: *          'C'. When TRANSR is 'N' the LDA is N+1 when N is even and N
050: *          is odd. See the Note below for more details.
051: *
052: *          On exit, the Hermitian inverse of the original matrix, in the
053: *          same storage format.
054: *
055: *  INFO    (output) INTEGER
056: *          = 0:  successful exit
057: *          < 0:  if INFO = -i, the i-th argument had an illegal value
058: *          > 0:  if INFO = i, the (i,i) element of the factor U or L is
059: *                zero, and the inverse could not be computed.
060: *
061: *  Note:
062: *  =====
063: *
064: *  We first consider Standard Packed Format when N is even.
065: *  We give an example where N = 6.
066: *
067: *      AP is Upper             AP is Lower
068: *
069: *   00 01 02 03 04 05       00
070: *      11 12 13 14 15       10 11
071: *         22 23 24 25       20 21 22
072: *            33 34 35       30 31 32 33
073: *               44 45       40 41 42 43 44
074: *                  55       50 51 52 53 54 55
075: *
076: *
077: *  Let TRANSR = 'N'. RFP holds AP as follows:
078: *  For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
079: *  three columns of AP upper. The lower triangle A(4:6,0:2) consists of
080: *  conjugate-transpose of the first three columns of AP upper.
081: *  For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
082: *  three columns of AP lower. The upper triangle A(0:2,0:2) consists of
083: *  conjugate-transpose of the last three columns of AP lower.
084: *  To denote conjugate we place -- above the element. This covers the
085: *  case N even and TRANSR = 'N'.
086: *
087: *         RFP A                   RFP A
088: *
089: *                                -- -- --
090: *        03 04 05                33 43 53
091: *                                   -- --
092: *        13 14 15                00 44 54
093: *                                      --
094: *        23 24 25                10 11 55
095: *
096: *        33 34 35                20 21 22
097: *        --
098: *        00 44 45                30 31 32
099: *        -- --
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
385: