/* nag_real_svd(f02wec) Example Program * * Copyright 1990 Numerical Algorithms Group. * * Mark 1, 1990. */ #include #include #include #include #define EX1_MMAX 20 #define EX1_NMAX 10 #define EX2_MMAX 10 #define EX2_NMAX 20 static void ex1(void), ex2(void); int main(void) { Vprintf("f02wec Example Program Results\n"); Vscanf(" %*[^\n]"); /* Skip heading in data file */ ex1(); ex2(); return EXIT_SUCCESS; } static void ex1(void) { Integer tda = EX1_NMAX; Integer tdpt = EX1_NMAX; double a[EX1_MMAX][EX1_NMAX], b[EX1_MMAX], e[EX1_NMAX-1]; double pt[EX1_NMAX][EX1_NMAX], sv[EX1_NMAX], dummy[1]; Integer i, j, m, n, iter, failinfo; Boolean wantp, wantq; static NagError fail; Vprintf("Example 1\n"); Vscanf(" %*[^\n]"); /* Skip Example 1 heading */ Vscanf(" %*[^\n]"); Vscanf("%ld%ld", &m, &n); if (m > EX1_MMAX || n > EX1_NMAX) { Vprintf("m or n is out of range.\n"); Vprintf("m = %2ld, n = %2ld\n", m, n); } else { 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) Vscanf("%lf", &b[i]); /* Find the SVD of A. */ wantq = TRUE; wantp = TRUE; fail.print = TRUE; f02wec(m, n, &a[0][0], tda, (Integer)1, b, (Integer)1, wantq, dummy, (Integer)1, sv, wantp, &pt[0][0], tdpt, &iter, e, &failinfo, &fail); if (fail.code != NE_NOERROR) exit(EXIT_FAILURE); 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"); } } static void ex2(void) { Integer tda = EX2_NMAX; Integer tdq = EX2_MMAX; double a[EX2_MMAX][EX2_NMAX], e[EX2_NMAX-1]; double q[EX2_MMAX][EX2_MMAX], sv[EX2_MMAX], dummy[1]; Integer i, j, m, n, iter, ncolb, failinfo; Boolean wantp, wantq; static NagError fail; Vprintf("\nExample 2\n"); Vscanf(" %*[^\n]"); /* Skip Example 2 heading */ Vscanf(" %*[^\n]"); Vscanf("%ld%ld", &m, &n); if (m > EX2_MMAX || n > EX2_NMAX) { Vprintf("m or n is out of range.\n"); Vprintf("m = %2ld, n = %2ld\n", m, n); } 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 = TRUE; wantp = TRUE; ncolb = 0; fail.print = TRUE; f02wec(m, n, &a[0][0], tda, ncolb, dummy, (Integer)1, wantq, &q[0][0], tdq, sv, wantp, dummy, (Integer)1, &iter, e, &failinfo, &fail); if (fail.code != NE_NOERROR) exit(EXIT_FAILURE); 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"); } } }