NAG Library Manual, Mark 28.6
Interfaces:  FL   CL   CPP   AD 

NAG CL Interface Introduction
Example description
/* nag_correg_corrmat_nearest (g02aac) Example Program.
 *
 * Copyright 2022 Numerical Algorithms Group.
 *
 * Mark 28.6, 2022.
 */
/* Pre-processor includes */
#include <math.h>
#include <nag.h>
#include <stdio.h>
#include <string.h>

#ifdef __cplusplus
extern "C"
{
#endif
  static double NAG_CALL bound(Nag_OrderType order,
                               Nag_UploType  uplo,
                               Integer       n,
                               const double  g[],
                               Integer       pdg,
                               double        a[],
                               Integer       pda,
                               double        offdiag[],
                               Integer       ipiv[],
                               NagError *    fail);
#ifdef __cplusplus
}
#endif

int main(void)
{
  /* Integer scalar and array declarations */
  Integer  exit_status = 0;
  Integer  feval, i, iter, j, maxit, maxits, n;
  Integer  pda, pdg, pdx;
  Integer *iwork = 0;

  /* Double scalar and array declarations */
  double        errtol, dist, nrmgrd;
  double *      a = 0, *g = 0, *work = 0, *x = 0;
  Nag_NormType  norm;
  Nag_OrderType order;
  Nag_UploType  uplo;
  NagError      fail;

  INIT_FAIL(fail);

#ifdef NAG_COLUMN_MAJOR
#define A(I, J) a[(J - 1) * pda + I - 1]
#define G(I, J) g[(J - 1) * pdg + I - 1]
#define X(I, J) x[(J - 1) * pdx + I - 1]
  order = Nag_ColMajor;
#else
#define A(I, J) a[(I - 1) * pda + J - 1]
#define G(I, J) g[(I - 1) * pdg + J - 1]
#define X(I, J) x[(I - 1) * pdx + J - 1]
  order = Nag_RowMajor;
#endif

  printf("nag_correg_corrmat_nearest (g02aac) Example Program Results\n\n");

  /* Skip heading in data file */
  scanf("%*[^\n] ");

  /* Read in the problem size */
  scanf("%" NAG_IFMT "%*[^\n] ", &n);

  pda = n;
  pdg = n;
  pdx = n;
  if (!(a = NAG_ALLOC(pda * n, double)) || !(g = NAG_ALLOC(pdg * n, double)) ||
      !(iwork = NAG_ALLOC(n, Integer)) || !(work = NAG_ALLOC(n, double)) ||
      !(x = NAG_ALLOC(pdx * n, double)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

  /* Read in the matrix g */
  for (i = 1; i <= n; i++)
    for (j = 1; j <= n; j++)
      scanf("%lf", &G(i, j));
  scanf("%*[^\n] ");

  uplo = Nag_Lower;

  /*
   * Copy G into A to compute a bound on the distance to the NCM
   * We assume G is symmetric and copy only the lower triangle
   */
  for (j = 1; j <= n; j++)
    {
      for (i = j; i <= n; i++)
        A(i, j) = G(i, j);
    }

  dist = bound(order, uplo, n, g, pdg, a, pda, work, iwork, &fail);
  if (fail.code != NE_NOERROR)
    {
      printf("Error computing the upper bound on distance to the NCM.\n%s\n",
             fail.message);
      exit_status = 1;
      goto END;
    }
  printf(" Upper bound on the distance to the NCM: %12.3f\n\n", dist);

  /*
   * A was overwritten by bound. In order to compute the actual
   * distance to the nearest correlation matrix, store another copy of G in A
   */
  for (j = 1; j <= n; j++)
    {
      for (i = j; i <= n; i++)
        A(i, j) = G(i, j);
    }

  /* Set up method parameters */
  errtol = 1.00e-7;
  maxits = 200;
  maxit  = 10;

  /*
   * nag_correg_corrmat_nearest (g02aac)
   * Computes the nearest correlation matrix to a real square matrix,
   * using the method of Qi and Sun
   */
  nag_correg_corrmat_nearest(order, g, pdg, n, errtol, maxits, maxit, x, pdx,
                             &iter, &feval, &nrmgrd, &fail);
  if (fail.code != NE_NOERROR)
    {
      printf("Error from nag_correg_corrmat_nearest (g02aac).\n%s\n",
             fail.message);
      exit_status = 2;
      goto END;
    }
  printf(" Nearest Correlation Matrix\n");
  for (i = 1; i <= n; i++)
    {
      for (j = 1; j <= n; j++)
        printf("%11.5f%s", X(i, j), j % 4 ? " " : "\n");
    }

  /*
   * Compute distance between the original input matrix and the nearest
   * correlation matrix
   */
  for (j = 1; j <= n; j++)
    {
      for (i = j; i <= n; i++)
        A(i, j) -= X(i, j);
    }

  norm = Nag_FrobeniusNorm;
  /* nag_dsy_norm (f16rcc).
   * 1-norm, infinity-norm, Frobenius norm, largest absolute
   * element, real symmetric matrix
   */
  f16rcc(order, norm, uplo, n, a, pda, &dist, &fail);
  if (fail.code != NE_NOERROR)
    {
      printf("Error from nag_correg_corrmat_nearest (g02aac).\n%s\n",
             fail.message);
      exit_status = 3;
      goto END;
    }

  printf("\n");
  printf("\n");
  printf(" Distance between G and X:               %12.3f\n", dist);
  printf(" Number of Newton steps taken:%19" NAG_IFMT "\n", iter);
  printf(" Number of function evaluations:%17" NAG_IFMT "\n", feval);
  if (nrmgrd > errtol)
    printf(" Norm of gradient of last Newton step: %12.3f\n", nrmgrd);
  printf("\n");

END:
  NAG_FREE(a);
  NAG_FREE(g);
  NAG_FREE(iwork);
  NAG_FREE(work);
  NAG_FREE(x);

  return exit_status;
}

static double NAG_CALL bound(Nag_OrderType order,
                             Nag_UploType  uplo,
                             Integer       n,
                             const double  g[],
                             Integer       pdg,
                             double        a[],
                             Integer       pda,
                             double        offdiag[],
                             Integer       ipiv[],
                             NagError *    fail)
{
  /* Scalars */
  Integer      i, j;
  double       delta, ubound;
  Nag_NormType norm;
  Nag_UploType uplo_c;

  delta  = 1.0E-7;
  ubound = 0.0;
  uplo_c = uplo;

  /* Modified cholesky routines expect column major storage */
  if (order == Nag_RowMajor)
    uplo_c = (uplo_c == Nag_Lower ? Nag_Upper : Nag_Lower);

  /* nag_matop_real_modified_cholesky (f01mdc)
   * Compute the modified Cholesky factorization of A
   */
  f01mdc(uplo_c, n, a, pda, offdiag, ipiv, delta, fail);
  if (fail->code != NE_NOERROR)
    return 0.0;

  /* Compute the perturbed matrix A + E
   * nag_matop_real_mod_chol_perturbed_a (f01mec)
   * Compute the modified Cholesky factorization of A
   */
  f01mec(uplo_c, n, a, pda, offdiag, ipiv, fail);
  if (fail->code != NE_NOERROR)
    return 0.0;

  /* Overwrite A with D^(-1/2)(G + E)D^(-1/2) - G
   * where D = diag(G + E)
   */
  for (j = 1; j <= n; j++)
    {
      for (i = j + 1; i <= n; i++)
        {
          A(i, j) = A(i, j) / sqrt(A(i, i) * A(j, j)) - G(i, j);
        }
      A(j, j) = 1.0E0 - G(j, j);
    }

  /* Bound is given by ||D^(-1/2)(G + E)D^(-1/2) - G|| */
  norm = Nag_FrobeniusNorm;
  /* nag_dsy_norm (f16rcc).
   * 1-norm, infinity-norm, Frobenius norm, largest absolute
   * element, real symmetric matrix
   */
  f16rcc(order, norm, uplo, n, a, pda, &ubound, fail);
  if (fail->code != NE_NOERROR)
    return 0.0;

  return ubound;
}