Test inversion by verifying identity product

#include "../include/algebraic.h"

dcomplex* matrix_mult(

		 dcomplex *a1, const magma_int_t r1, const magma_int_t c1,

		 dcomplex *a2, const magma_int_t r2, const magma_int_t c2

) {

  dcomplex *result = NULL;

  dcomplex value;

  if (c1 == r2) {

    result = new dcomplex[r1 * c2]();

    for (magma_int_t i = 0; i < r1; i++) {

      for (magma_int_t j = 0; j < c2; j++) {

	value = 0.0 + I * 0.0;

	for (magma_int_t k = 0; k < c1; k++) {

	  value += a1[i * c1 + k] * a2[k * c2 + j];

	}

	result[i * c2 + j] = value;

      }

    }

  }

  return result;

}

void print_matrix(dcomplex *a, const magma_int_t rows, const magma_int_t columns) {

  for (magma_int_t i = 0; i < rows; i++) {

    for (magma_int_t j = 0; j < columns; j++) {

      printf("\t(%.5le,%.5le)", real(a[i * columns + j]), imag(a[i * columns + j]));

    }

    printf("\n");

  }

}

int main() {

  int result = (int)SUCCESS;

#ifdef USE_MAGMA

  const double tolerance = 1.0e-6;

  dcomplex **test_matrix;

  dcomplex *vec_matrix = new dcomplex[nn];

  dcomplex *vec_inv_matrix = new dcomplex[nn];

  for (int ivi = 0; ivi < nn; ivi++) vec_matrix[ivi] = 1.0 + ivi + I * (1.0 + ivi);

  dcomplex *old_matrix = new dcomplex[nn];

  dcomplex *prod_matrix;

  for (int ivi = 0; ivi < nn; ivi++) {

    vec_matrix[ivi] = 2.0 + ivi - I * (2 + ivi) / 2.0;

    old_matrix[ivi] = vec_matrix[ivi];

  }

  test_matrix = &vec_matrix;

  vec_inv_matrix[0] = -1.0 + I * 1.0;

  vec_inv_matrix[1] = 0.5 - I * 0.5;

  vec_inv_matrix[2] = 0.75 - I * 0.75;

  vec_inv_matrix[3] = -0.25 + I * 0.25;

#ifdef USE_MAGMA

  if (sizeof(np_int) != sizeof(magma_int_t)) {

    printf("ERROR: sizeof(np_int) = %ld; sizeof(magma_int_t) = %ld\n",

  }

#endif

  if (result == 0) {

    dcomplex difference;

    invert_matrix(test_matrix, n, result, n);

    for (int tvi = 0; tvi < 4; tvi++) {

      dcomplex difference = vec_matrix[tvi] - vec_inv_matrix[tvi];

    prod_matrix = matrix_mult(vec_matrix, n, n, old_matrix, n, n);

    for (int tvi = 0; tvi < n; tvi++) {

      for (int tvj = 0; tvj < n; tvj++) {

	difference = (tvi == tvj) ? prod_matrix[n * tvi + tvj] - 1.0 : prod_matrix[n * tvi + tvj];

	if (real(difference) < -tolerance || real(difference) > tolerance) result = 1;

	if (imag(difference) < -tolerance || imag(difference) > tolerance) result = 1;

      }

    if (result != 0) printf("ERROR: failed matrix inversion test!\n");

    }

  delete[] vec_inv_matrix;

    if (result != 0) {

      printf("ERROR: failed matrix inversion test!\n");

    }

    printf("INFO: original matrix:\n");

    print_matrix(old_matrix, n, n);

    printf("INFO: inverted matrix:\n");

    print_matrix(vec_matrix, n, n);

    printf("INFO: product matrix:\n");

    print_matrix(prod_matrix, n, n);

    delete[] prod_matrix;

  }

  delete[] vec_matrix;

  delete[] old_matrix;

#ifdef USE_MAGMA

  magma_finalize();

#endif