LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dtftri()

subroutine dtftri ( character  transr,
character  uplo,
character  diag,
integer  n,
double precision, dimension( 0: * )  a,
integer  info 
)

DTFTRI

Download DTFTRI + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 DTFTRI computes the inverse of a triangular matrix A stored in RFP
 format.

 This is a Level 3 BLAS version of the algorithm.
Parameters
[in]TRANSR
          TRANSR is CHARACTER*1
          = 'N':  The Normal TRANSR of RFP A is stored;
          = 'T':  The Transpose TRANSR of RFP A is stored.
[in]UPLO
          UPLO is CHARACTER*1
          = 'U':  A is upper triangular;
          = 'L':  A is lower triangular.
[in]DIAG
          DIAG is CHARACTER*1
          = 'N':  A is non-unit triangular;
          = 'U':  A is unit triangular.
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in,out]A
          A is DOUBLE PRECISION array, dimension (0:nt-1);
          nt=N*(N+1)/2. On entry, the triangular factor of a Hermitian
          Positive Definite matrix A in RFP format. RFP format is
          described by TRANSR, UPLO, and N as follows: If TRANSR = 'N'
          then RFP A is (0:N,0:k-1) when N is even; k=N/2. RFP A is
          (0:N-1,0:k) when N is odd; k=N/2. IF TRANSR = 'T' then RFP is
          the transpose of RFP A as defined when
          TRANSR = 'N'. The contents of RFP A are defined by UPLO as
          follows: If UPLO = 'U' the RFP A contains the nt elements of
          upper packed A; If UPLO = 'L' the RFP A contains the nt
          elements of lower packed A. The LDA of RFP A is (N+1)/2 when
          TRANSR = 'T'. When TRANSR is 'N' the LDA is N+1 when N is
          even and N is odd. See the Note below for more details.

          On exit, the (triangular) inverse of the original matrix, in
          the same storage format.
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -i, the i-th argument had an illegal value
          > 0: if INFO = i, A(i,i) is exactly zero.  The triangular
               matrix is singular and its inverse can not be computed.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
  We first consider Rectangular Full Packed (RFP) Format when N is
  even. We give an example where N = 6.

      AP is Upper             AP is Lower

   00 01 02 03 04 05       00
      11 12 13 14 15       10 11
         22 23 24 25       20 21 22
            33 34 35       30 31 32 33
               44 45       40 41 42 43 44
                  55       50 51 52 53 54 55


  Let TRANSR = 'N'. RFP holds AP as follows:
  For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
  three columns of AP upper. The lower triangle A(4:6,0:2) consists of
  the transpose of the first three columns of AP upper.
  For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
  three columns of AP lower. The upper triangle A(0:2,0:2) consists of
  the transpose of the last three columns of AP lower.
  This covers the case N even and TRANSR = 'N'.

         RFP A                   RFP A

        03 04 05                33 43 53
        13 14 15                00 44 54
        23 24 25                10 11 55
        33 34 35                20 21 22
        00 44 45                30 31 32
        01 11 55                40 41 42
        02 12 22                50 51 52

  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
  transpose of RFP A above. One therefore gets:


           RFP A                   RFP A

     03 13 23 33 00 01 02    33 00 10 20 30 40 50
     04 14 24 34 44 11 12    43 44 11 21 31 41 51
     05 15 25 35 45 55 22    53 54 55 22 32 42 52


  We then consider Rectangular Full Packed (RFP) Format when N is
  odd. We give an example where N = 5.

     AP is Upper                 AP is Lower

   00 01 02 03 04              00
      11 12 13 14              10 11
         22 23 24              20 21 22
            33 34              30 31 32 33
               44              40 41 42 43 44


  Let TRANSR = 'N'. RFP holds AP as follows:
  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
  three columns of AP upper. The lower triangle A(3:4,0:1) consists of
  the transpose of the first two columns of AP upper.
  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
  three columns of AP lower. The upper triangle A(0:1,1:2) consists of
  the transpose of the last two columns of AP lower.
  This covers the case N odd and TRANSR = 'N'.

         RFP A                   RFP A

        02 03 04                00 33 43
        12 13 14                10 11 44
        22 23 24                20 21 22
        00 33 34                30 31 32
        01 11 44                40 41 42

  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
  transpose of RFP A above. One therefore gets:

           RFP A                   RFP A

     02 12 22 00 01             00 10 20 30 40 50
     03 13 23 33 11             33 11 21 31 41 51
     04 14 24 34 44             43 44 22 32 42 52

Definition at line 200 of file dtftri.f.

201*
202* -- LAPACK computational routine --
203* -- LAPACK is a software package provided by Univ. of Tennessee, --
204* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
205*
206* .. Scalar Arguments ..
207 CHARACTER TRANSR, UPLO, DIAG
208 INTEGER INFO, N
209* ..
210* .. Array Arguments ..
211 DOUBLE PRECISION A( 0: * )
212* ..
213*
214* =====================================================================
215*
216* .. Parameters ..
217 DOUBLE PRECISION ONE
218 parameter( one = 1.0d+0 )
219* ..
220* .. Local Scalars ..
221 LOGICAL LOWER, NISODD, NORMALTRANSR
222 INTEGER N1, N2, K
223* ..
224* .. External Functions ..
225 LOGICAL LSAME
226 EXTERNAL lsame
227* ..
228* .. External Subroutines ..
229 EXTERNAL xerbla, dtrmm, dtrtri
230* ..
231* .. Intrinsic Functions ..
232 INTRINSIC mod
233* ..
234* .. Executable Statements ..
235*
236* Test the input parameters.
237*
238 info = 0
239 normaltransr = lsame( transr, 'N' )
240 lower = lsame( uplo, 'L' )
241 IF( .NOT.normaltransr .AND. .NOT.lsame( transr, 'T' ) ) THEN
242 info = -1
243 ELSE IF( .NOT.lower .AND. .NOT.lsame( uplo, 'U' ) ) THEN
244 info = -2
245 ELSE IF( .NOT.lsame( diag, 'N' ) .AND. .NOT.lsame( diag, 'U' ) )
246 $ THEN
247 info = -3
248 ELSE IF( n.LT.0 ) THEN
249 info = -4
250 END IF
251 IF( info.NE.0 ) THEN
252 CALL xerbla( 'DTFTRI', -info )
253 RETURN
254 END IF
255*
256* Quick return if possible
257*
258 IF( n.EQ.0 )
259 $ RETURN
260*
261* If N is odd, set NISODD = .TRUE.
262* If N is even, set K = N/2 and NISODD = .FALSE.
263*
264 IF( mod( n, 2 ).EQ.0 ) THEN
265 k = n / 2
266 nisodd = .false.
267 ELSE
268 nisodd = .true.
269 END IF
270*
271* Set N1 and N2 depending on LOWER
272*
273 IF( lower ) THEN
274 n2 = n / 2
275 n1 = n - n2
276 ELSE
277 n1 = n / 2
278 n2 = n - n1
279 END IF
280*
281*
282* start execution: there are eight cases
283*
284 IF( nisodd ) THEN
285*
286* N is odd
287*
288 IF( normaltransr ) THEN
289*
290* N is odd and TRANSR = 'N'
291*
292 IF( lower ) THEN
293*
294* SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:n1-1) )
295* T1 -> a(0,0), T2 -> a(0,1), S -> a(n1,0)
296* T1 -> a(0), T2 -> a(n), S -> a(n1)
297*
298 CALL dtrtri( 'L', diag, n1, a( 0 ), n, info )
299 IF( info.GT.0 )
300 $ RETURN
301 CALL dtrmm( 'R', 'L', 'N', diag, n2, n1, -one, a( 0 ),
302 $ n, a( n1 ), n )
303 CALL dtrtri( 'U', diag, n2, a( n ), n, info )
304 IF( info.GT.0 )
305 $ info = info + n1
306 IF( info.GT.0 )
307 $ RETURN
308 CALL dtrmm( 'L', 'U', 'T', diag, n2, n1, one, a( n ), n,
309 $ a( n1 ), n )
310*
311 ELSE
312*
313* SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:n2-1)
314* T1 -> a(n1+1,0), T2 -> a(n1,0), S -> a(0,0)
315* T1 -> a(n2), T2 -> a(n1), S -> a(0)
316*
317 CALL dtrtri( 'L', diag, n1, a( n2 ), n, info )
318 IF( info.GT.0 )
319 $ RETURN
320 CALL dtrmm( 'L', 'L', 'T', diag, n1, n2, -one, a( n2 ),
321 $ n, a( 0 ), n )
322 CALL dtrtri( 'U', diag, n2, a( n1 ), n, info )
323 IF( info.GT.0 )
324 $ info = info + n1
325 IF( info.GT.0 )
326 $ RETURN
327 CALL dtrmm( 'R', 'U', 'N', diag, n1, n2, one, a( n1 ),
328 $ n, a( 0 ), n )
329*
330 END IF
331*
332 ELSE
333*
334* N is odd and TRANSR = 'T'
335*
336 IF( lower ) THEN
337*
338* SRPA for LOWER, TRANSPOSE and N is odd
339* T1 -> a(0), T2 -> a(1), S -> a(0+n1*n1)
340*
341 CALL dtrtri( 'U', diag, n1, a( 0 ), n1, info )
342 IF( info.GT.0 )
343 $ RETURN
344 CALL dtrmm( 'L', 'U', 'N', diag, n1, n2, -one, a( 0 ),
345 $ n1, a( n1*n1 ), n1 )
346 CALL dtrtri( 'L', diag, n2, a( 1 ), n1, info )
347 IF( info.GT.0 )
348 $ info = info + n1
349 IF( info.GT.0 )
350 $ RETURN
351 CALL dtrmm( 'R', 'L', 'T', diag, n1, n2, one, a( 1 ),
352 $ n1, a( n1*n1 ), n1 )
353*
354 ELSE
355*
356* SRPA for UPPER, TRANSPOSE and N is odd
357* T1 -> a(0+n2*n2), T2 -> a(0+n1*n2), S -> a(0)
358*
359 CALL dtrtri( 'U', diag, n1, a( n2*n2 ), n2, info )
360 IF( info.GT.0 )
361 $ RETURN
362 CALL dtrmm( 'R', 'U', 'T', diag, n2, n1, -one,
363 $ a( n2*n2 ), n2, a( 0 ), n2 )
364 CALL dtrtri( 'L', diag, n2, a( n1*n2 ), n2, info )
365 IF( info.GT.0 )
366 $ info = info + n1
367 IF( info.GT.0 )
368 $ RETURN
369 CALL dtrmm( 'L', 'L', 'N', diag, n2, n1, one,
370 $ a( n1*n2 ), n2, a( 0 ), n2 )
371 END IF
372*
373 END IF
374*
375 ELSE
376*
377* N is even
378*
379 IF( normaltransr ) THEN
380*
381* N is even and TRANSR = 'N'
382*
383 IF( lower ) THEN
384*
385* SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) )
386* T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0)
387* T1 -> a(1), T2 -> a(0), S -> a(k+1)
388*
389 CALL dtrtri( 'L', diag, k, a( 1 ), n+1, info )
390 IF( info.GT.0 )
391 $ RETURN
392 CALL dtrmm( 'R', 'L', 'N', diag, k, k, -one, a( 1 ),
393 $ n+1, a( k+1 ), n+1 )
394 CALL dtrtri( 'U', diag, k, a( 0 ), n+1, info )
395 IF( info.GT.0 )
396 $ info = info + k
397 IF( info.GT.0 )
398 $ RETURN
399 CALL dtrmm( 'L', 'U', 'T', diag, k, k, one, a( 0 ), n+1,
400 $ a( k+1 ), n+1 )
401*
402 ELSE
403*
404* SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) )
405* T1 -> a(k+1,0) , T2 -> a(k,0), S -> a(0,0)
406* T1 -> a(k+1), T2 -> a(k), S -> a(0)
407*
408 CALL dtrtri( 'L', diag, k, a( k+1 ), n+1, info )
409 IF( info.GT.0 )
410 $ RETURN
411 CALL dtrmm( 'L', 'L', 'T', diag, k, k, -one, a( k+1 ),
412 $ n+1, a( 0 ), n+1 )
413 CALL dtrtri( 'U', diag, k, a( k ), n+1, info )
414 IF( info.GT.0 )
415 $ info = info + k
416 IF( info.GT.0 )
417 $ RETURN
418 CALL dtrmm( 'R', 'U', 'N', diag, k, k, one, a( k ), n+1,
419 $ a( 0 ), n+1 )
420 END IF
421 ELSE
422*
423* N is even and TRANSR = 'T'
424*
425 IF( lower ) THEN
426*
427* SRPA for LOWER, TRANSPOSE and N is even (see paper)
428* T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1)
429* T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k
430*
431 CALL dtrtri( 'U', diag, k, a( k ), k, info )
432 IF( info.GT.0 )
433 $ RETURN
434 CALL dtrmm( 'L', 'U', 'N', diag, k, k, -one, a( k ), k,
435 $ a( k*( k+1 ) ), k )
436 CALL dtrtri( 'L', diag, k, a( 0 ), k, info )
437 IF( info.GT.0 )
438 $ info = info + k
439 IF( info.GT.0 )
440 $ RETURN
441 CALL dtrmm( 'R', 'L', 'T', diag, k, k, one, a( 0 ), k,
442 $ a( k*( k+1 ) ), k )
443 ELSE
444*
445* SRPA for UPPER, TRANSPOSE and N is even (see paper)
446* T1 -> B(0,k+1), T2 -> B(0,k), S -> B(0,0)
447* T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k
448*
449 CALL dtrtri( 'U', diag, k, a( k*( k+1 ) ), k, info )
450 IF( info.GT.0 )
451 $ RETURN
452 CALL dtrmm( 'R', 'U', 'T', diag, k, k, -one,
453 $ a( k*( k+1 ) ), k, a( 0 ), k )
454 CALL dtrtri( 'L', diag, k, a( k*k ), k, info )
455 IF( info.GT.0 )
456 $ info = info + k
457 IF( info.GT.0 )
458 $ RETURN
459 CALL dtrmm( 'L', 'L', 'N', diag, k, k, one, a( k*k ), k,
460 $ a( 0 ), k )
461 END IF
462 END IF
463 END IF
464*
465 RETURN
466*
467* End of DTFTRI
468*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
dtrmm
subroutine dtrmm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
DTRMM
Definition dtrmm.f:177
dtrtri
subroutine dtrtri(uplo, diag, n, a, lda, info)
DTRTRI
Definition dtrtri.f:109
Here is the call graph for this function:
Here is the caller graph for this function: