/* nag_complex_banded_eigensystem_solve (f12auc) Example Program.
*
* Copyright 2014 Numerical Algorithms Group.
*
* Mark 24, 2013.
*/
#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <naga02.h>
#include <nagf12.h>
#include <nagf16.h>
#include <nagx04.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_complex_banded_eigensystem_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("%ld%ld%ld%*[^\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_complex_banded_eigensystem_init (f12atc).
* But first get the required array lengths using arrays of length 1
* to store required lengths.
*/
licomm = -1;
lcomm = -1;
nag_complex_banded_eigensystem_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_complex_banded_eigensystem_init(n, nev, ncv, icomm, licomm, comm, lcomm,
&fail);
if (fail.code != NE_NOERROR)
{
printf("Error from nag_complex_sparse_eigensystem_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_complex_sparse_eigensystem_option (f12arc).
*/
nag_complex_sparse_eigensystem_option("SHIFTED INVERSE", icomm, comm, &fail);
nag_complex_sparse_eigensystem_option("GENERALIZED", icomm, comm, &fail);
if (fail.code != NE_NOERROR)
{
printf("Error from nag_complex_sparse_eigensystem_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 (a02bac). */
ch_sub = nag_complex(0.5 * h * rho - 1.0, 0.0);
ch_sup = nag_complex(-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 super-diagonal are (-1,0).
*/
for ( j=1; j<=n; j++)
{
AB(j,j) = nag_complex(4.0, 0.0);
MB(j,j) = nag_complex(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_complex_banded_eigensystem_solve (f12auc).
*/
nag_complex_banded_eigensystem_solve(kl, ku, ab, mb, sigma, &nconv, d,
v, resid, v, comm, icomm, &fail);
if (fail.code != NE_NOERROR)
{
printf("Error from nag_complex_banded_eigensystem_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_zgbmv (f16sbc).
* The scaled vector subtraction is performed using nag_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_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_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_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_zge_norm (f16uac)
* with the Frobenius norm on a complex general matrix with one column.
*/
nag_zge_norm(order, Nag_FrobeniusNorm, n, 1, ax, n, &anorm, &fail);
if (fail.code != NE_NOERROR)
{
printf("Error from nag_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_gen_real_mat_print (x04cac).
*/
printf(" Number of eigenvalues wanted = %ld\n",nev);
printf(" Number of eigenvalues converged = %ld\n",nconv);
printf(" Eigenvalues are closest to Sigma = ( %lf , %lf )\n\n",
sigma.re, sigma.im);
fflush(stdout);
nag_gen_real_mat_print(order, matrix, diag, nconv, ncol, d_print, nconv,
" Ritz values closest to sigma", 0, &fail);
if (fail.code != NE_NOERROR)
{
printf("Error from nag_gen_real_mat_print (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;
}