/* nag_sparse_nsym_precon_ssor_solve (f11ddc) Example Program. * * Copyright 2011, Numerical Algorithms Group. * * Mark 23, 2011. */ #include #include #include int main(void) { /* Scalars */ Integer exit_status = 0; double anorm, omega, sigmax, stplhs, stprhs, tol; Integer i, irevcm, iterm, itn, la, liwork, lwneed, lwork, m, maxitn, monit, n, nnz; /* Arrays */ char nag_enum_arg[100]; double *a = 0, *b = 0, *rdiag = 0, *wgt = 0, *work = 0, *x = 0; Integer *icol = 0, *irow = 0, *iwork = 0; /* NAG types */ Nag_SparseNsym_CheckData ckdd,ckxa; Nag_NormType norm; Nag_SparseNsym_PrecType precon; Nag_TransType trans; Nag_SparseNsym_Weight weight; Nag_SparseNsym_Method method; NagError fail, fail1; INIT_FAIL(fail); INIT_FAIL(fail1); printf("nag_sparse_nsym_precon_ssor_solve (f11ddc) Example Program Results"); printf("\n\n"); /* Skip heading in data file*/ scanf("%*[^\n]"); /* Read algorithmic parameters*/ scanf("%ld%ld%*[^\n]", &n, &m); scanf("%ld%*[^\n]", &nnz); la = 3 * nnz; lwork = MAX(n * (m + 3) + m * (m + 5) + 101, 7 * n + 100); liwork = 2 * n + 1; if ( !(a = NAG_ALLOC((la), double)) || !(b = NAG_ALLOC((n), double)) || !(rdiag = NAG_ALLOC((n), double)) || !(wgt = NAG_ALLOC((n), double)) || !(work = NAG_ALLOC((lwork), double)) || !(x = NAG_ALLOC((n), double)) || !(icol = NAG_ALLOC((la), Integer)) || !(irow = NAG_ALLOC((la), Integer)) || !(iwork = NAG_ALLOC((liwork), Integer)) ) { printf("Allocation failure\n"); exit_status = -1; goto END; } /* Read or initialize the parameters for the iterative solver*/ scanf("%s%*[^\n] ", nag_enum_arg); /* nag_enum_name_to_value(x04nac). * Converts NAG enum member name to value */ method = (Nag_SparseNsym_Method) nag_enum_name_to_value(nag_enum_arg); scanf("%s%*[^\n]", nag_enum_arg); precon = (Nag_SparseNsym_PrecType) nag_enum_name_to_value(nag_enum_arg); scanf("%s%*[^\n]", nag_enum_arg); norm = (Nag_NormType) nag_enum_name_to_value(nag_enum_arg); scanf("%ld%*[^\n]", &iterm); scanf("%lf%ld%*[^\n]", &tol, &maxitn); scanf("%lf%lf%*[^\n]", &anorm, &sigmax); scanf("%lf%*[^\n]", &omega); /* Read the non-zero elements of the matrix a*/ for (i = 0; i < nnz; i++) scanf("%lf%ld%ld%*[^\n] ", &a[i], &irow[i], &icol[i]); /* Read right-hand side vector b and initial approximate solution x*/ for (i = 0; i < n; i++) scanf("%lf", &b[i]); scanf("%*[^\n]"); for (i = 0; i < n; i++) scanf("%lf", &x[i]); weight = Nag_SparseNsym_UnWeighted; monit = 0; /* Call to initialize the solver * nag_sparse_nsym_basic_setup (f11bdc) * Real sparse nonsymmetric linear systems, setup routine */ nag_sparse_nsym_basic_setup(method, precon, norm, weight, iterm, n, m, tol, maxitn, anorm, sigmax, monit, &lwneed, work, lwork, &fail); /* Calculate reciprocal diagonal matrix elements if necessary*/ if (precon == Nag_SparseNsym_Prec) { for (i = 0; i < n; i++) iwork[i] = 0; for (i = 0; i < nnz; i++) { if (irow[i] == icol[i]) { iwork[irow[i]-1]++; if (a[i] == 0.0) { printf("Matrix has a zero diagonal element \n"); goto END; } rdiag[(irow[i]-1)] = 1.0/a[i]; } } for (i = 0; i < n; i++) { if (iwork[i] == 0) { printf("Matrix has a missing diagonal element \n"); goto END; } if (iwork[i] >= 2) { printf("Matrix has a multiple diagonal element \n"); goto END; } } } /* call solver repeatedly to solve the equations*/ irevcm = 0; ckxa = Nag_SparseNsym_Check; ckdd = Nag_SparseNsym_Check; while (irevcm!=4) { /* nag_sparse_nsym_basic_solver (f11bec) * Real sparse nonsymmetric linear systems, solver routine * preconditioned RGMRES, CGS, Bi-CGSTAB or TFQMR method */ nag_sparse_nsym_basic_solver(&irevcm, x, b, wgt, work, lwork, &fail); switch (irevcm) { case 1: /* Compute matrix-vector product using * nag_sparse_nsym_matvec (f11xac) * Real sparse nonsymmetric matrix vector multiply */ trans = Nag_NoTrans; nag_sparse_nsym_matvec(trans, n, nnz, a, irow, icol, ckxa, x, b, &fail1); ckxa = Nag_SparseNsym_NoCheck; break; case -1: /* Compute transposed matrix-vector product */ trans = Nag_Trans; nag_sparse_nsym_matvec(trans, n, nnz, a, irow, icol, ckxa, x, b, &fail1); ckxa = Nag_SparseNsym_NoCheck; break; case 2: /* SSOR preconditioning using * nag_sparse_nsym_precon_ssor_solve (f11ddc) * Solution of linear system involving preconditioning matrix generated * by applying SSOR to real sparse nonsymmetric matrix */ trans = Nag_NoTrans; nag_sparse_nsym_precon_ssor_solve(trans, n, nnz, a, irow, icol, rdiag, omega, ckdd, x, b,&fail1); ckdd = Nag_SparseNsym_NoCheck; break; } if (fail1.code != NE_NOERROR) irevcm = 6; } if (fail.code != NE_NOERROR) { printf("Error from nag_sparse_nsym_basic_solver (f11bec)\n%s\n", fail.message); exit_status = 1; goto END; } /* nag_sparse_nsym_basic_diagnostic (f11bfc) * Real sparse nonsymmetric linear systems, diagnostic */ nag_sparse_nsym_basic_diagnostic(&itn, &stplhs, &stprhs, &anorm, &sigmax, work, lwork, &fail); printf(" Converged in %11ld iterations\n", itn); printf(" Matrix norm = %9.3e\n", anorm); printf(" Final residual norm = %9.3e\n\n", stplhs); /* Output x*/ printf(" Solution of linear system\n"); for (i = 0; i < n; i++) printf("%16.4e\n", x[i]); END: if (a) NAG_FREE(a); if (b) NAG_FREE(b); if (rdiag) NAG_FREE(rdiag); if (wgt) NAG_FREE(wgt); if (work) NAG_FREE(work); if (x) NAG_FREE(x); if (icol) NAG_FREE(icol); if (irow) NAG_FREE(irow); if (iwork) NAG_FREE(iwork); return exit_status; }