/* nag_real_partial_svd (f02wgc) Example Program. * * Copyright 2009, Numerical Algorithms Group. * * Mark 9, 2009. */ /* Pre-processor includes */ #include #include #include #include #include #ifdef __cplusplus extern "C" { #endif static void NAG_CALL av(Integer *iflag, Integer m, Integer n, const double x[], double ax[], Nag_Comm *comm); #ifdef __cplusplus } #endif int main(void) { /*Integer scalar and array declarations */ Integer exit_status = 0; Integer i, m, n, nconv, ncv, nev; Integer pdu, pdv; Nag_Comm comm; NagError fail; /*Double scalar and array declarations */ static double ruser[1] = {-1.0}; double *resid = 0, *sigma = 0, *u = 0, *v = 0; Nag_OrderType order; INIT_FAIL(fail); printf("nag_real_partial_svd (f02wgc) Example Program Results\n\n"); /* For communication with user-supplied functions: */ comm.user = ruser; /* Skip heading in data file*/ scanf("%*[^\n] "); scanf("%ld%ld%ld%ld%*[^\n]", &m, &n, &nev, &ncv); #ifdef NAG_COLUMN_MAJOR order = Nag_ColMajor; pdu = m; pdv = n; #define U(I, J) u[(J-1)*pdu + I-1] #define V(I, J) v[(J-1)*pdv + I-1] #else order = Nag_RowMajor; pdu = ncv; pdv = ncv; #define U(I, J) u[(I-1)*pdu + J-1] #define V(I, J) v[(I-1)*pdv + J-1] #endif if (!(resid = NAG_ALLOC(m, double)) || !(sigma = NAG_ALLOC(ncv, double)) || !(u = NAG_ALLOC(m*ncv, double)) || !(v = NAG_ALLOC(n*ncv, double))) { printf("Allocation failure\n"); exit_status = -1; goto END; } /* * nag_real_partial_svd (f02wgc) * Computes leading terms in the singular value decomposition of * a real general matrix; also computes corresponding left and right * singular vectors. */ nag_real_partial_svd(order, m, n, nev, ncv, av, &nconv, sigma, u, pdu, v, pdv, resid, &comm, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_real_partial_svd (f02wgc).\n%s\n", fail.message); exit_status = 1; goto END; } /* Print computed residuals*/ printf("%s\n", " Singular Value Residual"); for (i = 0; i < nconv; i++) printf("%10.5f %16.2g\n", sigma[i], resid[i]); printf("\n"); END: NAG_FREE(resid); NAG_FREE(sigma); NAG_FREE(u); NAG_FREE(v); return exit_status; } static void NAG_CALL av(Integer *iflag, Integer m, Integer n, const double x[], double ax[], Nag_Comm *comm) { Integer i, j; double one = 1.0, zero = 0.0; double h, k, s, t; /* Matrix vector multiply subroutine Computing w <- A*x or w <- Trans(A)*x.*/ if (comm->user[0] == -1.0) { printf("(User-supplied callback av, first invocation.)\n"); comm->user[0] = 0.0; } h = one/(double)(m+1); k = one/(double)(n+1); if (*iflag == 1) { for (i = 0; i < m; i++) ax[i] = zero; t = zero; for (j = 0; j < n; j++) { t = t+k; s = zero; for (i = 0; i < MIN(j+1, m); i++) { s = s+h; ax[i] = ax[i]+k*s*(t-one)*x[j]; } for (i = j+1; i < m; i++) { s = s+h; ax[i] = ax[i]+k*t*(s-one)*x[j]; } } } else { for (i = 0; i < n; i++) ax[i] = zero; t = zero; for (j = 0; j < n; j++) { t = t+k; s = zero; for (i = 0; i < MIN(j+1, m); i++) { s = s+h; ax[j] = ax[j]+k*s*(t-one)*x[i]; } for (i = j+1; i < m; i++) { s = s+h; ax[j] = ax[j]+k*t*(s-one)*x[i]; } } } return; }