/* nag_real_banded_sparse_eigensystem_sol (f12agc) Example Program. * * Copyright 2005 Numerical Algorithms Group. * * Mark 8, 2005. */ #include #include #include #include #include #include #include int main(void) { /* Constants */ Integer licomm=140; double rho = 100.0; /* Scalars */ Complex res, eig; double alpha, h, resr, resi, sigmai, sigmar; Integer exit_status, i, j, k, kl, ksub, ksup, ku, lcomm; Integer ldab, n, nconv, ncv, nev, nx; /* Nag types */ Nag_Boolean first; NagError fail; /* Arrays */ double *ab=0, *ax=0, *comm=0, *eigvr=0, *eigvi=0, *eigest=0; double *mb=0, *mx=0, *resid=0, *v=0; Integer *icomm=0; exit_status = 0; INIT_FAIL(fail); Vprintf("nag_real_banded_sparse_eigensystem_sol (f12agc) Example Program" " Results\n"); /* Skip heading in data file */ Vscanf("%*[^\n] "); Vscanf("%ld%ld%ld%lf%lf%*[^\n] ", &nx, &nev, &ncv, &sigmar, &sigmai); n = nx * nx; /* ku, kl are number of superdiagonals and subdiagonals within the */ /* band of matrices A and M. */ kl = nx; ku = nx; /* ldab is the width of the band required to store the */ /* factorisationof the matrix A. */ ldab = 2*kl + ku + 1; lcomm = 60; /* Allocate memory. */ if ( !(ab = NAG_ALLOC(ldab * n, double)) || !(ax = NAG_ALLOC(n, double)) || !(comm = NAG_ALLOC(lcomm, double)) || !(eigvr = NAG_ALLOC(ncv, double)) || !(eigvi = NAG_ALLOC(ncv, double)) || !(eigest = NAG_ALLOC(ncv, double)) || !(mb = NAG_ALLOC(ldab * n, double)) || !(mx = NAG_ALLOC(n, double)) || !(resid = NAG_ALLOC(n, double)) || !(v = NAG_ALLOC(ncv * n, double)) || !(icomm = NAG_ALLOC(licomm, Integer)) ) { Vprintf("Allocation failure\n"); exit_status = -1; goto END; } /* Initialise communication arrays for problem using nag_real_banded_sparse_eigensystem_init (f12afc). */ nag_real_banded_sparse_eigensystem_init(n, nev, ncv, icomm, licomm, comm, lcomm, &fail); /* Select the required spectrum using nag_real_sparse_eigensystem_option (f12adc). */ nag_real_sparse_eigensystem_option("SHIFTED IMAGINARY", icomm, comm, &fail); /* Select the problem type using nag_real_sparse_eigensystem_option (f12adc). */ nag_real_sparse_eigensystem_option("GENERALIZED", icomm, comm, &fail); /* Construct the matrix A in banded form and store in AB. */ /* Zero out matrices first */ for (j = 0; j < n*ldab; ++j) { ab[j] = 0.0; mb[j] = 0.0; } /* Main diagonal of A. */ k = kl + ku; for (j = 0; j < n; ++j) { ab[k] = 4.; mb[k] = 4.; k = k + ldab; } /* First subdiagonal and superdiagonal of A. */ ksup = kl + ku - 1; ksub = kl + ku + 1; h = .5 * rho / (double) (nx + 1); for (i = 1; i <= nx; ++i) { ksub = ksub + ldab; for (j = 1; j <= nx - 1; ++j) { ab[ksub] = - 1. + h; ab[ksup] = - 1. - h; ksup = ksup + ldab; ksub = ksub + ldab; } ksup = ksup + ldab; } ksub = kl + ku + 1 + ldab; ksup = kl + ku - 1; for (j = 1; j <= n - 1; ++j) { mb[ksup] = 1.; mb[ksub] = 1.; ksup = ksup + 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.; ab[ksub] = - 1.; ksup = ksup + ldab; ksub = ksub + ldab; } } /* Find eigenvalues closest in value to SIGMA and corresponding eigenvectors using nag_real_banded_sparse_eigensystem_sol (f12agc). */ nag_real_banded_sparse_eigensystem_sol(kl, ku, ab, mb, sigmar, sigmai, &nconv, eigvr, eigvi, v, resid, v, comm, icomm, &fail); if (fail.code == NE_NOERROR) { /* Compute the residual norm ||A*x - lambda*Mx||. */ first = Nag_TRUE; k = 0; for (j = 0; j <= nconv-1; ++j) { if (eigvi[j] == 0.) { /* Eigenvalue is Real. */ /* 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., &ab[kl], ldab, &v[k], 1, 0., ax, 1, &fail); /* mx <- MV_k. */ nag_dgbmv(Nag_ColMajor, Nag_NoTrans, n, n, kl, ku, 1., &mb[kl], ldab, &v[k], 1, 0., mx, 1, &fail); /* ax <- ax - (lambda_j) * mx. */ alpha = -eigvr[j]; /* Compute vector update using nag_daxpby (f16ecc). */ nag_daxpby(n, alpha, mx, 1, 1., ax, 1, &fail); /* resr = norm(ax). */ /* 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(eigvr[j]); } else if (first) { /* First or a conjugate pair of eigenvalues. */ /* Real part of residual Re[Ax-lamdaMx]. */ /* ax <- AV_k, where V_k is real part of kth Ritz vector. */ /* Compute matrix-vector multiply using nag_dgbmv (f16pbc). */ nag_dgbmv(Nag_ColMajor, Nag_NoTrans, n, n, kl, ku, 1., &ab[kl], ldab, &v[k], 1, 0., ax, 1, &fail); /* mx <- MV_k */ nag_dgbmv(Nag_ColMajor, Nag_NoTrans, n, n, kl, ku, 1., &mb[kl], ldab, &v[k], 1, 0., mx, 1, &fail); /* ax <- ax - (lambda_j).re * mx. */ /* Compute vector update using nag_daxpby (f16ecc). */ alpha = -eigvr[j]; nag_daxpby(n, alpha, mx, 1, 1., ax, 1, &fail); /* mx <- MW_k where W_k is imaginary part kth Ritz vector. */ nag_dgbmv(Nag_ColMajor, Nag_NoTrans, n, n, kl, ku, 1., &mb[kl], ldab, &v[k+n], 1, 0., mx, 1, &fail); /* ax <- ax - (lambda_j).im * mx. */ alpha = eigvi[j]; nag_daxpby(n, alpha, mx, 1, 1., ax, 1, &fail); /* resr = norm(ax). */ /* Compute 2-norm of Ritz estimates using nag_dge_norm (f16rac).*/ nag_dge_norm(Nag_ColMajor, Nag_FrobeniusNorm, n, 1, ax, n, &resr, &fail); /* Imaginary part of residual Im[Ax-lamdaMx]. */ /* ax <- AW_k. */ /* Compute matrix-vector multiply using nag_dgbmv (f16pbc). */ nag_dgbmv(Nag_ColMajor, Nag_NoTrans, n, n, kl, ku, 1., &ab[kl], ldab, &v[k+n], 1, 0., ax, 1, &fail); /* ax <- ax - (lambda_j).re * mx. */ alpha = -eigvr[j]; /* Compute vector update using nag_daxpby (f16ecc). */ nag_daxpby(n, alpha, mx, 1, 1., ax, 1, &fail); /* mx <- MV_k. */ nag_dgbmv(Nag_ColMajor, Nag_NoTrans, n, n, kl, ku, 1., &mb[kl], ldab, &v[k], 1, 0., mx, 1, &fail); /* ax <- ax - (lambda_j).im * mx. */ alpha = -eigvi[j]; nag_daxpby(n, alpha, mx, 1, 1., ax, 1, &fail); /* resi = norm(ax). */ /* Compute 2-norm of Ritz estimates using nag_dge_norm (f16rac).*/ nag_dge_norm(Nag_ColMajor, Nag_FrobeniusNorm, n, 1, ax, n, &resi, &fail); /* Scale residual using Ritz value. */ /* Assign to Complex type using nag_complex (a02bac) */ res = nag_complex(resr,resi); eig = nag_complex(eigvr[j],eigvi[j]); eigest[j] = nag_complex_abs(res)/nag_complex_abs(eig); /* Set residual for second in conjugate pair. */ eigest[j+1] = eigest[j]; first = Nag_FALSE; } else { /* Second of complex conjugate pair; */ /* Already got residual from first in pair. */ first = Nag_TRUE; } k = k + n; } /* Print computed eigenvalues. */ Vprintf("\n The %4ld generalized Ritz values", nconv); Vprintf(" closest to \n"); Vprintf(" ( %7.4f, +-%7.4f ) are:\n\n", sigmar, sigmai); for (j = 0; j <= nconv-1; ++j) { if (eigest[j] <= 1.0e-10) { Vprintf("%8ld%5s( %7.4f, %7.4f ).\n", j+1, "",eigvr[j],eigvi[j]); } else { Vprintf("%8ld%5s( %7.4f, %7.4f )%5s%18.16f.\n", j+1, "",eigvr[j],eigvi[j]," *** ",eigest[j]); } } } else { Vprintf(" Error from nag_real_banded_sparse_eigensystem_sol (f12agc)." "\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 (eigvr) NAG_FREE(eigvr); if (eigvi) NAG_FREE(eigvi); if (eigest) NAG_FREE(eigest); if (mb) NAG_FREE(mb); if (mx) NAG_FREE(mx); if (resid) NAG_FREE(resid); if (v) NAG_FREE(v); if (icomm) NAG_FREE(icomm); return exit_status; }