MAGMA  1.2.0
MatrixAlgebraonGPUandMulticoreArchitectures
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups
dgetri_gpu.cpp
Go to the documentation of this file.
1 /*
2  -- MAGMA (version 1.2.0) --
3  Univ. of Tennessee, Knoxville
4  Univ. of California, Berkeley
5  Univ. of Colorado, Denver
6  May 2012
7 
8  @generated d Thu May 10 22:26:48 2012
9 
10 */
11 #include "common_magma.h"
12 
13 // === Define what BLAS to use ============================================
14 #define PRECISION_d
15 #if (defined(PRECISION_s) || defined(PRECISION_d))
16  #define magma_dgemm magmablas_dgemm
17  #define magma_dtrsm magmablas_dtrsm
18 #endif
19 
20 #if (GPUSHMEM >= 200)
21 #if (defined(PRECISION_s))
22  #undef magma_sgemm
23  #define magma_sgemm magmablas_sgemm_fermi80
24 #endif
25 #endif
26 // === End defining what BLAS to use ======================================
27 
28 extern "C" magma_int_t
29 magma_dgetri_gpu( magma_int_t n, double *dA, magma_int_t lda,
30  magma_int_t *ipiv, double *dwork, magma_int_t lwork,
31  magma_int_t *info )
32 {
33 /* -- MAGMA (version 1.2.0) --
34  Univ. of Tennessee, Knoxville
35  Univ. of California, Berkeley
36  Univ. of Colorado, Denver
37  May 2012
38 
39  Purpose
40  =======
41 
42  DGETRI computes the inverse of a matrix using the LU factorization
43  computed by DGETRF. This method inverts U and then computes inv(A) by
44  solving the system inv(A)*L = inv(U) for inv(A).
45 
46  Note that it is generally both faster and more accurate to use DGESV,
47  or DGETRF and DGETRS, to solve the system AX = B, rather than inverting
48  the matrix and multiplying to form X = inv(A)*B. Only in special
49  instances should an explicit inverse be computed with this routine.
50 
51  Arguments
52  =========
53 
54  N (input) INTEGER
55  The order of the matrix A. N >= 0.
56 
57  dA (input/output) DOUBLE_PRECISION array on the GPU, dimension (LDA,N)
58  On entry, the factors L and U from the factorization
59  A = P*L*U as computed by DGETRF_GPU.
60  On exit, if INFO = 0, the inverse of the original matrix A.
61 
62  LDA (input) INTEGER
63  The leading dimension of the array A. LDA >= max(1,N).
64 
65  IPIV (input) INTEGER array, dimension (N)
66  The pivot indices from DGETRF; for 1<=i<=N, row i of the
67  matrix was interchanged with row IPIV(i).
68 
69  DWORK (workspace/output) COMPLEX*16 array on the GPU, dimension (MAX(1,LWORK))
70 
71  LWORK (input) INTEGER
72  The dimension of the array DWORK. LWORK >= N*NB, where NB is
73  the optimal blocksize returned by magma_get_dgetri_nb(n).
74 
75  Unlike LAPACK, this version does not currently support a
76  workspace query, because the workspace is on the GPU.
77 
78  INFO (output) INTEGER
79  = 0: successful exit
80  < 0: if INFO = -i, the i-th argument had an illegal value
81  > 0: if INFO = i, U(i,i) is exactly zero; the matrix is
82  singular and its cannot be computed.
83 
84  ===================================================================== */
85 
86  /* Local variables */
87  double c_one = MAGMA_D_ONE;
88  double c_neg_one = MAGMA_D_NEG_ONE;
89  double *dL = dwork;
90  magma_int_t ldl = n;
91  magma_int_t nb = magma_get_dgetri_nb(n);
92  magma_int_t j, jmax, jb, jp;
93 
94  *info = 0;
95  if (n < 0)
96  *info = -1;
97  else if (lda < max(1,n))
98  *info = -3;
99  else if ( lwork < n*nb )
100  *info = -6;
101 
102  if (*info != 0) {
103  magma_xerbla( __func__, -(*info) );
104  return *info;
105  }
106 
107  /* Quick return if possible */
108  if ( n == 0 )
109  return *info;
110 
111  /* Invert the triangular factor U */
112  magma_dtrtri_gpu( MagmaUpper, MagmaNonUnit, n, dA, lda, info );
113  if ( *info != 0 )
114  return *info;
115 
116  jmax = ((n-1) / nb)*nb;
117  for( j = jmax; j >= 0; j -= nb ) {
118  jb = min( nb, n-j );
119 
120  // copy current block column of L to work space,
121  // then replace with zeros in A.
122  magmablas_dlacpy( MagmaUpperLower, n-j, jb,
123  &dA[j + j*lda], lda,
124  &dL[j ], ldl );
125  magmablas_dlaset( MagmaLower, n-j, jb, &dA[j + j*lda], lda );
126 
127  // compute current block column of Ainv
128  // Ainv(:, j:j+jb-1)
129  // = ( U(:, j:j+jb-1) - Ainv(:, j+jb:n) L(j+jb:n, j:j+jb-1) )
130  // * L(j:j+jb-1, j:j+jb-1)^{-1}
131  // where L(:, j:j+jb-1) is stored in dL.
132  if ( j+jb < n ) {
133  magma_dgemm( MagmaNoTrans, MagmaNoTrans, n, jb, n-j-jb,
134  c_neg_one, &dA[(j+jb)*lda], lda,
135  &dL[ j+jb ], ldl,
136  c_one, &dA[ j*lda], lda );
137  }
138  magma_dtrsm( MagmaRight, MagmaLower, MagmaNoTrans, MagmaUnit,
139  n, jb, c_one,
140  &dL[j ], ldl,
141  &dA[j*lda], lda );
142  }
143 
144  // Apply column interchanges
145  for( j = n-2; j >= 0; --j ) {
146  jp = ipiv[j] - 1;
147  if ( jp != j ) {
148  magmablas_dswap( n, &dA[ j*lda ], 1, &dA[ jp*lda ], 1 );
149  }
150  }
151 
152  return *info;
153 }