forked from HPAC/mr3smp
-
Notifications
You must be signed in to change notification settings - Fork 0
The MRRR algorithm for multi-core and SMP systems
License
SebastianAchilles/mr3smp
Â
Â
Folders and files
| Name | Name | Last commit message | Last commit date | |
|---|---|---|---|---|
 |  | |||
 |  | |||
 |  | |||
 |  | |||
 |  | |||
 |  | |||
 |  | |||
 |  | |||
 |  | |||
Repository files navigation
These files provide the routine 'mrrr'. Its purpose is to compute all
or a subset of eigenvalues and optionally eigenvectors of a symmetric
tridiagonal matrix. 'mrrr' is based on the algorithm of Multiple
Relatively Robust Representations (MRRR). The routine thereby targets
multi-core processors and SMP systems made out of them. The
implementation is based on LAPACK's routine 'dstemr'.
Building the archive libmrrr.a:
--------------------------------
1) In the folder INSTALL, depending on the compiler used,
replace the file 'make.inc' (which assumes GNU gcc)
with the corresponding file.
For example if you want to use Intel's compilers:
$cp ./INSTALL/make.inc.intel ./make.inc
If you do not find a file matching your compiler, then use any of
the files as a basis to edit.
2) Edit the file 'make.inc' according to your system. Especially, set
the variable INCLAPACK to zero if you plan to link to LAPACK and
BLAS, instead of compiling the necessary routines and including
them in the archive. (LAPACK 3.3 or equivalent is recommended.)
3) Run make.
Known issues:
-------------
- PROBLEM: Spinlocks are not supported
SOLUTION: Compile the library with SPINLOCK_SUPPORT=0 in 'make.inc'
- PROBLEM: The C99 feature of complex numbers is not supported
SOLUTION: Use a newer compiler or compile the library with
COMPLEX_SUPPORT=0 (default) in 'make.inc'. In the latter case the
dense Hermitian routines will not be available in the archive
Using libmrrr.a:
----------------
The folder EXAMPLES contains simple examples of how to use the routine
'mrrr' in C and Fortran code. Edit the 'Makefile' in these folders if
you do not use the GNU compilers and run 'make' to compile the
examples.
In general, the code that calls 'mrrr' needs to be linked to the
library 'libmrrr.a', which is created in the LIB folder.
Below are given the C and Fortran prototypes of the function 'mrrr'.
For more information please see 'INCLUDE/mrrr.h'.
######################################################################
# C function prototype: #
######################################################################
# #
# int mrrr(char *jobz, char *range, int *n, double *D, #
# double *E, double *vl, double *vu, int *il, int *iu, #
# int *tryrac, int *m, double *W, double *Z, int *ldz, #
# int *Zsupp); #
# #
# INPUTS: #
# ------- #
# jobz "N" or "n" - compute only eigenvalues #
# "V" or "v" - compute also eigenvectors #
# "C" or "c" - count the maximal number of #
# locally computed eigenvectors #
# range "A" or "a" - all #
# "V" or "v" - by interval: (VL,VU] #
# "I" or "i" - by index: IL-IU #
# n Matrix size #
# ldz Leading dimension of eigenvector matrix Z; #
# often equal to matrix size n #
# #
# INPUT + OUTPUT: #
# --------------- #
# D (double[n]) Diagonal elements of tridiagonal T. #
# (On output the array will be overwritten). #
# E (double[n]) Off-diagonal elements of tridiagonal T. #
# First n-1 elements contain off-diagonals, #
# the last element can have an abitrary value. #
# (On output the array will be overwritten.) #
# vl If range="V", lower bound of interval #
# (vl,vu], on output refined. #
# If range="A" or "I" not referenced as input. #
# On output the interval (vl,vu] contains ALL #
# the computed eigenvalues. #
# vu If range="V", upper bound of interval #
# (vl,vu], on output refined. #
# If range="A" or "I" not referenced as input. #
# On output the interval (vl,vu] contains ALL #
# the computed eigenvalues. #
# il If range="I", lower index (1-based indexing) of #
# the subset 'il' to 'iu'. #
# If range="A" or "V" not referenced as input. #
# On output the eigenvalues with index 'il' to 'iu'#
# are computed. #
# iu If range="I", upper index (1-based indexing) of #
# the subset 'il' to 'iu'. #
# If range="A" or "V" not referenced as input. #
# On output the eigenvalues with index 'il' to 'iu'#
# are computed. #
# tryrac 0 - do not try to achieve high relative accuracy.#
# 1 - relative accuracy will be attempted; #
# on output it is set to zero if high relative #
# accuracy is not achieved. #
# #
# OUTPUT: #
# ------- #
# m Number of eigenvalues and eigenvectors computed. #
# If jobz="C", 'm' will be set to the number of #
# eigenvalues/eigenvectors that will be computed. #
# W (double[n]) Eigenvalues #
# The first 'm' entries contain the eigenvalues. #
# Z Eigenvectors #
# (double[n*m]) Enough space must be provided to store the #
# vectors. 'm' should be bigger or equal #
# to 'n' for range="A" and 'iu-il+1' for range="I".#
# For range="V" the minimum 'm' can be queried #
# using jobz="C". #
# Zsupp Support of eigenvectors, which is given by #
# (double[2*n]) Z[2*i] to Z[2*i+1] for the i-th eigenvector #
# stored locally (1-based indexing). #
# #
######################################################################
######################################################################
# Fortran prototype: #
######################################################################
# #
# MRRR(JOBZ, RANGE, N, D, E, VL, VU, IL, IU, TRYRAC, M, W, #
# Z, LDZ, ZSUPP, INFO) #
# #
# CHARACTER JOBZ, RANGE #
# INTEGER N, IL, IU, TRYRAC, M, LDZ, ZSUPP(*), INFO #
# DOUBLE PRECISION D(*), D(*), VL, VU, W(*), Z(*,*) #
######################################################################
The number of threads created can be specified via the environment
variable PMR_NUM_THREADS. In case this variable is not
defined, then the routine will create as many threads as specified by
the variable DEFAULT_NUM_THREADS (which is set in 'mrrr.h').
ADDITIONAL ROUTINES FOR DENSE STANDARD AND GENERALIZED EIGENPROBLEMS
--------------------------------------------------------------------
For convinience there are also routines for the standard and
generalized symmetric and Hermitian eigenproblem added that make use
of the multi-threaded MRRR. The functionality is similar to LAPACK's
routines for this operations. The routines using packed storage should
be avoided whenever possible, since their performance is much worse.
The routines are (using a LAPACK like notation):
----------------------------------------------------------------------
Standard eigenvector problem A*x = lambda*x:
Computing all or a subset of eigenvalues, and optionally
eigenvectors
----------------------------------------------------------------------
dsyeig: A is symmetric
zheeig: A is Hermitian
dspeig: A is symmetric and stored in packed storage
zhpeig: A is Hermitian and stored in packed storage
----------------------------------------------------------------------
----------------------------------------------------------------------
Standard eigenvector problem A*x = lambda*B*x or A*B*x = lambda*x
or B*A*x = lambda*x:
Computing all or a subset of eigenvalues, and optionally
eigenvectors
----------------------------------------------------------------------
dsygeig: A and B are symmetric and B is definite
zhegeig: A and B are Hermitian and B is definite
dspgeig: A and B are symmetric, stored in packed storage,
and B is definite
zhpgeig: A and B are Hermitian, stored in packed storage,
and B is definite
----------------------------------------------------------------------
The dense standard eigenvalue problem:
--------------------------------------
For the reduction and back-transformation the routines call the
corresponding LAPACK routines. Therefore LAPACK is required to run
the routines. Simply code snippets of how to call the functions are
included in the examples folder.
The prototypes are similar to LAPACK's "dsyevx" and are as follows.
The symmetric case:
######################################################################
# C function prototype: #
######################################################################
# #
# int dsyeig(char *jobz, char *range, char *uplo, int *n, double *A, #
# int *lda, double *vl, double *vu, int *il, int *iu, #
# int *m, double *W, double *Z, int *ldz); #
# #
# Arguments: #
# ---------- #
# #
# INPUTS: #
# ------- #
# jobz "N" - compute only eigenvalues #
# "V" - compute also eigenvectors #
# range "A" - all #
# "V" - by interval: (VL,VU] #
# "I" - by index: IL-IU #
# uplo "L" - Upper triangle of A stored #
# "U" - Lower triangle of A stored #
# n Order of the matrix A #
# lda Leading dimension of matrix A; #
# often equal to matrix size n #
# ldz Leading dimension of eigenvector matrix Z; #
# often equal to matrix size n #
# #
# INPUT + OUTPUT: #
# --------------- #
# A On entry symmetric input matrix, stored #
# (double[lda*n]) in column major ordering. Depending on the #
# value of 'uplo' only the upper or lower #
# triangular part is referenced. #
# (On output the array will be overwritten). #
# vl If range="V", lower bound of interval #
# (vl,vu], on output refined. #
# If range="A" or "I" not referenced as input. #
# On output the interval (vl,vu] contains ALL #
# the computed eigenvalues. #
# vu If range="V", upper bound of interval #
# (vl,vu], on output refined. #
# If range="A" or "I" not referenced as input. #
# On output the interval (vl,vu] contains ALL #
# the computed eigenvalues. #
# il If range="I", lower index (1-based indexing) of #
# the subset 'il' to 'iu'. #
# If range="A" or "V" not referenced as input. #
# On output the eigenvalues with index il to iu #
# are computed. #
# iu If range="I", upper index (1-based indexing) of #
# the subset 'il' to 'iu'. #
# If range="A" or "V" not referenced as input. #
# On output the eigenvalues with index il to iu #
# are computed. #
# #
# OUTPUT: #
# ------- #
# m Number of eigenvalues and eigenvectors computed. #
# W (double[n]) Eigenvalues #
# The first 'm' entries contain the eigenvalues. #
# Z Eigenvectors. #
# (double[n*m]) Enough space must be provided to store the #
# vectors. 'm' should be bigger or equal #
# to 'n' for range="A" or "V" and 'iu-il+1' for #
# range="I". #
# #
# RETURN VALUE: #
# ------------- #
# 0 - Success #
# 1 - Wrong input parameter #
# 2 - Misc errors #
######################################################################
######################################################################
# Fortran prototype: #
######################################################################
# #
# DSYEIG(JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU, M, W, #
# Z, LDZ, INFO) #
# #
# CHARACTER JOBZ, RANGE, UPLO #
# INTEGER N, IL, IU, M, LDA, LDZ, INFO #
# DOUBLE PRECISION A(*,*), VL, VU, W(*), Z(*,*) #
# #
######################################################################
The symmetric case (packed storage):
######################################################################
# C function prototype: #
######################################################################
# #
# int dspeig(char *jobz, char *range, char *uplo, int *n, #
# double *AP, double *vl, double *vu, int *il, int *iu, #
# int *m, double *W, double *Z, int *ldz); #
# #
# Arguments: #
# ---------- #
# #
# INPUTS: #
# ------- #
# jobz "N" - compute only eigenvalues #
# "V" - compute also eigenvectors #
# range "A" - all #
# "V" - by interval: (VL,VU] #
# "I" - by index: IL-IU #
# uplo "L" - Upper triangle of A stored #
# "U" - Lower triangle of A stored #
# n Order of the matrix A #
# ldz Leading dimension of eigenvector matrix Z; #
# often equal to matrix size n #
# #
# INPUT + OUTPUT: #
# --------------- #
# AP (double[s]) On entry symmetric input matrix, stored #
# s = (n*(n+1))/2 in column major ordering. Depending on the #
# value of 'uplo' only the upper or lower #
# triangular part is stored. #
# (On output the array will be overwritten). #
# vl If range="V", lower bound of interval #
# (vl,vu], on output refined. #
# If range="A" or "I" not referenced as input. #
# On output the interval (vl,vu] contains ALL #
# the computed eigenvalues. #
# vu If range="V", upper bound of interval #
# (vl,vu], on output refined. #
# If range="A" or "I" not referenced as input. #
# On output the interval (vl,vu] contains ALL #
# the computed eigenvalues. #
# il If range="I", lower index (1-based indexing) of #
# the subset 'il' to 'iu'. #
# If range="A" or "V" not referenced as input. #
# On output the eigenvalues with index il to iu #
# are computed. #
# iu If range="I", upper index (1-based indexing) of #
# the subset 'il' to 'iu'. #
# If range="A" or "V" not referenced as input. #
# On output the eigenvalues with index il to iu #
# are computed. #
# #
# OUTPUT: #
# ------- #
# m Number of eigenvalues and eigenvectors computed #
# W (double[n]) Eigenvalues #
# The first 'm' entries contain the eigenvalues #
# Z Eigenvectors #
# (double[n*m]) Enough space must be provided to store the #
# vectors. 'm' should be bigger or equal #
# to 'n' for range="A" or "V" and 'iu-il+1' for #
# range="I" #
# #
# RETURN VALUE: #
# ------------- #
# 0 - Success #
# 1 - Wrong input parameter #
# 2 - Misc errors #
# #
######################################################################
######################################################################
# Fortran prototype: #
######################################################################
# #
# DSPEIG(JOBZ, RANGE, UPLO, N, AP, VL, VU, IL, IU, M, W, #
# Z, LDZ, INFO) #
# #
# CHARACTER JOBZ, RANGE, UPLO #
# INTEGER N, IL, IU, M, LDZ, INFO #
# DOUBLE PRECISION AP(*), VL, VU, W(*), Z(*,*) #
# #
######################################################################
The Hermitian case (only added to 'libmrrr.a' if specified in
'make.inc'):
######################################################################
# C function prototype: #
######################################################################
# #
# int zheeig(char *jobz, char *range, char *uplo, int *n, #
# double complex *A, int *lda, double *vl, double *vu, #
# int *il, int *iu, int *m, double *W, double complex *Z, #
# int *ldz); #
# #
# Arguments: #
# ---------- #
# #
# INPUTS: #
# ------- #
# jobz "N" - compute only eigenvalues #
# "V" - compute also eigenvectors #
# range "A" - all #
# "V" - by interval: (VL,VU] #
# "I" - by index: IL-IU #
# uplo "L" - Upper triangle of A stored #
# "U" - Lower triangle of A stored #
# n Order of the matrix A #
# lda Leading dimension of matrix A; #
# often equal to matrix size n #
# ldz Leading dimension of eigenvector matrix Z; #
# often equal to matrix size n #
# #
# INPUT + OUTPUT: #
# --------------- #
# A On entry Hermitian input matrix, stored #
# (double complex in column major ordering. Depending on the #
# [lda*n]) value of 'uplo' only the upper or lower #
# triangular part is referenced. #
# (On output the array will be overwritten). #
# vl If range="V", lower bound of interval #
# (vl,vu], on output refined. #
# If range="A" or "I" not referenced as input. #
# On output the interval (vl,vu] contains ALL #
# the computed eigenvalues. #
# vu If range="V", upper bound of interval #
# (vl,vu], on output refined. #
# If range="A" or "I" not referenced as input. #
# On output the interval (vl,vu] contains ALL #
# the computed eigenvalues. #
# il If range="I", lower index (1-based indexing) of #
# the subset 'il' to 'iu'. #
# If range="A" or "V" not referenced as input. #
# On output the eigenvalues with index il to iu #
# are computed. #
# iu If range="I", upper index (1-based indexing) of #
# the subset 'il' to 'iu'. #
# If range="A" or "V" not referenced as input. #
# On output the eigenvalues with index il to iu #
# are computed. #
# #
# OUTPUT: #
# ------- #
# m Number of eigenvalues and eigenvectors computed. #
# W (double[n]) Eigenvalues #
# The first 'm' entries contain the eigenvalues. #
# Z Eigenvectors. #
# (double complex Enough space must be provided to store the #
# [n*m]) vectors. 'm' should be bigger or equal #
# to 'n' for range="A" or "V" and 'iu-il+1' for #
# range="I". #
# #
# RETURN VALUE: #
# ------------- #
# 0 - Success #
# 1 - Wrong input parameter #
# 2 - Misc errors #
######################################################################
######################################################################
# Fortran prototype: #
######################################################################
# #
# DSYEIG(JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU, M, W, #
# Z, LDZ, INFO) #
# #
# CHARACTER JOBZ, RANGE, UPLO #
# INTEGER N, IL, IU, M, LDA, LDZ, INFO #
# DOUBLE PRECISION VL, VU, W(*) #
# COMPLEX*16 A(*,*), Z(*,*) #
# #
######################################################################
The Hermitian case (packed storage; only added to if specified in
'make.inc'):
######################################################################
# C function prototype: #
######################################################################
# #
# int zhpeig(char *jobz, char *range, char *uplo, int *n, #
# double complex *AP, double *vl, double *vu, #
# int *il, int *iu, int *m, double *W, #
# double complex *Z, int *ldz); #
# #
# Arguments: #
# ---------- #
# #
# INPUTS: #
# ------- #
# jobz "N" - compute only eigenvalues #
# "V" - compute also eigenvectors #
# range "A" - all #
# "V" - by interval: (VL,VU] #
# "I" - by index: IL-IU #
# uplo "L" - Upper triangle of A stored #
# "U" - Lower triangle of A stored #
# n Order of the matrix A #
# ldz Leading dimension of eigenvector matrix Z; #
# often equal to matrix size n #
# #
# INPUT + OUTPUT: #
# --------------- #
# AP (double On entry Hermitian input matrix, stored #
# complex[s]) in packed storage by columns. Depending on the #
# s = (n*(n+1))/2 value of 'uplo' only the upper or lower #
# triangular part is stored. #
# (On output the array will be overwritten). #
# vl If range="V", lower bound of interval #
# (vl,vu], on output refined. #
# If range="A" or "I" not referenced as input. #
# On output the interval (vl,vu] contains ALL #
# the computed eigenvalues. #
# vu If range="V", upper bound of interval #
# (vl,vu], on output refined. #
# If range="A" or "I" not referenced as input. #
# On output the interval (vl,vu] contains ALL #
# the computed eigenvalues. #
# il If range="I", lower index (1-based indexing) of #
# the subset 'il' to 'iu'. #
# If range="A" or "V" not referenced as input. #
# On output the eigenvalues with index il to iu #
# are computed. #
# iu If range="I", upper index (1-based indexing) of #
# the subset 'il' to 'iu'. #
# If range="A" or "V" not referenced as input. #
# On output the eigenvalues with index il to iu #
# are computed. #
# #
# OUTPUT: #
# ------- #
# m Number of eigenvalues and eigenvectors computed. #
# W (double[n]) Eigenvalues #
# The first 'm' entries contain the eigenvalues. #
# Z Eigenvectors #
# (double Enough space must be provided to store the #
# complex[n*m]) vectors. 'm' should be bigger or equal #
# to 'n' for range="A" or "V" and 'iu-il+1' for #
# range="I". #
# #
# RETURN VALUE: #
# ------------- #
# 0 - Success #
# 1 - Wrong input parameter #
# 2 - Misc errors #
# #
######################################################################
######################################################################
# Fortran prototype: #
######################################################################
# #
# ZHPEIG(JOBZ, RANGE, UPLO, N, AP, VL, VU, IL, IU, M, W, #
# Z, LDZ, INFO) #
# #
# CHARACTER JOBZ, RANGE, UPLO #
# INTEGER N, IL, IU, M, LDZ, INFO #
# DOUBLE PRECISION VL, VU, W(*) #
# COMPLEX*16 AP(*), Z(*,*) #
# #
######################################################################
The dense generalized definite eigenvalue problem:
--------------------------------------------------
Using LAPACK notation the routines are called DSYGEIG and ZHEGEIG.
Simply code snippets of how to call the functions are included in the
examples folder.
The problems that can be solved have the following form:
A*x = lambda*B*x or A*B*x = lambda*x or B*A*x = lambda*x,
where A and B are symmetric/Hermitian and B is definite.
Symmetric-definite case:
######################################################################
# C function prototype: #
######################################################################
# #
# int dsygeig(int *itype, char *jobz, char *range, char *uplo, #
# int *n, double *A, int *lda, double *B, int *ldb, #
# double *vl, double *vu, int *il, int *iu, int *m, #
# double *W); #
# #
# Arguments: #
# ---------- #
# #
# INPUTS: #
# ------- #
# itype 1 - A*x = lambda*B*x #
# 2 - A*B*x = lambda*x #
# 3 - B*A*x = lambda*x #
# jobz "N" - compute only eigenvalues #
# "V" - compute also eigenvectors #
# range "A" - all #
# "V" - by interval: (VL,VU] #
# "I" - by index: IL-IU #
# uplo "L" - Upper triangle of A and B stored #
# "U" - Lower triangle of A and B stored #
# n Order of the matrix A and B #
# lda Leading dimension of matrix A; #
# often equal to matrix size n #
# ldb Leading dimension of eigenvector matrix B; #
# often equal to matrix size n #
# #
# INPUT + OUTPUT: #
# --------------- #
# A On entry symmetric input matrix, stored #
# (double[lda*n]) in column major ordering. Depending on the #
# value of 'uplo' only the upper or lower #
# triangular part is referenced #
# On output the array will contain the #
# 'm' computed eigenvectors #
# B On entry symmetric definite input matrix, #
# (double[ldb*n]) stored in column major ordering. #
# Depending on the value of 'uplo' only the upper #
# or lower triangular part is referenced #
# On output overwritten #
# vl If range="V", lower bound of interval #
# (vl,vu], on output refined #
# If range="A" or "I" not referenced as input #
# On output the interval (vl,vu] contains ALL #
# the computed eigenvalues #
# vu If range="V", upper bound of interval #
# (vl,vu], on output refined #
# If range="A" or "I" not referenced as input. #
# On output the interval (vl,vu] contains ALL #
# the computed eigenvalues #
# il If range="I", lower index (1-based indexing) of #
# the subset 'il' to 'iu' #
# If range="A" or "V" not referenced as input. #
# On output the eigenvalues with index il to iu #
# are computed #
# iu If range="I", upper index (1-based indexing) of #
# the subset 'il' to 'iu' #
# If range="A" or "V" not referenced as input #
# On output the eigenvalues with index il to iu #
# are computed #
# #
# OUTPUT: #
# ------- #
# m Number of eigenvalues and eigenvectors computed #
# W (double[n]) Eigenvalues #
# The first 'm' entries contain the eigenvalues #
# #
# NOTICE: The routine will allocate work space of size #
# double[n*n] for range="A" or "V" and double[m*n] #
# for range="I" #
# #
# RETURN VALUE: #
# ------------- #
# 0 - Success #
# 1 - Wrong input parameter #
# 2 - Misc errors #
# #
######################################################################
######################################################################
# Fortran prototype: #
######################################################################
# #
# DSYGEIG(ITYPE, JOBZ, RANGE, UPLO, N, A, LDA, B, LDB, #
# VL, VU, IL, IU, M, W, INFO); #
# #
# CHARACTER JOBZ, RANGE, UPLO #
# INTEGER N, IL, IU, M, LDA, LDB, INFO #
# DOUBLE PRECISION A(*,*), B(*,*), VL, VU, W(*) #
# #
######################################################################
Symmetric-definite case (packed storage):
######################################################################
# C function prototype: #
######################################################################
# #
# int dspgeig(int *itype, char *jobz, char *range, char *uplo, #
# int *n, double *AP, double *BP, double *vl, #
# double *vu, int *il, int *iu, int *m, double *W, #
# double *Z, int *ldz); #
# #
# Arguments: #
# ---------- #
# #
# INPUTS: #
# ------- #
# itype 1 - A*x = lambda*B*x #
# 2 - A*B*x = lambda*x #
# 3 - B*A*x = lambda*x #
# jobz "N" - compute only eigenvalues #
# "V" - compute also eigenvectors #
# range "A" - all #
# "V" - by interval: (VL,VU] #
# "I" - by index: IL-IU #
# uplo "L" - Upper triangle of A and B stored #
# "U" - Lower triangle of A and B stored #
# n Order of the matrix A and B #
# ldz Leading dimension of matrix Z; #
# often equal to matrix size n #
# #
# INPUT + OUTPUT: #
# --------------- #
# AP On entry symmetric input matrix, stored in #
# (double[s]) packed format by columns. Depending on the #
# s = (n*(n+1))/2 value of 'uplo' only the upper or lower #
# triangular part is stored #
# On output the array will contain the #
# 'm' computed eigenvectors #
# BP On entry symmetric definite input matrix, stored #
# (double[s]) in packed format by columns. #
# s = (n*(n+1))/2 Depending on the value of 'uplo' only the upper #
# or lower triangular part is stored #
# On output overwritten #
# vl If range="V", lower bound of interval #
# (vl,vu], on output refined #
# If range="A" or "I" not referenced as input #
# On output the interval (vl,vu] contains ALL #
# the computed eigenvalues #
# vu If range="V", upper bound of interval #
# (vl,vu], on output refined #
# If range="A" or "I" not referenced as input. #
# On output the interval (vl,vu] contains ALL #
# the computed eigenvalues #
# il If range="I", lower index (1-based indexing) of #
# the subset 'il' to 'iu' #
# If range="A" or "V" not referenced as input. #
# On output the eigenvalues with index il to iu #
# are computed #
# iu If range="I", upper index (1-based indexing) of #
# the subset 'il' to 'iu' #
# If range="A" or "V" not referenced as input #
# On output the eigenvalues with index il to iu #
# are computed #
# #
# OUTPUT: #
# ------- #
# m Number of eigenvalues and eigenvectors computed #
# W (double[n]) Eigenvalues #
# The first 'm' entries contain the eigenvalues #
# Z (double[m*n]) Eigenvectors #
# Enough space must be provided to store the #
# vectors. 'm' should be bigger or equal #
# to 'n' for range="A" or "V" and 'iu-il+1' for #
# range="I". #
# #
# NOTICE: The routine will allocate work space of size #
# double[n*n] for range="A" or "V" and double[m*n] #
# for range="I" #
# #
# RETURN VALUE: #
# ------------- #
# 0 - Success #
# 1 - Wrong input parameter #
# 2 - Misc errors #
# #
######################################################################
######################################################################
# Fortran prototype: #
######################################################################
# #
# DSPGEIG(ITYPE, JOBZ, RANGE, UPLO, N, AP, BP, #
# VL, VU, IL, IU, M, W, Z, LDZ, INFO) #
# #
# CHARACTER JOBZ, RANGE, UPLO #
# INTEGER N, IL, IU, M, LDZ, INFO #
# DOUBLE PRECISION AP(*), BP(*), VL, VU, W(*), Z(*,*) #
# #
######################################################################
Hermitian-definite case:
######################################################################
# C function prototype: #
######################################################################
# #
# int zhegeig(int *itype, char *jobz, char *range, char *uplo, #
# int *n, double complex *A, int *lda, #
# double complex *B, int *ldb, double *vl, double *vu, #
# int *il, int *iu, int *m, double *W); #
# #
# Arguments: #
# ---------- #
# #
# INPUTS: #
# ------- #
# itype 1 - A*x = lambda*B*x #
# 2 - A*B*x = lambda*x #
# 3 - B*A*x = lambda*x #
# jobz "N" - compute only eigenvalues #
# "V" - compute also eigenvectors #
# range "A" - all #
# "V" - by interval: (VL,VU] #
# "I" - by index: IL-IU #
# uplo "L" - Upper triangle of A and B stored #
# "U" - Lower triangle of A and B stored #
# n Order of the matrix A and B #
# lda Leading dimension of matrix A; #
# often equal to matrix size n #
# ldb Leading dimension of eigenvector matrix B; #
# often equal to matrix size n #
# #
# INPUT + OUTPUT: #
# --------------- #
# A On entry symmetric input matrix, stored #
# (double in column major ordering. Depending on the #
# complex[lda*n]) value of 'uplo' only the upper or lower #
# triangular part is referenced #
# On output the array will contain the #
# 'm' computed eigenvectors #
# B On entry symmetric definite input matrix, #
# (double stored in column major ordering. #
# complex[ldb*n]) Depending on the value of 'uplo' only the upper #
# or lower triangular part is referenced #
# On output overwritten #
# vl If range="V", lower bound of interval #
# (vl,vu], on output refined #
# If range="A" or "I" not referenced as input #
# On output the interval (vl,vu] contains ALL #
# the computed eigenvalues #
# vu If range="V", upper bound of interval #
# (vl,vu], on output refined #
# If range="A" or "I" not referenced as input. #
# On output the interval (vl,vu] contains ALL #
# the computed eigenvalues #
# il If range="I", lower index (1-based indexing) of #
# the subset 'il' to 'iu' #
# If range="A" or "V" not referenced as input. #
# On output the eigenvalues with index il to iu #
# are computed #
# iu If range="I", upper index (1-based indexing) of #
# the subset 'il' to 'iu' #
# If range="A" or "V" not referenced as input #
# On output the eigenvalues with index il to iu #
# are computed #
# #
# OUTPUT: #
# ------- #
# m Number of eigenvalues and eigenvectors computed #
# W (double[n]) Eigenvalues #
# The first 'm' entries contain the eigenvalues #
# #
# NOTICE: The routine will allocate work space of size #
# double[n*n] for range="A" or "V" and double[m*n] #
# for range="I" #
# #
# RETURN VALUE: #
# ------------- #
# 0 - Success #
# 1 - Wrong input parameter #
# 2 - Misc errors #
# #
######################################################################
######################################################################
# Fortran prototype: #
######################################################################
# #
# ZHEGEIG(ITYPE, JOBZ, RANGE, UPLO, N, A, LDA, B, LDB, #
# VL, VU, IL, IU, M, W, INFO); #
# #
# CHARACTER JOBZ, RANGE, UPLO #
# INTEGER N, IL, IU, M, LDA, LDB, INFO #
# DOUBLE PRECISION VL, VU, W(*) #
# COMPLEX*16 A(*,*), B(*,*) #
# #
######################################################################
Hermitian-definite case (packed storage):
######################################################################
# C function prototype: #
######################################################################
# #
# int zhpgeig(int *itype, char *jobz, char *range, char *uplo, #
# int *n, double complex *AP, double complex *BP, #
# double *vl, double *vu, int *il, int *iu, int *m, #
# double *W, double complex *Z, int *ldz); #
# #
# Arguments: #
# ---------- #
# #
# INPUTS: #
# ------- #
# itype 1 - A*x = lambda*B*x #
# 2 - A*B*x = lambda*x #
# 3 - B*A*x = lambda*x #
# jobz "N" - compute only eigenvalues #
# "V" - compute also eigenvectors #
# range "A" - all #
# "V" - by interval: (VL,VU] #
# "I" - by index: IL-IU #
# uplo "L" - Upper triangle of A and B stored #
# "U" - Lower triangle of A and B stored #
# n Order of the matrix A and B #
# ldz Leading dimension of matrix Z; #
# often equal to matrix size n #
# #
# INPUT + OUTPUT: #
# --------------- #
# AP On entry Hermitian input matrix, stored in #
# (double packed format by columns. Depending on the #
# complex[s]) value of 'uplo' only the upper or lower #
# s = (n*(n+1))/2 triangular part is stored #
# On output the array will contain the #
# 'm' computed eigenvectors #
# BP On entry Hermitian definite input matrix, stored #
# (double in packed format by columns. #
# complex[s]) Depending on the value of 'uplo' only the upper #
# s = (n*(n+1))/2 or lower triangular part is stored #
# On output overwritten #
# vl If range="V", lower bound of interval #
# (vl,vu], on output refined #
# If range="A" or "I" not referenced as input #
# On output the interval (vl,vu] contains ALL #
# the computed eigenvalues #
# vu If range="V", upper bound of interval #
# (vl,vu], on output refined #
# If range="A" or "I" not referenced as input. #
# On output the interval (vl,vu] contains ALL #
# the computed eigenvalues #
# il If range="I", lower index (1-based indexing) of #
# the subset 'il' to 'iu' #
# If range="A" or "V" not referenced as input. #
# On output the eigenvalues with index il to iu #
# are computed #
# iu If range="I", upper index (1-based indexing) of #
# the subset 'il' to 'iu' #
# If range="A" or "V" not referenced as input #
# On output the eigenvalues with index il to iu #
# are computed #
# #
# OUTPUT: #
# ------- #
# m Number of eigenvalues and eigenvectors computed #
# W (double[n]) Eigenvalues #
# The first 'm' entries contain the eigenvalues #
# Z (double Eigenvectors #
# complex[m*n]) Enough space must be provided to store the #
# vectors. 'm' should be bigger or equal #
# to 'n' for range="A" or "V" and 'iu-il+1' for #
# range="I". #
# #
# NOTICE: The routine will allocate work space of size #
# double[n*n] for range="A" or "V" and double[m*n] #
# for range="I" #
# #
# RETURN VALUE: #
# ------------- #
# 0 - Success #
# 1 - Wrong input parameter #
# 2 - Misc errors #
# #
######################################################################
######################################################################
# Fortran prototype: #
######################################################################
# #
# ZHPGEIG(ITYPE, JOBZ, RANGE, UPLO, N, AP, BP, #
# VL, VU, IL, IU, M, W, Z, LDZ, INFO) #
# #
# CHARACTER JOBZ, RANGE, UPLO #
# INTEGER N, IL, IU, M, LDZ, INFO #
# DOUBLE PRECISION VL, VU, W(*) #
# COMPLEX*16 AP(*), BP(*), Z(*,*) #
# #
######################################################################
COMMENTS AND BUGS:
petschow@aices.rwth-aachen.de
About
The MRRR algorithm for multi-core and SMP systems
Resources
License
Stars
Watchers
Forks
Releases
No releases published
Packages 0
No packages published
Languages
- Fortran 58.6%
- C 38.5%
- Makefile 2.2%
- Other 0.7%