arpack_general
Computes some of the eigenvalues and eigenvectors of the generalized nonsymmetric eigenvalue problem Ax = λBx using an implicitly restarted Arnoldi method (IRAM). The algorithm can also be used for the standard case B = I. The matrices A, B are real, but eigenvalues may be complex and occur in conjugate pairs.
NOTE: Function arpack_general is available in double precision only. |
Synopsis
#include <imsl.h>
d_complex *imsl_d_arpack_general (void fcn(), int n, int nev, …, 0)
Required Arguments
void fcn (int n, double x[], int task, double y[]) (Input)
User-supplied function to return matrix-vector operations or solutions of linear systems.
int n (Input)
The dimension of the problem.
double x[] (Input)
An array of size n containing the vector to which the operator will be applied.
int task (Input)
An enumeration type which specifies the operation to be performed. Variable task is an enumerated integer value associated with enum type Imsl_arpack_task. Table 15 describes the possible values.
task | Description |
|---|---|
IMSL_ARPACK_PREPARE | Take initial steps to prepare for the operations to follow. These steps can include defining data for the matrices, factorizations for upcoming linear system solves or recording the vectors used in the operations. |
IMSL_ARPACK_A_X | y = Ax |
IMSL_ARPACK_B_X | y = Bx |
IMSL_ARPACK_INV_SHIFT_X | y = (A ‑ σB)-1x, σ real, or y = Re{(A ‑ σB)-1x}, σ complex, or y = Im{(A ‑ σB)-1x}, σ complex, B general or = I |
IMSL_ARPACK_INV_B_X | y = B-1x |
double y[] (Output)
An array of size n containing the result of a matrix-vector operation or the solution of a linear system.
int n (Input)
The dimension of the problem.
int nev (Input)
The number of eigenvalues to be computed. It is required that 0 < nev < n - 1.
Return Value
A pointer to the nev eigenvalues of the general eigenvalue problem. Complex conjugate eigenvalues are stored consecutively, with the eigenvalue with positive imaginary part in the first place. To release this space, use imsl_free. If no value can be computed, then NULL is returned.
Synopsis with Optional Arguments
#include <imsl.h>
d_complex *imsl_d_arpack_general (void fcn(), int n, int nev,
IMSL_XGUESS, double xguess[],
IMSL_ITMAX, int itmax,
IMSL_TOLERANCE, double tol,
IMSL_SHIFT, d_complex shift,
IMSL_EIGVAL_LOCATION, Imsl_arpack_eigval_location eigval_loc,
IMSL_EIG_PROBLEM_TYPE, Imsl_arpack_problem_type problem_type,
IMSL_EIG_SOLVE_MODE, Imsl_arpack_solve_mode mode,
IMSL_NUM_ARNOLDI_VECTORS, int ncv,
IMSL_NUM_ACCURATE_EIGVALS, int *n_acc,
IMSL_VECTORS, double **evec,
IMSL_VECTORS_USER, double evecu[],
IMSL_EVECU_COL_DIM, int evecu_col_dim,
IMSL_RETURN_USER, d_complex evalu[],
IMSL_FCN_W_DATA, void fcn(), void *data,
0)
Optional Arguments
IMSL_XGUESS, double xguess[] (Input)
A non-zero vector of size n containing the starting vector for the implicitly restarted Arnoldi method.
A non-zero vector of size n containing the starting vector for the implicitly restarted Arnoldi method.
By default, a random starting vector is computed internally.
IMSL_ITMAX, int itmax (Input)
The maximum number of Arnoldi iterations.
Default: itmax = 1000.
The maximum number of Arnoldi iterations.
Default: itmax = 1000.
IMSL_TOLERANCE, double tol (Input)
Tolerance value used in the criterion for the acceptability of the relative accuracy of the Ritz values.
Default: tol = imsl_f_machine(3).
Tolerance value used in the criterion for the acceptability of the relative accuracy of the Ritz values.
Default: tol = imsl_f_machine(3).
IMSL_SHIFT, d_complex shift (Input)
The shift value used in the shift-invert spectral transformations.
Default: shift = {0, 0}.
The shift value used in the shift-invert spectral transformations.
Default: shift = {0, 0}.
IMSL_EIGVAL_LOCATION, Imsl_arpack_eigval_location eigval_loc (Input)
An enumeration type which specifies the location of the eigenvalues to compute.
An enumeration type which specifies the location of the eigenvalues to compute.
eigval_loc | Description |
|---|---|
IMSL_ARPACK_LARGEST_MAGNITUDE | Compute eigenvalues of largest magnitude. |
IMSL_ARPACK_SMALLEST_MAGNITUDE | Compute eigenvalues of smallest magnitude. |
IMSL_ARPACK_LARGEST_REAL_PART | Compute eigenvalues of largest algebraic real part. |
IMSL_ARPACK_SMALLEST_REAL_PART | Compute eigenvalues of smallest algebraic real part. |
IMSL_ARPACK_LARGEST_IMAG_PART | Compute eigenvalues of largest imaginary part magnitude. |
IMSL_ARPACK_SMALLEST_IMAG_PART | Compute eigenvalues of smallest imaginary part magnitude. |
For computational modes that use a spectral transformation the eigenvalue location refers to the transformed—not the original—problem. See the Description section for an example.
Default: eigval_loc = IMSL_ARPACK_LARGEST_MAGNITUDE.
Default: eigval_loc = IMSL_ARPACK_LARGEST_MAGNITUDE.
IMSL_EIG_PROBLEM_TYPE, Imsl_arpack_problem_type problem_type (Input)
An enumeration type that indicates if a standard or generalized eigenvalue problem is to be solved.
An enumeration type that indicates if a standard or generalized eigenvalue problem is to be solved.
problem_type | Description |
|---|---|
IMSL_ARPACK_STANDARD | Solve standard problem, Ax = λx. |
IMSL_ARPACK_GENERALIZED | Solve generalized problem, Ax = λBx. |
Default: problem_type = IMSL_ARPACK_STANDARD.
IMSL_EIG_SOLVE_MODE, Imsl_arpack_solve_mode mode (Input)
An enumeration type indicating which computational mode is used for the eigenvalue computation. Variables problem_type and mode together define the tasks that must be provided in the user-supplied function. The following table describes the values variable mode can take, the feasible combinations with variable problem_type and the related tasks:
An enumeration type indicating which computational mode is used for the eigenvalue computation. Variables problem_type and mode together define the tasks that must be provided in the user-supplied function. The following table describes the values variable mode can take, the feasible combinations with variable problem_type and the related tasks:
mode | problem_type | Required tasks |
|---|---|---|
IMSL_ARPACK_REGULAR | IMSL_ARPACK_STANDARD | y = Ax |
IMSL_ARPACK_REGULAR_INVERSE | IMSL_ARPACK_GENERALIZED | y = Ax, y = Bx, y = B-1x |
IMSL_ARPACK_SHIFT_INVERT | IMSL_ARPACK_STANDARD | y = (A ‑ σI)-1x |
IMSL_ARPACK_SHIFT_INVERT | IMSL_ARPACK_GENERALIZED | y = Bx, y = (A ‑ σB)-1x |
IMSL_ARPACK_REAL_SHIFT_INVERT | IMSL_ARPACK_GENERALIZED | y = Ax, y = Bx, y = Re{(A ‑ σB)-1x} |
IMSL_ARPACK_IMAG_SHIFT_INVERT | IMSL_ARPACK_GENERALIZED | y = Ax, y = Bx, y = Im{(A ‑ σB)-1x} |
Default: mode = IMSL_ARPACK_REGULAR.
IMSL_NUM_ARNOLDI_VECTORS, int ncv (Input)
The number of Arnoldi vectors generated in each iteration of the Arnoldi method. It is required that nev + 2 <= ncv <= n. A value ncv >= min(2*nev + 1, n) is recommended.
The number of Arnoldi vectors generated in each iteration of the Arnoldi method. It is required that nev + 2 <= ncv <= n. A value ncv >= min(2*nev + 1, n) is recommended.
Default: ncv = min(2*nev + 1, n).
IMSL_NUM_ACCURATE_EIGVALS, int *n_acc (Output)
The number of eigenvalues that the algorithm was able to compute accurately. This number can be smaller than nev.
The number of eigenvalues that the algorithm was able to compute accurately. This number can be smaller than nev.
IMSL_VECTORS, double **evec (Output)
The address of a pointer to an array of size n × (nev+1) containing the B-orthonormalized eigenvectors corresponding to the n_acc converged eigenvalues. Typically, double *evec is declared, and &evec is used as an argument. For a closer description of the array content, see optional argument IMSL_VECTORS_USER.
The address of a pointer to an array of size n × (nev+1) containing the B-orthonormalized eigenvectors corresponding to the n_acc converged eigenvalues. Typically, double *evec is declared, and &evec is used as an argument. For a closer description of the array content, see optional argument IMSL_VECTORS_USER.
IMSL_VECTORS_USER, double evecu[] (Output)
A user-defined array of size n × (nev+1) containing the n_acc B-orthonormalized eigenvectors of the eigenvalue problem in compact form. The eigenvectors are stored column-wise in the same order as the eigenvalues. An eigenvector corresponding to a real eigenvalue is real and represented by a single column. For a complex conjugate pair of eigenvalues, the real and imaginary parts of the eigenvector related to the eigenvalue with positive imaginary part are stored in two consecutive columns of array evecu. If an eigenvalue is complex and has no complex conjugate counterpart due to the choice of nev, then the corresponding eigenvector is stored in two consecutive columns of evecu.
A user-defined array of size n × (nev+1) containing the n_acc B-orthonormalized eigenvectors of the eigenvalue problem in compact form. The eigenvectors are stored column-wise in the same order as the eigenvalues. An eigenvector corresponding to a real eigenvalue is real and represented by a single column. For a complex conjugate pair of eigenvalues, the real and imaginary parts of the eigenvector related to the eigenvalue with positive imaginary part are stored in two consecutive columns of array evecu. If an eigenvalue is complex and has no complex conjugate counterpart due to the choice of nev, then the corresponding eigenvector is stored in two consecutive columns of evecu.
IMSL_EVECU_COL_DIM, int evecu_col_dim (Input)
The column dimension of evecu.
Default: evecu_col_dim = nev + 1
The column dimension of evecu.
Default: evecu_col_dim = nev + 1
IMSL_RETURN_USER, d_complex evalu[] (Output)
An array of size nev containing the accurately computed eigenvalues in the first n_acc locations. Complex conjugate pairs of eigenvalues are stored consecutively in evalu.
An array of size nev containing the accurately computed eigenvalues in the first n_acc locations. Complex conjugate pairs of eigenvalues are stored consecutively in evalu.
IMSL_FCN_W_DATA, void fcn (int n, double x[], int task, double y[]), void *data, (Input/Output)
User-supplied function to return matrix-vector operations or solutions of linear systems, which also accepts a pointer to data that is supplied by the user. data is a pointer to the data to be passed to the user-supplied function.
User-supplied function to return matrix-vector operations or solutions of linear systems, which also accepts a pointer to data that is supplied by the user. data is a pointer to the data to be passed to the user-supplied function.
NOTE: The possibility to supply user-data via IMSL_FCN_W_DATA is an important feature of arpack_general. It allows the user to transfer problem-specific data to the algorithm without the need to define global data. See Passing Data to User-Supplied Functions in the introduction to this manual for more details. |
Description
Function imsl_d_arpack_general, which is based on ARPACK subroutines DNAUPD and DNEUPD (see the ARPACK Users' Guide, Lehoucq et al. (1998)), computes selected eigenvalue-eigenvector pairs for generalized nonsymmetric eigenvalue problems of the form
Ax = λBx.
Here, A is a real general and B a positive semi-definite matrix. For B = I, the generalized problem reduces to the standard nonsymmetric eigenvalue problem.
The ARPACK routine DNAUPD implements a variant of the Arnoldi method and uses reverse communication to obtain the required matrix-vector products or solutions of linear systems for the iterations. Responses to these requests are made by calling the user-defined function fcn. User data can be made available for the evaluations by optional argument IMSL_FCN_W_DATA.
For a given problem, the requested responses depend on the settings of optional arguments IMSL_EIG_PROBLEM_TYPE and IMSL_EIG_SOLVE_MODE. For each response, a corresponding task must be defined in the user-defined function fcn. The Mode/problem type combinations table under optional argument IMSL_EIG_SOLVE_MODE shows which tasks have to be defined for a certain problem.
The following code snippet shows the complete list of tasks available for fcn and their meaning:
void fcn(int n, double x[], int itask, double y[])
{
switch (itask) {
/*
* Define responses to different tasks for the generalized
* eigenvalue problem
* A*x = lambda * B * x,
* which includes the ordinary case B = I.
*/
case IMSL_ARPACK_PREPARE:
/*
* Take initial steps to prepare for the operations
* that follow. Note that arpack_general internally
* always calls fcn with this enum value, even if it is
* not required by the user.
*/
break;
case IMSL_ARPACK_A_X:
/*
* Compute matrix-vector product y = A * x
*/
break;
case IMSL_ARPACK_B_X:
/*
* Compute matrix-vector product y = B * x
*/
break;
case IMSL_ARPACK_INV_SHIFT_X:
/*
* Compute the matrix-vector product
* z = inv(A - sigma * B) * x,
* and return
* y = z, if mode = IMSL_ARPACK_SHIFT_INVERT, sigma real,
* y = Re{z}, if mode = IMSL_ARPACK_REAL_SHIFT_INVERT, sigma complex,
* y = Im{z}, if mode = IMSL_ARPACK_IMAG_SHIFT_INVERT, sigma complex.
*
* Usually, matrix A - sigma * B is not directly inverted.
* Instead, a factorization of A - sigma * B is determined,
* and the factors are used to compute z via backsolves.
*
* Example:
* If an LU factorization of A - sigma * B exists, then
* A - sigma * B = P * L * U,
* P a permutation matrix. Vector z can then be determined
* as solution of the linear system
* L * U * z = trans(P) * x.
* The LU factorization only has to be computed once, for
* example outside of fcn or within IMSL_ARPACK_PREPARE.
*/
break;
case IMSL_ARPACK_INV_B_X:
/*
* Compute matrix-vector product
* y = inv(B) * x.
* Usually, matrix B is not directly inverted.
* Instead, a factorization of B is determined, and the
* factors are used to compute y via backsolves.
*
* Example:
* If matrix B is positive definite, then a Cholesky
* factorization B = L * trans(L) exists. Vector y can then
* be determined by solving the linear system
* L * trans(L) * y = x.
* The Cholesky factorization only has to be computed once,
* for example outside of fcn or within IMSL_ARPACK_PREPARE.
*/
break;
default:
/*
* Define error conditions, if necessary.
*/
break;
}
}
Internally, imsl_d_arpack_general first determines the eigenvalues for the problem specified by optional arguments IMSL_EIG_SOLVE_MODE and IMSL_EIG_PROBLEM_TYPE.
Table 19 shows the matrices whose eigenvalues are determined for a given combination of these optional arguments.
mode | problem_type | Matrix |
|---|---|---|
IMSL_ARPACK_REGULAR | IMSL_ARPACK_STANDARD | A |
IMSL_ARPACK_REGULAR_INVERSE | IMSL_ARPACK_GENERALIZED | B-1A |
IMSL_ARPACK_SHIFT_INVERT | IMSL_ARPACK_STANDARD | (A ‑ σI)-1, σ real |
IMSL_ARPACK_SHIFT_INVERT | IMSL_ARPACK_GENERALIZED | (A ‑ σB)-1B, σ real |
IMSL_ARPACK_REAL_SHIFT_INVERT | IMSL_ARPACK_GENERALIZED | Re{(A ‑ σB)-1B},σ complex |
IMSL_ARPACK_IMAG_SHIFT_INVERT | IMSL_ARPACK_GENERALIZED | Im{(A ‑ σB)-1B}, σ complex |
Note that the eigenvalue location defined by optional argument IMSL_EIGVAL_LOCATION always refers to the matrices of Table 19.
For example, for mode=IMSL_ARPACK_SHIFT_INVERT, problem_type=IMSL_ARPACK_STANDARD, and eigval_loc=IMSL_ARPACK_LARGEST_MAGNITUDE, the eigenvalues of largest magnitude of the shift-inverted matrix (A ‑ σI)-1 are computed. Because of the relationship
these eigenvalues correspond to the eigenvalues of the original problem Ax = λx that are closest to the shift σ in absolute value.
In a second step, implemented via ARPACK routine DNEUPD, imsl_d_arpack_general internally transforms the eigenvalues back to the eigenvalues of the original problem Ax = λBx or Ax = λx and computes eigenvectors, if required.
Besides matrix A being real and general, the modes for the generalized eigenproblem require some additional properties of matrix B that are summarized in Table 20:
mode | Matrix properties |
|---|---|
IMSL_ARPACK_REGULAR_INVERSE | B positive definite |
IMSL_ARPACK_SHIFT_INVERT | B positive semi-definite |
IMSL_ARPACK_REAL_SHIFT_INVERT | B positive semi-definite |
IMSL_ARPACK_IMAG_SHIFT_INVERT | B positive semi-definite |
If the nonsymmetric problem has complex eigenvalues in conjugate pairs, the eigenvectors are returned in a compact representation: If the eigenvalue λj has a positive imaginary part, the complex eigenvector is constructed from the relation wj = vj + ivj + 1. The real vectors vj, vj + 1 are consecutive columns of the arrays evec or evecu. The eigenvalue-eigenvector relationship is Awj = λjwj. Since A is real, is also an eigenvalue: . For purposes of checking results, the complex residual should be small in norm relative to the norm of A. Since the norms of and are identical, a check of the alternate relationship is not necessary. In the case of a real eigenvalue, the associated eigenvector is real and represented by a single column in evec or evecu.
Copyright notice for ARPACK
Copyright (c) 1996-2008 Rice University. Developed by D.C. Sorensen, R.B. Lehoucq, C. Yang, and K. Maschhoff. All rights reserved.
Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer listed in this license in the documentation and/or other materials provided with the distribution.
- Neither the name of the copyright holders nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
Examples
Example 1
The generalized eigenvalue problem Ax = λBx is solved using shift-invert strategies. The matrix A is tri-diagonal with the values 2 on the diagonal, -2 on the sub-diagonal and 3 on the super-diagonal. The matrix B is tri-diagonal positive definite with the values 4 on the diagonal and 1 on the off-diagonals. A complex shift σ = 0.4 + 0.6i is used. Two strategies of shift-invert are illustrated, y = Re{(A - σB)-1Bx} and y = Im{(A‑σB)-1Bx}. In each case, nev=6 eigenvalues are obtained, each with 3 pairs of complex conjugate values.
#include <math.h>
#include <imsl.h>
#include <omp.h>
#include <stdio.h>
static void ax(int n, double x[], double y[]);
static void bx(int n, double x[], double y[]);
static void fcn_w_data(int n, double x[], int itask, double y[],
void *data);
static void compute_residuals_fcn(
void(*fcn) (int n, double x[], int itask, double y[], void *data),
int problem_type, int n, int nev, int nconv, d_complex eigvals[],
double evecu[], double ax[], void *data);
typedef struct {
d_complex *band_matrix;
int *ipvt;
d_complex *factor;
d_complex *work;
d_complex sigma;
int shift_strategy;
} imsl_arpack_data;
int main() {
int n, nev, ncv, n_acc, i, j;
int *ipvt = NULL;
d_complex sigma = { 0.4, 0.6 };
d_complex *a_matrix = NULL, *factor = NULL, *eigvals = NULL;
d_complex *work = NULL;
double *evecu = NULL, *rwork = NULL;
imsl_arpack_data usr_data;
Imsl_arpack_solve_mode mode[] = {
IMSL_ARPACK_REAL_SHIFT_INVERT, IMSL_ARPACK_IMAG_SHIFT_INVERT
};
n = 100;
nev = 6;
ncv = 30;
eigvals = (d_complex *)malloc(nev * sizeof(d_complex));
evecu = (double *)malloc(n * (nev + 1) * sizeof(double));
rwork = (double *)malloc((nev + 3 * n) * sizeof(double));
/* Allocate arrays needed in the LU factorization */
ipvt = (int *)malloc(n * sizeof(int));
a_matrix = (d_complex *)malloc(3 * n * sizeof(d_complex));
factor = (d_complex *)malloc(4 * n * sizeof(d_complex));
work = (d_complex *)malloc(2 * n * sizeof(d_complex));
if (!eigvals || !evecu || !rwork || !ipvt || !a_matrix || !factor
|| !work) {
printf("Memory allocation error\n");
goto FREE_SPACE;
}
usr_data.band_matrix = a_matrix;
usr_data.ipvt = ipvt;
usr_data.factor = factor;
usr_data.sigma = sigma;
usr_data.work = work;
for (i = 0; i <= 1; i++) {
usr_data.shift_strategy = i;
imsl_d_arpack_general(NULL, n, nev,
IMSL_EIG_PROBLEM_TYPE, IMSL_ARPACK_GENERALIZED,
IMSL_EIG_SOLVE_MODE, mode[i],
IMSL_SHIFT, sigma,
IMSL_NUM_ARNOLDI_VECTORS, ncv,
IMSL_NUM_ACCURATE_EIGVALS, &n_acc,
IMSL_FCN_W_DATA, fcn_w_data, &usr_data,
IMSL_VECTORS_USER, evecu,
IMSL_RETURN_USER, eigvals,
0);
printf("\nNumber of requested eigenvalues : %d\n", nev);
printf("Number of accurate (converged) eigenvalues : %d\n\n",
n_acc);
compute_residuals_fcn(fcn_w_data, 1, n, nev, n_acc, eigvals,
evecu, rwork, &usr_data);
/*
* Display eigenvalues and corresponding residuals
* || A * x - lambda * B * x || / |lambda|
*/
if (i == 0) {
printf(" Largest magnitude eigenvalues, real shift\n");
printf(" =========================================\n");
}
else {
printf(" Largest magnitude eigenvalues, imaginary shift"
"\n");
printf(" =============================================="
"\n");
}
printf("%30s%26s\n", "Eigenvalues (Real, Imag)",
"Relative residuals");
for (j = 0; j < n_acc; j++) {
printf("(%14.8lf, %14.8lf)%20.8lf\n", eigvals[j].re,
eigvals[j].im, rwork[j]);
}
}
FREE_SPACE:
if (eigvals)
free(eigvals);
if (ipvt)
free(ipvt);
if (a_matrix)
free(a_matrix);
if (factor)
free(factor);
if (work)
free(work);
if (rwork)
free(rwork);
if (evecu)
free(evecu);
}
static void fcn_w_data(int n, double x[], int itask, double y[],
void *data)
{
int j, shift_strategy;
int *ipvt = NULL;
d_complex *c_matrix = NULL, *factor = NULL, *work = NULL;
d_complex cl, cdiag, cu;
d_complex sigma;
imsl_arpack_data *usr_data = (imsl_arpack_data *)data;
c_matrix = usr_data->band_matrix;
ipvt = usr_data->ipvt;
factor = usr_data->factor;
sigma = usr_data->sigma;
work = usr_data->work;
shift_strategy = usr_data->shift_strategy;
switch (itask) {
case IMSL_ARPACK_PREPARE:
/*
* Create tridiagonal matrix
* C := A - shift * B
* in complex arithmetic.
*/
cl = imsl_zd_convert(-2.0 - sigma.re, -sigma.im);
cdiag = imsl_zd_convert(2.0 - 4.0 * sigma.re, -4.0 * sigma.im);
cu = imsl_zd_convert(3.0 - sigma.re, -sigma.im);
for (j = 1; j <= n; j++) {
c_matrix[j - 1] = cu;
c_matrix[n + j - 1] = cdiag;
c_matrix[2 * n + j - 1] = cl;
}
/* Compute LU factorization of tridiagonal matrix */
imsl_z_lin_sol_gen_band(n, c_matrix, 1, 1, NULL,
IMSL_FACTOR_USER, ipvt, factor,
IMSL_FACTOR_ONLY,
0);
if (imsl_error_type() != 0) {
imsl_set_user_fcn_return_flag(1);
}
break;
case IMSL_ARPACK_A_X:
ax(n, x, y);
break;
case IMSL_ARPACK_B_X:
bx(n, x, y);
break;
case IMSL_ARPACK_INV_SHIFT_X:
/*
* Solve (A - sigma * M) * y = x in complex arithmetic
*/
for (j = 0; j < n; j++) {
work[j] = imsl_zd_convert(x[j], 0.0);
}
imsl_z_lin_sol_gen_band(n, NULL, 1, 1, work,
IMSL_FACTOR_USER, ipvt, factor,
IMSL_RETURN_USER, &work[n],
IMSL_SOLVE_ONLY,
0);
if (imsl_error_type() != 0) {
imsl_set_user_fcn_return_flag(2);
}
if (shift_strategy == 0) {
for (j = 0; j < n; j++) {
y[j] = work[n + j].re;
}
}
else if (shift_strategy == 1) {
for (j = 0; j < n; j++) {
y[j] = work[n + j].im;
}
}
break;
default:
imsl_set_user_fcn_return_flag(3);
break;
}
}
/*
* Matrix-vector multiplication function
*
* Computes the matrix vector multiplication y <- M*x, where M is
* an n by n symmetric tridiagonal matrix with 4 on the diagonal, 1
* on the subdiagonal and superdiagonal.
*/
static void bx(int n, double x[], double y[]) {
int j;
y[0] = 4.0 * x[0] + x[1];
for (j = 2; j <= n - 1; j++) {
y[j - 1] = x[j - 2] + 4.0 * x[j - 1] + x[j];
}
y[n - 1] = x[n - 2] + 4.0 * x[n - 1];
}
/*
* Matrix-vector multiplication function
*
* Compute the matrix vector multiplication y <- A*x where A is an
* n by n symmetric tridiagonal matrix with 2 on the diagonal, -2
* on the subdiagonal and 3 on the superdiagonal.
*/
static void ax(int n, double x[], double y[]) {
int j;
y[0] = 2.0 * x[0] + 3.0 * x[1];
for (j = 2; j <= n - 1; j++) {
y[j - 1] = -2.0 * x[j - 2] + 2.0 * x[j - 1] + 3.0 * x[j];
}
y[n - 1] = -2.0 * x[n - 2] + 2.0 * x[n - 1];
}
/*
* Compute residuals
* || A * x - lambda * B * x || / |lambda|,
* including the case B = I.
*
* problem_type = 0 (standard) --> size(rwork) >= nev + 2 * n
* problem_type = 1 (generalized) --> size(rwork) >= nev + 3 * n
*/
static void compute_residuals_fcn(
void(*fcn) (int n, double x[], int itask, double y[], void *data),
int problem_type, int n, int nev, int nconv, d_complex eigvals[],
double evecu[], double rwork[], void *data) {
int first, i, j;
double temp;
/*
* The following computations assume that complex conjugate
* pairs of eigenvalues are stored consecutively and that the
* imaginary part of the first eigenvalue is > 0, as
* guaranteed by arpack_general.
* The computed residuals are stored in rwork[0:nconv-1].
* This example actually uses the general case only, but
* contains also the standard case if the user wants to
* compute residuals for his own standard problems.
*/
first = 1;
if (problem_type == 0) { /* standard problem */
/*
* Compute the residual norm
*
* || A * x - lambda * x ||
*
* for the n_acc accurately computed eigenvalues and
* eigenvectors.
*/
for (i = 0; i < nconv; i++) {
if (eigvals[i].im == 0.0) {
/*
* Ritz value is real.
*/
/* Copy eigenvectors into ax */
for (j = 0; j < n; j++) {
rwork[nev + j] = evecu[j * (nev + 1) + i];
}
fcn(n, &rwork[nev], IMSL_ARPACK_A_X, &rwork[nev + n],
data);
for (j = 0; j < n; j++) {
rwork[nev + n + j] -= eigvals[i].re * rwork[nev + j];
}
rwork[i] = imsl_d_vector_norm(n, &rwork[nev + n], 0);
if (fabs(eigvals[i].re) != 0.0) {
rwork[i] /= fabs(eigvals[i].re);
}
}
else if (first) {
/*
* Compute real part of A * x - lambda * x,
*
* A * x_re - lambda_re * x_re + lambda_im * x_im
*/
for (j = 0; j < n; j++) {
rwork[nev + j] = evecu[j * (nev + 1) + i];
}
fcn(n, &rwork[nev], IMSL_ARPACK_A_X, &rwork[nev + n],
data);
for (j = 0; j < n; j++) {
rwork[nev + n + j] -= eigvals[i].re *
evecu[j * (nev + 1) + i];
rwork[nev + n + j] += eigvals[i].im *
evecu[j * (nev + 1) + i + 1];
}
/*
* Compute
* || A * x_re - lambda_re * x_re + lambda_im * x_im ||
*/
rwork[i] = imsl_d_vector_norm(n, &rwork[nev + n], 0);
/*
* Compute imaginary part of A * x - lambda * x,
*
* A * x_im - lambda_im * x_re - lambda_re * x_im
*/
for (j = 0; j < n; j++) {
rwork[nev + j] = evecu[j * (nev + 1) + i + 1];
}
fcn(n, &rwork[nev], IMSL_ARPACK_A_X, &rwork[nev + n],
data);
for (j = 0; j < n; j++) {
rwork[nev + n + j] -= eigvals[i].im *
evecu[j * (nev + 1) + i];
rwork[nev + n + j] -= eigvals[i].re *
evecu[j * (nev + 1) + i + 1];
}
/*
* Compute || A * x - lambda * x ||
*/
rwork[i] = hypot(rwork[i], imsl_d_vector_norm(n,
&rwork[nev + n], 0));
/*
* Compute res := || A*x - lambda * x || / || lambda ||
*/
temp = hypot(eigvals[i].re, eigvals[i].im);
if (temp != 0.0) {
rwork[i] /= temp;
}
/*
* Take into account that
* || A * x - lambda * x || =
* || conj(A * x - lambda * x) ||
*/
rwork[i + 1] = rwork[i];
first = 0;
}
else {
first = 1;
}
}
}
else { /* generalized problem */
/*
* Compute the residual norm
*
* || A * x - lambda * B * x || / | lambda |
*
* for the n_acc accurately computed eigenvalues and
* eigenvectors.
*/
for (i = 0; i < nconv; i++) {
if (eigvals[i].im == 0.0) {
/*
* Ritz value is real.
*/
/* Copy eigenvectors into rwork */
for (j = 0; j < n; j++) {
rwork[nev + j] = evecu[j * (nev + 1) + i];
}
fcn(n, &rwork[nev], IMSL_ARPACK_A_X, &rwork[nev + n],
data);
fcn(n, &rwork[nev], IMSL_ARPACK_B_X, &rwork[nev + 2 * n],
data);
for (j = 0; j < n; j++) {
rwork[nev + n + j] -= eigvals[i].re *
rwork[nev + 2 * n + j];
}
rwork[i] = imsl_d_vector_norm(n, &rwork[nev + n], 0);
if (fabs(eigvals[i].re) != 0.0) {
rwork[i] /= fabs(eigvals[i].re);
}
}
else if (first) {
/*
* Ritz value is complex.
*/
/*
* Compute real part of A * x - lambda * B * x,
*
* A * x_re - lambda_re * B * x_re +
* lambda_im * B * x_im
*/
for (j = 0; j < n; j++) {
rwork[nev + j] = evecu[j * (nev + 1) + i];
}
fcn(n, &rwork[nev], IMSL_ARPACK_A_X, &rwork[nev + n],
data);
fcn(n, &rwork[nev], IMSL_ARPACK_B_X, &rwork[nev + 2 * n],
data);
for (j = 0; j < n; j++) {
rwork[nev + n + j] -= eigvals[i].re *
rwork[nev + 2 * n + j];
}
for (j = 0; j < n; j++) {
rwork[nev + j] = evecu[j * (nev + 1) + i + 1];
}
fcn(n, &rwork[nev], IMSL_ARPACK_B_X, &rwork[nev + 2 * n],
data);
for (j = 0; j < n; j++) {
rwork[nev + n + j] += eigvals[i].im *
rwork[nev + 2 * n + j];
}
/*
* Compute
* || A * x_re - lambda_re * B * x_re +
* lambda_im * B * x_im ||
*/
rwork[i] = imsl_d_vector_norm(n, &rwork[nev + n], 0);
/*
* Compute imaginary part of A*x - lambda * B * x,
*
* A * x_im - lambda_im * B * x_re
* - lambda_re * B * x_im
*/
for (j = 0; j < n; j++) {
rwork[nev + j] = evecu[j * (nev + 1) + i + 1];
}
fcn(n, &rwork[nev], IMSL_ARPACK_A_X, &rwork[nev + n],
data);
fcn(n, &rwork[nev], IMSL_ARPACK_B_X, &rwork[nev + 2 * n],
data);
for (j = 0; j < n; j++) {
rwork[nev + n + j] -= eigvals[i].re *
rwork[nev + 2 * n + j];
}
for (j = 0; j < n; j++) {
rwork[nev + j] = evecu[j * (nev + 1) + i];
}
fcn(n, &rwork[nev], IMSL_ARPACK_B_X, &rwork[nev + 2 * n],
data);
for (j = 0; j < n; j++) {
rwork[nev + n + j] -= eigvals[i].im *
rwork[nev + 2 * n + j];
}
/*
* Compute || A * x - lambda * x ||
*/
rwork[i] = hypot(rwork[i], imsl_d_vector_norm(n,
&rwork[nev + n], 0));
/*
* Compute res := || A*x - lambda * x || / || lambda ||
*/
temp = hypot(eigvals[i].re, eigvals[i].im);
if (temp != 0.0) {
rwork[i] /= temp;
}
/*
* Take into account that
* || A * x - lambda * x || =
* || conj(A * x - lambda * x) ||
*/
rwork[i + 1] = rwork[i];
first = 0;
}
else {
first = 1;
}
}
}
}
Output
Number of requested eigenvalues : 6
Number of accurate (converged) eigenvalues : 6
Largest magnitude eigenvalues, real shift
=========================================
Eigenvalues (Real, Imag) Relative residuals
( 0.50000000, 0.59581177) 0.00000000
( 0.50000000, -0.59581177) 0.00000000
( 0.50000000, 0.63311769) 0.00000000
( 0.50000000, -0.63311769) 0.00000000
( 0.50000000, 0.55827553) 0.00000000
( 0.50000000, -0.55827553) 0.00000000
Number of requested eigenvalues : 6
Number of accurate (converged) eigenvalues : 6
Largest magnitude eigenvalues, imaginary shift
==============================================
Eigenvalues (Real, Imag) Relative residuals
( 0.50000000, 0.59581177) 0.00000000
( 0.50000000, -0.59581177) 0.00000000
( 0.50000000, 0.63311769) 0.00000000
( 0.50000000, -0.63311769) 0.00000000
( 0.50000000, 0.55827553) 0.00000000
( 0.50000000, -0.55827553) 0.00000000
Warning Errors
IMSL_ARPACK_MAX_ITER_REACHED | The maximum number of iterations has been reached. All possible eigenvalues have been found. Variable "n_acc" returns the number of wanted converged Ritz values. |
IMSL_ARPACK_NO_SHIFTS_APPLIED | No shifts could be applied during a cycle of the implicitly restarted Arnoldi iteration. One possibility is to increase the size of "ncv" = # relative to "nev" = #. |
Fatal Errors
IMSL_START_VECTOR_ZERO | The starting vector "xguess" is zero. Use a non-zero vector instead. |
IMSL_UNABLE_TO_BUILD_ARNOLDI | The algorithm was not able to build an Arnoldi factorization. The size of the current Arnoldi factorization is #. Use of a different starting vector "xguess" may help. |
IMSL_QR_CONVERGENCE_ERROR | The iteration for an eigenvalue failed to converge during the processing of the implicitly restarted Arnoldi method. |
IMSL_QR_CONVERGENCE_ERROR_1 | The iteration for an eigenvalue failed to converge during the postprocessing phase of the implicitly restarted Arnoldi method. |
IMSL_SCHUR_FORM_REORD_ERROR | Reordering of the Schur form failed because some eigenvalues are too close to separate (the problem is very ill-conditioned). |
IMSL_EIGVEC_COMPUTATION_ERROR | The eigenvector computation during the postprocessing phase of the implicitly restarted Arnoldi method failed. |
IMSL_SYMM_TRIDIAG_QL_QR_ERROR | The eigenvalue calculation via the symmetric tridiagonal QL or QR algorithm during the post-processing phase of the implicitly restarted Arnoldi method failed. |
IMSL_ARPACK_NO_EIGVALS_FOUND | The implicitly restarted Arnoldi method did not find any eigenvalues to sufficient accuracy. Use of a different starting vector "xguess", a larger iteration number "itmax", a different number "ncv" of Arnoldi vectors or a different problem type and/or solve mode may help. |
IMSL_DIFF_N_CONV_RITZ_VALUES | The number of converged Ritz values computed by the iteratively restarted Arnoldi method differs from the number of converged Ritz values determined during the post-processing phase. |
IMSL_RAYLEIGH_DENOM_ZERO | The denominator of the Rayleigh quotient of the generalized eigenvector w, numbered #, is equal to zero. More specifically, the B-norm of the eigenvector, sqrt(ctrans(w)*B*w), is zero. |
