arpack_symmetric
Computes some of the eigenvalues and eigenvectors of the generalized real symmetric eigenvalue problem Ax = λBx using an implicitly restarted Arnoldi method (IRAM). The algorithm can also be used for the standard case B = I.
NOTE: Function arpack_symmetric is available in double precision only. |
Synopsis
#include <imsl.h>
double *imsl_d_arpack_symmetric (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 9 lists the following 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, 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.
Return Value
A pointer to the nev eigenvalues of the symmetric eigenvalue problem. To release this space, use imsl_free. If no value can be computed, then NULL is returned.
Synopsis with Optional Arguments
#include <imsl.h>
double *imsl_d_arpack_symmetric (void fcn(), int n, int nev,
IMSL_XGUESS, double xguess[],
IMSL_ITMAX, int itmax,
IMSL_TOLERANCE, double tol,
IMSL_SHIFT, double 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_LANCZOS_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, double 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, double shift (Input)
The shift value used in the shift-invert spectral transformations.
Default: shift = 0.
The shift value used in the shift-invert spectral transformations.
Default: shift = 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_ALGEBRAIC | Compute algebraically largest eigenvalues. |
IMSL_ARPACK_SMALLEST_ALGEBRAIC | Compute algebraically smallest eigenvalues. |
IMSL_ARPACK_LARGEST_MAGNITUDE | Compute eigenvalues of largest magnitude. |
IMSL_ARPACK_SMALLEST_MAGNITUDE | Compute eigenvalues of smallest magnitude. |
IMSL_ARPACK_BOTH_ENDS | Compute eigenvalues from both ends of the spectrum. |
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_ALGEBRAIC.
Default: eigval_loc = IMSL_ARPACK_LARGEST_ALGEBRAIC.
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_BUCKLING | IMSL_ARPACK_GENERALIZED | y = Ax, y = (A ‑ σB)-1x |
IMSL_ARPACK_CAYLEY | IMSL_ARPACK_GENERALIZED | y = Ax, y = Bx, y = (A ‑ σB)-1x |
Default: mode = IMSL_ARPACK_REGULAR.
IMSL_NUM_LANCZOS_VECTORS, int ncv (Input)
The number of Lanczos vectors generated in each iteration of the Arnoldi method. It is required that nev +1 <= ncv <= n. A value ncv >= min(2*nev, n) is recommended.
The number of Lanczos vectors generated in each iteration of the Arnoldi method. It is required that nev +1 <= ncv <= n. A value ncv >= min(2*nev, n) is recommended.
Default: ncv = min(2*nev, 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 containing the B-orthonormalized eigenvectors of the eigenvalue problem in the first n_acc columns. Typically, double *evec is declared, and &evec is used as an argument.
The address of a pointer to an array of size n × nev containing the B-orthonormalized eigenvectors of the eigenvalue problem in the first n_acc columns. Typically, double *evec is declared, and &evec is used as an argument.
IMSL_VECTORS_USER, double evecu[] (Output)
A user-defined array of size n × nev containing the B-orthonormalized eigenvectors of the eigenvalue problem in the first n_acc columns.
A user-defined array of size n × nev containing the B-orthonormalized eigenvectors of the eigenvalue problem in the first n_acc columns.
IMSL_EVECU_COL_DIM, int evecu_col_dim (Input)
The column dimension of evecu.
Default: evecu_col_dim = nev
The column dimension of evecu.
Default: evecu_col_dim = nev
IMSL_RETURN_USER, double evalu[] (Output)
An array of size nev containing the accurately computed eigenvalues in the first n_acc locations.
An array of size nev containing the accurately computed eigenvalues in the first n_acc locations.
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. See Passing Data to User-Supplied Functions in the introduction to this manual for more details.
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. See Passing Data to User-Supplied Functions in the introduction to this manual for more details.
NOTE: The possibility to supply user-data via IMSL_FCN_W_DATA is an important feature of arpack_symmetric. It allows the user to transfer problem-specific data to the algorithm without the need to define global data. See the documentation examples (Example 2, Example 3, and Example 4) for use cases. |
Description
Function imsl_d_arpack_symmetric, which is based on ARPACK subroutines DSAUPD and DSEUPD (see the ARPACK Users' Guide, Lehoucq et al. (1998)), computes selected eigenvalue-eigenvector pairs for generalized symmetric eigenvalue problems of the form
Ax = λBx.
Here, A and B are symmetric matrices. For B = I, the generalized problem reduces to the standard symmetric eigenvalue problem.
The ARPACK routine DSAUPD implements a variant of the Lanczos 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_symmetric 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 matrix-vector product
* y = inv(A - sigma * B) * x.
* 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 y via backsolves.
*
* Example:
* If an LU factorization of A - sigma * B exists, then
* A - sigma * B = P * L * U,
* P a permutation matrix. Vector y can then be determined
* as solution of the linear system
* L * U * y = 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_symmetric first determines the eigenvalues for the problem specified by optional arguments IMSL_EIG_SOLVE_MODE and IMSL_EIG_PROBLEM_TYPE.
Table 13 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 |
IMSL_ARPACK_SHIFT_INVERT | IMSL_ARPACK_GENERALIZED | (A ‑ σB)-1B |
IMSL_ARPACK_BUCKLING | IMSL_ARPACK_GENERALIZED | (A ‑ σB)-1A |
IMSL_ARPACK_CAYLEY | IMSL_ARPACK_GENERALIZED | (A ‑ σB)-1(A + σB) |
Note that the eigenvalue location defined by optional argument IMSL_EIGVAL_LOCATION always refers to the matrices of Table 13.
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, imsl_d_arpack_symmetric internally transforms the eigenvalues back to the eigenvalues of the original problem Ax = λBx or Ax = λx and computes eigenvectors, if required.
Besides matrices A and B always being symmetric, the modes for the generalized eigenproblem require some additional matrix properties summarized in Table 14:
mode | Matrix properties |
|---|---|
IMSL_ARPACK_REGULAR_INVERSE | B positive definite |
IMSL_ARPACK_SHIFT_INVERT | B positive semi-definite |
IMSL_ARPACK_BUCKLING | A positive semi-definite, B indefinite |
IMSL_ARPACK_CAYLEY | B positive semi-definite |
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
Eigenvalues and eigenfunctions of the Laplacian operator
defined by
‑Δu = λu
on the unit square [0,1] × [0,1] with zero Dirichlet boundary values, are approximated.
The full set of eigenvalues and their eigenfunctions are given by the sequence
where m,n are positive integers.
This provides a check on the accuracy of the numerical results. Equally spaced divided differences on the unit square are used to approximate ‑ Δu. A few eigenvalues of smallest magnitude and their eigenvectors are requested. The problem reduces to a symmetric matrix eigenvalue computation. The user function code provides the second order divided difference operator applied to an input vector under task IMSL_ARPACK_A_X.
#include <stdio.h>
#include <stdlib.h>
#include <imsl.h>
static void fcn(int n, double x[], int iact, double result[]);
static void tx(int nx, double x[], double y[]);
static void ax(int nx, double v[], double w[]);
static void daxpy(int n, double alpha, double x[], double y[]);
int main() {
int n, nx, nev, n_acc, i, j;
double *eigvals = NULL, *evecu = NULL, *res = NULL;
nx = 10;
n = nx * nx;
nev = 5;
/* Allocate memory for eigenvectors */
evecu = (double *)malloc(n * nev * sizeof(double));
/* Allocate memory for auxiliary array */
res = (double *)malloc((2 * n + nev) * sizeof(double));
eigvals = imsl_d_arpack_symmetric(fcn, n, nev,
IMSL_EIGVAL_LOCATION, IMSL_ARPACK_SMALLEST_MAGNITUDE,
IMSL_NUM_ACCURATE_EIGVALS, &n_acc,
IMSL_VECTORS_USER, evecu,
0);
printf("Number of requested eigenvalues : %d\n", nev);
printf("Number of accurate (converged) eigenvalues : %d\n", n_acc);
for (i = 0; i < n_acc; i++) {
/*
* Compute the residual norm || A * x - lambda * x || for the
* n_acc accurately computed eigenvalues and eigenvectors.
*/
/* Compute A * x - lambda * x */
for (j = 0; j < n; j++) {
res[nev + j] = evecu[j * nev + i];
}
ax(nx, &res[nev], &res[nev + n]);
for (j = 0; j < n; j++) {
res[nev + n + j] -= eigvals[i] * res[nev + j];
}
/* Compute relative residuals */
res[i] = imsl_d_vector_norm(n, &res[nev + n], 0);
if (eigvals[i] != 0.0) {
res[i] /= eigvals[i];
}
}
/*
* Display eigenvalues and residuals
*/
printf("\n Smallest Laplacian eigenvalues\n");
printf("%14s%25s\n", "Eigenvalues", "Relative residuals");
for (i = 0; i < n_acc; i++) {
printf("%14.8lf%20.8lf\n", eigvals[i], res[i]);
}
/* Print first 2D Laplacian eigenfunction at Grid Points */
for (j = 0; j < n; j++) {
res[j] = evecu[j * nev];
}
imsl_d_write_matrix("First 2D Laplacian Eigenfunction at Grid Points",
nx, nx, res, 0);
if (eigvals)
imsl_free(eigvals);
if (evecu)
free(evecu);
if (res)
free(res);
}
static void fcn(int n, double x[], int itask, double y[])
{
int nx = 10; /* n = nx * nx */
switch (itask) {
case IMSL_ARPACK_PREPARE:
/* Nothing to prepare, but formally handle this case */
break;
case IMSL_ARPACK_A_X:
ax(nx, x, y);
break;
default:
imsl_set_user_fcn_return_flag(1);
break;
}
}
/*
* Matrix-vector function
*
* The matrix used is the 2 dimensional discrete Laplacian on the
* unit square with zero Dirichlet boundary condition.
*
* Computes y <- A * x, where A is the nx*nx by nx*nx block
* tridiagonal matrix
*
* | T -I |
* | -I T -I |
* A = | -I T |
* | ... -I |
* | -I T |
*
* Function tx() is called to compute y <- T * x.
*/
static void ax(int nx, double x[], double y[]) {
int i, j, lo, n2;
double h2;
tx(nx, x, y);
daxpy(nx, -1.0, &x[nx], y);
for (j = 2; j <= nx - 1; j++) {
lo = (j - 1)*nx;
tx(nx, &x[lo], &y[lo]);
daxpy(nx, -1.0, &x[lo - nx], &y[lo]);
daxpy(nx, -1.0, &x[lo + nx], &y[lo]);
}
lo = (nx - 1) * nx;
tx(nx, &x[lo], &y[lo]);
daxpy(nx, -1.0, &x[lo - nx], &y[lo]);
/*
* Scale the vector w by (1 / (h * h)), where h is the
* mesh size.
*/
n2 = nx * nx;
h2 = 1.0 / ((double)((nx + 1)*(nx + 1)));
for (i = 0; i < n2; i++) {
y[i] /= h2;
}
}
static void tx(int nx, double x[], double y[]) {
int j;
double dd, dl, du;
/*
* Compute the matrix vector multiplication y <- T * x
* where T is an nx by nx tridiagonal matrix with dd on the
* diagonal, dl on the subdiagonal, and du on the superdiagonal.
*/
dd = 4.0;
dl = -1.0;
du = -1.0;
y[0] = dd * x[0] + du * x[1];
for (j = 2; j <= nx - 1; j++) {
y[j - 1] = dl * x[j - 2] + dd * x[j - 1] + du * x[j];
}
y[nx - 1] = dl * x[nx - 2] + dd * x[nx - 1];
}
static void daxpy(int n, double alpha, double x[], double y[]) {
int i;
/*
* Compute y <- alpha * x + y
*/
for (i = 0; i < n; i++) {
y[i] += alpha * x[i];
}
}
Output
Number of requested eigenvalues : 5
Number of accurate (converged) eigenvalues : 5
Smallest Laplacian eigenvalues
Eigenvalues Relative residuals
19.60540077 0.00000000
48.21934544 0.00000000
48.21934544 0.00000000
76.83329011 0.00000000
93.32640277 0.00000000
First 2D Laplacian Eigenfunction at Grid Points
1 2 3 4 5
1 0.0144 0.0277 0.0387 0.0466 0.0507
2 0.0277 0.0531 0.0743 0.0894 0.0973
3 0.0387 0.0743 0.1038 0.1250 0.1360
4 0.0466 0.0894 0.1250 0.1504 0.1637
5 0.0507 0.0973 0.1360 0.1637 0.1781
6 0.0507 0.0973 0.1360 0.1637 0.1781
7 0.0466 0.0894 0.1250 0.1504 0.1637
8 0.0387 0.0743 0.1038 0.1250 0.1360
9 0.0277 0.0531 0.0743 0.0894 0.0973
10 0.0144 0.0277 0.0387 0.0466 0.0507
6 7 8 9 10
1 0.0507 0.0466 0.0387 0.0277 0.0144
2 0.0973 0.0894 0.0743 0.0531 0.0277
3 0.1360 0.1250 0.1038 0.0743 0.0387
4 0.1637 0.1504 0.1250 0.0894 0.0466
5 0.1781 0.1637 0.1360 0.0973 0.0507
6 0.1781 0.1637 0.1360 0.0973 0.0507
7 0.1637 0.1504 0.1250 0.0894 0.0466
8 0.1360 0.1250 0.1038 0.0743 0.0387
9 0.0973 0.0894 0.0743 0.0531 0.0277
10 0.0507 0.0466 0.0387 0.0277 0.0144
Example 2
In this example, the eigenvalues and eigenfunctions of the 1D Laplacian operator
on the unit interval [0,1] with boundary conditions u(0) =u(1) = 0 are approximated. Equally spaced divided differences are used for the operator, which yields a symmetric tri-diagonal matrix. The eigenvalues and (normed) eigenfunctions are
This example shows that using inverse iteration for approximating the largest reciprocals of eigenvalues is more efficient than directly computing the smallest magnitude eigenvalues by matrix-vector products.
By using mode IMSL_ARPACK_SHIFT_INVERT, the algorithm first computes the largest eigenvalues of the shift-inverse matrix (A ‑ σI)-1, here with σ = 0. These eigenvalues are then transformed back to the smallest eigenvalues of A ‑ σI = A, a positive definite matrix. When user-defined function fcn is entered with task IMSL_ARPACK_PREPARE, the LU factorization of the shifted matrix A ‑ σI is computed. When fcn is later entered with task IMSL_ARPACK_INV_SHIFT_X, the LU factorization is available to efficiently compute y = (A ‑ σI)-1x = A-1x via LUy = x.
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <imsl.h>
static void ax(int nx, double x[], double y[]);
static void fcn_w_data(int n, double x[], int itask, double y[],
void *data);
typedef struct {
double shift;
double *band_matrix;
int *ipvt;
double *factor;
} imsl_arpack_data;
int main() {
int n, nev, ncv, n_acc, i, j;
int *ipvt = NULL;
double shift = 0.0;
double *a_matrix = NULL, *factor = NULL, *eigvals = NULL;
double *evecu = NULL, *res = NULL;
imsl_arpack_data usr_data;
n = 100;
nev = 4;
ncv = 10;
shift = 0.0;
/* Allocate memory for eigenvectors */
evecu = (double *)malloc(n * nev * sizeof(double));
/* Allocate arrays needed in the LU factorization */
ipvt = (int *)malloc(n * sizeof(int));
a_matrix = (double *)malloc(3 * n * sizeof(double));
factor = (double *)malloc(4 * n * sizeof(double));
/* Allocate memory for auxiliary array */
res = (double *)malloc((2 * n + nev) * sizeof(double));
if (!evecu || !ipvt || !a_matrix || !factor || !res) {
printf("Memory allocation error\n");
goto FREE_SPACE;
}
usr_data.band_matrix = a_matrix;
usr_data.ipvt = ipvt;
usr_data.factor = factor;
usr_data.shift = shift;
eigvals = imsl_d_arpack_symmetric(NULL, n, nev,
IMSL_EIG_SOLVE_MODE, IMSL_ARPACK_SHIFT_INVERT,
IMSL_NUM_LANCZOS_VECTORS, ncv,
IMSL_NUM_ACCURATE_EIGVALS, &n_acc,
IMSL_FCN_W_DATA, fcn_w_data, &usr_data,
IMSL_VECTORS_USER, evecu,
IMSL_SHIFT, shift,
0);
printf("Number of requested eigenvalues : %d\n", nev);
printf("Number of accurate (converged) eigenvalues : %d\n", n_acc);
for (i = 0; i < n_acc; i++) {
/*
* Compute the residual norm || A * x - lambda * x || for the
* n_acc accurately computed eigenvalues and eigenvectors.
*/
/* Compute A * x - lambda * x */
for (j = 0; j < n; j++) {
res[nev + j] = evecu[j * nev + i];
}
ax(n, &res[nev], &res[nev + n]);
for (j = 0; j < n; j++) {
res[nev + n + j] -= eigvals[i] * res[nev + j];
}
/* Compute relative residuals */
res[i] = imsl_d_vector_norm(n, &res[nev + n], 0);
if (fabs(eigvals[i]) != 0.0) {
res[i] /= fabs(eigvals[i]);
}
}
/*
* Display eigenvalues and residuals
*/
printf("\n Largest Laplacian eigenvalues near zero shift\n");
printf("%14s%25s\n", "Eigenvalues", "Relative residuals");
for (i = 0; i < n_acc; i++) {
printf("%14.8lf%20.8lf\n", eigvals[i], res[i]);
}
FREE_SPACE:
if (eigvals)
imsl_free(eigvals);
if (ipvt)
free(ipvt);
if (a_matrix)
free(a_matrix);
if (factor)
free(factor);
if (res)
free(res);
if (evecu)
free(evecu);
}
static void fcn_w_data(int n, double x[], int itask, double y[],
void *data)
{
int j;
int *ipvt = NULL;
double shift, h2;
double *a_matrix = NULL, *factor = NULL;
imsl_arpack_data *usr_data = (imsl_arpack_data *)data;
shift = usr_data->shift;
a_matrix = usr_data->band_matrix;
ipvt = usr_data->ipvt;
factor = usr_data->factor;
switch (itask) {
case IMSL_ARPACK_PREPARE:
/* Create symmetric tridiagonal matrix in band storage format */
h2 = 1.0 / (((double)(n + 1)) * (n + 1));
for (j = 1; j <= n; j++) {
a_matrix[j - 1] = -1.0 / h2;
a_matrix[n + j - 1] = 2.0 / h2 - shift;
a_matrix[2 * n + j - 1] = a_matrix[j - 1];
}
/* Compute LU factorization of tridiagonal matrix */
imsl_d_lin_sol_gen_band(n, a_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_INV_SHIFT_X:
/* Solve (A - shift * I) y = x */
imsl_d_lin_sol_gen_band(n, NULL, 1, 1, x,
IMSL_FACTOR_USER, ipvt, factor,
IMSL_RETURN_USER, y,
IMSL_SOLVE_ONLY,
0);
if (imsl_error_type() != 0) {
imsl_set_user_fcn_return_flag(2);
}
break;
default:
imsl_set_user_fcn_return_flag(3);
break;
}
}
/*
* Matrix-vector function
* The matrix is the 1 dimensional discrete Laplacian on the
* interval [0, 1] with zero Dirichlet boundary condition.
*/
static void ax(int n, double x[], double y[]) {
int i;
double h2;
y[0] = 2.0 * x[0] - x[1];
for (i = 1; i <= n - 2; i++) {
y[i] = -x[i - 1] + 2.0 * x[i] - x[i + 1];
}
y[n - 1] = -x[n - 2] + 2.0 * x[n - 1];
/*
* Scale the vector y by (1 / h ^ 2).
*/
h2 = 1.0 / (((double)(n + 1)) * (n + 1));
for (i = 0; i < n; i++) {
y[i] /= h2;
}
}
Output
Number of requested eigenvalues : 4
Number of accurate (converged) eigenvalues : 4
Largest Laplacian eigenvalues near zero shift
Eigenvalues Relative residuals
9.86880868 0.00000000
39.46568728 0.00000000
88.76200274 0.00000000
157.71006404 0.00000000
Example 3
In this example, a generalized problem is solved using the regular inverse mode. The problem comes from using equally spaced linear finite element test functions to approximate eigenvalues and eigenfunctions of the 1D Laplacian operator ‑Δ ≡ ‑d2/dx2, defined by
‑Δu = λu,
on the unit interval [0,1] with boundary conditions u(0) =u(1)=0. This is Example 2 but solved using finite elements and the eigenvalues scaled by 1/π2 so that λn = n2, n = 1,2,… . In matrix notation, we have the matrix problem Ax = λBx, with both A and B tri-diagonal and symmetric. The matrix B is non-singular.
The user function fcn requires the solution of a tri-diagonal system of linear equations with the input vector x as right-hand side, By = x. When fcn is entered with task IMSL_ARPACK_PREPARE, the LU factorization of matrix B is computed. When fcn is later called with task IMSL_ARPACK_INV_SHIFT_X, the LU factorization is available to efficiently compute y = B-1x.
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <imsl.h>
static void ax(int nx, 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);
typedef struct {
double *band_matrix;
int *ipvt;
double *factor;
} imsl_arpack_data;
int main() {
int n, nev, ncv, n_acc, i, j;
int *ipvt = NULL;
double *a_matrix = NULL, *factor = NULL, *eigvals = NULL;
double *evecu = NULL, *res = NULL;
imsl_arpack_data usr_data;
n = 100;
nev = 4;
ncv = 10;
/* Allocate memory for eigenvectors */
evecu = (double *)malloc(n * nev * sizeof(double));
/* Allocate arrays needed in the LU factorization */
ipvt = (int *)malloc(n * sizeof(int));
a_matrix = (double *)malloc(3 * n * sizeof(double));
factor = (double *)malloc(4 * n * sizeof(double));
/* Allocate memory for auxiliary array */
res = (double *)malloc((nev + 3 * n) * sizeof(double));
if (!evecu || !ipvt || !a_matrix || !factor || !res) {
printf("Memory allocation error\n");
goto FREE_SPACE;
}
usr_data.band_matrix = a_matrix;
usr_data.ipvt = ipvt;
usr_data.factor = factor;
eigvals = imsl_d_arpack_symmetric(NULL, n, nev,
IMSL_EIGVAL_LOCATION, IMSL_ARPACK_SMALLEST_MAGNITUDE,
IMSL_EIG_PROBLEM_TYPE, IMSL_ARPACK_GENERALIZED,
IMSL_EIG_SOLVE_MODE, IMSL_ARPACK_REGULAR_INVERSE,
IMSL_NUM_LANCZOS_VECTORS, ncv,
IMSL_NUM_ACCURATE_EIGVALS, &n_acc,
IMSL_FCN_W_DATA, fcn_w_data, &usr_data,
IMSL_VECTORS_USER, evecu,
0);
printf("Number of requested eigenvalues : %d\n", nev);
printf("Number of accurate (converged) eigenvalues : %d\n", n_acc);
for (i = 0; i < n_acc; i++) {
/*
* Compute the residual norm || A * x - lambda * B * x ||
* for the n_acc accurately computed eigenvalues and
* eigenvectors.
*/
/* Compute A * x - lambda * B * x */
for (j = 0; j < n; j++) {
res[j + nev] = evecu[j * nev + i];
}
ax(n, &res[nev], &res[nev + n]);
bx(n, &res[nev], &res[nev + 2 * n]);
for (j = 0; j < n; j++) {
res[nev + n + j] -= eigvals[i] * res[nev + 2 * n + j];
}
/* Compute relative residuals */
res[i] = imsl_d_vector_norm(n, &res[nev + n], 0);
if (fabs(eigvals[i]) != 0.0) {
res[i] /= fabs(eigvals[i]);
}
}
/*
* Display eigenvalues and residuals
*/
printf("\n Smallest Laplacian eigenvalues\n");
printf("%14s%25s\n", "Eigenvalues", "Relative residuals");
for (i = 0; i < n_acc; i++) {
printf("%14.8lf%20.8lf\n", eigvals[i], res[i]);
}
FREE_SPACE:
if (eigvals)
imsl_free(eigvals);
if (ipvt)
free(ipvt);
if (a_matrix)
free(a_matrix);
if (factor)
free(factor);
if (evecu)
free(evecu);
if (res)
free(res);
}
static void fcn_w_data(int n, double x[], int itask, double y[],
void *data)
{
int j;
int *ipvt = NULL;
double h, r1, r2;
double *b_matrix = NULL, *factor = NULL;
double pi = 3.1415926535897932384;
imsl_arpack_data *usr_data = (imsl_arpack_data *)data;
b_matrix = usr_data->band_matrix;
ipvt = usr_data->ipvt;
factor = usr_data->factor;
switch (itask) {
case IMSL_ARPACK_PREPARE:
/* Create symmetric tridiagonal matrix B */
h = pi / (n + 1);
r1 = (2.0 / 3.0) * h;
r2 = (1.0 / 6.0) * h;
for (j = 0; j < n; j++) {
b_matrix[j] = r2;
b_matrix[n + j] = r1;
b_matrix[2 * n + j] = r2;
}
/* Compute LU factorization of tridiagonal matrix B */
imsl_d_lin_sol_gen_band(n, b_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_B_X:
/* Solve B * y = x */
imsl_d_lin_sol_gen_band(n, NULL, 1, 1, x,
IMSL_FACTOR_USER, ipvt, factor,
IMSL_RETURN_USER, y,
IMSL_SOLVE_ONLY,
0);
if (imsl_error_type() != 0) {
imsl_set_user_fcn_return_flag(2);
}
break;
default:
imsl_set_user_fcn_return_flag(3);
break;
}
}
/*
* Matrix-vector function
*
* The matrix is the 1 - dimensional mass matrix
* on the interval [0, 1].
*/
static void bx(int n, double x[], double y[]) {
int j;
double h;
double pi = 3.1415926535897932384;
y[0] = 4.0 * x[0] + x[1];
for (j = 1; j < n - 1; j++) {
y[j] = x[j - 1] + 4.0 * x[j] + x[j + 1];
}
y[n - 1] = x[n - 2] + 4.0 * x[n - 1];
/*
* Scale the vector w by h.
*/
h = pi / ((n + 1) * 6.0);
for (j = 0; j < n; j++) {
y[j] *= h;
}
}
/*
* Matrix-vector function
*
* The matrix used is the stiffness matrix obtained from the finite
* element discretization of the 1 - dimensional discrete Laplacian
* on the interval [0, 1] with zero Dirichlet boundary condition
* using piecewise linear elements.
*/
static void ax(int n, double x[], double y[]) {
int j;
double h;
double pi = 3.1415926535897932384;
y[0] = 2.0 * x[0] - x[1];
for (j = 1; j < n - 1; j++) {
y[j] = -x[j - 1] + 2.0 * x[j] - x[j + 1];
}
y[n - 1] = -x[n - 2] + 2.0 * x[n - 1];
/*
* Scale the vector w by (1 / h).
*/
h = pi / (n + 1);
for (j = 0; j < n; j++) {
y[j] /= h;
}
}
Output
Number of requested eigenvalues : 4
Number of accurate (converged) eigenvalues : 4
Smallest Laplacian eigenvalues
Eigenvalues Relative residuals
1.00008063 0.00000000
4.00129018 0.00000000
9.00653261 0.00000000
16.02065092 0.00000000
Example 4
In this example, a generalized problem Ax = λBx, with both A and B tri-diagonal and symmetric and B non-singular, is solved using the shift-invert spectral transformation mode. The problem stems from using equally spaced linear finite element test functions to approximate eigenvalues and eigenfunctions of the 1D Laplacian operator ‑Δ≡-d2/dx2, defined by
-Δu = λu,
on the unit interval [0,1] with boundary conditions u(0) =u(1)=0. This is Example 2, but solved using finite elements.
The algorithm iteratively computes the eigenvalues v of largest magnitude of the transformed system
(A-σB)-1 Bx = vx
where
and the shift parameter σ = 0.
These eigenvalues are then transformed back to the eigenvalues of smallest magnitude of the original system, λ = 1/v, and associated eigenvectors are determined.The user function fcn requires the solution of a tri-diagonal system of linear equations with the input vector x as right-hand side, (A - σB) y=x, where σ = 0. When fcn is entered with task IMSL_ARPACK_PREPARE, the LU factorization of matrix A is computed. When fcn is later called with task IMSL_ARPACK_INV_SHIFT_X, the LU factorization is available to efficiently compute y=(A - σB)-1x = A-1x.
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <imsl.h>
static void ax(int nx, 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);
typedef struct {
double *band_matrix;
int *ipvt;
double *factor;
double shift;
} imsl_arpack_data;
int main() {
int n, nev, ncv, n_acc, i, j;
int *ipvt = NULL;
double shift = 0.0;
double *a_matrix = NULL, *factor = NULL, *eigvals = NULL;
double *evecu = NULL, *res = NULL;
imsl_arpack_data usr_data;
n = 100;
nev = 4;
ncv = 10;
/* Allocate memory for eigenvectors */
evecu = (double *)malloc(n * nev * sizeof(double));
/* Allocate arrays needed in the LU factorization */
ipvt = (int *)malloc(n * sizeof(int));
a_matrix = (double *)malloc(3 * n * sizeof(double));
factor = (double *)malloc(4 * n * sizeof(double));
/* Allocate memory for auxiliary array */
res = (double *)malloc((nev + 3 * n) * sizeof(double));
if (!evecu || !ipvt || !a_matrix || !factor || !res) {
printf("Memory allocation error\n");
goto FREE_SPACE;
}
usr_data.band_matrix = a_matrix;
usr_data.ipvt = ipvt;
usr_data.factor = factor;
usr_data.shift = shift;
eigvals = imsl_d_arpack_symmetric(NULL, n, nev,
IMSL_EIGVAL_LOCATION, IMSL_ARPACK_LARGEST_MAGNITUDE,
IMSL_EIG_PROBLEM_TYPE, IMSL_ARPACK_GENERALIZED,
IMSL_EIG_SOLVE_MODE, IMSL_ARPACK_SHIFT_INVERT,
IMSL_SHIFT, shift,
IMSL_NUM_LANCZOS_VECTORS, ncv,
IMSL_NUM_ACCURATE_EIGVALS, &n_acc,
IMSL_FCN_W_DATA, fcn_w_data, &usr_data,
IMSL_VECTORS_USER, evecu,
0);
printf("Number of requested eigenvalues : %d\n", nev);
printf("Number of accurate (converged) eigenvalues : %d\n", n_acc);
for (i = 0; i < n_acc; i++) {
/*
* Compute the residual norm || A * x - lambda * B * x ||
* for the n_acc accurately computed eigenvalues and
* eigenvectors.
*/
/* Compute A * x - lambda * B * x */
for (j = 0; j < n; j++) {
res[j + nev] = evecu[j * nev + i];
}
ax(n, &res[nev], &res[nev + n]);
bx(n, &res[nev], &res[nev + 2 * n]);
for (j = 0; j < n; j++) {
res[nev + n + j] -= eigvals[i] * res[nev + 2 * n + j];
}
/* Compute relative residuals */
res[i] = imsl_d_vector_norm(n, &res[nev + n], 0);
if (fabs(eigvals[i]) != 0.0) {
res[i] /= fabs(eigvals[i]);
}
}
/*
* Display eigenvalues and residuals
*/
printf("\n Smallest Laplacian eigenvalues\n");
printf("%14s%25s\n", "Eigenvalues", "Relative residuals");
for (i = 0; i < n_acc; i++) {
printf("%14.8lf%20.8lf\n", eigvals[i], res[i]);
}
FREE_SPACE:
if (eigvals)
imsl_free(eigvals);
if (ipvt)
free(ipvt);
if (a_matrix)
free(a_matrix);
if (factor)
free(factor);
if (evecu)
free(evecu);
if (res)
free(res);
}
static void fcn_w_data(int n, double x[], int itask, double y[],
void *data)
{
int j;
int *ipvt = NULL;
double h, r1, r2, shift;
double *b_matrix = NULL, *factor = NULL;
imsl_arpack_data *usr_data = (imsl_arpack_data *)data;
b_matrix = usr_data->band_matrix;
ipvt = usr_data->ipvt;
factor = usr_data->factor;
shift = usr_data->shift;
switch (itask) {
case IMSL_ARPACK_PREPARE:
/* Create symmetric tridiagonal matrix (A - shift * B) */
h = 1.0 / (n + 1);
r1 = (2.0 / 3.0) * h;
r2 = (1.0 / 6.0) * h;
for (j = 1; j <= n; j++) {
b_matrix[j - 1] = -1.0 / h - shift * r2;
b_matrix[n + j - 1] = 2.0 / h - shift * r1;
b_matrix[2 * n + j - 1] = b_matrix[j - 1];
}
/* Compute LU factorization of tridiagonal matrix */
imsl_d_lin_sol_gen_band(n, b_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_B_X:
bx(n, x, y);
break;
case IMSL_ARPACK_INV_SHIFT_X:
/* Solve (A - shift * B) * y = x */
imsl_d_lin_sol_gen_band(n, NULL, 1, 1, x,
IMSL_FACTOR_USER, ipvt, factor,
IMSL_RETURN_USER, y,
IMSL_SOLVE_ONLY,
0);
if (imsl_error_type() != 0) {
imsl_set_user_fcn_return_flag(2);
}
break;
default:
imsl_set_user_fcn_return_flag(3);
break;
}
}
/*
* Matrix-vector function B*x
* The matrix used is the 1 - dimensional mass matrix
* on the interval [0, 1].
*/
static void bx(int n, double x[], double y[]) {
int j;
double h;
y[0] = 4.0 * x[0] + x[1];
for (j = 1; j < n - 1; j++) {
y[j] = x[j - 1] + 4.0 * x[j] + x[j + 1];
}
y[n - 1] = x[n - 2] + 4.0 * x[n - 1];
/*
* Scale the vector w by h.
*/
h = 1.0 / (6.0 * (n + 1));
for (j = 0; j < n; j++) {
y[j] *= h;
}
}
/*
* Matrix-vector function A*x
*
* The matrix is the finite element discretization of the
* 1 - dimensional discrete Laplacian on [0, 1] with zero
* Dirichlet boundary condition using piecewise linear
* elements.
*/
static void ax(int n, double x[], double y[]) {
int j;
double h;
y[0] = 2.0 * x[0] - x[1];
for (j = 1; j < n - 1; j++) {
y[j] = -x[j - 1] + 2.0 * x[j] - x[j + 1];
}
y[n - 1] = -x[n - 2] + 2.0 * x[n - 1];
/*
* Scale the vector w by (1 / h)
*/
h = 1.0 / ((double)(n + 1));
for (j = 0; j < n; j++) {
y[j] /= h;
}
}
Output
Number of requested eigenvalues : 4
Number of accurate (converged) eigenvalues : 4
Smallest Laplacian eigenvalues
Eigenvalues Relative residuals
9.87040017 0.00000000
39.49115121 0.00000000
88.89091388 0.00000000
158.11748683 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_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. |
