/* nag_dgesvj (f08kjc) 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 i, j, lwork, m, mv, n, n_vrows, n_vcols, pda, pdv, ranka; /* Arrays */ double *a = 0, *rcondu = 0, *rcondv = 0, *s = 0, *v = 0, *work = 0; char nag_enum_arg[40]; /* Nag Types */ Nag_OrderType order; Nag_MatrixType joba; Nag_LeftVecsType jobu; Nag_RightVecsType jobv; NagError fail; #ifdef NAG_COLUMN_MAJOR #define A(I, J) a[(J-1)*pda + I-1] #define V(I, J) v[(J-1)*pdv + I-1] order = Nag_ColMajor; #else #define A(I, J) a[(I-1)*pda + J-1] #define V(I, J) v[(I-1)*pdv + J-1] order = Nag_RowMajor; #endif INIT_FAIL(fail); printf("nag_dgesvj (f08kjc) Example Program Results\n\n"); /* Skip heading in data file*/ scanf("%*[^\n]"); scanf("%ld%ld%*[^\n]", &m, &n); if (n < 0 || m < n) { printf("Invalid m or n\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_MatrixType) 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); n_vcols = n; n_vrows = n; mv = 0; if (jobv==Nag_RightVecsMV) { scanf("%"NAG_IFMT, &mv); n_vrows = mv; } else if (jobv==Nag_NotRightVecs) { n_vrows = 1; n_vcols = 1; } scanf("%*[^\n]"); #ifdef NAG_COLUMN_MAJOR pda = m; pdv = n_vrows; #else pda = n; pdv = n_vcols; #endif lwork = n + m; if (!(a = NAG_ALLOC(m*n, double)) || !(rcondu = NAG_ALLOC(m, double)) || !(rcondv = NAG_ALLOC(m, double)) || !(s = NAG_ALLOC(n, double)) || !(v = NAG_ALLOC(n_vrows*n_vcols, double)) || !(work = NAG_ALLOC(lwork, double)) ) { printf("Allocation failure\n"); exit_status = -1; goto END; } /* Read the m by n matrix A from data file*/ if (joba == Nag_GeneralMatrix) { for (i = 1; i <= m; i++) for (j = 1; j <= n; j++) scanf("%lf", &A(i, j)); } else if (joba == Nag_UpperMatrix) { for (i = 1; i <= m; i++) for (j = i; j <= n; j++) scanf("%lf", &A(i, j)); } else { for (i = 1; i <= m; i++) for (j = 1; j <= i; j++) scanf("%lf", &A(i, j)); } scanf("%*[^\n]"); /* jobv==Nag_RightVecsMV means that the first mv rows of v must be set. */ if (jobv==Nag_RightVecsMV) { for (i = 1; i <= mv; i++) for (j = 1; j <= n; j++) scanf("%lf", &V(i, j)); scanf("%*[^\n]"); } /* nag_dgesvj (f08kjc) * Compute the singular values and left and right singular vectors * of A (A = U*S*V, m>=n). */ nag_dgesvj(order, joba, jobu, jobv, m, n, a, pda, s, mv, v, pdv, work, lwork, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_dgesvj (f08kjc).\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 solution*/ printf("Singular values\n "); for (j = 0; j < n; j++) printf("%8.4f", s[j]); printf("\n\n"); if (fabs(work[0] - 1.0) > eps) printf("Values need scaling by factor = %13.5e\n\n", work[0]); ranka = (Integer) work[1]; printf("Rank of A = %5ld\n\n",ranka); if (jobu != Nag_NotLeftVecs) { /* nag_gen_real_mat_print (x04cac) * Print real general matrix (easy-to-use) */ fflush(stdout); nag_gen_real_mat_print(order, Nag_GeneralMatrix, Nag_NonUnitDiag, m, ranka, a, pda, "Left spanning 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) { printf("\n"); fflush(stdout); nag_gen_real_mat_print(order, Nag_GeneralMatrix, Nag_NonUnitDiag, n, n, v, pdv, "Right singular vectors", 0, &fail); } else if (jobv == Nag_RightVecsMV) { printf("\n"); fflush(stdout); nag_gen_real_mat_print(order, Nag_GeneralMatrix, Nag_NonUnitDiag, mv, n, v, pdv, "Right singular vectors applied to input V", 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); 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; } /* Print the approximate error bounds for the singular values and vectors. */ printf("\nError estimate for the singular values\n"); printf("%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 (v) NAG_FREE(v); if (work) NAG_FREE(work); return exit_status; }