/* nag_dgejsv (f08khc) Example Program. * * Copyright 2011 Numerical Algorithms Group. * * Mark 23, 2011. */ #include #include #include #include #include #include #include int main(void) { /* Scalars */ double eps, serrbd; Integer exit_status = 0; Integer pda, pdu, pdv; Integer i, j, m, n, n_uvecs, n_vvecs; /* Arrays */ double *a = 0, *rcondu = 0, *rcondv = 0, *s = 0, *u = 0, *v = 0; double work[7]; Integer iwork[3]; char nag_enum_arg[40]; /* Nag Types */ Nag_OrderType order; Nag_Preprocess joba; Nag_LeftVecsType jobu; Nag_RightVecsType jobv; Nag_ZeroCols jobr; Nag_TransType jobt; Nag_Perturb jobp; NagError fail; #ifdef NAG_COLUMN_MAJOR #define A(I, J) a[(J-1)*pda + I-1] order = Nag_ColMajor; #else #define A(I, J) a[(I-1)*pda + J-1] order = Nag_RowMajor; #endif INIT_FAIL(fail); printf("nag_dgejsv (f08khc) Example Program Results\n\n"); jobu = Nag_LeftSpan; jobv = Nag_RightVecs; jobr = Nag_ZeroColsRestrict; jobt = Nag_NoTrans; jobp = Nag_PerturbOff; /* Skip heading in data file*/ scanf("%*[^\n]"); scanf("%ld%ld%*[^\n]", &m, &n); if (n < 0 || m < n) { printf("Invalid n or nrhs\n"); exit_status = 1; goto END;; } /* Read Nag type arguments by name and convert to value */ scanf(" %s%*[^\n]", nag_enum_arg); /* nag_enum_name_to_value(x04nac). * Converts NAG enum member name to value */ joba = (Nag_Preprocess) nag_enum_name_to_value(nag_enum_arg); scanf(" %s%*[^\n]", nag_enum_arg); jobu = (Nag_LeftVecsType) nag_enum_name_to_value(nag_enum_arg); scanf(" %s%*[^\n]", nag_enum_arg); jobv = (Nag_RightVecsType) nag_enum_name_to_value(nag_enum_arg); scanf(" %s%*[^\n]", nag_enum_arg); jobr = (Nag_ZeroCols) nag_enum_name_to_value(nag_enum_arg); scanf(" %s%*[^\n]", nag_enum_arg); jobt = (Nag_TransType) nag_enum_name_to_value(nag_enum_arg); scanf(" %s%*[^\n]", nag_enum_arg); jobp = (Nag_Perturb) nag_enum_name_to_value(nag_enum_arg); /* Size of u and v depends on some of the above Nag type arguments. */ n_uvecs = 1; if (jobu==Nag_LeftVecs) { n_uvecs = m; } else if (jobu==Nag_LeftSpan) { n_uvecs = n; } else if (jobu==Nag_NotLeftWork && jobv==Nag_RightVecs && jobt==Nag_Trans && m==n) { n_uvecs = m; } if (jobv==Nag_NotRightVecs) { n_vvecs = 1; } else { n_vvecs = n; } #ifdef NAG_COLUMN_MAJOR pda = m; pdu = m; pdv = n; #else pda = n; pdu = n_uvecs; pdv = n_vvecs; #endif if (!(a = NAG_ALLOC(m*n, double)) || !(rcondu = NAG_ALLOC(m, double)) || !(rcondv = NAG_ALLOC(m, double)) || !(s = NAG_ALLOC(n, double)) || !(u = NAG_ALLOC(m*n_uvecs, double)) || !(v = NAG_ALLOC(n_vvecs*n_vvecs, double))) { printf("Allocation failure\n"); exit_status = -1; goto END; } /* Read the m by n matrix A from data file*/ for (i = 1; i <= m; i++) for (j = 1; j <= n; j++) scanf("%lf", &A(i, j)); scanf("%*[^\n]"); /* nag_dgejsv (f08khc) * Compute the singular values and left and right singular vectors * of A (A = U*S*V^T, m>=n). */ nag_dgejsv(order, joba, jobu, jobv, jobr, jobt, jobp, m, n, a, pda, s, u, pdu, v, pdv, work, iwork, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_dgejsv (f08khc).\n%s\n", fail.message); exit_status = 1; goto END; } /* Get the machine precision, eps and compute the approximate * error bound for the computed singular values. Note that for * the 2-norm, s[0] = norm(A). */ eps = nag_machine_precision; serrbd = eps * s[0]; /* Print (possibly scaled) singular values. */ if (fabs(work[0] - work[1]) < 2.0 * eps) { /* No scaling required*/ printf("Singular values\n"); for (j = 0; j < n; j++) printf("%8.4f", s[j]); } else { printf("Scaled singular values\n"); for (j = 0; j < n; j++) printf("%8.4f", s[j]); printf("\nFor true singular values, multiply by a/b,\n"); printf("where a = %f and b = %f", work[0], work[1]); } printf("\n\n"); /* Print left and right (spanning) singular vectors, if requested. using * nag_gen_real_mat_print (x04cac) * Print real general matrix (easy-to-use) */ if (jobu==Nag_LeftVecs || jobu==Nag_LeftSpan) { fflush(stdout); nag_gen_real_mat_print(order, Nag_GeneralMatrix, Nag_NonUnitDiag, m, n, u, pdu, "Left singular vectors", 0, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_gen_real_mat_print (x04cac).\n%s\n", fail.message); exit_status = 1; goto END; } } if (jobv==Nag_RightVecs || jobv==Nag_RightVecsJRots) { printf("\n"); fflush(stdout); nag_gen_real_mat_print(order, Nag_GeneralMatrix, Nag_NonUnitDiag, n, n, v, pdv, "Right singular vectors", 0, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_gen_real_mat_print (x04cac).\n%s\n", fail.message); exit_status = 1; goto END; } } /* nag_ddisna (f08flc) * Estimate reciprocal condition numbers for the singular vectors. */ nag_ddisna(Nag_LeftSingVecs, m, n, s, rcondu, &fail); if (fail.code == NE_NOERROR) nag_ddisna(Nag_RightSingVecs, m, n, s, rcondv, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_ddisna (f08flc).\n%s\n", fail.message); exit_status = 1; goto END; } if (joba==Nag_ColpivRrankCond || joba==Nag_FullpivRrankCond) { printf("\n\nEstimate of the condition number of column equilibrated A\n"); printf("%11.1e", work[2]); } /* Print the approximate error bounds for the singular values and vectors. */ printf("\n\nError estimate for the singular values\n%11.1e", serrbd); printf("\n\nError estimates for left singular vectors\n"); for (i = 0; i < n; i++) printf("%11.1e", serrbd/rcondu[i]); printf("\n\nError estimates for right singular vectors\n"); for (i = 0; i < n; i++) printf("%11.1e", serrbd/rcondv[i]); printf("\n"); END: if (a) NAG_FREE(a); if (rcondu) NAG_FREE(rcondu); if (rcondv) NAG_FREE(rcondv); if (s) NAG_FREE(s); if (u) NAG_FREE(u); if (v) NAG_FREE(v); return exit_status; }