/* 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;
}