/* nag_mv_rot_procrustes (g03bcc) Example Program.
*
* Copyright 2019 Numerical Algorithms Group.
*
* Mark 27.0, 2019.
*
*/
#include <nag.h>
#include <stdio.h>
#define R(I, J) r[(I) *tdr + J]
#define X(I, J) x[(I) *tdx + J]
#define Y(I, J) y[(I) *tdy + J]
#define YHAT(I, J) yhat[(I) *tdy + J]
int main(void)
{
Integer exit_status = 0, i, j, m, n, tdr, tdx, tdy;
char nag_enum_arg[40];
double alpha, *r = 0, *res = 0, rss, *x = 0, *y = 0, *yhat = 0;
Nag_RotationScale scale;
Nag_TransNorm stand;
NagError fail;
INIT_FAIL(fail);
printf("nag_mv_rot_procrustes (g03bcc) Example Program Results\n\n");
/* Skip heading in data file */
scanf("%*[^\n]");
scanf("%" NAG_IFMT "", &n);
scanf("%" NAG_IFMT "", &m);
scanf("%39s", nag_enum_arg);
/* nag_enum_name_to_value (x04nac).
* Converts NAG enum member name to value
*/
stand = (Nag_TransNorm) nag_enum_name_to_value(nag_enum_arg);
scanf("%39s", nag_enum_arg);
scale = (Nag_RotationScale) nag_enum_name_to_value(nag_enum_arg);
if (m >= 1 && n >= m) {
if (!(r = NAG_ALLOC(m * m, double)) ||
!(res = NAG_ALLOC(n, double)) ||
!(x = NAG_ALLOC(n * m, double)) ||
!(y = NAG_ALLOC(n * m, double)) || !(yhat = NAG_ALLOC(n * m, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
tdr = m;
tdx = m;
tdy = m;
}
else {
printf("Invalid m or n.\n");
exit_status = 1;
return exit_status;
}
for (i = 0; i < n; ++i) {
for (j = 0; j < m; ++j)
scanf("%lf", &X(i, j));
}
for (i = 0; i < n; ++i) {
for (j = 0; j < m; ++j)
scanf("%lf", &Y(i, j));
}
/* nag_mv_rot_procrustes (g03bcc).
* Procrustes rotations
*/
nag_mv_rot_procrustes(stand, scale, n, m, x, tdx, y, tdy,
yhat, r, tdr, &alpha, &rss, res, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_mv_rot_procrustes (g03bcc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
printf("\n Rotation Matrix\n\n");
for (i = 0; i < m; ++i) {
for (j = 0; j < m; ++j)
printf(" %7.3f ", R(i, j));
printf("\n");
}
if (scale == Nag_LsqScale) {
printf("\n%s%10.3f\n", " Scale factor = ", alpha);
}
printf("\n Target Matrix \n\n");
for (i = 0; i < n; ++i) {
for (j = 0; j < m; ++j)
printf(" %7.3f ", Y(i, j));
printf("\n");
}
printf("\n Fitted Matrix\n\n");
for (i = 0; i < n; ++i) {
for (j = 0; j < m; ++j)
printf(" %7.3f ", YHAT(i, j));
printf("\n");
}
printf("\n%s%10.3f\n", "RSS = ", rss);
END:
NAG_FREE(r);
NAG_FREE(res);
NAG_FREE(x);
NAG_FREE(y);
NAG_FREE(yhat);
return exit_status;
}