/* nag_complex_sparse_eigensystem_iter (f12apc) Example Program. * * Copyright 2005 Numerical Algorithms Group. * * Mark 8, 2005. */ #include #include #include #include #include #include #include #include static void my_zgttrf(Integer, Complex *, Complex *, Complex *, Complex *, Integer *, Integer *); static void my_zgttrs(Integer, Complex *, Complex *, Complex *, Complex *, Integer *, Complex *, Complex *); int main(void) { /* Constants */ Integer licomm=140, imon=0; /* Scalars */ Complex rho, s1, s2, s3, sigma; double estnrm, hr1, hr2, sr; Integer exit_status, info, irevcm, j, lcomm, n, nconv; Integer ncv, nev, niter, nshift, nx; /* Nag types */ NagError fail; /* Arrays */ Complex *comm=0, *eigv=0, *eigest=0, *dd=0, *dl=0, *du=0; Complex *du2=0, *resid=0, *v=0; Integer *icomm=0, *ipiv=0; /* Ponters */ Complex *mx=0, *x=0, *y=0; exit_status = 0; INIT_FAIL(fail); Vprintf("nag_complex_sparse_eigensystem_iter (f12apc) Example Program " "Results\n"); /* Skip heading in data file */ Vscanf("%*[^\n] "); Vscanf("%ld%ld%ld%*[^\n] ", &nx, &nev, &ncv); n = nx * nx; lcomm = 3*n + 3*ncv*ncv + 5*ncv + 60; /* Allocate memory */ if ( !(comm = NAG_ALLOC(lcomm, Complex)) || !(eigv = NAG_ALLOC(ncv, Complex)) || !(eigest = NAG_ALLOC(ncv, Complex)) || !(dd = NAG_ALLOC(n, Complex)) || !(dl = NAG_ALLOC(n, Complex)) || !(du = NAG_ALLOC(n, Complex)) || !(du2 = NAG_ALLOC(n, Complex)) || !(resid = NAG_ALLOC(n, Complex)) || !(v = NAG_ALLOC(n * ncv, Complex)) || !(icomm = NAG_ALLOC(licomm, Integer)) || !(ipiv = NAG_ALLOC(n, Integer)) ) { Vprintf("Allocation failure\n"); exit_status = -1; goto END; } /* Initialise communication arrays for problem using nag_complex_sparse_eigensystem_init (f12anc). */ nag_complex_sparse_eigensystem_init(n, nev, ncv, icomm, licomm, comm, lcomm, &fail); /* Select the required mode using nag_complex_sparse_eigensystem_option (f12arc). */ nag_complex_sparse_eigensystem_option("SHIFTED INVERSE", icomm, comm, &fail); /* Set values for sigma and rho */ /* Assign to Complex type using nag_complex (a02bac) */ sigma = nag_complex(0.0, 0.0); rho = nag_complex(10.0, 0.0); hr1 = (double) (n+1); /* 1/h */ hr2 = hr1*hr1; /* 1/(h*h) */ sr = 0.5*hr1*rho.re; /* s/h */ /* Assign to Complex type using nag_complex (a02bac) */ s1 = nag_complex(-hr2-sr,0.0); /* -1/(h*h) - s/h */ s3 = nag_complex(-hr2+sr,0.0); /* -1/(h*h) + s/h */ s2 = nag_complex(2.0*hr2,0.0); /* Compute Complex subtraction using nag_complex_subtract (a02cbc). */ s2 = nag_complex_subtract(s2,sigma); /* two/(h*h) - sigma */ for (j = 0; j <= n - 2; ++j) { dl[j] = s1; dd[j] = s2; du[j] = s3; } dd[n - 1] = s2; my_zgttrf(n, dl, dd, du, du2, ipiv, &info); irevcm = 0; REVCOMLOOP: /* repeated calls to reverse communication routine nag_complex_sparse_eigensystem_iter (f12apc). */ nag_complex_sparse_eigensystem_iter(&irevcm, resid, v, &x, &y, &mx, &nshift, comm, icomm, &fail); if (irevcm != 5) { if (irevcm == -1 || irevcm == 1) { /* Perform x <--- OP*x = inv[A-sigma*I]*x */ my_zgttrs(n, dl, dd, du, du2, ipiv, x, y); } else if (irevcm == 4 && imon == 1) { /* If imon=1, get monitoring information using nag_complex_sparse_eigensystem_monit (f12asc). */ nag_complex_sparse_eigensystem_monit(&niter, &nconv, eigv, eigest, icomm, comm); /* Compute 2-norm of Ritz estimates using nag_zge_norm (f16uac). */ nag_zge_norm(Nag_ColMajor, Nag_FrobeniusNorm, nev, 1, eigest, nev, &estnrm, &fail); Vprintf("Iteration %3ld, ", niter); Vprintf(" No. converged = %3ld,", nconv); Vprintf(" norm of estimates = %16.8e\n\n", estnrm); } goto REVCOMLOOP; } if (fail.code == NE_NOERROR) { /* Post-Process using nag_complex_sparse_eigensystem_sol (f12aqc) to compute eigenvalues/vectors. */ nag_complex_sparse_eigensystem_sol(&nconv, eigv, v, sigma, resid, v, comm, icomm, &fail); Vprintf("\n"); Vprintf(" The %4ld Ritz values of smallest magnitude are:\n\n", nconv); for (j = 0; j <= nconv-1; ++j) { Vprintf("%8ld%5s( %12.4f , %12.4f )\n", j+1, "", eigv[j].re, eigv[j].im); } } else { Vprintf(" Error from nag_complex_sparse_eigensystem_iter (f12apc)." "\n%s\n", fail.message); exit_status = 1; goto END; } END: if (comm) NAG_FREE(comm); if (eigv) NAG_FREE(eigv); if (eigest) NAG_FREE(eigest); if (dd) NAG_FREE(dd); if (dl) NAG_FREE(dl); if (du) NAG_FREE(du); if (du2) NAG_FREE(du2); if (resid) NAG_FREE(resid); if (v) NAG_FREE(v); if (icomm) NAG_FREE(icomm); if (ipiv) NAG_FREE(ipiv); return exit_status; } static void my_zgttrf(Integer n, Complex dl[], Complex d[], Complex du[], Complex du2[], Integer ipiv[], Integer *info) { /* A simple C version of the Lapack routine zgttrf with argument checking removed */ /* Scalars */ Complex temp, fact, z1; Integer i; /* Function Body */ *info = 0; for (i = 0; i < n; ++i) { ipiv[i] = i; } for (i = 0; i < n - 2; ++i) { du2[i] = nag_complex(0.0,0.0); } for (i = 0; i < n - 2; ++i) { if (fabs(d[i].re)+fabs(d[i].im) >= fabs(dl[i].re)+fabs(dl[i].im)) { /* No row interchange required, eliminate dl[i]. */ if (fabs(d[i].re)+fabs(d[i].im) != 0.0) { /* Compute Complex division using nag_complex_divide (a02cdc). */ fact = nag_complex_divide(dl[i],d[i]); dl[i] = fact; /* Compute Complex multiply using nag_complex_multiply (a02ccc). */ fact = nag_complex_multiply(fact,du[i]); /* Compute Complex subtraction using nag_complex_subtract (a02cbc). */ d[i+1] = nag_complex_subtract(d[i+1],fact); } } else { /* Interchange rows I and I+1, eliminate dl[I] */ /* Compute Complex division using nag_complex_divide (a02cdc). */ fact = nag_complex_divide(d[i],dl[i]); d[i] = dl[i]; dl[i] = fact; temp = du[i]; du[i] = d[i+1]; /* Compute Complex multiply using nag_complex_multiply (a02ccc). */ z1 = nag_complex_multiply(fact,d[i+1]); /* Compute Complex subtraction using nag_complex_subtract (a02cbc). */ d[i+1] = nag_complex_subtract(temp,z1); du2[i] = du[i+1]; /* Compute Complex multiply using nag_complex_multiply (a02ccc). */ du[i+1] = nag_complex_multiply(fact,du[i+1]); /* Perform Complex negation using nag_complex_negate (a02cec). */ du[i+1] = nag_complex_negate(du[i+1]); ipiv[i] = i + 1; } } if (n > 1) { i = n - 2; if (fabs(d[i].re)+fabs(d[i].im) >= fabs(dl[i].re)+fabs(dl[i].im)) { if (fabs(d[i].re)+fabs(d[i].im) != 0.0) { /* Compute Complex division using nag_complex_divide (a02cdc). */ fact = nag_complex_divide(dl[i],d[i]); dl[i] = fact; /* Compute Complex multiply using nag_complex_multiply (a02ccc). */ fact = nag_complex_multiply(fact,du[i]); /* Compute Complex subtraction using nag_complex_subtract (a02cbc). */ d[i+1] = nag_complex_subtract(d[i+1],fact); } } else { /* Compute Complex division using nag_complex_divide (a02cdc). */ fact = nag_complex_divide(d[i],dl[i]); d[i] = dl[i]; dl[i] = fact; temp = du[i]; du[i] = d[i+1]; /* Compute Complex multiply using nag_complex_multiply (a02ccc). */ z1 = nag_complex_multiply(fact,d[i+1]); /* Compute Complex subtraction using nag_complex_subtract (a02cbc). */ d[i+1] = nag_complex_subtract(temp,z1); ipiv[i] = i + 1; } } /* Check for a zero on the diagonal of U. */ for (i = 0; i < n; ++i) { if (fabs(d[i].re)+fabs(d[i].im) == 0.0) { *info = i; goto END; } } END: return; } static void my_zgttrs(Integer n, Complex dl[], Complex d[], Complex du[], Complex du2[], Integer ipiv[], Complex b[], Complex y[]) { /* A simple C version of the Lapack routine zgttrs with argument checking removed, the number of right-hand-sides=1, Trans='N' */ /* Scalars */ Complex temp, z1; Integer i; /* Solve L*x = b. */ for (i = 0; i < n ; ++i) { y[i] = b[i]; } for (i = 0; i < n - 1; ++i) { if (ipiv[i] == i) { /* y[i+1] = y[i+1] - dl[i]*y[i] */ /* Compute Complex multiply using nag_complex_multiply (a02ccc). */ temp = nag_complex_multiply(dl[i],y[i]); /* Compute Complex subtraction using nag_complex_subtract (a02cbc). */ y[i+1] = nag_complex_subtract(y[i+1],temp); } else { temp = y[i]; y[i] = y[i+1]; /* Compute Complex multiply using nag_complex_multiply (a02ccc). */ z1 = nag_complex_multiply(dl[i],y[i]); /* Compute Complex subtraction using nag_complex_subtract (a02cbc). */ y[i+1] = nag_complex_subtract(temp,z1); } } /* Solve U*x = b. */ /* Compute Complex division using nag_complex_divide (a02cdc). */ y[n-1] = nag_complex_divide(y[n-1],d[n-1]); if (n > 1) { /* Compute Complex multiply using nag_complex_multiply (a02ccc). */ temp = nag_complex_multiply(du[n-2],y[n-1]); /* Compute Complex subtraction using nag_complex_subtract (a02cbc). */ z1 = nag_complex_subtract(y[n-2],temp); /* Compute Complex division using nag_complex_divide (a02cdc). */ y[n-2] = nag_complex_divide(z1,d[n-2]); } for (i = n - 3; i >= 0; --i) { /* y[i] = (y[i]-du[i]*y[i+1]-du2[i]*y[i+2])/d[i]; */ /* Compute Complex multiply using nag_complex_multiply (a02ccc). */ temp = nag_complex_multiply(du[i],y[i+1]); z1 = nag_complex_multiply(du2[i],y[i+2]); /* Compute Complex addition using nag_complex_add (a02cac). */ temp = nag_complex_add(temp,z1); /* Compute Complex subtraction using nag_complex_subtract (a02cbc). */ z1 = nag_complex_subtract(y[i],temp); /* Compute Complex division using nag_complex_divide (a02cdc). */ y[i] = nag_complex_divide(z1,d[i]); } return; }