/* 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; 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); Vprintf("nag_real_symm_banded_sparse_eigensystem_sol (f12fgc) Example " "Program Results\n\n"); /* Skip heading in data file */ Vscanf("%*[^\n] "); Vscanf("%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)) ) { Vprintf("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); /* 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. */ /* Main diagonal of A. */ h2 = 1.0 / ((nx + 1) * (nx + 1)); k = kl + ku; for (j = 0; j <= n-1; ++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 = 1; i <= nx; ++i) { ksup = ksup + ldab; for (j = 1; 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 = 1; i <= nx - 1; ++i) { for (j = 1; 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. */ Vprintf("\n The %4ld Ritz values are:\n\n", nconv); for (j = 0; j <= nconv-1; ++j) { if (eigest[j] <= 1.0e-10) { Vprintf("%8ld%5s%7.4f\n", j+1, "", eigv[j]); } else { Vprintf("%8ld%5s%7.4f%5s%18.16f\n", j+1, "", eigv[j], " *** ",eigest[j]); } } } else { Vprintf(" 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; }