/* nag_mv_rot_procrustes (g03bcc) Example Program.
*
* Copyright 2024 Numerical Algorithms Group.
*
* Mark 30.1, 2024.
*
*/
#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;
}