MAGMA  magma-1.4.0 Matrix Algebra on GPU and Multicore Architectures
dgetri_gpu.cpp
Go to the documentation of this file.
1 /*
2  -- MAGMA (version 1.4.0) --
3  Univ. of Tennessee, Knoxville
4  Univ. of California, Berkeley
6  August 2013
7
8  @generated d Tue Aug 13 16:44:11 2013
9
10 */
11 #include "common_magma.h"
12
13 #define PRECISION_d
14
15 extern "C" magma_int_t
17  magma_int_t *ipiv, double *dwork, magma_int_t lwork,
18  magma_int_t *info )
19 {
20 /* -- MAGMA (version 1.4.0) --
21  Univ. of Tennessee, Knoxville
22  Univ. of California, Berkeley
24  August 2013
25
26  Purpose
27  =======
28  DGETRI computes the inverse of a matrix using the LU factorization
29  computed by DGETRF. This method inverts U and then computes inv(A) by
30  solving the system inv(A)*L = inv(U) for inv(A).
31
32  Note that it is generally both faster and more accurate to use DGESV,
33  or DGETRF and DGETRS, to solve the system AX = B, rather than inverting
34  the matrix and multiplying to form X = inv(A)*B. Only in special
35  instances should an explicit inverse be computed with this routine.
36
37  Arguments
38  =========
39  N (input) INTEGER
40  The order of the matrix A. N >= 0.
41
42  dA (input/output) DOUBLE_PRECISION array on the GPU, dimension (LDDA,N)
43  On entry, the factors L and U from the factorization
44  A = P*L*U as computed by DGETRF_GPU.
45  On exit, if INFO = 0, the inverse of the original matrix A.
46
47  LDDA (input) INTEGER
48  The leading dimension of the array A. LDDA >= max(1,N).
49
50  IPIV (input) INTEGER array, dimension (N)
51  The pivot indices from DGETRF; for 1<=i<=N, row i of the
52  matrix was interchanged with row IPIV(i).
53
54  DWORK (workspace/output) DOUBLE_PRECISION array on the GPU, dimension (MAX(1,LWORK))
55
56  LWORK (input) INTEGER
57  The dimension of the array DWORK. LWORK >= N*NB, where NB is
58  the optimal blocksize returned by magma_get_dgetri_nb(n).
59
60  Unlike LAPACK, this version does not currently support a
61  workspace query, because the workspace is on the GPU.
62
63  INFO (output) INTEGER
64  = 0: successful exit
65  < 0: if INFO = -i, the i-th argument had an illegal value
66  > 0: if INFO = i, U(i,i) is exactly zero; the matrix is
67  singular and its cannot be computed.
68
69  ===================================================================== */
70
71  #define dA(i, j) (dA + (i) + (j)*ldda)
72  #define dL(i, j) (dL + (i) + (j)*lddl)
73
74  /* Local variables */
75  double c_one = MAGMA_D_ONE;
76  double c_neg_one = MAGMA_D_NEG_ONE;
77  double *dL = dwork;
78  magma_int_t lddl = n;
80  magma_int_t j, jmax, jb, jp;
81
82  *info = 0;
83  if (n < 0)
84  *info = -1;
85  else if (ldda < max(1,n))
86  *info = -3;
87  else if ( lwork < n*nb )
88  *info = -6;
89
90  if (*info != 0) {
91  magma_xerbla( __func__, -(*info) );
92  return *info;
93  }
94
95  /* Quick return if possible */
96  if ( n == 0 )
97  return *info;
98
99  /* Invert the triangular factor U */
100  magma_dtrtri_gpu( MagmaUpper, MagmaNonUnit, n, dA, ldda, info );
101  if ( *info != 0 )
102  return *info;
103
104  jmax = ((n-1) / nb)*nb;
105  for( j = jmax; j >= 0; j -= nb ) {
106  jb = min( nb, n-j );
107
108  // copy current block column of L to work space,
109  // then replace with zeros in A.
111  dA(j,j), ldda,
112  dL(j,0), lddl );
113  magmablas_dlaset( MagmaLower, n-j, jb, dA(j,j), ldda );
114
115  // compute current block column of Ainv
116  // Ainv(:, j:j+jb-1)
117  // = ( U(:, j:j+jb-1) - Ainv(:, j+jb:n) L(j+jb:n, j:j+jb-1) )
118  // * L(j:j+jb-1, j:j+jb-1)^{-1}
119  // where L(:, j:j+jb-1) is stored in dL.
120  if ( j+jb < n ) {
121  magma_dgemm( MagmaNoTrans, MagmaNoTrans, n, jb, n-j-jb,
122  c_neg_one, dA(0,j+jb), ldda,
123  dL(j+jb,0), lddl,
124  c_one, dA(0,j), ldda );
125  }
127  n, jb, c_one,
128  dL(j,0), lddl,
129  dA(0,j), ldda );
130  }
131
132  // Apply column interchanges
133  for( j = n-2; j >= 0; --j ) {
134  jp = ipiv[j] - 1;
135  if ( jp != j ) {
136  magmablas_dswap( n, dA(0,j), 1, dA(0,jp), 1 );
137  }
138  }
139
140  return *info;
141 }
#define MagmaUpperLower
Definition: magma.h:63
#define min(a, b)
Definition: common_magma.h:86
#define MAGMA_D_ONE
Definition: magma.h:176
magma_int_t magma_dtrtri_gpu(char uplo, char diag, magma_int_t n, double *dA, magma_int_t ldda, magma_int_t *info)
Definition: dtrtri_gpu.cpp:16
#define __func__
Definition: common_magma.h:65
#define MagmaUpper
Definition: magma.h:61
#define dL(i, j)
magma_int_t ldda
void magma_dtrsm(magma_side_t side, magma_uplo_t uplo, magma_trans_t trans, magma_diag_t diag, magma_int_t m, magma_int_t n, double alpha, magmaDouble_const_ptr dA, magma_int_t ldda, magmaDouble_ptr dB, magma_int_t lddb)
void magmablas_dlaset(magma_uplo_t uplo, magma_int_t m, magma_int_t n, magmaDouble_ptr dA, magma_int_t ldda)
int magma_int_t
Definition: magmablas.h:12
#define dA(i, j)
#define dwork(dev, i, j)
void magmablas_dlacpy(magma_uplo_t uplo, magma_int_t m, magma_int_t n, magmaDouble_const_ptr dA, magma_int_t ldda, magmaDouble_ptr dB, magma_int_t lddb)
#define MagmaLower
Definition: magma.h:62
void magma_xerbla(const char *srname, magma_int_t info)
Definition: xerbla.cpp:8
magma_int_t magma_dgetri_gpu(magma_int_t n, double *dA, magma_int_t ldda, magma_int_t *ipiv, double *dwork, magma_int_t lwork, magma_int_t *info)
Definition: dgetri_gpu.cpp:16
#define MagmaNonUnit
Definition: magma.h:65
void magma_dgemm(magma_trans_t transA, magma_trans_t transB, magma_int_t m, magma_int_t n, magma_int_t k, double alpha, magmaDouble_const_ptr dA, magma_int_t ldda, magmaDouble_const_ptr dB, magma_int_t lddb, double beta, magmaDouble_ptr dC, magma_int_t lddc)
#define MagmaRight
Definition: magma.h:69
#define MagmaUnit
Definition: magma.h:66
#define MAGMA_D_NEG_ONE
Definition: magma.h:178
#define MagmaNoTrans
Definition: magma.h:57
#define max(a, b)
Definition: common_magma.h:82
magma_int_t magma_get_dgetri_nb(magma_int_t m)
Definition: get_nb.cpp:574