/* Some examples of how to use lapacke library in C */

/* sudo apt install liblapacke-dev */

/* Compile with: gcc lapacke_examples.c -o lapacke_examples -llapacke*/

#include <stdio.h>		//for printf

// lapacke include
#include <lapacke.h>


int main() {
  int i,j;
  
  
  double complex A[2][2];
  A[0][0] = 1.0; A[0][1] = 0.0;
  A[1][0] = 1.0; A[1][1] = 3.0;
  
  printf("--- Matriz A:\n");
  for(i=0; i < 2; i++){
    for(j=0; j < 2; j++){
      printf("\t %f",creal(A[i][j]));
    }
    printf("\n");
  }
  printf("\n");
  printf("\n");
  
  /* Here we'll store the complex eigvenvalues */
  double complex complex_eigval[2];
  /* Here we'll store the complex eigvenvectors */
  double complex complex_eigvec[2];
  
  LAPACKE_zgeev(LAPACK_ROW_MAJOR, 'N', 'V', 2, (lapack_complex_double*)A, 2, (lapack_complex_double*)complex_eigval, NULL, 2, (lapack_complex_double*)complex_eigvec, 2);
  
  printf("--- Eigenvalores de A:\n");
  for(i=0; i < 2; i++){
    printf("\t %f",creal(complex_eigval[i]));
  }
  printf("\n");
  
  printf("--- Eigenvectores de A:\n");
  for(i=0; i < 2*2; i+=2){
    j=0;
    while (j < 2){
      printf("\t %f",creal(complex_eigvec[i+j]));
      j++;
    }
    printf("\n");
  }
  printf("\n");
  printf("\n");
  
  
  /* because of the ordering by ROWS (LAPACK_ROW_MAJOR), the following is necessary to get one eigenvector */
  int eig_i;
  double complex this_eigval;
  double complex this_eigvec[2];
  
  printf("--- Primer eigenvector y eigenvalor de A:\n");
  /* setting which eigenvalue/eigenvector to extract*/
  eig_i = 0;
  
  /* obtaining eigenvalue*/
  this_eigval = complex_eigval[eig_i];
  /* obtaining eigenvector*/
  for(i=eig_i,j=0; i < 2*2; i+=2,j++){
    this_eigvec[j] = complex_eigvec[i];
  }
  
  /* printing eigenvalue and eigenvector */
  printf("Eigenvalor: %f \n",creal(this_eigval));
  printf("Eigenvector: \n");
  for(j=0; j < 2; j++){
    printf("\t %f",creal(this_eigvec[j]));
    printf("\n");
  }
  
  printf("--- Segundo eigenvector y eigenvalor de A:\n");
  /* setting which eigenvalue/eigenvector to extract*/
  eig_i = 1;
  
  /* obtaining eigenvalue*/
  this_eigval = complex_eigval[eig_i];
  /* obtaining eigenvector*/
  for(i=eig_i,j=0; i < 2*2; i+=2,j++){
    this_eigvec[j] = complex_eigvec[i];
  }
  
  /* printing eigenvalue and eigenvector */
  printf("Eigenvalor: %f \n",creal(this_eigval));
  printf("Eigenvector: \n");
  for(j=0; j < 2; j++){
    printf("\t %f",creal(this_eigvec[j]));
    printf("\n");
  }
  
  printf("\n");
  printf("\n");
  
  printf("--- Inversa de A:\n");
  /* resetting the original values of A,
  since LAPACKE_zgeev modifies its input
  for some reason*/
  A[0][0] = 1.0; A[0][1] = 0.0;
  A[1][0] = 1.0; A[1][1] = 3.0;
  
  int pivot[2];
  
  double complex A_inv[2*2];
  for(i=0; i < 2; i++){
    for(j=0; j < 2; j++){
      A_inv[(2*i)+j] = A[i][j];
    }
  }
  
  LAPACKE_zgetrf(LAPACK_ROW_MAJOR, 2, 2, (lapack_complex_double*)A_inv, 2, pivot); 
  LAPACKE_zgetri(LAPACK_ROW_MAJOR, 2, (lapack_complex_double*)A_inv, 2, pivot);
  
  for(i=0; i < 2*2; i+=2){
    j=0;
    while (j < 2){
      printf("\t %f",creal(A_inv[i+j]));
      j++;
    }
    printf("\n");
  }
  printf("\n");
  
  return 0;
}