/* nag_superlu_refine_lu (f11mhc) Example Program. * * Copyright 2005 Numerical Algorithms Group. * * Mark 8, 2005. */ #include #include #include #include #include /* Table of constant values */ static Integer c__1 = 1; static Integer c__80 = 80; static Integer c__0 = 0; int main(int argc, char *argv[]) { FILE *fpin, *fpout; char *outfile = 0; double flop, thresh; Integer exit_status = 0, i, j; Integer n, nnz, nnzl, nnzu, nrhs, nzlmx, nzlumx, nzumx; double *a = 0, *b = 0, *berr = 0, *ferr = 0, *lval = 0; double *uval = 0, *x = 0; Integer *icolzp = 0, *il = 0, *iprm = 0, *irowix = 0; Integer *iu = 0; /* Nag types */ Nag_OrderType order = Nag_ColMajor; Nag_MatrixType matrix = Nag_GeneralMatrix; Nag_DiagType diag = Nag_NonUnitDiag; Nag_ColumnPermutationType ispec; Nag_TransType trans; NagError fail; INIT_FAIL(fail); /* Check for command-line IO options */ fpin = nag_example_file_io(argc, argv, "-data", NULL); fpout = nag_example_file_io(argc, argv, "-results", NULL); (void) nag_example_file_io(argc, argv, "-nag_write", &outfile); fprintf(fpout, "nag_superlu_refine_lu (f11mhc) Example Program Results\n\n"); /* Skip heading in data file */ fscanf(fpin, "%*[^\n] "); /* Read order of matrix and number of right hand sides */ fscanf(fpin, "%ld%ld%*[^\n] ", &n, &nrhs); /* Read the matrix A */ if (!(icolzp = NAG_ALLOC(n+1, Integer))) { fprintf(fpout, "Allocation failure\n"); exit_status = -1; goto END; } for (i = 1; i <= n + 1; ++i) fscanf(fpin, "%ld%*[^\n] ", &icolzp[i - 1]); nnz = icolzp[n] - 1; /* Allocate memory */ if (!(irowix = NAG_ALLOC(nnz, Integer)) || !(a = NAG_ALLOC(nnz, double)) || !(il = NAG_ALLOC(7*n+8*nnz+4, Integer)) || !(iu = NAG_ALLOC(2*n+8*nnz+1, Integer)) || !(uval = NAG_ALLOC(8*nnz, double)) || !(lval = NAG_ALLOC(8*nnz, double)) || !(b = NAG_ALLOC(n * nrhs, double)) || !(x = NAG_ALLOC(n * nrhs, double)) || !(berr = NAG_ALLOC(nrhs, double)) || !(ferr = NAG_ALLOC(nrhs, double)) || !(iprm = NAG_ALLOC(7*n, Integer))) { fprintf(fpout, "Allocation failure\n"); exit_status = -1; goto END; } for (i = 1; i <= nnz; ++i) fscanf(fpin, "%lf%ld%*[^\n] ", &a[i - 1], &irowix[i - 1]); /* Read the right hand sides */ for (j = 1; j <= nrhs; ++j) { for (i = 1; i <= n; ++i) { fscanf(fpin, "%lf", &x[j*n + i - n - 1]); b[j*n + i - n - 1] = x[j*n + i - n - 1]; } fscanf(fpin, "%*[^\n] "); } /* Calculate COLAMD permutation */ ispec = Nag_Sparse_Colamd; /* nag_superlu_column_permutation (f11mdc). * Real sparse nonsymmetric linear systems, setup for * nag_superlu_lu_factorize (f11mec) */ nag_superlu_column_permutation(ispec, n, icolzp, irowix, iprm, &fail); if (fail.code != NE_NOERROR) { fprintf(fpout, "Error from nag_superlu_column_permutation (f11mdc).\n%s\n", fail.message); exit_status = 1; goto END; } /* Factorise */ thresh = 1.; nzlmx = 8*nnz; nzlumx = 8*nnz; nzumx = 8*nnz; /* nag_superlu_lu_factorize (f11mec). * LU factorization of real sparse matrix */ nag_superlu_lu_factorize(n, irowix, a, iprm, thresh, nzlmx, &nzlumx, nzumx, il, lval, iu, uval, &nnzl, &nnzu, &flop, &fail); if (fail.code != NE_NOERROR) { fprintf(fpout, "Error from nag_superlu_lu_factorize (f11mec).\n%s\n", fail.message); exit_status = 1; goto END; } /* Compute solution in array X */ trans = Nag_NoTrans; /* nag_superlu_solve_lu (f11mfc). * Solution of real sparse simultaneous linear equations * (coefficient matrix already factorized) */ nag_superlu_solve_lu(order, trans, n, iprm, il, lval, iu, uval, nrhs, x, n, &fail); if (fail.code != NE_NOERROR) { fprintf(fpout, "Error from nag_superlu_solve_lu (f11mfc).\n%s\n", fail.message); exit_status = 1; goto END; } /* Improve solution, and compute backward errors and estimated */ /* bounds on the forward errors */ /* nag_superlu_refine_lu (f11mhc). * Refined solution with error bounds of real system of * linear equations, multiple right-hand sides */ nag_superlu_refine_lu(order, trans, n, icolzp, irowix, a, iprm, il, lval, iu, uval, nrhs, b, n, x, n, ferr, berr, &fail); if (fail.code != NE_NOERROR) { fprintf(fpout, "Error from nag_superlu_refine_lu (f11mhc).\n%s\n", fail.message); exit_status = 1; goto END; } /* Print solution */ fprintf(fpout, "\n"); /* nag_gen_real_mat_print (x04cac). * Print real general matrix (easy-to-use) */ if (outfile) fclose(fpout); nag_gen_real_mat_print(order, matrix, diag, n, nrhs, x, n, "Solutions", outfile, &fail); if (outfile && !(fpout = fopen(outfile, "a"))) { exit_status = 2; goto END; } if (fail.code != NE_NOERROR) { fprintf(fpout, "Error from nag_gen_real_mat_print (x04cac).\n%s\n", fail.message); exit_status = 1; goto END; } /* nag_gen_real_mat_print_comp (x04cbc). * Print real general matrix (comprehensive) */ if (outfile) fclose(fpout); nag_gen_real_mat_print_comp(order, matrix, diag, nrhs, c__1, ferr, nrhs, "%8.2g", "Estimated Forward Error", Nag_NoLabels, NULL, Nag_NoLabels, NULL, c__80, c__0, outfile, &fail); if (outfile && !(fpout = fopen(outfile, "a"))) { exit_status = 2; goto END; } if (fail.code != NE_NOERROR) { fprintf(fpout, "Error from nag_gen_real_mat_print_comp (x04cbc).\n%s\n", fail.message); exit_status = 1; goto END; } /* nag_gen_real_mat_print_comp (x04cbc), see above. */ if (outfile) fclose(fpout); nag_gen_real_mat_print_comp(order, matrix, diag, nrhs, c__1, berr, nrhs, "%8.2g", "Backward Error", Nag_NoLabels, NULL, Nag_NoLabels, NULL, c__80, c__0, outfile, &fail); if (outfile && !(fpout = fopen(outfile, "a"))) { exit_status = 2; goto END; } if (fail.code != NE_NOERROR) { fprintf(fpout, "Error from nag_gen_real_mat_print_comp (x04cbc).\n%s\n", fail.message); exit_status = 1; goto END; } END: if (fpin != stdin) fclose(fpin); if (fpout != stdout) fclose(fpout); if (a) NAG_FREE(a); if (b) NAG_FREE(b); if (berr) NAG_FREE(berr); if (ferr) NAG_FREE(ferr); if (lval) NAG_FREE(lval); if (uval) NAG_FREE(uval); if (x) NAG_FREE(x); if (icolzp) NAG_FREE(icolzp); if (il) NAG_FREE(il); if (iprm) NAG_FREE(iprm); if (irowix) NAG_FREE(irowix); if (iu) NAG_FREE(iu); return exit_status; }