/* nag_eigen_real_symm_sparse_arnoldi (f02fkc) Example Program.
*
* Copyright 2024 Numerical Algorithms Group.
*
* Mark 30.2, 2024.
*/
#include <math.h>
#include <nag.h>
/* User-defined Functions */
#ifdef __cplusplus
extern "C" {
#endif
static void NAG_CALL myoption(Integer icomm[], double com[], Integer *istat,
Nag_Comm *comm);
static void NAG_CALL mymonit(Integer ncv, Integer niter, Integer nconv,
const double w[], const double rzest[],
Integer *istat, Nag_Comm *comm);
#ifdef __cplusplus
}
#endif
int main(void) {
/* Scalars */
double h2, sigma;
Integer exit_status = 0;
Integer fileid, fmode, i, imon, j, k, lo, maxit, mode;
Integer n, nconv, ncv, nev, nnz, nx, prtlvl, tdv;
/* Local Arrays */
double *w = 0, *a = 0, *resid = 0, *v = 0;
double user[1];
Integer *icol = 0, *irow = 0;
Integer iuser[5];
const char *filename = "f02fkce.monit";
/* Nag Types */
Nag_Comm comm;
NagError fail;
INIT_FAIL(fail);
comm.user = user;
comm.iuser = iuser;
user[0] = 0.0;
iuser[0] = 0;
/* Output preamble */
printf(" nag_eigen_real_symm_sparse_arnoldi (f02fkc) ");
printf("Example Program Results\n\n");
fflush(stdout);
/* Skip heading in data file */
scanf("%*[^\n] ");
/* Read in problem size and parameters */
scanf("%" NAG_IFMT "%*[^\n]%" NAG_IFMT "%*[^\n]%" NAG_IFMT "", &nx, &nev,
&ncv);
scanf("%*[^\n]%lf%*[^\n]", &sigma);
n = nx * nx;
nnz = 3 * n - 2 * nx;
tdv = n;
if (!(resid = NAG_ALLOC((ncv), double)) || !(a = NAG_ALLOC((nnz), double)) ||
!(icol = NAG_ALLOC((nnz), Integer)) ||
!(irow = NAG_ALLOC((nnz), Integer)) || !(w = NAG_ALLOC((ncv), double)) ||
!(v = NAG_ALLOC((tdv) * (ncv), double))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Construct A in sparse (SCS) format where:
* A_{i,i} = 4/(h*h)
* A_{i+1,i) = -1/(h*h)
* A_{i+nx,i} = -1/(h*h)
*/
h2 = 1.0 / (double)((nx + 1) * (nx + 1));
/* Main Diagonal of A */
k = 0;
for (i = 1; i <= n; i++) {
irow[k] = i;
icol[k] = i;
a[k] = 4.0 / h2;
k++;
}
/* First subdiagonal of A. */
for (i = 1; i <= nx; i++) {
lo = (i - 1) * nx;
for (j = lo + 1; j <= lo + nx - 1; j++) {
irow[k] = j + 1;
icol[k] = j;
a[k] = -1.0 / h2;
k++;
}
}
/* nx-th subdiagonal of A. */
for (i = 1; i < nx; i++) {
lo = (i - 1) * nx;
for (j = lo + 1; j <= lo + nx; j++) {
irow[k] = j + nx;
icol[k] = j;
a[k] = -1.0 / h2;
k++;
}
}
/* Set some options via iuser array and routine argument OPTION.
* iuser[0] = print level, iuser[1] = iteration limit,
* iuser[2]>0 means shifted-invert mode
* iuser[3]>0 means print monitoring info.
*/
scanf("%" NAG_IFMT "%*[^\n]%" NAG_IFMT "%*[^\n]", &prtlvl, &maxit);
scanf("%" NAG_IFMT "%*[^\n]%" NAG_IFMT "%*[^\n]", &mode, &imon);
if (imon > 0) {
/* Open the monitoring file for writing using
* nag_file_open (x04acc).
* If prtlvl >=10 internal monitoring information is also written.
*/
fmode = 1;
nag_file_open(filename, fmode, &fileid, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_file_open (x04acc) %s\n", fail.message);
exit_status = 1;
goto END;
}
iuser[4] = fileid;
}
iuser[0] = prtlvl;
iuser[1] = maxit;
iuser[2] = mode;
iuser[3] = imon;
/* Compute eigenvalues and eigenvectors using
* nag_eigen_real_symm_sparse_arnoldi (f02fkc).
* selected eigenvalues of real general matrix (driver).
*/
nag_eigen_real_symm_sparse_arnoldi(n, nnz, a, irow, icol, nev, ncv, sigma,
mymonit, myoption, &nconv, w, v, tdv,
resid, &comm, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_eigen_real_symm_sparse_arnoldi (f02fkc)\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
if (imon > 0) {
/* Close the monitoring file using
* nag_file_close (x04adc).
*/
nag_file_close(fileid, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_file_close (x04adc) %s\n", fail.message);
exit_status = 1;
goto END;
}
}
printf(" The %4" NAG_IFMT " ", nconv);
printf(" Ritz values of closest to %13.5e are \n", sigma);
for (i = 0; i < nconv; i++) {
/* nag_machine_precision (x02ajc) */
if (resid[i] > (double)(100 * n) * nag_machine_precision) {
printf("%7" NAG_IFMT " %13.5e %13.5e\n", i + 1, w[i], resid[i]);
} else {
printf("%8" NAG_IFMT " %13.5e\n", i + 1, w[i]);
}
}
END:
NAG_FREE(w);
NAG_FREE(a);
NAG_FREE(v);
NAG_FREE(resid);
NAG_FREE(icol);
NAG_FREE(irow);
return exit_status;
}
static void NAG_CALL myoption(Integer icomm[], double com[], Integer *istat,
Nag_Comm *comm) {
NagError fail1;
char rec[26];
INIT_FAIL(fail1);
/* Set options using
* nag_sparseig_real_symm_option (f12fdc).
*/
if (comm->iuser[0] > 0) {
sprintf(rec, "Print Level=%5" NAG_IFMT, comm->iuser[0]);
fail1.code = 1;
nag_sparseig_real_symm_option(rec, icomm, com, &fail1);
*istat = MAX(*istat, fail1.code);
}
if (comm->iuser[1] > 100) {
sprintf(rec, "Iteration Limit=%5" NAG_IFMT, comm->iuser[1]);
fail1.code = 1;
nag_sparseig_real_symm_option(rec, icomm, com, &fail1);
*istat = MAX(*istat, fail1.code);
}
if (comm->iuser[2] > 0) {
fail1.code = 1;
nag_sparseig_real_symm_option("Shifted Inverse", icomm, com, &fail1);
*istat = MAX(*istat, fail1.code);
}
if (comm->iuser[3] > 0) {
fail1.code = 1;
sprintf(rec, "Monitoring=%5" NAG_IFMT, comm->iuser[4]);
nag_sparseig_real_symm_option(rec, icomm, com, &fail1);
*istat = MAX(*istat, fail1.code);
}
}
static void NAG_CALL mymonit(Integer ncv, Integer niter, Integer nconv,
const double w[], const double rzest[],
Integer *istat, Nag_Comm *comm) {
Integer i;
char line[100];
if (comm->iuser[3] > 0) {
/* Write lines to the file we opened for monitoring using
* nag_file_line_write (x04bac).
*/
if (niter == 1 && comm->iuser[2] > 0) {
sprintf(line, " Arnoldi basis vectors used: %4" NAG_IFMT "\n", ncv);
nag_file_line_write(comm->iuser[4], line);
sprintf(line, " The following Ritz values (mu) are related to the\n");
nag_file_line_write(comm->iuser[4], line);
sprintf(line, " true eigenvalues (lambda) by lambda = sigma + 1/mu\n");
nag_file_line_write(comm->iuser[4], line);
}
sprintf(line, "\n Iteration number %4" NAG_IFMT "\n", niter);
nag_file_line_write(comm->iuser[4], line);
sprintf(line,
" Ritz values converged so far (%4" NAG_IFMT ") and their Ritz "
"estimates:\n",
nconv);
nag_file_line_write(comm->iuser[4], line);
for (i = 0; i < nconv; i++) {
sprintf(line, " %4" NAG_IFMT " %13.5e %13.5e\n", i + 1, w[i], rzest[i]);
nag_file_line_write(comm->iuser[4], line);
}
sprintf(line, " Next (unconverged) Ritz value:\n");
nag_file_line_write(comm->iuser[4], line);
sprintf(line, " %4" NAG_IFMT " %13.5e\n", nconv + 1, w[nconv]);
nag_file_line_write(comm->iuser[4], line);
}
*istat = 0;
}