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