/* 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 */ 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"); /* Skip heading in data file*/ scanf("%*[^\n] "); scanf("%ld%ld%ld%ld%*[^\n]", &m, &n, &nev, &ncv); #ifdef NAG_COLUMN_MAJOR pdu = m; #define U(I, J) u[(J-1)*pdu + I-1] pdv = n; #define V(I, J) v[(J-1)*pdv + I-1] order = Nag_ColMajor; #else pdu = ncv; #define U(I, J) u[(I-1)*pdu + J-1] pdv = ncv; #define V(I, J) v[(I-1)*pdv + J-1] order = Nag_RowMajor; #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; } /* Initialize for problem.*/ /* * 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: if (resid) NAG_FREE(resid); if (sigma) NAG_FREE(sigma); if (u) NAG_FREE(u); if (v) NAG_FREE(v); return exit_status; } /* Pre-processor defines */ #define ONE 1.00e+0 #define ZERO 0.00e+0 static void NAG_CALL av(Integer *iflag, Integer m, Integer n, const double x[], double ax[], Nag_Comm *comm) { Integer i, j; double one, zero; double h, k, s, t; /* Assign parameter values to corresponding variables */ one = ONE; zero = ZERO; /* Matrix vector subroutines*/ /* Computes w <- A*x or w <- Trans(A)*x.*/ 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; }