/* nag_sparseig_real_symm_band_solve (f12fgc) Example Program.
*
* Copyright 2024 Numerical Algorithms Group.
*
* Mark 30.1, 2024.
*/
#include <math.h>
#include <nag.h>
#include <stdio.h>
int main(void) {
/* Constants */
Integer licomm = 140;
/* Scalars */
double alpha, h2, resr, sigma = 0.0;
Integer exit_status, i, j, k, kl, ku, ksub, ksup, lcomm;
Integer ldab, n, nconv, ncv, nev, nx;
/* Nag types */
NagError fail;
/* Arrays */
double *ab = 0, *ax = 0, *comm = 0, *eigv = 0, *eigest = 0, *mb = 0;
double *resid = 0, *v = 0;
Integer *icomm = 0;
exit_status = 0;
INIT_FAIL(fail);
printf("nag_sparseig_real_symm_band_solve (f12fgc)"
" Example Program Results\n\n");
/* Skip heading in data file */
scanf("%*[^\n] ");
scanf("%" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT "%*[^\n] ", &nx, &nev, &ncv);
n = nx * nx;
kl = nx;
ku = nx;
ldab = 2 * kl + ku + 1;
lcomm = 3 * n + ncv * ncv + 8 * ncv + 60;
/* Allocate memory */
if (!(ab = NAG_ALLOC(ldab * n, double)) || !(ax = NAG_ALLOC(n, double)) ||
!(comm = NAG_ALLOC(lcomm, double)) || !(eigv = NAG_ALLOC(ncv, double)) ||
!(eigest = NAG_ALLOC(ncv, double)) || !(mb = NAG_ALLOC(1, double)) ||
!(resid = NAG_ALLOC(n, double)) || !(v = NAG_ALLOC(n * ncv, double)) ||
!(icomm = NAG_ALLOC(licomm, Integer))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Initialize communication arrays for problem using
nag_sparseig_real_symm_band_init (f12ffc). */
nag_sparseig_real_symm_band_init(n, nev, ncv, icomm, licomm, comm, lcomm,
&fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_sparseig_real_symm_band_init"
" (f12ffc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
/* Construct the matrix A in banded form and store in AB */
/* ku, kl are number of superdiagonals and subdiagonals within the
band of matrices A and M. */
for (j = 0; j < n * ldab; ++j) {
ab[j] = 0.0;
}
/* Main diagonal of A. */
h2 = 1.0 / ((nx + 1) * (nx + 1));
k = kl + ku;
for (j = 0; j < n; ++j) {
ab[k] = 4.0 / h2;
k = k + ldab;
}
/* First subdiagonal and superdiagonal of A. */
ksup = kl + ku - 1;
ksub = kl + ku + 1;
for (i = 0; i < nx; ++i) {
ksup = ksup + ldab;
for (j = 0; j < nx - 1; ++j) {
ab[ksub] = -1.0 / h2;
ab[ksup] = -1.0 / h2;
ksup = ksup + ldab;
ksub = ksub + ldab;
}
ksub = ksub + ldab;
}
/* kl-th subdiagonal and ku-th superdiagonal. */
ksup = kl + nx * ldab;
ksub = 2 * kl + ku;
for (i = 0; i < nx - 1; ++i) {
for (j = 0; j < nx; ++j) {
ab[ksup] = -1.0 / h2;
ab[ksub] = -1.0 / h2;
ksup = ksup + ldab;
ksub = ksub + ldab;
}
}
/* Find eigenvalues of largest magnitude and the corresponding
* eigenvectors using nag_sparseig_real_symm_band_solve (f12fgc).
*/
nag_sparseig_real_symm_band_solve(kl, ku, ab, mb, sigma, &nconv, eigv, v,
resid, v, comm, icomm, &fail);
if (fail.code == NE_NOERROR) {
/* Compute the residual norm ||A*x - lambda*x||. */
k = 0;
for (j = 0; j <= nconv - 1; ++j) {
/* ax <- AV_k, where V_k is k-th Ritz vector. */
/* Compute matrix-vector multiply using nag_blast_dgbmv (f16pbc). */
nag_blast_dgbmv(Nag_ColMajor, Nag_NoTrans, n, n, kl, ku, 1.0, &ab[kl],
ldab, &v[k], 1, 0.0, ax, 1, &fail);
/* ax <- ax - (lambda_j) * V_k. */
alpha = -eigv[j];
/* Compute vector update using nag_blast_daxpby (f16ecc). */
nag_blast_daxpby(n, alpha, &v[k], 1, 1.0, ax, 1, &fail);
/* Compute 2-norm of Ritz estimates using nag_blast_dge_norm (f16rac). */
nag_blast_dge_norm(Nag_ColMajor, Nag_FrobeniusNorm, n, 1, ax, n, &resr,
&fail);
/* Scale. */
eigest[j] = resr / fabs(eigv[j]);
k = k + n;
}
/* Print computed eigenvalues. */
printf("\n The %4" NAG_IFMT " Ritz values are:\n\n", nconv);
for (j = 0; j <= nconv - 1; ++j) {
if (eigest[j] <= 1.0e-10) {
printf("%8" NAG_IFMT "%5s%7.4f\n", j + 1, "", eigv[j]);
} else {
printf("%8" NAG_IFMT "%5s%7.4f%5s%18.16f\n", j + 1, "", eigv[j],
" *** ", eigest[j]);
}
}
} else {
printf(" Error from "
"nag_sparseig_real_symm_band_solve (f12fgc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
END:
NAG_FREE(ab);
NAG_FREE(ax);
NAG_FREE(comm);
NAG_FREE(eigv);
NAG_FREE(eigest);
NAG_FREE(mb);
NAG_FREE(resid);
NAG_FREE(v);
NAG_FREE(icomm);
return exit_status;
}