/* 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; NagError fail; Nag_RotationScale scale; Nag_TransNorm stand; char char_scale[2], char_stand[2]; double alpha, *r=0, *res=0, rss, *x=0, *y=0, *yhat=0; INIT_FAIL(fail); Vprintf("nag_mv_procustes (g03bcc) Example Program Results\n\n"); /* Skip heading in data file */ Vscanf("%*[^\n]"); Vscanf("%ld",&n); Vscanf("%ld",&m); Vscanf("%s",char_stand); Vscanf("%s",char_scale); 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)) ) { Vprintf("Allocation failure\n"); exit_status = -1; goto END; } tdr = m; tdx = m; tdy = m; } else { Vprintf("Invalid m or n.\n"); exit_status = 1; return exit_status; } for (i = 0; i < n; ++i) { for (j = 0; j < m; ++j) Vscanf("%lf",&X(i,j)); } for (i = 0; i < n; ++i) { for (j = 0; j < m; ++j) Vscanf("%lf",&Y(i,j)); } if (*char_stand == 'N') { stand = Nag_NoTransNorm; } else if (*char_stand == 'Z') { stand = Nag_Orig; } else if (*char_stand == 'C') { stand = Nag_OrigCentroid; } else if (*char_stand == 'U') { stand = Nag_Norm; } else if (*char_stand == 'S') { stand = Nag_OrigNorm; } else if (*char_stand == 'M') { stand = Nag_OrigNormCentroid; } if (*char_scale == 'S') { scale = Nag_LsqScale; } else if (*char_scale == 'U') { scale = Nag_NotLsqScale; } /* 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) { Vprintf("Error from nag_mv_procustes (g03bcc).\n%s\n", fail.message); exit_status = 1; goto END; } Vprintf("\n Rotation Matrix\n\n"); for (i = 0; i < m; ++i) { for (j = 0; j < m; ++j) Vprintf(" %7.3f ",R(i,j)); Vprintf("\n"); } if (*char_scale == 'S' || *char_scale == 's') { Vprintf("\n%s%10.3f\n"," Scale factor = ",alpha); } Vprintf("\n Target Matrix \n\n"); for (i = 0; i < n; ++i) { for (j = 0; j < m; ++j) Vprintf(" %7.3f ",Y(i,j)); Vprintf("\n"); } Vprintf("\n Fitted Matrix\n\n"); for (i = 0; i < n; ++i) { for (j = 0; j < m; ++j) Vprintf(" %7.3f ",YHAT(i,j)); Vprintf("\n"); } Vprintf("\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; }