/* nag_multi_normal_pdf_vector (g01lbc) Example Program. * * Copyright 2011, Numerical Algorithms Group. * * Mark 23, 2011. */ #include #include #include #include #include #define X(I, J) x[(J) *pdx + I] #define SIG(I, J) sig[(J) *pdsig + I] int main(void) { /* Integer scalar and array declarations */ Integer i, j, k, n, pdx, pdsig, rank; Integer exit_status = 0; /* NAG structures and types */ NagError fail; Nag_MatrixType iuld; Nag_Boolean ilog; /* Double scalar and array declarations */ double *x = 0, *xmu = 0, *sig = 0, *pdf = 0; /* Character scalar and array declarations */ char ciuld[40], cilog[40]; /* Initialise the error structure */ INIT_FAIL(fail); printf("nag_multi_normal_pdf_vector (g01lbc) Example Program Results\n\n"); /* Skip heading in data file */ scanf("%*[^\n] "); /* Read in the problem size */ scanf("%ld%ld%39s%39s%*[^\n] ",&k,&n,ciuld,cilog); ilog = (Nag_Boolean) nag_enum_name_to_value(cilog); iuld = (Nag_MatrixType) nag_enum_name_to_value(ciuld); pdx = n; pdsig = (iuld==Nag_DiagonalMatrix) ? 1 : n; if (!(x = NAG_ALLOC(pdx*k, double)) || !(xmu = NAG_ALLOC(n,double)) || !(sig = NAG_ALLOC(pdsig*n, double)) || !(pdf = NAG_ALLOC(k,double))) { printf("Allocation failure\n"); exit_status = -1; goto END; } /* Read in and echo the vector of means */ for (i = 0; i < n; i++) scanf("%lf", &xmu[i]); scanf("%*[^\n] "); printf("Vector of Means:\n"); for (i = 0; i < n; i++) printf("%8.4f ", xmu[i]); printf("\n\n"); /* Read in and echo the covariance matrix */ if (iuld == Nag_DiagonalMatrix) { for (i = 0; i < n; i++) scanf("%lf", &SIG(1,i)); scanf("%*[^\n] "); printf("Diagonal Elements of Covariance Matrix:\n"); for (i = 0; i < n; i++) printf("%8.4f ", SIG(1,i)); printf("\n\n"); } else if (iuld == Nag_LowerMatrix || iuld == Nag_LowerFactored) { for (i = 0; i < n; i++) { for (j = 0; j <= i; j++) scanf("%lf", &SIG(i,j)); scanf("%*[^\n] "); } /* Print details of Covariance matrix using * nag_gen_real_mat_print (x04cac) */ if (iuld == Nag_LowerMatrix) nag_gen_real_mat_print(Nag_ColMajor, Nag_LowerMatrix, Nag_NonUnitDiag, n, n, sig, pdsig, "Covariance Matrix:", NULL, &fail); else nag_gen_real_mat_print(Nag_ColMajor, Nag_LowerMatrix, Nag_NonUnitDiag, n, n, sig, pdsig, "Lower Triangular Cholesky Factor " "of Covariance Matrix:", NULL, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_gen_real_mat_print (x04cac).\n%s\n", fail.message); exit_status = 1; goto END; } } else { /* Upper triangle of matrix or factor is stored */ for (i = 0; i < n; i++) { for (j = i; j < n; j++) scanf("%lf", &SIG(i,j)); scanf("%*[^\n] "); } /* Print details of Covariance matrix using * nag_gen_real_mat_print (x04cac) */ if (iuld == Nag_UpperMatrix) nag_gen_real_mat_print(Nag_ColMajor, Nag_UpperMatrix, Nag_NonUnitDiag, n, n, sig, pdsig, "Covariance Matrix:", NULL, &fail); else nag_gen_real_mat_print(Nag_ColMajor, Nag_UpperMatrix, Nag_NonUnitDiag, n, n, sig, pdsig, "Upper Triangular Cholesky Factor " "of Covariance Matrix:", NULL, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_gen_real_mat_print (x04cac).\n%s\n", fail.message); exit_status = 1; goto END; } } /* Read in the points at which to evaluate the PDF */ for (j = 0; j < k; j++) for (i = 0; i < n; i++) scanf("%lf", &X(i,j)); /* nag_multi_normal_pdf_vector (g01lbc): evaluate the PDF */ nag_multi_normal_pdf_vector(ilog, k, n, x, pdx, xmu, iuld, sig, pdsig, pdf, &rank, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_multi_normal_pdf_vector (g01lbc).\n%s\n", fail.message); exit_status = 1; goto END; } /* Display results */ printf("\nRank of the covariance matrix: %ld\n\n",rank); if (ilog) printf(" log(PDF) X\n"); else printf(" PDF X\n"); printf(" ------------------------------------------------\n"); for (j = 0; j < k; j++) { printf(" %13.4e", pdf[j]); for (i = 0; i < n; i++) printf(" %8.4f", X(i,j)); printf("\n"); } END: NAG_FREE(x); NAG_FREE(xmu); NAG_FREE(sig); NAG_FREE(pdf); return exit_status; }