/* nag_eigen_real_gen_sparse_arnoldi (f02ekc) Example Program.
*
* NAGPRODCODE Version.
*
* Copyright 2016 Numerical Algorithms Group.
*
* Mark 26, 2016.
*/
#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <naga02.h>
#include <nagf02.h>
#include <nagf12.h>
#include <nagx02.h>
#include <nagx04.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 Complex w[], const double rzest[],
Integer *istat, Nag_Comm *comm);
#ifdef __cplusplus
}
#endif
int main(void)
{
/* Scalars */
double one = 1.0, two = 2.0, three = 3.0;
double h, rho, s, sigma;
Integer exit_status = 0;
Integer fileid, fmode, i, imon, k, maxit, mode;
Integer n, nconv, ncv, nev, nnz, nx, prtlvl, tdv;
/* Local Arrays */
Complex *w = 0;
double *a = 0, *resid = 0, *v = 0;
double user[1];
Integer *icolzp = 0, *irowix = 0;
Integer iuser[5];
const char *filename = "f02ekce.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_gen_sparse_arnoldi (f02ekc) ");
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]%lf%*[^\n]", &rho, &sigma);
n = nx * nx;
nnz = 3 * n - 2;
tdv = n;
if (!(resid = NAG_ALLOC((ncv), double)) ||
!(a = NAG_ALLOC((nnz), double)) ||
!(icolzp = NAG_ALLOC((n + 1), Integer)) ||
!(irowix = NAG_ALLOC((nnz), Integer)) ||
!(w = NAG_ALLOC((ncv), Complex)) ||
!(v = NAG_ALLOC((tdv) * (ncv), double))
)
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Construct A in compressed column storage (CCS) format where:
* A_{i,i} = 2 + i
* A_{i+1,i) = 3
* A_{i,i+1} = rho/(2n+2) - 1
*/
h = one / (double) (n + 1);
s = rho * h / two - one;
a[0] = two + one;
a[1] = three;
icolzp[0] = 1;
irowix[0] = 1;
irowix[1] = 2;
k = 3;
for (i = 2; i <= n - 1; i++) {
icolzp[i - 1] = k;
irowix[k - 1] = i - 1;
irowix[k] = i;
irowix[k + 1] = i + 1;
a[k - 1] = s;
a[k] = two + (double) (i);
a[k + 1] = three;
k = k + 3;
}
icolzp[n - 1] = k;
icolzp[n] = k + 2;
irowix[k - 1] = n - 1;
irowix[k] = n;
a[k - 1] = s;
a[k] = two + (double) (n);
/* 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_open_file (x04acc):
* open unit number for reading, writing or appending, and
* associate unit with named file
*/
/* If prtlvl >=10 internal monitoring information in addition to whatever is
* written to fileid using mymonit.
*/
fmode = 1;
nag_open_file(filename, fmode, &fileid, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_open_file (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_gen_sparse_arnoldi (f02ekc):
* selected eigenvalues of real general matrix driver.
*/
nag_eigen_real_gen_sparse_arnoldi(n, nnz, a, icolzp, irowix, nev, ncv,
sigma, mymonit, myoption, &nconv, w, v,
tdv, resid, &comm, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_eigen_real_gen_sparse_arnoldi (f02ekc)\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
if (imon > 0) {
/* Close the monitoring file using
* nag_close_file (x04adc):
* close file associated with given unit number
*/
nag_close_file(fileid, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_close_file (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 = 1; i <= nconv; i++) {
/* nag_machine_precision (x02ajc) */
if (resid[i - 1] > (double) (100 * n) * nag_machine_precision) {
printf("%7" NAG_IFMT " ( %13.5e, %13.5e) ", i, w[i - 1].re,
w[i - 1].im);
printf("%13.5e\n", resid[i - 1]);
}
else {
printf("%8" NAG_IFMT " ( %13.5e, %13.5e) \n", i, w[i - 1].re,
w[i - 1].im);
}
}
END:
NAG_FREE(w);
NAG_FREE(a);
NAG_FREE(v);
NAG_FREE(resid);
NAG_FREE(icolzp);
NAG_FREE(irowix);
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);
if (comm->iuser[0] > 0) {
sprintf(rec, "Print Level=%5" NAG_IFMT, comm->iuser[0]);
fail1.code = 1;
/* Set print level using
* nag_real_sparse_eigensystem_option (f12adc)
* Set a single option from a string.
*/
nag_real_sparse_eigensystem_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;
/* Set iteration limit using
* nag_real_sparse_eigensystem_option (f12adc)
* Set a single option from a string.
*/
nag_real_sparse_eigensystem_option(rec, icomm, com, &fail1);
*istat = MAX(*istat, fail1.code);
}
if (comm->iuser[2] > 0) {
fail1.code = 1;
/* Set computational mode to shifted inverse real. */
nag_real_sparse_eigensystem_option("Shifted Inverse Real", icomm, com,
&fail1);
*istat = MAX(*istat, fail1.code);
}
if (comm->iuser[3] > 0) {
fail1.code = 1;
/* Switch monitoring on and use the fileid stored in iuser[4]. */
sprintf(rec, "Monitoring=%5" NAG_IFMT, comm->iuser[4]);
nag_real_sparse_eigensystem_option(rec, icomm, com, &fail1);
*istat = MAX(*istat, fail1.code);
}
}
static void NAG_CALL mymonit(Integer ncv, Integer niter, Integer nconv,
const Complex 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_write_line (x04bac):
* write formatted record to external file.
*/
if (niter == 1 && comm->iuser[2] > 0) {
sprintf(line, " Arnoldi basis vectors used: %4" NAG_IFMT "\n", ncv);
nag_write_line(comm->iuser[4], line);
sprintf(line, " The following Ritz values (mu) are related to the\n");
nag_write_line(comm->iuser[4], line);
sprintf(line, " true eigenvalues (lambda) by lambda = sigma + 1/mu\n");
nag_write_line(comm->iuser[4], line);
}
sprintf(line, "\n Iteration number %4" NAG_IFMT "\n", niter);
nag_write_line(comm->iuser[4], line);
sprintf(line,
" Ritz values converged so far (%4" NAG_IFMT ") and their Ritz "
"estimates:\n", nconv);
nag_write_line(comm->iuser[4], line);
for (i = 1; i <= nconv; i++) {
sprintf(line,
" %4" NAG_IFMT " (%13.5e,%13.5e) %13.5e\n",
i, w[i - 1].re, w[i - 1].im, rzest[i - 1]);
nag_write_line(comm->iuser[4], line);
}
sprintf(line, " Next (unconverged) Ritz value:\n");
nag_write_line(comm->iuser[4], line);
sprintf(line,
" %4" NAG_IFMT " (%13.5e,%13.5e)\n",
nconv + 1, w[nconv].re, w[nconv].im);
nag_write_line(comm->iuser[4], line);
}
*istat = 0;
}