001:       SUBROUTINE ZTFTRI( TRANSR, UPLO, DIAG, 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, DIAG
013:       INTEGER            INFO, N
014: *     ..
015: *     .. Array Arguments ..
016:       COMPLEX*16         A( 0: * )
017: *     ..
018: *
019: *  Purpose
020: *  =======
021: *
022: *  ZTFTRI computes the inverse of a triangular matrix A stored in RFP
023: *  format.
024: *
025: *  This is a Level 3 BLAS version of the algorithm.
026: *
027: *  Arguments
028: *  =========
029: *
030: *  TRANSR    (input) CHARACTER
031: *          = 'N':  The Normal TRANSR of RFP A is stored;
032: *          = 'C':  The Conjugate-transpose TRANSR of RFP A is stored.
033: *
034: *  UPLO    (input) CHARACTER
035: *          = 'U':  A is upper triangular;
036: *          = 'L':  A is lower triangular.
037: *
038: *  DIAG    (input) CHARACTER
039: *          = 'N':  A is non-unit triangular;
040: *          = 'U':  A is unit triangular.
041: *
042: *  N       (input) INTEGER
043: *          The order of the matrix A.  N >= 0.
044: *
045: *  A       (input/output) COMPLEX*16 array, dimension ( N*(N+1)/2 );
046: *          On entry, the triangular matrix A in RFP format. RFP format
047: *          is described by TRANSR, UPLO, and N as follows: If TRANSR =
048: *          'N' then RFP A is (0:N,0:k-1) when N is even; k=N/2. RFP A is
049: *          (0:N-1,0:k) when N is odd; k=N/2. IF TRANSR = 'C' then RFP is
050: *          the Conjugate-transpose of RFP A as defined when
051: *          TRANSR = 'N'. The contents of RFP A are defined by UPLO as
052: *          follows: If UPLO = 'U' the RFP A contains the nt elements of
053: *          upper packed A; If UPLO = 'L' the RFP A contains the nt
054: *          elements of lower packed A. The LDA of RFP A is (N+1)/2 when
055: *          TRANSR = 'C'. When TRANSR is 'N' the LDA is N+1 when N is
056: *          even and N is odd. See the Note below for more details.
057: *
058: *          On exit, the (triangular) inverse of the original matrix, in
059: *          the same storage format.
060: *
061: *  INFO    (output) INTEGER
062: *          = 0: successful exit
063: *          < 0: if INFO = -i, the i-th argument had an illegal value
064: *          > 0: if INFO = i, A(i,i) is exactly zero.  The triangular
065: *               matrix is singular and its inverse can not be computed.
066: *
067: *  Notes:
068: *  ======
069: *
070: *  We first consider Standard Packed Format when N is even.
071: *  We give an example where N = 6.
072: *
073: *      AP is Upper             AP is Lower
074: *
075: *   00 01 02 03 04 05       00
076: *      11 12 13 14 15       10 11
077: *         22 23 24 25       20 21 22
078: *            33 34 35       30 31 32 33
079: *               44 45       40 41 42 43 44
080: *                  55       50 51 52 53 54 55
081: *
082: *
083: *  Let TRANSR = 'N'. RFP holds AP as follows:
084: *  For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
085: *  three columns of AP upper. The lower triangle A(4:6,0:2) consists of
086: *  conjugate-transpose of the first three columns of AP upper.
087: *  For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
088: *  three columns of AP lower. The upper triangle A(0:2,0:2) consists of
089: *  conjugate-transpose of the last three columns of AP lower.
090: *  To denote conjugate we place -- above the element. This covers the
091: *  case N even and TRANSR = 'N'.
092: *
093: *         RFP A                   RFP A
094: *
095: *                                -- -- --
096: *        03 04 05                33 43 53
097: *                                   -- --
098: *        13 14 15                00 44 54
099: *                                      --
100: *        23 24 25                10 11 55
101: *
102: *        33 34 35                20 21 22
103: *        --
104: *        00 44 45                30 31 32
105: *        -- --
106: *        01 11 55                40 41 42
107: *        -- -- --
108: *        02 12 22                50 51 52
109: *
110: *  Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
111: *  transpose of RFP A above. One therefore gets:
112: *
113: *
114: *           RFP A                   RFP A
115: *
116: *     -- -- -- --                -- -- -- -- -- --
117: *     03 13 23 33 00 01 02    33 00 10 20 30 40 50
118: *     -- -- -- -- --                -- -- -- -- --
119: *     04 14 24 34 44 11 12    43 44 11 21 31 41 51
120: *     -- -- -- -- -- --                -- -- -- --
121: *     05 15 25 35 45 55 22    53 54 55 22 32 42 52
122: *
123: *
124: *  We next  consider Standard Packed Format when N is odd.
125: *  We give an example where N = 5.
126: *
127: *     AP is Upper                 AP is Lower
128: *
129: *   00 01 02 03 04              00
130: *      11 12 13 14              10 11
131: *         22 23 24              20 21 22
132: *            33 34              30 31 32 33
133: *               44              40 41 42 43 44
134: *
135: *
136: *  Let TRANSR = 'N'. RFP holds AP as follows:
137: *  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
138: *  three columns of AP upper. The lower triangle A(3:4,0:1) consists of
139: *  conjugate-transpose of the first two   columns of AP upper.
140: *  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
141: *  three columns of AP lower. The upper triangle A(0:1,1:2) consists of
142: *  conjugate-transpose of the last two   columns of AP lower.
143: *  To denote conjugate we place -- above the element. This covers the
144: *  case N odd  and TRANSR = 'N'.
145: *
146: *         RFP A                   RFP A
147: *
148: *                                   -- --
149: *        02 03 04                00 33 43
150: *                                      --
151: *        12 13 14                10 11 44
152: *
153: *        22 23 24                20 21 22
154: *        --
155: *        00 33 34                30 31 32
156: *        -- --
157: *        01 11 44                40 41 42
158: *
159: *  Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
160: *  transpose of RFP A above. One therefore gets:
161: *
162: *
163: *           RFP A                   RFP A
164: *
165: *     -- -- --                   -- -- -- -- -- --
166: *     02 12 22 00 01             00 10 20 30 40 50
167: *     -- -- -- --                   -- -- -- -- --
168: *     03 13 23 33 11             33 11 21 31 41 51
169: *     -- -- -- -- --                   -- -- -- --
170: *     04 14 24 34 44             43 44 22 32 42 52
171: *
172: *  =====================================================================
173: *
174: *     .. Parameters ..
175:       COMPLEX*16         CONE
176:       PARAMETER          ( CONE = ( 1.0D+0, 0.0D+0 ) )
177: *     ..
178: *     .. Local Scalars ..
179:       LOGICAL            LOWER, NISODD, NORMALTRANSR
180:       INTEGER            N1, N2, K
181: *     ..
182: *     .. External Functions ..
183:       LOGICAL            LSAME
184:       EXTERNAL           LSAME
185: *     ..
186: *     .. External Subroutines ..
187:       EXTERNAL           XERBLA, ZTRMM, ZTRTRI
188: *     ..
189: *     .. Intrinsic Functions ..
190:       INTRINSIC          MOD
191: *     ..
192: *     .. Executable Statements ..
193: *
194: *     Test the input parameters.
195: *
196:       INFO = 0
197:       NORMALTRANSR = LSAME( TRANSR, 'N' )
198:       LOWER = LSAME( UPLO, 'L' )
199:       IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'C' ) ) THEN
200:          INFO = -1
201:       ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
202:          INFO = -2
203:       ELSE IF( .NOT.LSAME( DIAG, 'N' ) .AND. .NOT.LSAME( DIAG, 'U' ) )
204:      +         THEN
205:          INFO = -3
206:       ELSE IF( N.LT.0 ) THEN
207:          INFO = -4
208:       END IF
209:       IF( INFO.NE.0 ) THEN
210:          CALL XERBLA( 'ZTFTRI', -INFO )
211:          RETURN
212:       END IF
213: *
214: *     Quick return if possible
215: *
216:       IF( N.EQ.0 )
217:      +   RETURN
218: *
219: *     If N is odd, set NISODD = .TRUE.
220: *     If N is even, set K = N/2 and NISODD = .FALSE.
221: *
222:       IF( MOD( N, 2 ).EQ.0 ) THEN
223:          K = N / 2
224:          NISODD = .FALSE.
225:       ELSE
226:          NISODD = .TRUE.
227:       END IF
228: *
229: *     Set N1 and N2 depending on LOWER
230: *
231:       IF( LOWER ) THEN
232:          N2 = N / 2
233:          N1 = N - N2
234:       ELSE
235:          N1 = N / 2
236:          N2 = N - N1
237:       END IF
238: *
239: *
240: *     start execution: there are eight cases
241: *
242:       IF( NISODD ) THEN
243: *
244: *        N is odd
245: *
246:          IF( NORMALTRANSR ) THEN
247: *
248: *           N is odd and TRANSR = 'N'
249: *
250:             IF( LOWER ) THEN
251: *
252: *             SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:n1-1) )
253: *             T1 -> a(0,0), T2 -> a(0,1), S -> a(n1,0)
254: *             T1 -> a(0), T2 -> a(n), S -> a(n1)
255: *
256:                CALL ZTRTRI( 'L', DIAG, N1, A( 0 ), N, INFO )
257:                IF( INFO.GT.0 )
258:      +            RETURN
259:                CALL ZTRMM( 'R', 'L', 'N', DIAG, N2, N1, -CONE, A( 0 ),
260:      +                     N, A( N1 ), N )
261:                CALL ZTRTRI( 'U', DIAG, N2, A( N ), N, INFO )
262:                IF( INFO.GT.0 )
263:      +            INFO = INFO + N1
264:                IF( INFO.GT.0 )
265:      +            RETURN
266:                CALL ZTRMM( 'L', 'U', 'C', DIAG, N2, N1, CONE, A( N ), N,
267:      +                     A( N1 ), N )
268: *
269:             ELSE
270: *
271: *             SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:n2-1)
272: *             T1 -> a(n1+1,0), T2 -> a(n1,0), S -> a(0,0)
273: *             T1 -> a(n2), T2 -> a(n1), S -> a(0)
274: *
275:                CALL ZTRTRI( 'L', DIAG, N1, A( N2 ), N, INFO )
276:                IF( INFO.GT.0 )
277:      +            RETURN
278:                CALL ZTRMM( 'L', 'L', 'C', DIAG, N1, N2, -CONE, A( N2 ),
279:      +                     N, A( 0 ), N )
280:                CALL ZTRTRI( 'U', DIAG, N2, A( N1 ), N, INFO )
281:                IF( INFO.GT.0 )
282:      +            INFO = INFO + N1
283:                IF( INFO.GT.0 )
284:      +            RETURN
285:                CALL ZTRMM( 'R', 'U', 'N', DIAG, N1, N2, CONE, A( N1 ),
286:      +                     N, A( 0 ), N )
287: *
288:             END IF
289: *
290:          ELSE
291: *
292: *           N is odd and TRANSR = 'C'
293: *
294:             IF( LOWER ) THEN
295: *
296: *              SRPA for LOWER, TRANSPOSE and N is odd
297: *              T1 -> a(0), T2 -> a(1), S -> a(0+n1*n1)
298: *
299:                CALL ZTRTRI( 'U', DIAG, N1, A( 0 ), N1, INFO )
300:                IF( INFO.GT.0 )
301:      +            RETURN
302:                CALL ZTRMM( 'L', 'U', 'N', DIAG, N1, N2, -CONE, A( 0 ),
303:      +                     N1, A( N1*N1 ), N1 )
304:                CALL ZTRTRI( 'L', DIAG, N2, A( 1 ), N1, INFO )
305:                IF( INFO.GT.0 )
306:      +            INFO = INFO + N1
307:                IF( INFO.GT.0 )
308:      +            RETURN
309:                CALL ZTRMM( 'R', 'L', 'C', DIAG, N1, N2, CONE, A( 1 ),
310:      +                     N1, A( N1*N1 ), N1 )
311: *
312:             ELSE
313: *
314: *              SRPA for UPPER, TRANSPOSE and N is odd
315: *              T1 -> a(0+n2*n2), T2 -> a(0+n1*n2), S -> a(0)
316: *
317:                CALL ZTRTRI( 'U', DIAG, N1, A( N2*N2 ), N2, INFO )
318:                IF( INFO.GT.0 )
319:      +            RETURN
320:                CALL ZTRMM( 'R', 'U', 'C', DIAG, N2, N1, -CONE,
321:      +                     A( N2*N2 ), N2, A( 0 ), N2 )
322:                CALL ZTRTRI( 'L', DIAG, N2, A( N1*N2 ), N2, INFO )
323:                IF( INFO.GT.0 )
324:      +            INFO = INFO + N1
325:                IF( INFO.GT.0 )
326:      +            RETURN
327:                CALL ZTRMM( 'L', 'L', 'N', DIAG, N2, N1, CONE,
328:      +                     A( N1*N2 ), N2, A( 0 ), N2 )
329:             END IF
330: *
331:          END IF
332: *
333:       ELSE
334: *
335: *        N is even
336: *
337:          IF( NORMALTRANSR ) THEN
338: *
339: *           N is even and TRANSR = 'N'
340: *
341:             IF( LOWER ) THEN
342: *
343: *              SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) )
344: *              T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0)
345: *              T1 -> a(1), T2 -> a(0), S -> a(k+1)
346: *
347:                CALL ZTRTRI( 'L', DIAG, K, A( 1 ), N+1, INFO )
348:                IF( INFO.GT.0 )
349:      +            RETURN
350:                CALL ZTRMM( 'R', 'L', 'N', DIAG, K, K, -CONE, A( 1 ),
351:      +                     N+1, A( K+1 ), N+1 )
352:                CALL ZTRTRI( 'U', DIAG, K, A( 0 ), N+1, INFO )
353:                IF( INFO.GT.0 )
354:      +            INFO = INFO + K
355:                IF( INFO.GT.0 )
356:      +            RETURN
357:                CALL ZTRMM( 'L', 'U', 'C', DIAG, K, K, CONE, A( 0 ), N+1,
358:      +                     A( K+1 ), N+1 )
359: *
360:             ELSE
361: *
362: *              SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) )
363: *              T1 -> a(k+1,0) ,  T2 -> a(k,0),   S -> a(0,0)
364: *              T1 -> a(k+1), T2 -> a(k), S -> a(0)
365: *
366:                CALL ZTRTRI( 'L', DIAG, K, A( K+1 ), N+1, INFO )
367:                IF( INFO.GT.0 )
368:      +            RETURN
369:                CALL ZTRMM( 'L', 'L', 'C', DIAG, K, K, -CONE, A( K+1 ),
370:      +                     N+1, A( 0 ), N+1 )
371:                CALL ZTRTRI( 'U', DIAG, K, A( K ), N+1, INFO )
372:                IF( INFO.GT.0 )
373:      +            INFO = INFO + K
374:                IF( INFO.GT.0 )
375:      +            RETURN
376:                CALL ZTRMM( 'R', 'U', 'N', DIAG, K, K, CONE, A( K ), N+1,
377:      +                     A( 0 ), N+1 )
378:             END IF
379:          ELSE
380: *
381: *           N is even and TRANSR = 'C'
382: *
383:             IF( LOWER ) THEN
384: *
385: *              SRPA for LOWER, TRANSPOSE and N is even (see paper)
386: *              T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1)
387: *              T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k
388: *
389:                CALL ZTRTRI( 'U', DIAG, K, A( K ), K, INFO )
390:                IF( INFO.GT.0 )
391:      +            RETURN
392:                CALL ZTRMM( 'L', 'U', 'N', DIAG, K, K, -CONE, A( K ), K,
393:      +                     A( K*( K+1 ) ), K )
394:                CALL ZTRTRI( 'L', DIAG, K, A( 0 ), K, INFO )
395:                IF( INFO.GT.0 )
396:      +            INFO = INFO + K
397:                IF( INFO.GT.0 )
398:      +            RETURN
399:                CALL ZTRMM( 'R', 'L', 'C', DIAG, K, K, CONE, A( 0 ), K,
400:      +                     A( K*( K+1 ) ), K )
401:             ELSE
402: *
403: *              SRPA for UPPER, TRANSPOSE and N is even (see paper)
404: *              T1 -> B(0,k+1),     T2 -> B(0,k),   S -> B(0,0)
405: *              T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k
406: *
407:                CALL ZTRTRI( 'U', DIAG, K, A( K*( K+1 ) ), K, INFO )
408:                IF( INFO.GT.0 )
409:      +            RETURN
410:                CALL ZTRMM( 'R', 'U', 'C', DIAG, K, K, -CONE,
411:      +                     A( K*( K+1 ) ), K, A( 0 ), K )
412:                CALL ZTRTRI( 'L', DIAG, K, A( K*K ), K, INFO )
413:                IF( INFO.GT.0 )
414:      +            INFO = INFO + K
415:                IF( INFO.GT.0 )
416:      +            RETURN
417:                CALL ZTRMM( 'L', 'L', 'N', DIAG, K, K, CONE, A( K*K ), K,
418:      +                     A( 0 ), K )
419:             END IF
420:          END IF
421:       END IF
422: *
423:       RETURN
424: *
425: *     End of ZTFTRI
426: *
427:       END
428: