/* nag_complex_banded_eigensystem_solve (f12auc) Example Program. * * Copyright 2013 Numerical Algorithms Group. * * Mark 24, 2013. */ #include #include #include #include #include #include #include #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"); /* 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; jku) AB(j-ku,j) = minusone; } for (j=2; j