1*> \brief \b CTFTRI
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CTFTRI + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctftri.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctftri.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctftri.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE CTFTRI( TRANSR, UPLO, DIAG, N, A, INFO )
22*
23*       .. Scalar Arguments ..
24*       CHARACTER          TRANSR, UPLO, DIAG
25*       INTEGER            INFO, N
26*       ..
27*       .. Array Arguments ..
28*       COMPLEX            A( 0: * )
29*       ..
30*
31*
32*> \par Purpose:
33*  =============
34*>
35*> \verbatim
36*>
37*> CTFTRI computes the inverse of a triangular matrix A stored in RFP
38*> format.
39*>
40*> This is a Level 3 BLAS version of the algorithm.
41*> \endverbatim
42*
43*  Arguments:
44*  ==========
45*
46*> \param[in] TRANSR
47*> \verbatim
48*>          TRANSR is CHARACTER*1
49*>          = 'N':  The Normal TRANSR of RFP A is stored;
50*>          = 'C':  The Conjugate-transpose TRANSR of RFP A is stored.
51*> \endverbatim
52*>
53*> \param[in] UPLO
54*> \verbatim
55*>          UPLO is CHARACTER*1
56*>          = 'U':  A is upper triangular;
57*>          = 'L':  A is lower triangular.
58*> \endverbatim
59*>
60*> \param[in] DIAG
61*> \verbatim
62*>          DIAG is CHARACTER*1
63*>          = 'N':  A is non-unit triangular;
64*>          = 'U':  A is unit triangular.
65*> \endverbatim
66*>
67*> \param[in] N
68*> \verbatim
69*>          N is INTEGER
70*>          The order of the matrix A.  N >= 0.
71*> \endverbatim
72*>
73*> \param[in,out] A
74*> \verbatim
75*>          A is COMPLEX array, dimension ( N*(N+1)/2 );
76*>          On entry, the triangular matrix A in RFP format. RFP format
77*>          is described by TRANSR, UPLO, and N as follows: If TRANSR =
78*>          'N' then RFP A is (0:N,0:k-1) when N is even; k=N/2. RFP A is
79*>          (0:N-1,0:k) when N is odd; k=N/2. IF TRANSR = 'C' then RFP is
80*>          the Conjugate-transpose of RFP A as defined when
81*>          TRANSR = 'N'. The contents of RFP A are defined by UPLO as
82*>          follows: If UPLO = 'U' the RFP A contains the nt elements of
83*>          upper packed A; If UPLO = 'L' the RFP A contains the nt
84*>          elements of lower packed A. The LDA of RFP A is (N+1)/2 when
85*>          TRANSR = 'C'. When TRANSR is 'N' the LDA is N+1 when N is
86*>          even and N is odd. See the Note below for more details.
87*>
88*>          On exit, the (triangular) inverse of the original matrix, in
89*>          the same storage format.
90*> \endverbatim
91*>
92*> \param[out] INFO
93*> \verbatim
94*>          INFO is INTEGER
95*>          = 0: successful exit
96*>          < 0: if INFO = -i, the i-th argument had an illegal value
97*>          > 0: if INFO = i, A(i,i) is exactly zero.  The triangular
98*>               matrix is singular and its inverse can not be computed.
99*> \endverbatim
100*
101*  Authors:
102*  ========
103*
104*> \author Univ. of Tennessee
105*> \author Univ. of California Berkeley
106*> \author Univ. of Colorado Denver
107*> \author NAG Ltd.
108*
109*> \date December 2016
110*
111*> \ingroup complexOTHERcomputational
112*
113*> \par Further Details:
114*  =====================
115*>
116*> \verbatim
117*>
118*>  We first consider Standard Packed Format when N is even.
119*>  We give an example where N = 6.
120*>
121*>      AP is Upper             AP is Lower
122*>
123*>   00 01 02 03 04 05       00
124*>      11 12 13 14 15       10 11
125*>         22 23 24 25       20 21 22
126*>            33 34 35       30 31 32 33
127*>               44 45       40 41 42 43 44
128*>                  55       50 51 52 53 54 55
129*>
130*>
131*>  Let TRANSR = 'N'. RFP holds AP as follows:
132*>  For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
133*>  three columns of AP upper. The lower triangle A(4:6,0:2) consists of
134*>  conjugate-transpose of the first three columns of AP upper.
135*>  For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
136*>  three columns of AP lower. The upper triangle A(0:2,0:2) consists of
137*>  conjugate-transpose of the last three columns of AP lower.
138*>  To denote conjugate we place -- above the element. This covers the
139*>  case N even and TRANSR = 'N'.
140*>
141*>         RFP A                   RFP A
142*>
143*>                                -- -- --
144*>        03 04 05                33 43 53
145*>                                   -- --
146*>        13 14 15                00 44 54
147*>                                      --
148*>        23 24 25                10 11 55
149*>
150*>        33 34 35                20 21 22
151*>        --
152*>        00 44 45                30 31 32
153*>        -- --
154*>        01 11 55                40 41 42
155*>        -- -- --
156*>        02 12 22                50 51 52
157*>
158*>  Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
159*>  transpose of RFP A above. One therefore gets:
160*>
161*>
162*>           RFP A                   RFP A
163*>
164*>     -- -- -- --                -- -- -- -- -- --
165*>     03 13 23 33 00 01 02    33 00 10 20 30 40 50
166*>     -- -- -- -- --                -- -- -- -- --
167*>     04 14 24 34 44 11 12    43 44 11 21 31 41 51
168*>     -- -- -- -- -- --                -- -- -- --
169*>     05 15 25 35 45 55 22    53 54 55 22 32 42 52
170*>
171*>
172*>  We next  consider Standard Packed Format when N is odd.
173*>  We give an example where N = 5.
174*>
175*>     AP is Upper                 AP is Lower
176*>
177*>   00 01 02 03 04              00
178*>      11 12 13 14              10 11
179*>         22 23 24              20 21 22
180*>            33 34              30 31 32 33
181*>               44              40 41 42 43 44
182*>
183*>
184*>  Let TRANSR = 'N'. RFP holds AP as follows:
185*>  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
186*>  three columns of AP upper. The lower triangle A(3:4,0:1) consists of
187*>  conjugate-transpose of the first two   columns of AP upper.
188*>  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
189*>  three columns of AP lower. The upper triangle A(0:1,1:2) consists of
190*>  conjugate-transpose of the last two   columns of AP lower.
191*>  To denote conjugate we place -- above the element. This covers the
192*>  case N odd  and TRANSR = 'N'.
193*>
194*>         RFP A                   RFP A
195*>
196*>                                   -- --
197*>        02 03 04                00 33 43
198*>                                      --
199*>        12 13 14                10 11 44
200*>
201*>        22 23 24                20 21 22
202*>        --
203*>        00 33 34                30 31 32
204*>        -- --
205*>        01 11 44                40 41 42
206*>
207*>  Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
208*>  transpose of RFP A above. One therefore gets:
209*>
210*>
211*>           RFP A                   RFP A
212*>
213*>     -- -- --                   -- -- -- -- -- --
214*>     02 12 22 00 01             00 10 20 30 40 50
215*>     -- -- -- --                   -- -- -- -- --
216*>     03 13 23 33 11             33 11 21 31 41 51
217*>     -- -- -- -- --                   -- -- -- --
218*>     04 14 24 34 44             43 44 22 32 42 52
219*> \endverbatim
220*>
221*  =====================================================================
222      SUBROUTINE CTFTRI( TRANSR, UPLO, DIAG, N, A, INFO )
223*
224*  -- LAPACK computational routine (version 3.7.0) --
225*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
226*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
227*     December 2016
228*
229*     .. Scalar Arguments ..
230      CHARACTER          TRANSR, UPLO, DIAG
231      INTEGER            INFO, N
232*     ..
233*     .. Array Arguments ..
234      COMPLEX            A( 0: * )
235*     ..
236*
237*  =====================================================================
238*
239*     .. Parameters ..
240      COMPLEX            CONE
241      PARAMETER          ( CONE = ( 1.0E+0, 0.0E+0 ) )
242*     ..
243*     .. Local Scalars ..
244      LOGICAL            LOWER, NISODD, NORMALTRANSR
245      INTEGER            N1, N2, K
246*     ..
247*     .. External Functions ..
248      LOGICAL            LSAME
249      EXTERNAL           LSAME
250*     ..
251*     .. External Subroutines ..
252      EXTERNAL           XERBLA, CTRMM, CTRTRI
253*     ..
254*     .. Intrinsic Functions ..
255      INTRINSIC          MOD
256*     ..
257*     .. Executable Statements ..
258*
259*     Test the input parameters.
260*
261      INFO = 0
262      NORMALTRANSR = LSAME( TRANSR, 'N' )
263      LOWER = LSAME( UPLO, 'L' )
264      IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'C' ) ) THEN
265         INFO = -1
266      ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
267         INFO = -2
268      ELSE IF( .NOT.LSAME( DIAG, 'N' ) .AND. .NOT.LSAME( DIAG, 'U' ) )
269     $         THEN
270         INFO = -3
271      ELSE IF( N.LT.0 ) THEN
272         INFO = -4
273      END IF
274      IF( INFO.NE.0 ) THEN
275         CALL XERBLA( 'CTFTRI', -INFO )
276         RETURN
277      END IF
278*
279*     Quick return if possible
280*
281      IF( N.EQ.0 )
282     $   RETURN
283*
284*     If N is odd, set NISODD = .TRUE.
285*     If N is even, set K = N/2 and NISODD = .FALSE.
286*
287      IF( MOD( N, 2 ).EQ.0 ) THEN
288         K = N / 2
289         NISODD = .FALSE.
290      ELSE
291         NISODD = .TRUE.
292      END IF
293*
294*     Set N1 and N2 depending on LOWER
295*
296      IF( LOWER ) THEN
297         N2 = N / 2
298         N1 = N - N2
299      ELSE
300         N1 = N / 2
301         N2 = N - N1
302      END IF
303*
304*
305*     start execution: there are eight cases
306*
307      IF( NISODD ) THEN
308*
309*        N is odd
310*
311         IF( NORMALTRANSR ) THEN
312*
313*           N is odd and TRANSR = 'N'
314*
315            IF( LOWER ) THEN
316*
317*             SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:n1-1) )
318*             T1 -> a(0,0), T2 -> a(0,1), S -> a(n1,0)
319*             T1 -> a(0), T2 -> a(n), S -> a(n1)
320*
321               CALL CTRTRI( 'L', DIAG, N1, A( 0 ), N, INFO )
322               IF( INFO.GT.0 )
323     $            RETURN
324               CALL CTRMM( 'R', 'L', 'N', DIAG, N2, N1, -CONE, A( 0 ),
325     $                     N, A( N1 ), N )
326               CALL CTRTRI( 'U', DIAG, N2, A( N ), N, INFO )
327               IF( INFO.GT.0 )
328     $            INFO = INFO + N1
329               IF( INFO.GT.0 )
330     $            RETURN
331               CALL CTRMM( 'L', 'U', 'C', DIAG, N2, N1, CONE, A( N ), N,
332     $                     A( N1 ), N )
333*
334            ELSE
335*
336*             SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:n2-1)
337*             T1 -> a(n1+1,0), T2 -> a(n1,0), S -> a(0,0)
338*             T1 -> a(n2), T2 -> a(n1), S -> a(0)
339*
340               CALL CTRTRI( 'L', DIAG, N1, A( N2 ), N, INFO )
341               IF( INFO.GT.0 )
342     $            RETURN
343               CALL CTRMM( 'L', 'L', 'C', DIAG, N1, N2, -CONE, A( N2 ),
344     $                     N, A( 0 ), N )
345               CALL CTRTRI( 'U', DIAG, N2, A( N1 ), N, INFO )
346               IF( INFO.GT.0 )
347     $            INFO = INFO + N1
348               IF( INFO.GT.0 )
349     $            RETURN
350               CALL CTRMM( 'R', 'U', 'N', DIAG, N1, N2, CONE, A( N1 ),
351     $                     N, A( 0 ), N )
352*
353            END IF
354*
355         ELSE
356*
357*           N is odd and TRANSR = 'C'
358*
359            IF( LOWER ) THEN
360*
361*              SRPA for LOWER, TRANSPOSE and N is odd
362*              T1 -> a(0), T2 -> a(1), S -> a(0+n1*n1)
363*
364               CALL CTRTRI( 'U', DIAG, N1, A( 0 ), N1, INFO )
365               IF( INFO.GT.0 )
366     $            RETURN
367               CALL CTRMM( 'L', 'U', 'N', DIAG, N1, N2, -CONE, A( 0 ),
368     $                     N1, A( N1*N1 ), N1 )
369               CALL CTRTRI( 'L', DIAG, N2, A( 1 ), N1, INFO )
370               IF( INFO.GT.0 )
371     $            INFO = INFO + N1
372               IF( INFO.GT.0 )
373     $            RETURN
374               CALL CTRMM( 'R', 'L', 'C', DIAG, N1, N2, CONE, A( 1 ),
375     $                     N1, A( N1*N1 ), N1 )
376*
377            ELSE
378*
379*              SRPA for UPPER, TRANSPOSE and N is odd
380*              T1 -> a(0+n2*n2), T2 -> a(0+n1*n2), S -> a(0)
381*
382               CALL CTRTRI( 'U', DIAG, N1, A( N2*N2 ), N2, INFO )
383               IF( INFO.GT.0 )
384     $            RETURN
385               CALL CTRMM( 'R', 'U', 'C', DIAG, N2, N1, -CONE,
386     $                     A( N2*N2 ), N2, A( 0 ), N2 )
387               CALL CTRTRI( 'L', DIAG, N2, A( N1*N2 ), N2, INFO )
388               IF( INFO.GT.0 )
389     $            INFO = INFO + N1
390               IF( INFO.GT.0 )
391     $            RETURN
392               CALL CTRMM( 'L', 'L', 'N', DIAG, N2, N1, CONE,
393     $                     A( N1*N2 ), N2, A( 0 ), N2 )
394            END IF
395*
396         END IF
397*
398      ELSE
399*
400*        N is even
401*
402         IF( NORMALTRANSR ) THEN
403*
404*           N is even and TRANSR = 'N'
405*
406            IF( LOWER ) THEN
407*
408*              SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) )
409*              T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0)
410*              T1 -> a(1), T2 -> a(0), S -> a(k+1)
411*
412               CALL CTRTRI( 'L', DIAG, K, A( 1 ), N+1, INFO )
413               IF( INFO.GT.0 )
414     $            RETURN
415               CALL CTRMM( 'R', 'L', 'N', DIAG, K, K, -CONE, A( 1 ),
416     $                     N+1, A( K+1 ), N+1 )
417               CALL CTRTRI( 'U', DIAG, K, A( 0 ), N+1, INFO )
418               IF( INFO.GT.0 )
419     $            INFO = INFO + K
420               IF( INFO.GT.0 )
421     $            RETURN
422               CALL CTRMM( 'L', 'U', 'C', DIAG, K, K, CONE, A( 0 ), N+1,
423     $                     A( K+1 ), N+1 )
424*
425            ELSE
426*
427*              SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) )
428*              T1 -> a(k+1,0) ,  T2 -> a(k,0),   S -> a(0,0)
429*              T1 -> a(k+1), T2 -> a(k), S -> a(0)
430*
431               CALL CTRTRI( 'L', DIAG, K, A( K+1 ), N+1, INFO )
432               IF( INFO.GT.0 )
433     $            RETURN
434               CALL CTRMM( 'L', 'L', 'C', DIAG, K, K, -CONE, A( K+1 ),
435     $                     N+1, A( 0 ), N+1 )
436               CALL CTRTRI( 'U', DIAG, K, A( K ), N+1, INFO )
437               IF( INFO.GT.0 )
438     $            INFO = INFO + K
439               IF( INFO.GT.0 )
440     $            RETURN
441               CALL CTRMM( 'R', 'U', 'N', DIAG, K, K, CONE, A( K ), N+1,
442     $                     A( 0 ), N+1 )
443            END IF
444         ELSE
445*
446*           N is even and TRANSR = 'C'
447*
448            IF( LOWER ) THEN
449*
450*              SRPA for LOWER, TRANSPOSE and N is even (see paper)
451*              T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1)
452*              T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k
453*
454               CALL CTRTRI( 'U', DIAG, K, A( K ), K, INFO )
455               IF( INFO.GT.0 )
456     $            RETURN
457               CALL CTRMM( 'L', 'U', 'N', DIAG, K, K, -CONE, A( K ), K,
458     $                     A( K*( K+1 ) ), K )
459               CALL CTRTRI( 'L', DIAG, K, A( 0 ), K, INFO )
460               IF( INFO.GT.0 )
461     $            INFO = INFO + K
462               IF( INFO.GT.0 )
463     $            RETURN
464               CALL CTRMM( 'R', 'L', 'C', DIAG, K, K, CONE, A( 0 ), K,
465     $                     A( K*( K+1 ) ), K )
466            ELSE
467*
468*              SRPA for UPPER, TRANSPOSE and N is even (see paper)
469*              T1 -> B(0,k+1),     T2 -> B(0,k),   S -> B(0,0)
470*              T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k
471*
472               CALL CTRTRI( 'U', DIAG, K, A( K*( K+1 ) ), K, INFO )
473               IF( INFO.GT.0 )
474     $            RETURN
475               CALL CTRMM( 'R', 'U', 'C', DIAG, K, K, -CONE,
476     $                     A( K*( K+1 ) ), K, A( 0 ), K )
477               CALL CTRTRI( 'L', DIAG, K, A( K*K ), K, INFO )
478               IF( INFO.GT.0 )
479     $            INFO = INFO + K
480               IF( INFO.GT.0 )
481     $            RETURN
482               CALL CTRMM( 'L', 'L', 'N', DIAG, K, K, CONE, A( K*K ), K,
483     $                     A( 0 ), K )
484            END IF
485         END IF
486      END IF
487*
488      RETURN
489*
490*     End of CTFTRI
491*
492      END
493