/* nag_sum_sqs_combine (g02bzc) Example Program. * * Copyright 2013 Numerical Algorithms Group. * * Mark 24, 2013. */ #include #include #include #include #include #define X(I,J) x[(order == Nag_ColMajor) ? (J)*pdx + (I) : (I)*pdx + (J)] int main(void) { /* Integer scalar and array declarations */ Integer b, i, j, ierr, lc, pdx, m, n, iwt; Integer exit_status = 0; /* NAG structures and types */ NagError fail; Nag_SumSquare mean; Nag_OrderType order = Nag_ColMajor; /* Double scalar and array declarations */ double alpha, xsw, ysw; double *wt = 0, *x = 0, *xc = 0, *xmean = 0, *yc = 0, *ymean = 0; /* Character scalar and array declarations */ char cmean[40]; /* Initialise the error structure */ INIT_FAIL(fail); printf("nag_sum_sqs_combine (g02bzc) Example Program Results\n\n"); /* Skip heading in data file */ scanf("%*[^\n] "); /* Read in the problem defining variables */ scanf("%39s%ld%*[^\n] ",cmean,&m); mean = (Nag_SumSquare) nag_enum_name_to_value(cmean); /* Allocate memory for output arrays */ lc = (m*m+m)/2; if (!(xmean = NAG_ALLOC(m, double)) || !(ymean = NAG_ALLOC(m, double)) || !(xc = NAG_ALLOC(lc, double)) || !(yc = NAG_ALLOC(lc, double))) { printf("Allocation failure\n"); exit_status = -1; goto END; } /* Loop over each block of data */ for (b = 0;;) { /* Read in the number of observations in this block and a flag indicating * whether weights have been supplied (iwt = 1) or not (iwt = 0). */ ierr = scanf("%ld%ld",&n,&iwt); if (ierr == EOF || ierr < 2) break; scanf("%*[^\n] "); /* Keep a running total of the number of blocks of data */ b++; /* Reallocate X to the required size */ NAG_FREE(x); pdx = n; if (!(x = NAG_ALLOC(pdx*m, double))) { printf("Allocation failure\n"); exit_status = -1; goto END; } /* Read in the data for this block */ if (iwt) { /* Weights supplied, so reallocate X to the required size */ NAG_FREE(wt); if (!(wt = NAG_ALLOC(n, double))) { printf("Allocation failure\n"); exit_status = -1; goto END; } for (i = 0; i < n; i++) { for (j = 0; j < m; j++) scanf("%lf",&X(i,j)); scanf("%lf",&wt[i]); } } else { /* No weights */ NAG_FREE(wt); wt = 0; for (i = 0; i < n; i++) for (j = 0; j < m; j++) scanf("%lf",&X(i,j)); } scanf("%*[^\n] "); /* Call nag_sum_sqs (g02buc) to summarise this block of data */ if (b == 1) { /* This is the first block of data, so summarise the data into * xmean and xc. */ nag_sum_sqs(order,mean,n,m,x,pdx,wt,&xsw,xmean,xc,&fail); if (fail.code != NE_NOERROR) { printf("Error from nag_sum_sqs (g02buc).\n%s\n", fail.message); exit_status = 1; goto END; } } else { /* This is not the first block of data, so summarise the data into * ymean and yc. */ nag_sum_sqs(order,mean,n,m,x,pdx,wt,&ysw,ymean,yc,&fail); if (fail.code != NE_NOERROR) { printf("Error from nag_sum_sqs (g02buc).\n%s\n", fail.message); exit_status = 1; goto END; } /* Call nag_sum_sqs_combine (g02bzc) to update the running summaries */ nag_sum_sqs_combine(mean,m,&xsw,xmean,xc,ysw,ymean,yc,&fail); if (fail.code != NE_NOERROR) { printf("Error from nag_sum_sqs_combine (g02bzc).\n%s\n", fail.message); exit_status = 1; goto END; } } } /* Display results */ printf(" Means\n "); for (i = 0; i < m; i++) printf("%14.4f",xmean[i]); printf("\n\n"); /* Call nag_pack_real_mat_print (x04ccc) to print the sums of squares */ nag_pack_real_mat_print(Nag_ColMajor,Nag_Upper,Nag_NonUnitDiag, m, xc, "Sums of squares and cross-products", NULL, &fail); if (xsw>1.0 && mean==Nag_AboutMean && fail.code == NE_NOERROR) { /* Convert the sums of squares and cross-products to a covariance matrix */ alpha = 1.0/(xsw-1.0); for (i = 0; i < lc; i++) xc[i] *= alpha; printf("\n"); nag_pack_real_mat_print(Nag_ColMajor,Nag_Upper,Nag_NonUnitDiag, m, xc, "Covariance matrix", NULL, &fail); } if (fail.code != NE_NOERROR) { printf("Error from nag_pack_real_mat_print (x04ccc).\n%s\n", fail.message); exit_status = 1; goto END; } END: NAG_FREE(x); NAG_FREE(wt); NAG_FREE(xc); NAG_FREE(xmean); NAG_FREE(yc); NAG_FREE(ymean); return(exit_status); }