/* nag_sparse_sym_precon_ssor_solve (f11jdc) Example Program. * * Copyright 2011, Numerical Algorithms Group. * * Mark 23, 2011. */ #include #include #include int main(void) { /* Scalars */ Integer exit_status = 0; double anorm, omega, sigerr, sigmax, sigtol, stplhs, stprhs, tol; Integer i, irevcm, iterm, itn, its, j, listr, lcneed, lcomm, maxitn, maxits, monit, n, nnz, nnz1; /* Arrays */ char nag_enum_arg[100]; double *a = 0, *b = 0, *rdiag = 0, *wgt = 0, *commarray = 0, *x = 0; Integer *icol = 0, *irow = 0, *istr = 0; /* NAG types */ Nag_NormType norm; Nag_SparseNsym_Method method; Nag_SparseNsym_PrecType precon; Nag_SparseSym_Bisection sigcmp; Nag_SparseSym_CheckData ckjd, ckxe; Nag_SparseSym_Dups dup; Nag_SparseSym_Weight weight; Nag_SparseSym_Zeros zero; NagError fail, fail1; INIT_FAIL(fail); printf("nag_sparse_sym_precon_ssor_solve (f11jdc) Example Program Results"); printf("\n\n"); /* Skip heading in data file*/ scanf("%*[^\n]"); /* Read algorithmic parameters*/ scanf("%ld%*[^\n]", &n); scanf("%ld%*[^\n]", &nnz); /* Allocate memory */ listr = n + 1; lcomm = 6 * n + 120; if ( !(a = NAG_ALLOC(nnz, double)) || !(b = NAG_ALLOC(n, double)) || !(rdiag = NAG_ALLOC(n, double)) || !(wgt = NAG_ALLOC(n, double)) || !(commarray = NAG_ALLOC(lcomm, double)) || !(x = NAG_ALLOC(n, double)) || !(icol = NAG_ALLOC(nnz, Integer)) || !(irow = NAG_ALLOC(nnz, Integer)) || !(istr = NAG_ALLOC(listr, Integer)) ) { printf("Allocation failure\n"); exit_status = -1; goto END; } 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); sigcmp = (Nag_SparseSym_Bisection) 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%ld%*[^\n]", &sigtol, &maxits); scanf("%lf%*[^\n]", &omega); /* Read the matrix a */ for (i = 0; i <= nnz-1; i++) scanf("%lf%"NAG_IFMT "%"NAG_IFMT "%*[^\n] ", &a[i], &irow[i], &icol[i]); /* Sort matrix a removing zero or duplicate elements using * nag_sparse_sym_sort (f11zbc). */ nnz1 = nnz; dup = Nag_SparseSym_RemoveDups; zero = Nag_SparseSym_RemoveZeros; nag_sparse_sym_sort (n, &nnz1, a, irow, icol, dup, zero, istr, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_sparse_sym_sort (f11zbc).\n%s\n", fail.message); exit_status = 1; goto END; } if (nnz != nnz1) { printf("Warning, input Matrix has zero or duplicate elements\n"); printf(" nnz has been reduced from %ld to %ld\n", nnz,nnz1); nnz = nnz1; } /* Check for zero diagonal matrix elements and calculate reciprocals.*/ for (i = 0; i < n; i++) { /* j points to last element in row i */ j = istr[i+1]-2; if (irow[j] == icol[j]) rdiag[irow[j]-1] = 1.0/a[j]; else { printf("Matrix has a missing element for diagonal %ld\n",i); goto END; } } /* Read right-hand side vector b and initial approximate solution x*/ for (i = 0; i <= n-1; i++) scanf("%lf", &b[i]); scanf("%*[^\n]"); for (i = 0; i <= n-1; i++) scanf("%lf", &x[i]); /* Initialize the basic symmteric solver (f11gec) using * nag_sparse_sym_basic_setup (f11gdc) */ weight = Nag_SparseSym_UnWeighted; monit = 0; nag_sparse_sym_basic_setup(method, precon, sigcmp, norm, weight, iterm, n, tol, maxitn, anorm, sigmax, sigtol, maxits, monit, &lcneed, commarray, lcomm, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_sparse_sym_setup (f11gdc).\n%s\n", fail.message); exit_status = 2; goto END; } /* call solver repeatedly to solve the equations */ irevcm = 0; ckxe = Nag_SparseSym_Check; ckjd = Nag_SparseSym_Check; while (1) { /* nag_sparse_sym_basic_solver(f11gec). * Real sparse symmetric linear systems, preconditioned conjugate gradient * or Lanczos method. */ nag_sparse_sym_basic_solver(&irevcm, x, b, wgt, commarray, lcomm, &fail); if (irevcm != 4) { INIT_FAIL(fail1); switch (irevcm) { case 1: /* Compute sparse symmetric matrix vector product using * nag_sparse_sym_matvec (f11xec). */ nag_sparse_sym_matvec(n, nnz, a, irow, icol, ckxe, x, b, &fail1); ckxe = Nag_SparseSym_NoCheck; break; case 2: /* SSOR preconditioning * nag_sparse_sym_precon_ssor_solve (f11jdc). * Solution of linear system involving preconditioning matrix * generated by applying SSOR to real sparse symmetric matrix */ nag_sparse_sym_precon_ssor_solve(n, nnz, a, irow, icol, rdiag, omega, ckjd, x, b, &fail1); ckjd = Nag_SparseSym_NoCheck; } if (fail1.code != NE_NOERROR) irevcm = 6; } else if (fail.code != NE_NOERROR) { printf("Error from nag_sparse_sym_basic_solver (f11gec).\n%s\n", fail.message); exit_status = 3; goto END; } else goto END_LOOP; } END_LOOP: /* Obtain and print diagnostic statistics using * nag_sparse_sym_basic_diagnostic (f11gfc). */ nag_sparse_sym_basic_diagnostic(&itn, &stplhs, &stprhs, &anorm, &sigmax, &its, &sigerr, commarray, lcomm, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_sparse_sym_basic_diagnostic (f11gfc).\n%s\n", fail.message); exit_status = 4; goto END; } printf("Converged in %10ld iterations \n", itn); printf("Final residual norm = %11.3e\n\n", stplhs); /* Output solution */ printf("%16s\n","Solution"); for (i = 0; i <= n - 1; 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 (commarray) NAG_FREE(commarray); if (x) NAG_FREE(x); if (icol) NAG_FREE(icol); if (irow) NAG_FREE(irow); if (istr) NAG_FREE(istr); return exit_status; }