/* nag_real_svd (f02wec) Example Program. * * Copyright 1990 Numerical Algorithms Group. * * Mark 1, 1990. * Mark 8 revised, 2004. */ #include #include #include #include #define EX1_MMAX 20 #define EX1_NMAX 10 static int ex1(), ex2(); int main(void) { Integer exit_status_ex1 = 0; Integer exit_status_ex2 = 0; printf("nag_real_svd (f02wec) Example Program Results\n"); scanf(" %*[^\n]"); /* Skip heading in data file */ exit_status_ex1 = ex1(); exit_status_ex2 = ex2(); return (exit_status_ex1 == 0 && exit_status_ex2 == 0) ? 0 : 1; } #define A(I, J) a[(I) *tda + J] #define B(I, J) b[(I) *tdb + J] #define PT(I, J) pt[(I) *tdpt + J] #define Q(I, J) q[(I) *tdq + J] static int ex1(void) { Nag_Boolean wantp, wantq; Integer exit_status = 0, failinfo, i, iter, j, m, n, ncolb, tda, tdb, tdpt; NagError fail; double *a = 0, *b = 0, *dummy = 0, *e = 0, *pt = 0, *sv = 0; INIT_FAIL(fail); printf("Example 1\n"); scanf(" %*[^\n]"); /* Skip Example 1 heading */ scanf(" %*[^\n]"); scanf("%ld%ld", &m, &n); if (m >= 0 && n >= 0) { ncolb = 1; if (!(a = NAG_ALLOC(m*n, double)) || !(b = NAG_ALLOC(m*ncolb, double)) || !(e = NAG_ALLOC(MIN(m, n)-1, double)) || !(pt = NAG_ALLOC(n*n, double)) || !(sv = NAG_ALLOC(MIN(m, n), double)) || !(dummy = NAG_ALLOC(1, double))) { printf("Allocation failure\n"); exit_status = -1; goto END; } tda = n; tdb = ncolb; tdpt = n; } else { printf("Invalid m or n.\n"); exit_status = 1; return exit_status; } scanf(" %*[^\n]"); for (i = 0; i < m; ++i) for (j = 0; j < n; ++j) scanf("%lf", &A(i, j)); scanf(" %*[^\n]"); for (i = 0; i < m; ++i) for (j = 0; j < ncolb; ++j) scanf("%lf", &B(i, j)); /* Find the SVD of A. */ wantq = Nag_TRUE; wantp = Nag_TRUE; /* nag_real_svd (f02wec). * SVD of real matrix */ nag_real_svd(m, n, a, tda, ncolb, b, tdb, wantq, dummy, (Integer) 1, sv, wantp, pt, tdpt, &iter, e, &failinfo, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_real_svd (f02wec).\n%s\n", fail.message); exit_status = 1; goto END; } printf("Singular value decomposition of A\n\n"); printf("Singular values\n"); for (i = 0; i < n; ++i) printf(" %8.4f", sv[i]); printf("\n\n"); printf("Left-hand singular vectors, by column\n"); for (i = 0; i < m; ++i) { for (j = 0; j < n; ++j) printf(" %8.4f", A(i, j)); printf("\n"); } printf("\n"); printf("Right-hand singular vectors, by column\n"); for (i = 0; i < n; ++i) { for (j = 0; j < n; ++j) printf(" %8.4f", PT(j, i)); printf("\n"); } printf("\n"); printf("Vector Q'*B\n"); for (i = 0; i < m; ++i) printf(" %8.4f", b[i]); printf("\n\n"); END: if (a) NAG_FREE(a); if (b) NAG_FREE(b); if (e) NAG_FREE(e); if (pt) NAG_FREE(pt); if (sv) NAG_FREE(sv); if (dummy) NAG_FREE(dummy); return exit_status; } static int ex2(void) { Nag_Boolean wantp, wantq; Integer exit_status = 0, failinfo, i, iter, j, m, n, ncolb, tda, tdq; NagError fail; double *a = 0, *dummy = 0, *e = 0, *q = 0, *sv = 0; INIT_FAIL(fail); printf("\nExample 2\n"); scanf(" %*[^\n]"); /* Skip Example 2 heading */ scanf(" %*[^\n]"); scanf("%ld%ld", &m, &n); if (m >= 0 && n >= 0) { if (!(a = NAG_ALLOC(m*n, double)) || !(q = NAG_ALLOC(m*m, double)) || !(e = NAG_ALLOC(MIN(m, n)-1, double)) || !(sv = NAG_ALLOC(MIN(m, n), double)) || !(dummy = NAG_ALLOC(1, double))) { printf("Allocation failure\n"); exit_status = -1; goto END; } tda = n; tdq = m; } else { } scanf(" %*[^\n]"); for (i = 0; i < m; ++i) for (j = 0; j < n; ++j) scanf("%lf", &A(i, j)); /* Find the SVD of A. */ wantq = Nag_TRUE; wantp = Nag_TRUE; ncolb = 0; /* nag_real_svd (f02wec), see above. */ nag_real_svd(m, n, a, tda, ncolb, dummy, (Integer) 1, wantq, q, tdq, sv, wantp, dummy, (Integer) 1, &iter, e, &failinfo, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_real_svd (f02wec).\n%s\n", fail.message); exit_status = 1; goto END; } printf("Singular value decomposition of A\n\n\n"); printf("Singular values\n\n"); for (i = 0; i < m; ++i) printf(" %8.4f", sv[i]); printf("\n\n"); printf("Left-hand singular vectors, by column\n\n"); for (i = 0; i < m; ++i) { for (j = 0; j < m; ++j) printf(" %8.4f", Q(i, j)); printf("\n"); } printf("Right-hand singular vectors, by column\n\n"); for (i = 0; i < n; ++i) { for (j = 0; j < m; ++j) printf(" %8.4f", A(j, i)); printf("\n"); } END: if (a) NAG_FREE(a); if (q) NAG_FREE(q); if (e) NAG_FREE(e); if (sv) NAG_FREE(sv); if (dummy) NAG_FREE(dummy); return exit_status; }