/* nag_real_symm_banded_sparse_eigensystem_sol (f12fgc) Example Program. * * Copyright 2005 Numerical Algorithms Group. * * Mark 8, 2005. */ #include #include #include #include #include #include int main(void) { /* Constants */ Integer licomm = 140; /* Scalars */ double alpha, h2, resr, sigma = 0.0; Integer exit_status, i, j, k, kl, ku, ksub, ksup, lcomm; Integer ldab, n, nconv, ncv, nev, nx; /* Nag types */ NagError fail; /* Arrays */ double *ab = 0, *ax = 0, *comm = 0, *eigv = 0, *eigest = 0, *mb = 0; double *resid = 0, *v = 0; Integer *icomm = 0; exit_status = 0; INIT_FAIL(fail); printf("nag_real_symm_banded_sparse_eigensystem_sol (f12fgc)" " Example Program Results\n\n"); /* Skip heading in data file */ scanf("%*[^\n] "); scanf("%ld%ld%ld%*[^\n] ", &nx, &nev, &ncv); n = nx * nx; kl = nx; ku = nx; ldab = 2*kl + ku + 1; lcomm = 3*n + ncv*ncv + 8*ncv + 60; /* Allocate memory */ if (!(ab = NAG_ALLOC(ldab * n, double)) || !(ax = NAG_ALLOC(n, double)) || !(comm = NAG_ALLOC(lcomm, double)) || !(eigv = NAG_ALLOC(ncv, double)) || !(eigest = NAG_ALLOC(ncv, double)) || !(mb = NAG_ALLOC(1, double)) || !(resid = NAG_ALLOC(n, double)) || !(v = NAG_ALLOC(n * ncv, double)) || !(icomm = NAG_ALLOC(licomm, Integer))) { printf("Allocation failure\n"); exit_status = -1; goto END; } /* Initialise communication arrays for problem using nag_real_symm_banded_sparse_eigensystem_init (f12ffc). */ nag_real_symm_banded_sparse_eigensystem_init(n, nev, ncv, icomm, licomm, comm, lcomm, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_real_symm_banded_sparse_eigensystem_init" " (f12ffc).\n%s\n", fail.message); exit_status = 1; 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. */ for (j = 0; j < n*ldab; ++j) { ab[j] = 0.0; } /* Main diagonal of A. */ h2 = 1.0 / ((nx + 1) * (nx + 1)); k = kl + ku; for (j = 0; j < n; ++j) { ab[k] = 4.0 / h2; k = k + ldab; } /* First subdiagonal and superdiagonal of A. */ ksup = kl + ku - 1; ksub = kl + ku + 1; for (i = 0; i < nx; ++i) { ksup = ksup + ldab; for (j = 0; j < nx-1; ++j) { ab[ksub] = -1.0 / h2; ab[ksup] = -1.0 / h2; ksup = ksup + ldab; ksub = ksub + ldab; } ksub = ksub + ldab; } /* kl-th subdiagonal and ku-th super-diagonal. */ ksup = kl + nx*ldab; ksub = 2*kl + ku; for (i = 0; i < nx-1; ++i) { for (j = 0; j < nx; ++j) { ab[ksup] = -1.0 / h2; ab[ksub] = -1.0 / h2; ksup = ksup + ldab; ksub = ksub + ldab; } } /* Find eigenvalues of largest magnitude and the corresponding * eigenvectors using nag_real_symm_banded_sparse_eigensystem_sol (f12fgc). */ nag_real_symm_banded_sparse_eigensystem_sol(kl, ku, ab, mb, sigma, &nconv, eigv, v, resid, v, comm, icomm, &fail); if (fail.code == NE_NOERROR) { /* Compute the residual norm ||A*x - lambda*x||. */ k = 0; for (j = 0; j <= nconv-1; ++j) { /* ax <- AV_k, where V_k is kth Ritz vector. */ /* Compute matrix-vector multiply using nag_dgbmv (f16pbc). */ nag_dgbmv(Nag_ColMajor, Nag_NoTrans, n, n, kl, ku, 1.0, &ab[kl], ldab, &v[k], 1, 0.0, ax, 1, &fail); /* ax <- ax - (lambda_j) * V_k. */ alpha = -eigv[j]; /* Compute vector update using nag_daxpby (f16ecc). */ nag_daxpby(n, alpha, &v[k], 1, 1.0, ax, 1, &fail); /* Compute 2-norm of Ritz estimates using nag_dge_norm (f16rac).*/ nag_dge_norm(Nag_ColMajor, Nag_FrobeniusNorm, n, 1, ax, n, &resr, &fail); /* Scale. */ eigest[j] = resr / fabs(eigv[j]); k = k + n; } /* Print computed eigenvalues. */ printf("\n The %4ld Ritz values are:\n\n", nconv); for (j = 0; j <= nconv-1; ++j) { if (eigest[j] <= 1.0e-10) { printf("%8ld%5s%7.4f\n", j+1, "", eigv[j]); } else { printf("%8ld%5s%7.4f%5s%18.16f\n", j+1, "", eigv[j], " *** ", eigest[j]); } } } else { printf(" Error from " "nag_real_symm_banded_sparse_eigensystem_sol (f12fgc).\n%s\n", fail.message); exit_status = 1; goto END; } END: if (ab) NAG_FREE(ab); if (ax) NAG_FREE(ax); if (comm) NAG_FREE(comm); if (eigv) NAG_FREE(eigv); if (eigest) NAG_FREE(eigest); if (mb) NAG_FREE(mb); if (resid) NAG_FREE(resid); if (v) NAG_FREE(v); if (icomm) NAG_FREE(icomm); return exit_status; }