/* Copyright 2004 INAF - Osservatorio Astronomico di Cagliari

   Licensed under the Apache License, Version 2.0 (the "License");
   you may not use this file except in compliance with the License.
   You may obtain a copy of the License at

       http://www.apache.org/licenses/LICENSE-2.0

   Unless required by applicable law or agreed to in writing, software
   distributed under the License is distributed on an "AS IS" BASIS,
   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   See the License for the specific language governing permissions and
   limitations under the License.
 */

/*! \file magma_calls.cpp
 *
 * \brief Implementation of the interface with MAGMA libraries.
 */
#ifndef INCLUDE_TYPES_H_
#include "../include/types.h"
#endif

#ifdef USE_MAGMA
#ifndef INCLUDE_MAGMA_CALLS_H_
#include "../include/magma_calls.h"
#endif

void magma_zinvert(dcomplex **mat, np_int n, int &jer) {
  // magma_int_t result = magma_init();
  magma_int_t err = MAGMA_SUCCESS;
  magma_queue_t queue = NULL;
  magma_device_t dev = 0;
  magma_queue_create(dev, &queue);
  magmaDoubleComplex *dwork; // dwork - workspace
  magma_int_t ldwork; // size of dwork
  magma_int_t *piv , info; // piv - array of indices of inter -
  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
  ldwork = m * magma_get_zgetri_nb(m); // optimal block size
  // allocate matrices
  err = magma_zmalloc(&d_a, mm); // device memory for a
  err = magma_zmalloc(&dwork, ldwork); // dev. mem. for ldwork
  piv = new magma_int_t[m]; // host mem.
  magma_zsetmatrix(m, m, a, m, d_a , m, queue); // copy a -> d_a
  
  magma_zgetrf_gpu(m, m, d_a, m, piv, &info);
  magma_zgetri_gpu(m, d_a, m, piv, dwork, ldwork, &info);
  
  magma_zgetmatrix( m, m, d_a , m, a, m, queue); // copy d_a -> a
  delete[] piv; // free host memory
  magma_free(d_a); // free device memory
  magma_queue_destroy(queue); // destroy queue
  // result = magma_finalize();
  jer = (int)err;
}
#endif
