1
/* This software was developed by Bruce Hendrickson and Robert Leland *
2
* at Sandia National Laboratories under US Department of Energy *
3
* contract DE-AC04-76DP00789 and is copyrighted by Sandia Corporation. */
7
#include "Gmsh_printf.h"
9
/* Finds first extended eigenpair of system corresponding to
10
tridiagonal T using using Rafael's bisection technique. */
12
void get_extval(alpha, beta, j, ritzval, s, eigtol, wnorm_g, sigma, extval, v, work1, work2)
13
double *alpha; /* j-vector of Lanczos scalars (using elements 1 to j) */
14
double *beta; /* (j+1)-vector of " " (has 0 element but using 1 to j-1) */
15
int j; /* number of Lanczos iterations taken */
16
double ritzval; /* Ritz value */
17
double *s; /* Ritz vector (length n, re-computed in this routine) */
18
double eigtol; /* tolerance on eigenpair */
19
double wnorm_g; /* W-norm of n-vector g, the rhs in the extended eig. problem */
20
double sigma; /* the norm constraint on the extended eigenvector */
21
double *extval; /* the extended eigenvalue this routine computes */
22
double *v; /* the j-vector solving the extended eig problem in T */
23
double *work1; /* j-vector of workspace */
24
double *work2; /* j-vector of workspace */
26
extern int DEBUG_EVECS; /* debug flag for eigen computation */
27
double lambda_low; /* lower bound on extended eval */
28
double lambda_high; /* upper bound on extended eval */
29
double tol; /* bisection tolerance */
30
double norm_v; /* norm of the extended T eigenvector v */
31
double lambda; /* the parameter that iterates to extval */
32
int cnt; /* debug iteration counter */
33
double diff; /* distance between lambda limits */
34
double norm(), Tevec();
35
void tri_solve(), cpvec();
37
/* Compute the Ritz vector */
38
Tevec(alpha, beta - 1, j, ritzval, s);
40
/* Shouldn't happen, but just in case ... */
44
if (DEBUG_EVECS > 0) {
45
printf("Degenerate extended eigenvector problem (g = 0).\n");
48
/* ... not really an extended eigenproblem; just return Ritz pair */
51
/* Set up the bisection parameters */
52
lambda_low = ritzval - wnorm_g / sigma;
53
lambda_high = ritzval - (wnorm_g / sigma) * s[1];
54
lambda = 0.5 * (lambda_low + lambda_high);
55
tol = eigtol * eigtol * (1 + fabs(lambda_low) + fabs(lambda_high));
57
if (DEBUG_EVECS > 2) {
58
printf("Computing extended eigenpairs of T\n");
59
printf(" target norm_v (= sigma) %g\n", sigma);
60
printf(" bisection tolerance %g\n", tol);
62
if (DEBUG_EVECS > 3) {
63
printf(" lambda iterates to the extended eigenvalue\n");
64
printf(" lambda_low lambda lambda_high norm_v\n");
67
/* Bisection loop - iterate until norm constraint is satisfied */
71
lambda = 0.5 * (lambda_low + lambda_high);
72
tri_solve(alpha, beta, j, lambda, v, wnorm_g, work1, work2);
73
norm_v = norm(v, 1, j);
74
if (DEBUG_EVECS > 3) {
75
printf("%2i %18.16f %18.16f %18.16f %g\n",
76
cnt++, lambda_low, lambda, lambda_high, norm_v);
82
diff = lambda_high - lambda_low;
85
/* Return the extended eigenvalue (eigvec is automatically returned) */
added
removed
contrib/Chaco/eigen/get_extval.c
