Enable possibility to send matrix inversion to different GPUs

src/include/algebraic.h

@@ -36,7 +36,8 @@
  * \param ier: `int &` Reference to an integer variable for returning a result flag.
  * \param max_size: `np_int` The maximum expected size (required by some call-backs,
  * optional, defaults to 0).
+ * \param target_device: `int` ID of target GPU, if available (defaults to 0).
  */
-void invert_matrix(dcomplex **mat, np_int size, int &ier, np_int max_size=0);
+void invert_matrix(dcomplex **mat, np_int size, int &ier, np_int max_size=0, int target_device=0);
 
 #endif

src/include/magma_calls.h

@@ -31,7 +31,8 @@
  * \param mat: Matrix of complex. The matrix to be inverted.
  * \param n: `np_int` The number of rows and columns of the [n x n] matrix.
  * \param jer: `int &` Reference to an integer return flag.
+ * \param device_id: `int` ID of the device for matrix inversion offloading.
  */
-void magma_zinvert(dcomplex **mat, np_int n, int &jer);
+void magma_zinvert(dcomplex **mat, np_int n, int &jer, int device_id=0);
 
 #endif

src/libnptm/algebraic.cpp

@@ -44,10 +44,10 @@ extern void lucin(dcomplex **mat, np_int max_size, np_int size, int &ier);
 
 using namespace std;
 
-void invert_matrix(dcomplex **mat, np_int size, int &ier, np_int max_size) {
+void invert_matrix(dcomplex **mat, np_int size, int &ier, np_int max_size, int target_device) {
   ier = 0;
 #ifdef USE_MAGMA
-  magma_zinvert(mat, size, ier);
+  magma_zinvert(mat, size, ier, target_device);
 #elif defined USE_LAPACK
   zinvert(mat, size, ier);
 #else

src/libnptm/magma_calls.cpp

@@ -27,19 +27,19 @@
 #include "../include/magma_calls.h"
 #endif
 
-void magma_zinvert(dcomplex **mat, np_int n, int &jer) {
+void magma_zinvert(dcomplex **mat, np_int n, int &jer, int device_id) {
   // magma_int_t result = magma_init();
   magma_int_t err = MAGMA_SUCCESS;
   magma_queue_t queue = NULL;
-  magma_device_t dev = 0;
+  magma_device_t dev = (magma_device_t)device_id;
   magma_queue_create(dev, &queue);
-  magmaDoubleComplex *dwork; // dwork - workspace
+  magmaDoubleComplex *dwork; // workspace
   magma_int_t ldwork; // size of dwork
-  magma_int_t *piv , info; // piv - array of indices of inter -
+  magma_int_t *piv , info; // array of pivot indices
   magma_int_t m = (magma_int_t)n; // changed rows; a - mxm matrix
-  magma_int_t mm = m * m; // size of a, r, c
-  magmaDoubleComplex *a = (magmaDoubleComplex *)&(mat[0][0]); // a - mxm matrix on the host
-  magmaDoubleComplex *d_a; // d_a - mxm matrix a on the device
+  magma_int_t mm = m * m; // size of a
+  magmaDoubleComplex *a = (magmaDoubleComplex *)&(mat[0][0]); // pointer to first element on host
+  magmaDoubleComplex *d_a; // pointer to first element on device
   ldwork = m * magma_get_zgetri_nb(m); // optimal block size
   // allocate matrices
   err = magma_zmalloc(&d_a, mm); // device memory for a