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