djacobi_rci_c.c

/*
********************************************************************************
*   Copyright(C) 2004-2013 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.
*   
*   *Third Party trademarks are the property of their respective owners.
*   
*   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 programm computes the Jacobi matrix of the function on the basis of RCI
  using the central difference.
*/

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

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

int
main ()
{
  /* user's objective function */
  extern void extendet_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 *) malloc (sizeof (double) * n * m)) == NULL)
    {
      printf ("\n#fail: jacobi matrix allocation failed\n");
    }
  if ((x = (double *) malloc (sizeof (double) * n)) == NULL)
    {
      printf ("\n#fail: vector X allocation failed\n");
      free (a);
      return 1;
    }
  if ((f1 = (double *) malloc (sizeof (double) * n)) == NULL)
    {
      printf ("\n#fail: vector F1 allocation failed\n");
      free (x);
      free (a);
      return 1;
    }
  if ((f2 = (double *) malloc (sizeof (double) * n)) == NULL)
    {
      printf ("\n#fail: vector F2 allocation failed\n");
      free (f1);
      free (x);
      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) */
          extendet_powell (&m, &n, x, f1);
        }
      else if (rci_request == 2)
        {
          /* calculate function value f1 = f(x-eps) */
          extendet_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 ();
  free (f2);
  free (f1);
  free (x);
  free (a);
  printf ("#pass\n");
  fflush (0);
  return 0;
}

/* 
    routine for extendet 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
extendet_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;
}
For more complete information about compiler optimizations, see our Optimization Notice.