/* nag_superlu_condition_number_lu (f11mgc) Example Program. * * Copyright 2005 Numerical Algorithms Group. * * Mark 8, 2005. */ #include #include #include #include int main(void) { double anorm, d, flop, rcond, thresh; Integer exit_status = 0, i, n, nnz, nnzl, nnzu, nzlmx, nzlumx, nzumx; double *a = 0, *lval = 0, *uval = 0; Integer *icolzp = 0, *il = 0, *iprm = 0, *irowix = 0; Integer *iu = 0; /* Nag types */ NagError fail; Nag_ColumnPermutationType ispec; Nag_NormType norm; INIT_FAIL(fail); printf("nag_superlu_condition_number_lu (f11mgc) Example Program " "Results\n\n"); /* Skip heading in data file */ scanf("%*[^\n] "); /* Read order of matrix */ scanf("%ld%*[^\n] ", &n); /* Read the matrix A */ if (!(icolzp = NAG_ALLOC(n+1, Integer))) { printf("Allocation failure\n"); exit_status = -1; goto END; } for (i = 1; i <= n + 1; ++i) scanf("%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)) || !(iprm = NAG_ALLOC(7*n, Integer))) { printf("Allocation failure\n"); exit_status = -1; goto END; } for (i = 1; i <= nnz; ++i) scanf("%lf%ld%*[^\n] ", &a[i - 1], &irowix[i - 1]); /* 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) { printf( "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) { printf("Error from nag_superlu_lu_factorize (f11mec).\n%s\n", fail.message); exit_status = 1; goto END; } /* Calculate norm */ norm = Nag_RealOneNorm; /* nag_superlu_matrix_norm (f11mlc). * 1-norm, infinity-norm, largest absolute element, real * general matrix */ nag_superlu_matrix_norm(norm, &anorm, n, icolzp, irowix, a, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_superlu_matrix_norm (f11mlc).\n%s\n", fail.message); exit_status = 1; goto END; } /* Calculate condition number */ /* nag_superlu_condition_number_lu (f11mgc). * Estimate condition number of real matrix, matrix already * factorized by nag_superlu_lu_factorize (f11mec) */ nag_superlu_condition_number_lu(norm, n, il, lval, iu, uval, anorm, &rcond, &fail); if (fail.code != NE_NOERROR) { printf( "Error from nag_superlu_condition_number_lu (f11mgc).\n%s\n", fail.message); exit_status = 1; goto END; } /* Output result */ d = 1 / rcond; printf("\n Norm, Condition number\n %7.3f,%7.3f\n\n", anorm, d); END: if (a) NAG_FREE(a); if (lval) NAG_FREE(lval); if (uval) NAG_FREE(uval); 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; }