/* nag_mv_procustes (g03bcc) Example Program. * * Copyright 1998 Numerical Algorithms Group. * * Mark 5, 1998. * Mark 8 revised, 2004. * */ #include #include #include #include #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_procustes (g03bcc) Example Program Results\n\n"); /* Skip heading in data file */ scanf("%*[^\n]"); scanf("%ld", &n); scanf("%ld", &m); scanf("%s", 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("%s", 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_procustes (g03bcc). * Procrustes rotations */ nag_mv_procustes(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_procustes (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: if (r) NAG_FREE(r); if (res) NAG_FREE(res); if (x) NAG_FREE(x); if (y) NAG_FREE(y); if (yhat) NAG_FREE(yhat); return exit_status; }