/* nag_sparseig_complex_band_solve (f12auc) Example Program.
*
* Copyright 2024 Numerical Algorithms Group.
*
* Mark 30.3, 2024.
*/
#include <math.h>
#include <nag.h>
#define AB(I, J) ab[(J - 1) * pdab + kl + ku + I - J]
#define MB(I, J) mb[(J - 1) * pdab + kl + ku + I - J]
int main(void) {
Integer exit_status = 0;
Complex one = {1.0, 0.0};
Complex minusone = {-1.0, 0.0};
Complex zero = {0.0, 0.0};
Nag_Boolean print_res = Nag_FALSE;
/* Scalars */
Complex ch_sub, ch_sup, alpha, sigma;
double h, rho, anorm;
Integer j, kl, ku, lcomm, licomm, n, ncol, nconv, ncv, nev, nx;
Integer pdab, pdmb, pdv;
/* Arrays */
Complex commd[1];
Integer icommd[1];
Complex *ab = 0, *ax = 0, *comm = 0, *d = 0, *mb = 0;
Complex *mx = 0, *resid = 0, *v = 0;
double *d_print = 0;
Integer *icomm = 0;
/* NAG types */
NagError fail;
Nag_OrderType order = Nag_ColMajor;
Nag_MatrixType matrix = Nag_GeneralMatrix;
Nag_DiagType diag = Nag_NonUnitDiag;
Nag_TransType trans = Nag_NoTrans;
INIT_FAIL(fail);
printf("nag_sparseig_complex_band_solve (f12auc) Example Program "
"Results\n\n");
fflush(stdout);
/* Skip heading in data file */
scanf("%*[^\n] ");
/* Read mesh size, number of eigenvalues wanted and size of subspace. */
scanf("%" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT "%*[^\n] ", &nx, &nev, &ncv);
/* Read complex value close to which eigenvalues are sought. */
scanf(" ( %lf , %lf ) %*[^\n] ", &sigma.re, &sigma.im);
n = nx * nx;
/* Initialize for the solution of a complex banded eigenproblem using
* nag_sparseig_complex_band_init (f12atc).
* But first get the required array lengths using arrays of length 1
* to store required lengths.
*/
licomm = -1;
lcomm = -1;
nag_sparseig_complex_band_init(n, nev, ncv, icommd, licomm, commd, lcomm,
&fail);
licomm = icommd[0];
lcomm = (Integer)(commd[0].re);
if (!(icomm = NAG_ALLOC((MAX(1, licomm)), Integer)) ||
!(comm = NAG_ALLOC((MAX(1, lcomm)), Complex))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
nag_sparseig_complex_band_init(n, nev, ncv, icomm, licomm, comm, lcomm,
&fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_sparseig_complex_init (f12anc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
/* Set options to show that the problem is a generalized eigenproblem and
* that eigenvalues are to be computed using shifted inverses using
* nag_sparseig_complex_option (f12arc).
*/
nag_sparseig_complex_option("SHIFTED INVERSE", icomm, comm, &fail);
nag_sparseig_complex_option("GENERALIZED", icomm, comm, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_sparseig_complex_option (f12arc).\n%s\n",
fail.message);
exit_status = 2;
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.
*/
kl = nx;
ku = nx;
pdab = 2 * kl + ku + 1;
pdmb = pdab;
if (!(ab = NAG_ALLOC(pdab * n, Complex)) ||
!(mb = NAG_ALLOC(pdmb * n, Complex))) {
printf("Allocation failure\n");
exit_status = -2;
goto END;
}
/* Zero out AB and MB. */
for (j = 0; j < pdab * n; j++) {
ab[j] = zero;
mb[j] = zero;
}
rho = 100.0;
h = 1.0 / (double)(nx + 1);
/* assign complex values using nag_complex_create (a02bac). */
ch_sub = nag_complex_create(0.5 * h * rho - 1.0, 0.0);
ch_sup = nag_complex_create(-0.5 * h * rho - 1.0, 0.0);
/* Set main diagonal of A and B to (4,0);
* first subdiagonal and superdiagonal of A are (-1 +/- rho/2h,0)..
* Outermost subdiagonal and superdiagonal are (-1,0).
*/
for (j = 1; j <= n; j++) {
AB(j, j) = nag_complex_create(4.0, 0.0);
MB(j, j) = nag_complex_create(4.0, 0.0);
if ((j - 1) % nx != 0 && j != n) {
AB(j + 1, j) = ch_sub;
}
if (j % nx != 0 && j != 1) {
AB(j - 1, j) = ch_sup;
}
if (j <= n - kl)
AB(j + kl, j) = minusone;
if (j > ku)
AB(j - ku, j) = minusone;
}
for (j = 2; j < n; j++) {
MB(j + 1, j) = one;
}
for (j = 2; j < n; j++) {
MB(j - 1, j) = one;
}
/* Allocate v to hold eigenvectors and d to hold eigenvalues. */
pdv = n;
if (!(d = NAG_ALLOC((nev), Complex)) ||
!(v = NAG_ALLOC((pdv) * (ncv), Complex)) ||
!(resid = NAG_ALLOC((n), Complex))) {
printf("Allocation failure\n");
exit_status = -3;
goto END;
}
/* Find eigenvalues closest in value to sigma and corresponding
* eigenvectors using nag_sparseig_complex_band_solve (f12auc).
*/
nag_sparseig_complex_band_solve(kl, ku, ab, mb, sigma, &nconv, d, v, resid, v,
comm, icomm, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_sparseig_complex_band_solve (f12auc).\n%s\n",
fail.message);
exit_status = 3;
goto END;
}
/* d_print stores eigenvalues (and possibly residuals) as real array. */
if (!(d_print = NAG_ALLOC(nconv * 3, double))) {
printf("Allocation failure\n");
exit_status = -4;
goto END;
}
for (j = 0; j < nconv; j++) {
d_print[j] = (d[j]).re;
d_print[j + nconv] = (d[j]).im;
}
/* By default this example will not calculate and print residuals.
* It will, however if the following line is uncommented.
*/
/* print_res = Nag_TRUE; */
if (print_res) {
ncol = 3;
/* Compute the residual norm ||A*x - lambda*B*x||.
* Matrix vector products Ax and Bx, for complex banded A and B, are formed
* using nag_blast_zgbmv (f16sbc).
* The scaled vector subtraction is performed using nag_blast_zaxpby
* (f16gcc).
*/
if (!(ax = NAG_ALLOC(n, Complex)) || !(mx = NAG_ALLOC(n, Complex))) {
printf("Allocation failure\n");
exit_status = -4;
goto END;
}
for (j = 0; j < nconv; j++) {
nag_blast_zgbmv(order, trans, n, n, kl, ku, one, &ab[kl], pdab,
&v[pdv * j], 1, zero, ax, 1, &fail);
if (fail.code == NE_NOERROR) {
nag_blast_zgbmv(order, trans, n, n, kl, ku, one, &mb[kl], pdmb,
&v[pdv * j], 1, zero, mx, 1, &fail);
if (fail.code == NE_NOERROR) {
/* alpha holds -d[j] using nag_complex_negate (a02cec). */
alpha = nag_complex_negate(d[j]);
nag_blast_zaxpby(n, alpha, mx, 1, one, ax, 1, &fail);
}
}
if (fail.code != NE_NOERROR) {
printf("Error in calculating residual norm.\n%s\n", fail.message);
exit_status = 4;
goto END;
}
/* Get ||ax|| = ||Ax - lambda Bx|| using nag_blast_zge_norm (f16uac)
* with the Frobenius norm on a complex general matrix with one column.
*/
nag_blast_zge_norm(order, Nag_FrobeniusNorm, n, 1, ax, n, &anorm, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_blast_zge_norm (f16uac).\n%s\n", fail.message);
exit_status = 5;
goto END;
}
/* nag_complex_abs (a02dbc) returns the modulus of the eigenvalue */
d_print[j + 2 * nconv] = anorm / nag_complex_abs(d[j]);
}
} else {
ncol = 2;
}
printf("\n");
/* Print the eigenvalues (and possibly the residuals) using the matrix
* printer, nag_file_print_matrix_real_gen (x04cac).
*/
printf(" Number of eigenvalues wanted = %" NAG_IFMT "\n", nev);
printf(" Number of eigenvalues converged = %" NAG_IFMT "\n", nconv);
printf(" Eigenvalues are closest to Sigma = ( %lf , %lf )\n\n", sigma.re,
sigma.im);
fflush(stdout);
nag_file_print_matrix_real_gen(order, matrix, diag, nconv, ncol, d_print,
nconv, " Ritz values closest to sigma", 0,
&fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_file_print_matrix_real_gen (x04cac).\n%s\n",
fail.message);
exit_status = 6;
}
END:
NAG_FREE(ab);
NAG_FREE(ax);
NAG_FREE(comm);
NAG_FREE(d);
NAG_FREE(mb);
NAG_FREE(mx);
NAG_FREE(resid);
NAG_FREE(v);
NAG_FREE(d_print);
NAG_FREE(icomm);
return exit_status;
}