/*******************************************************************************
* Copyright 2004-2018 Intel Corporation All Rights Reserved.
*
* The source code,  information  and material  ("Material") contained  herein is
* owned by Intel Corporation or its  suppliers or licensors,  and  title to such
* Material remains with Intel  Corporation or its  suppliers or  licensors.  The
* Material  contains  proprietary  information  of  Intel or  its suppliers  and
* licensors.  The Material is protected by  worldwide copyright  laws and treaty
* provisions.  No part  of  the  Material   may  be  used,  copied,  reproduced,
* modified, published,  uploaded, posted, transmitted,  distributed or disclosed
* in any way without Intel's prior express written permission.  No license under
* any patent,  copyright or other  intellectual property rights  in the Material
* is granted to  or  conferred  upon  you,  either   expressly,  by implication,
* inducement,  estoppel  or  otherwise.  Any  license   under such  intellectual
* property rights must be express and approved by Intel in writing.
*
* Unless otherwise agreed by Intel in writing,  you may not remove or alter this
* notice or  any  other  notice   embedded  in  Materials  by  Intel  or Intel's
* suppliers or licensors in any way.
*******************************************************************************/

/*
*   Content : DJACOBI RCI example
*
******************************************************************************** 
*/
/* 
  The program computes the Jacobi matrix of the function on the basis of RCI
  using the central difference.
*/

#include <stdlib.h>
#include <stdio.h>
#include "mkl_service.h"

#include "mkl_rci.h"
#include "mkl_types.h"

int main ()
{
  /* user’s objective function */
  extern void extended_powell (MKL_INT *, MKL_INT *, double *, double *);
  /* n - number of function variables
     m - dimension of function value */
  MKL_INT n = 4, m = 4;
  /* jacobi matrix */
  double *a;
  /* solution vector. contains values x for f(x)
     temporary arrays f1 & f2 that contain f1 = f(x+eps) | f2 = f(x-eps) */
  double *x, *f1, *f2;
  /* precisions for jacobi_matrix calculation */
  double eps = 0.000001;
  /* jacobi-matrix solver handle */
  _JACOBIMATRIX_HANDLE_t handle;
  /* controls of rci cycle */
  MKL_INT successful, rci_request, i;

  if ((a = (double *) mkl_malloc(sizeof (double) * n * m, 64)) == NULL)
    {
      printf ("\n#fail: jacobi matrix allocation failed\n");
    }
  if ((x = (double *) mkl_malloc(sizeof (double) * n, 64)) == NULL)
    {
      printf ("\n#fail: vector X allocation failed\n");
      mkl_free (a);
      return 1;
    }
  if ((f1 = (double *) mkl_malloc(sizeof (double) * n, 64)) == NULL)
    {
      printf ("\n#fail: vector F1 allocation failed\n");
      mkl_free (x);
      mkl_free (a);
      return 1;
    }
  if ((f2 = (double *) mkl_malloc(sizeof (double) * n, 64)) == NULL)
    {
      printf ("\n#fail: vector F2 allocation failed\n");
      mkl_free (f1);
      mkl_free (x);
      mkl_free (a);
      return 1;
    }

  /* set the x values */
  for (i = 0; i < n; i++)
    x[i] = 10.0;
  /* initialize solver (allocate memory, set initial values) */
  if (djacobi_init (&handle, &n, &m, x, a, &eps) != TR_SUCCESS)
    {
      /* if function does not complete successfully then print error message */
      printf ("\n#fail: error in djacobi_init\n");
      fflush (0);
      MKL_Free_Buffers ();
      return 1;
    }
  /* set initial rci cycle variables */
  rci_request = 0;
  successful = 0;
  /* rci cycle */
  while (successful == 0)
    {
      /* call solver */
      if (djacobi_solve (&handle, f1, f2, &rci_request) != TR_SUCCESS)
        {
      /* if function does not complete successfully then print error message */
          printf ("\n#fail: error in djacobi_solve\n");
          fflush (0);
          MKL_Free_Buffers ();
          return 1;
        }
      if (rci_request == 1)
        {
          /* calculate function value f1 = f(x+eps) */
          extended_powell (&m, &n, x, f1);
        }
      else if (rci_request == 2)
        {
          /* calculate function value f1 = f(x-eps) */
          extended_powell (&m, &n, x, f2);
        }
      else if (rci_request == 0)
        /* exit rci cycle */
        successful = 1;
    } /* rci cycle */
  /* free handle memory */
  if (djacobi_delete (&handle) != TR_SUCCESS)
    {
      /* if function does not complete successfully then print error message */
      printf ("\n#fail: error in djacobi_delete\n");
      fflush (0);
      MKL_Free_Buffers ();
      return 1;
    }
  MKL_Free_Buffers ();
  mkl_free (f2);
  mkl_free (f1);
  mkl_free (x);
  mkl_free (a);
  printf ("#pass\n");
  fflush (0);
  return 0;
}

/* 
    routine for extended Powell function calculation
    m in : dimension of function value
    n in : number of function variables
    x in : vector for function calculation
    f out: function value f(x)
*/
void extended_powell (MKL_INT * m, MKL_INT * n, double *x, double *f)
{
  MKL_INT i;
  for (i = 0; i < (*n) / 4; i++)
    {
      f[4 * i] = x[4 * i] + 10.0 * x[4 * i + 1];
      f[4 * i + 1] = 2.2360679774 * (x[4 * i + 2] - x[4 * i + 3]);
      f[4 * i + 2] = (x[4 * i + 1] - 2.0 * x[4 * i + 2]) * 
                     (x[4 * i + 1] - 2.0 * x[4 * i + 2]);
      f[4 * i + 3] = 3.1622776601 * (x[4 * i] - x[4 * i + 3]) * 
                                    (x[4 * i] - x[4 * i + 3]);
    }
  return;
}