/* 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(void), ex2(void); int main(void) { Integer exit_status_ex1=0; Integer exit_status_ex2=0; Vprintf("nag_real_svd (f02wec) Example Program Results\n"); Vscanf(" %*[^\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); Vprintf("Example 1\n"); Vscanf(" %*[^\n]"); /* Skip Example 1 heading */ Vscanf(" %*[^\n]"); Vscanf("%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)) ) { Vprintf("Allocation failure\n"); exit_status = -1; goto END; } tda = n; tdb = ncolb; tdpt = n; } else { Vprintf("Invalid m or n.\n"); exit_status = 1; return exit_status; } Vscanf(" %*[^\n]"); for (i = 0; i < m; ++i) for (j = 0; j < n; ++j) Vscanf("%lf", &A(i,j)); Vscanf(" %*[^\n]"); for (i = 0; i < m; ++i) for (j = 0; j < ncolb; ++j) Vscanf("%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) { Vprintf("Error from nag_real_svd (f02wec).\n%s\n", fail.message); exit_status = 1; goto END; } Vprintf("Singular value decomposition of A\n\n"); Vprintf("Singular values\n"); for (i = 0; i < n; ++i) Vprintf(" %8.4f", sv[i]); Vprintf("\n\n"); Vprintf("Left-hand singular vectors, by column\n"); for (i = 0; i < m; ++i) { for (j = 0; j < n; ++j) Vprintf(" %8.4f", A(i,j)); Vprintf("\n"); } Vprintf("\n"); Vprintf("Right-hand singular vectors, by column\n"); for (i = 0; i < n; ++i) { for (j = 0; j < n; ++j) Vprintf(" %8.4f", PT(j,i)); Vprintf("\n"); } Vprintf("\n"); Vprintf("Vector Q'*B\n"); for (i = 0; i < m; ++i) Vprintf(" %8.4f", b[i]); Vprintf("\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); Vprintf("\nExample 2\n"); Vscanf(" %*[^\n]"); /* Skip Example 2 heading */ Vscanf(" %*[^\n]"); Vscanf("%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)) ) { Vprintf("Allocation failure\n"); exit_status = -1; goto END; } tda = n; tdq = m; } else { } Vscanf(" %*[^\n]"); for (i = 0; i < m; ++i) for (j = 0; j < n; ++j) Vscanf("%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) { Vprintf("Error from nag_real_svd (f02wec).\n%s\n", fail.message); exit_status = 1; goto END; } Vprintf("Singular value decomposition of A\n\n\n"); Vprintf("Singular values\n\n"); for (i = 0; i < m; ++i) Vprintf(" %8.4f", sv[i]); Vprintf("\n\n"); Vprintf("Left-hand singular vectors, by column\n\n"); for (i = 0; i < m; ++i) { for (j = 0; j < m; ++j) Vprintf(" %8.4f", Q(i,j)); Vprintf("\n"); } Vprintf("Right-hand singular vectors, by column\n\n"); for (i = 0; i < n; ++i) { for (j = 0; j < m; ++j) Vprintf(" %8.4f", A(j,i)); Vprintf("\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; }