/* nag_sparse_real_symm_precon_ssor_solve (f11jdc) Example Program.
*
* Copyright 2022 Numerical Algorithms Group.
*
* Mark 28.3, 2022.
*/
#include <nag.h>
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_SparseSym_Method method;
Nag_SparseSym_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_real_symm_precon_ssor_solve (f11jdc) Example Program Results");
printf("\n\n");
/* Skip heading in data file */
scanf("%*[^\n]");
/* Read algorithmic parameters */
scanf("%" NAG_IFMT "%*[^\n]", &n);
scanf("%" NAG_IFMT "%*[^\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("%99s%*[^\n] ", nag_enum_arg);
/* nag_enum_name_to_value (x04nac).
* Converts NAG enum member name to value
*/
method = (Nag_SparseSym_Method)nag_enum_name_to_value(nag_enum_arg);
scanf("%99s%*[^\n] ", nag_enum_arg);
precon = (Nag_SparseSym_PrecType)nag_enum_name_to_value(nag_enum_arg);
scanf("%99s%*[^\n] ", nag_enum_arg);
sigcmp = (Nag_SparseSym_Bisection)nag_enum_name_to_value(nag_enum_arg);
scanf("%99s%*[^\n] ", nag_enum_arg);
norm = (Nag_NormType)nag_enum_name_to_value(nag_enum_arg);
scanf("%" NAG_IFMT "%*[^\n] ", &iterm);
scanf("%lf%" NAG_IFMT "%*[^\n]", &tol, &maxitn);
scanf("%lf%lf%*[^\n]", &anorm, &sigmax);
scanf("%lf%" NAG_IFMT "%*[^\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_real_symm_sort (f11zbc).
*/
nnz1 = nnz;
dup = Nag_SparseSym_RemoveDups;
zero = Nag_SparseSym_RemoveZeros;
nag_sparse_real_symm_sort(n, &nnz1, a, irow, icol, dup, zero, istr, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_sparse_real_symm_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 %" NAG_IFMT " to %" NAG_IFMT
"\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 %" NAG_IFMT "\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_real_symm_basic_setup (f11gdc)
*/
weight = Nag_SparseSym_UnWeighted;
monit = 0;
nag_sparse_real_symm_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_real_symm_basic_solver (f11gec).
* Real sparse symmetric linear systems, preconditioned conjugate gradient
* or Lanczos method.
*/
nag_sparse_real_symm_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_real_symm_matvec (f11xec).
*/
nag_sparse_real_symm_matvec(n, nnz, a, irow, icol, ckxe, x, b, &fail1);
ckxe = Nag_SparseSym_NoCheck;
break;
case 2:
/* SSOR preconditioning
* nag_sparse_real_symm_precon_ssor_solve (f11jdc).
* Solution of linear system involving preconditioning matrix
* generated by applying SSOR to real sparse symmetric matrix
*/
nag_sparse_real_symm_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_real_symm_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_real_symm_basic_diag (f11gfc).
*/
nag_sparse_real_symm_basic_diag(&itn, &stplhs, &stprhs, &anorm, &sigmax, &its,
&sigerr, commarray, lcomm, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_sparse_real_symm_basic_diag (f11gfc).\n%s\n",
fail.message);
exit_status = 4;
goto END;
}
printf("Converged in %10" NAG_IFMT " 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:
NAG_FREE(a);
NAG_FREE(b);
NAG_FREE(rdiag);
NAG_FREE(wgt);
NAG_FREE(commarray);
NAG_FREE(x);
NAG_FREE(icol);
NAG_FREE(irow);
NAG_FREE(istr);
return exit_status;
}