MAGMA  magma-1.4.0 Matrix Algebra on GPU and Multicore Architectures
strtri.cpp File Reference
#include "common_magma.h"
Include dependency graph for strtri.cpp:

Go to the source code of this file.

## Macros

#define A(i, j)   ( A + (i) + (j)*lda )

#define dA(i, j)   (dA + (i) + (j)*ldda)

## Functions

magma_int_t magma_strtri (char uplo, char diag, magma_int_t n, float *A, magma_int_t lda, magma_int_t *info)

## Macro Definition Documentation

 #define A ( i, j ) ( A + (i) + (j)*lda )
 #define dA ( i, j ) (dA + (i) + (j)*ldda)

## Function Documentation

 magma_int_t magma_strtri ( char uplo, char diag, magma_int_t n, float * A, magma_int_t lda, magma_int_t * info )

Definition at line 14 of file strtri.cpp.

16 {
17 /* -- MAGMA (version 1.4.0) --
18  Univ. of Tennessee, Knoxville
19  Univ. of California, Berkeley
20  Univ. of Colorado, Denver
21  August 2013
22
23  Purpose
24  =======
25  STRTRI computes the inverse of a real upper or lower triangular
26  matrix A.
27
28  This is the Level 3 BLAS version of the algorithm.
29
30  Arguments
31  =========
32  UPLO (input) CHARACTER*1
33  = 'U': A is upper triangular;
34  = 'L': A is lower triangular.
35
36  DIAG (input) CHARACTER*1
37  = 'N': A is non-unit triangular;
38  = 'U': A is unit triangular.
39
40  N (input) INTEGER
41  The order of the matrix A. N >= 0.
42
43  A (input/output) REAL array, dimension (LDA,N)
44  On entry, the triangular matrix A. If UPLO = 'U', the
45  leading N-by-N upper triangular part of the array A contains
46  the upper triangular matrix, and the strictly lower
47  triangular part of A is not referenced. If UPLO = 'L', the
48  leading N-by-N lower triangular part of the array A contains
49  the lower triangular matrix, and the strictly upper
50  triangular part of A is not referenced. If DIAG = 'U', the
51  diagonal elements of A are also not referenced and are
52  assumed to be 1.
53  On exit, the (triangular) inverse of the original matrix, in
54  the same storage format.
55
56  LDA (input) INTEGER
57  The leading dimension of the array A. LDA >= max(1,N).
58
59  INFO (output) INTEGER
60  = 0: successful exit
61  < 0: if INFO = -i, the i-th argument had an illegal value
62  > 0: if INFO = i, A(i,i) is exactly zero. The triangular
63  matrix is singular and its inverse cannot be computed.
64
65  ===================================================================== */
66
67  #define A(i, j) ( A + (i) + (j)*lda )
68  #define dA(i, j) (dA + (i) + (j)*ldda)
69
70  /* Local variables */
71  char uplo_[2] = {uplo, 0};
72  char diag_[2] = {diag, 0};
73  magma_int_t ldda, nb, nn, j, jb;
74  float c_zero = MAGMA_S_ZERO;
75  float c_one = MAGMA_S_ONE;
76  float c_neg_one = MAGMA_S_NEG_ONE;
77  float *dA;
78
79  int upper = lapackf77_lsame(uplo_, "U");
80  int nounit = lapackf77_lsame(diag_, "N");
81
82  *info = 0;
83
84  if ((! upper) && (! lapackf77_lsame(uplo_, "L")))
85  *info = -1;
86  else if ((! nounit) && (! lapackf77_lsame(diag_, "U")))
87  *info = -2;
88  else if (n < 0)
89  *info = -3;
90  else if (lda < max(1,n))
91  *info = -5;
92
93  if (*info != 0) {
94  magma_xerbla( __func__, -(*info) );
95  return *info;
96  }
97
98  /* Quick return */
99  if ( n == 0 )
100  return *info;
101
102  /* Check for singularity if non-unit */
103  if (nounit) {
104  for ( j=0; j<n; ++j ) {
105  if ( MAGMA_S_EQUAL( *A(j,j), c_zero )) {
106  *info = j+1; // Fortran index
107  return *info;
108  }
109  }
110  }
111
112  /* Determine the block size for this environment */
113  nb = magma_get_spotrf_nb(n);
114
115  ldda = ((n+31)/32)*32;
116  if (MAGMA_SUCCESS != magma_smalloc( &dA, (n)*ldda )) {
117  *info = MAGMA_ERR_DEVICE_ALLOC;
118  return *info;
119  }
120
121  magma_queue_t stream[2];
122  magma_queue_create( &stream[0] );
123  magma_queue_create( &stream[1] );
124
125  if (nb <= 1 || nb >= n)
126  lapackf77_strtri(uplo_, diag_, &n, A, &lda, info);
127  else {
128  if (upper) {
129  /* Compute inverse of upper triangular matrix */
130  for (j=0; j<n; j=j+nb) {
131  jb = min(nb, (n-j));
132  magma_ssetmatrix( jb, (n-j),
133  A(j, j), lda,
134  dA(j, j), ldda );
135
136  /* Compute rows 1:j-1 of current block column */
138  MagmaNoTrans, MagmaNonUnit, j, jb,
139  c_one, dA(0,0), ldda, dA(0, j),ldda);
140
142  MagmaNoTrans, MagmaNonUnit, j, jb,
143  c_neg_one, dA(j,j), ldda, dA(0, j),ldda);
144
145  magma_sgetmatrix_async( jb, jb,
146  dA(j, j), ldda,
147  A(j, j), lda, stream[1] );
148
149  magma_sgetmatrix_async( j, jb,
150  dA(0, j), ldda,
151  A(0, j), lda, stream[0] );
152
153  magma_queue_sync( stream[1] );
154
155  /* Compute inverse of current diagonal block */
156  lapackf77_strtri(MagmaUpperStr, diag_, &jb, A(j,j), &lda, info);
157
158  magma_ssetmatrix( jb, jb,
159  A(j, j), lda,
160  dA(j, j), ldda );
161  }
162  }
163  else {
164  /* Compute inverse of lower triangular matrix */
165  nn=((n-1)/nb)*nb+1;
166
167  for(j=nn-1; j>=0; j=j-nb) {
168  jb=min(nb,(n-j));
169
170  if((j+jb) < n) {
171  magma_ssetmatrix( (n-j), jb,
172  A(j, j), lda,
173  dA(j, j), ldda );
174
175  /* Compute rows j+jb:n of current block column */
177  MagmaNoTrans, MagmaNonUnit, (n-j-jb), jb,
178  c_one, dA(j+jb,j+jb), ldda, dA(j+jb, j), ldda );
179
181  MagmaNoTrans, MagmaNonUnit, (n-j-jb), jb,
182  c_neg_one, dA(j,j), ldda, dA(j+jb, j), ldda );
183
184  magma_sgetmatrix_async( n-j-jb, jb,
185  dA(j+jb, j), ldda,
186  A(j+jb, j), lda, stream[1] );
187
188  magma_sgetmatrix_async( jb, jb,
189  dA(j,j), ldda,
190  A(j,j), lda, stream[0] );
191
192  magma_queue_sync( stream[0] );
193  }
194
195  /* Compute inverse of current diagonal block */
196  lapackf77_strtri(MagmaLowerStr, diag_, &jb, A(j,j), &lda, info);
197
198  magma_ssetmatrix( jb, jb,
199  A(j, j), lda,
200  dA(j, j), ldda );
201  }
202  }
203  }
204
205  magma_queue_destroy( stream[0] );
206  magma_queue_destroy( stream[1] );
207  magma_free( dA );
208
209  return *info;
210 }
#define MAGMA_S_EQUAL(u, v)
Definition: magma.h:183
void magma_strmm(magma_side_t side, magma_uplo_t uplo, magma_trans_t trans, magma_diag_t diag, magma_int_t m, magma_int_t n, float alpha, magmaFloat_const_ptr dA, magma_int_t ldda, magmaFloat_ptr dB, magma_int_t lddb)
#define min(a, b)
Definition: common_magma.h:86
#define MagmaLeft
Definition: magma.h:68
#define magma_queue_create(queuePtr)
Definition: magma.h:113
#define __func__
Definition: common_magma.h:65
#define MagmaUpper
Definition: magma.h:61
#define MAGMA_ERR_DEVICE_ALLOC
Definition: magma_types.h:276
#define magma_free(ptr)
Definition: magma.h:57
#define MAGMA_S_NEG_ONE
Definition: magma.h:200
int magma_int_t
Definition: magmablas.h:12
#define magma_sgetmatrix_async(m, n, dA_src, ldda, hB_dst, ldb, queue)
Definition: magmablas_s.h:714
void lapackf77_strtri(const char *uplo, const char *diag, const magma_int_t *n, float *A, const magma_int_t *lda, magma_int_t *info)
#define magma_queue_destroy(queue)
Definition: magma.h:116
#define magma_ssetmatrix(m, n, hA_src, lda, dB_dst, lddb)
Definition: magmablas_s.h:702
magma_int_t ldda
#define MagmaLower
Definition: magma.h:62
#define MAGMA_S_ONE
Definition: magma.h:198
void magma_xerbla(const char *srname, magma_int_t info)
Definition: xerbla.cpp:8
#define lapackf77_lsame
Definition: magma_lapack.h:23
#define MAGMA_S_ZERO
Definition: magma.h:197
void magma_strsm(magma_side_t side, magma_uplo_t uplo, magma_trans_t trans, magma_diag_t diag, magma_int_t m, magma_int_t n, float alpha, magmaFloat_const_ptr dA, magma_int_t ldda, magmaFloat_ptr dB, magma_int_t lddb)
#define MagmaLowerStr
Definition: magma.h:85
#define MagmaNonUnit
Definition: magma.h:65
#define MAGMA_SUCCESS
Definition: magma.h:106
#define A(i, j)
static magma_err_t magma_smalloc(magmaFloat_ptr *ptrPtr, size_t n)
Definition: magma.h:77
magma_int_t magma_get_spotrf_nb(magma_int_t m)
Definition: get_nb.cpp:29
#define MagmaRight
Definition: magma.h:69
#define MagmaUpperStr
Definition: magma.h:84
#define dA(i, j)
#define MagmaNoTrans
Definition: magma.h:57
#define max(a, b)
Definition: common_magma.h:82
#define magma_queue_sync(queue)
Definition: magma.h:119

Here is the call graph for this function:

Here is the caller graph for this function: