/* nag_sparseig_real_band_solve (f12agc) Example Program.
*
* Copyright 2024 Numerical Algorithms Group.
*
* Mark 30.2, 2024.
*/
#include <math.h>
#include <nag.h>
#include <stdio.h>
int main(void) {
/* Constants */
double rho = 100.0;
/* Scalars */
Complex res, eig;
double alpha, h, resr, resi, sigmai, sigmar;
Integer exit_status, i, j, k, kl, ksub, ksup, ku, lcomm, licomm;
Integer ldab, n, nconv, ncv, nev, nx;
/* Nag types */
Nag_Boolean first;
NagError fail;
/* Arrays */
double *ab = 0, *ax = 0, *comm = 0, *eigvr = 0, *eigvi = 0, *eigest = 0;
double *mb = 0, *mx = 0, *resid = 0, *v = 0;
Integer *icomm = 0;
exit_status = 0;
INIT_FAIL(fail);
printf("nag_sparseig_real_band_solve (f12agc) Example "
"Program Results\n");
/* Skip heading in data file */
scanf("%*[^\n] ");
scanf("%" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT "%lf%lf%*[^\n] ", &nx, &nev,
&ncv, &sigmar, &sigmai);
n = nx * nx;
/* ku, kl are number of superdiagonals and subdiagonals within the */
/* band of matrices A and M. */
kl = nx;
ku = nx;
/* ldab is the width of the band required to store the */
/* factorization of the matrix A. */
ldab = 2 * kl + ku + 1;
/* Allocate memory. */
if (!(ab = NAG_ALLOC(ldab * n, double)) || !(ax = NAG_ALLOC(n, double)) ||
!(eigvr = NAG_ALLOC(ncv, double)) || !(eigvi = NAG_ALLOC(ncv, double)) ||
!(eigest = NAG_ALLOC(ncv, double)) ||
!(mb = NAG_ALLOC(ldab * n, double)) || !(mx = NAG_ALLOC(n, double)) ||
!(resid = NAG_ALLOC(n, double)) || !(v = NAG_ALLOC(ncv * n, double))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Initialize communication arrays for problem using
nag_sparseig_real_band_init (f12afc).
The first call sets lcomm = licomm = -1 to perform a workspace
query. */
lcomm = licomm = -1;
if (!(comm = NAG_ALLOC(1, double)) || !(icomm = NAG_ALLOC(1, Integer))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
nag_sparseig_real_band_init(n, nev, ncv, icomm, licomm, comm, lcomm, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_sparseig_real_band_init "
"(f12afc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
lcomm = (Integer)comm[0];
licomm = icomm[0];
NAG_FREE(comm);
NAG_FREE(icomm);
if (!(comm = NAG_ALLOC(lcomm, double)) ||
!(icomm = NAG_ALLOC(licomm, Integer))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
nag_sparseig_real_band_init(n, nev, ncv, icomm, licomm, comm, lcomm, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_sparseig_real_band_init "
"(f12afc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
/* Select the required spectrum using
nag_sparseig_real_option (f12adc). */
nag_sparseig_real_option("SHIFTED IMAGINARY", icomm, comm, &fail);
/* Select the problem type using
nag_sparseig_real_option (f12adc). */
nag_sparseig_real_option("GENERALIZED", icomm, comm, &fail);
/* Construct the matrix A in banded form and store in AB. */
/* Zero out matrices first */
for (j = 0; j < n * ldab; ++j) {
ab[j] = 0.0;
mb[j] = 0.0;
}
/* Main diagonal of A. */
k = kl + ku;
for (j = 0; j < n; ++j) {
ab[k] = 4.;
mb[k] = 4.;
k = k + ldab;
}
/* First subdiagonal and superdiagonal of A. */
ksup = kl + ku - 1;
ksub = kl + ku + 1;
h = .5 * rho / (double)(nx + 1);
for (i = 1; i <= nx; ++i) {
ksub = ksub + ldab;
for (j = 1; j <= nx - 1; ++j) {
ab[ksub] = -1. + h;
ab[ksup] = -1. - h;
ksup = ksup + ldab;
ksub = ksub + ldab;
}
ksup = ksup + ldab;
}
ksub = kl + ku + 1 + ldab;
ksup = kl + ku - 1;
for (j = 1; j <= n - 1; ++j) {
mb[ksup] = 1.;
mb[ksub] = 1.;
ksup = ksup + ldab;
ksub = ksub + ldab;
}
/* KL-th subdiagonal and KU-th superdiagonal. */
ksup = kl + nx * ldab;
ksub = 2 * kl + ku;
for (i = 1; i <= nx - 1; ++i) {
for (j = 1; j <= nx; ++j) {
ab[ksup] = -1.;
ab[ksub] = -1.;
ksup = ksup + ldab;
ksub = ksub + ldab;
}
}
/* Find eigenvalues closest in value to SIGMA and corresponding
eigenvectors using nag_sparseig_real_band_solve (f12agc). */
nag_sparseig_real_band_solve(kl, ku, ab, mb, sigmar, sigmai, &nconv, eigvr,
eigvi, v, resid, v, comm, icomm, &fail);
if (fail.code == NE_NOERROR) {
/* Compute the residual norm ||A*x - lambda*Mx||. */
first = Nag_TRUE;
k = 0;
for (j = 0; j <= nconv - 1; ++j) {
if (eigvi[j] == 0.) {
/* Eigenvalue is Real. */
/* 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., &ab[kl],
ldab, &v[k], 1, 0., ax, 1, &fail);
/* mx <- MV_k. */
nag_blast_dgbmv(Nag_ColMajor, Nag_NoTrans, n, n, kl, ku, 1., &mb[kl],
ldab, &v[k], 1, 0., mx, 1, &fail);
/* ax <- ax - (lambda_j) * mx. */
alpha = -eigvr[j];
/* Compute vector update using nag_blast_daxpby (f16ecc). */
nag_blast_daxpby(n, alpha, mx, 1, 1., ax, 1, &fail);
/* resr = norm(ax). */
/* 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(eigvr[j]);
} else if (first) {
/* First or a conjugate pair of eigenvalues. */
/* Real part of residual Re[Ax-lamdaMx]. */
/* ax <- AV_k, where V_k is real part of 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., &ab[kl],
ldab, &v[k], 1, 0., ax, 1, &fail);
/* mx <- MV_k */
nag_blast_dgbmv(Nag_ColMajor, Nag_NoTrans, n, n, kl, ku, 1., &mb[kl],
ldab, &v[k], 1, 0., mx, 1, &fail);
/* ax <- ax - (lambda_j).re * mx. */
/* Compute vector update using nag_blast_daxpby (f16ecc). */
alpha = -eigvr[j];
nag_blast_daxpby(n, alpha, mx, 1, 1., ax, 1, &fail);
/* mx <- MW_k where W_k is imaginary part k-th Ritz vector. */
nag_blast_dgbmv(Nag_ColMajor, Nag_NoTrans, n, n, kl, ku, 1., &mb[kl],
ldab, &v[k + n], 1, 0., mx, 1, &fail);
/* ax <- ax - (lambda_j).im * mx. */
alpha = eigvi[j];
nag_blast_daxpby(n, alpha, mx, 1, 1., ax, 1, &fail);
/* resr = norm(ax). */
/* 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);
/* Imaginary part of residual Im[Ax-lamdaMx]. */
/* ax <- AW_k. */
/* Compute matrix-vector multiply using nag_blast_dgbmv
(f16pbc). */
nag_blast_dgbmv(Nag_ColMajor, Nag_NoTrans, n, n, kl, ku, 1., &ab[kl],
ldab, &v[k + n], 1, 0., ax, 1, &fail);
/* ax <- ax - (lambda_j).re * mx. */
alpha = -eigvr[j];
/* Compute vector update using nag_blast_daxpby (f16ecc). */
nag_blast_daxpby(n, alpha, mx, 1, 1., ax, 1, &fail);
/* mx <- MV_k. */
nag_blast_dgbmv(Nag_ColMajor, Nag_NoTrans, n, n, kl, ku, 1., &mb[kl],
ldab, &v[k], 1, 0., mx, 1, &fail);
/* ax <- ax - (lambda_j).im * mx. */
alpha = -eigvi[j];
nag_blast_daxpby(n, alpha, mx, 1, 1., ax, 1, &fail);
/* resi = norm(ax). */
/* Compute 2-norm of Ritz estimates using nag_blast_dge_norm
(f16rac). */
nag_blast_dge_norm(Nag_ColMajor, Nag_FrobeniusNorm, n, 1, ax, n, &resi,
&fail);
/* Scale residual using Ritz value. */
/* Assign to Complex type using nag_complex_create (a02bac) */
res = nag_complex_create(resr, resi);
eig = nag_complex_create(eigvr[j], eigvi[j]);
eigest[j] = nag_complex_abs(res) / nag_complex_abs(eig);
/* Set residual for second in conjugate pair. */
eigest[j + 1] = eigest[j];
first = Nag_FALSE;
} else {
/* Second of complex conjugate pair; */
/* Already got residual from first in pair. */
first = Nag_TRUE;
}
k = k + n;
}
/* Print computed eigenvalues. */
printf("\n The %4" NAG_IFMT " generalized Ritz values", nconv);
printf(" closest to \n");
printf(" ( %7.4f, +-%7.4f ) are:\n\n", sigmar, sigmai);
for (j = 0; j <= nconv - 1; ++j) {
if (eigest[j] <= 1.0e-10) {
printf("%8" NAG_IFMT "%5s( %7.4f, %7.4f ).\n", j + 1, "", eigvr[j],
eigvi[j]);
} else {
printf("%8" NAG_IFMT "%5s( %7.4f, %7.4f )%5s%18.16f.\n", j + 1, "",
eigvr[j], eigvi[j], " *** ", eigest[j]);
}
}
} else {
printf("Error from "
"nag_sparseig_real_band_solve (f12agc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
END:
NAG_FREE(ab);
NAG_FREE(ax);
NAG_FREE(comm);
NAG_FREE(eigvr);
NAG_FREE(eigvi);
NAG_FREE(eigest);
NAG_FREE(mb);
NAG_FREE(mx);
NAG_FREE(resid);
NAG_FREE(v);
NAG_FREE(icomm);
return exit_status;
}