在Microsoft Visual Studio 2010中配置使用LAPACK/ARPACK - Hellraider - 博客园
和是一个求解大规模稠密/稀疏矩阵问题的库，最近在做特征值问题时用到。ARPACK这库相当古老，最早是RICE的一帮人弄的。LAPACK也差不多，貌似是美帝某个.gov发起的。这俩源代码是Fortran，最近在哥大求到的一个项目需要移植到Windows平台下。方法是C语言的就用VS，Fortran的用的gFortran编译之，然后打成dll和lib。如果想要更高级点，还有配置。下面谈谈具体怎么搞：
0. 参考资料
LAPACK比较简单，讲的很是详细，做了下补充，关于冲突的问题可以看。但我觉得改成/MT就可以了，具体见下面的步骤。
关于ARPACK如何编译成C++库，最最原始的资料来自USC的，但因为版本问题有点过时了。
做了一些修改，但过程比较简单。
还有的步骤，内容也大同小异，还有这个也是。
此外有，目测是从中文翻译来的。
上给了一篇编译成现代版Fortran库的步骤，含ARKACK和LAPACK，但人直接无视Makefile，仅供参考。
关于ARPACK++，可以参考和。
我在上述基础上改了改，在VS2010中搞定。
：下载左边链接的Visual Studio版本即可。
：下载arpack96和对应的patch。ARPACK自带BLAS和LAPACK，没有依赖问题。
：&下载一个installer。比如我下的是mingw-get-inst-.exe。
：下载ARPACK++ beta version。ARPACK++依赖于ARPACK，LAPACK，SUPERLU。
：这里需要用到的是99年的2.0版。再早没有，再晚的也不行。
2. 编译、使用LAPACK
解压LAPACK，它是自带某版本的Visual Studio工程的。
先转换到VS2010，然后把\LIB清空，我们自己编译新的。按F2cBlas, tmglib, blas, CLAPACK的顺序依次编译，得到的lib就是我们需要的。
在使用LAPACK时，在：
项目属性-&配置属性-&C/C++-&常规-&附加包含目录中添加LAPACK\INCLUDE目录位置；
项目属性-&配置属性-&链接器-&常规-&附加库目录中添加LAPACK\LIB文件夹目录位置；
项目属性-&配置属性-&链接器-&输入-&附加依赖项中添加libf2cd.lib BLASd.lib clapackd.lib tmglibd.lib（release模式就不带d）；
项目属性-&配置属性-&C/C++-&代码生成-&运行库选成/MTd（或/MT）；
最后用一个程序测试一下：
lapacktest.cpp
#include &stdio.h&
#include &process.h&
//因为程序是C++，而CLAPACK是f2c程序转换的C语言版本，所以在此处用extern关键字调用
#include &f2c.h&
#include &clapack.h&
#define SIZE 4
int main()
char JOBU;
char JOBVT;
//数据类型integer是fortran里的。这里在C++下可以使用的原因是f2c.h文件中已经作了定义
integer M = SIZE;
integer N = SIZE;
integer LDA = M;
integer LDU = M;
integer LDVT = N;
integer LWORK;
integer INFO;
integer mn = min( M, N );
integer MN = max( M, N );
double a[SIZE*SIZE] = { 16.0, 5.0, 9.0 , 4.0, 2.0, 11.0, 7.0 , 14.0, 3.0, 10.0, 6.0, 15.0, 13.0, 8.0, 12.0, 1.0};
double s[SIZE];
double wk[201];
double uu[SIZE*SIZE];
double vt[SIZE*SIZE];
JOBU = 'A';
JOBVT = 'A';
LWORK = 201;
/* Subroutine int dgesvd_(char *jobu, char *jobvt, integer *m, integer *n,
doublereal *a, integer *lda, doublereal *s, doublereal *u, integer *
ldu, doublereal *vt, integer *ldvt, doublereal *work, integer *lwork,
integer *info)
dgesvd_( &JOBU, &JOBVT, &M, &N, a, &LDA, s, uu, &LDU, vt, &LDVT, wk, &LWORK, &INFO);
printf("INFO=%d \n", INFO );
for ( i= 0; i& SIZE; i++ ) {
printf("s[ %d ] = %f\n", i, s[ i ] );
system("pause");
成功执行，结果非零。完工~
只是用LAPACK的童鞋，到这里就完成了~
3. 安装MinGW
执行MinGW安装程序选择目录，比如我的安装目录是C:\MinGW\。记得至少得勾上C++，gfortran编译器和MSYS。
装完以后别忘了附加环境变量PATH：C:\MinGW\bin。
4. 解压并改动ARPACK源码
先把arpack96.tar.gz解压，比如我的ARPACK目录是E:\ARPACK\。
将patch.tar.gz解压至ARPACK目录并覆盖。
然后纠结的来了，这里有三个文件要修改，不改动会导致现代gfortran编译不过。
第一个是ARPACK\UTIL\second.f，里面的ETIME遍不过。整个子过程可以替换成：
1 SUBROUTINE SECOND( T )
CALL CPU_TIME(T)
第二个是ARPACK\BLAS\dnrm2.f，里面的ASSIGN遍不过。不过我们可以在netlib上找到BLAS里。下载并直接覆盖之。
第三个是ARPACK\BLAS\snrm2.f，问题同上，找到覆盖。
5. 改动ARPACK\ARmake.inc
里面需要改动的有3处：
home = e:/ARPACK （改成自己的ARPACK目录）
FC = gfortran
FFLAGS&= -O2
6. 编译&打包&VS2010配置
打开MinGW Shell，先
$ cd ARPACK/
(自己的ARPACK目录)
$ make lib
$ dllwrap --export-all-symbols BLAS/*.o LAPACK/*.o SRC/*.o UTIL/*.o -lgfortran --output-def arpack_win32.def -o arpack_win32.dll
出现warning无视。然后打开VS2010程序目录下Visual Studio Tools里面的Visual Studio命令提示(2010)，
& cd ARPACK\
(自己的ARPACK目录)
& lib /machine:i386 /def:arpack_win32.def
我们需要的就是ARPACK目录下生成的.dll和.lib。
VS2010配置就简单了：
项目属性-&配置属性-&链接器-&常规-&附加库目录中添加lib文件夹目录；
项目属性-&配置属性-&链接器-&输入-&附加依赖项中添加arpack_win32.lib；
项目属性-&配置属性-&C/C++-&代码生成-&运行库选成/MTd（或/MT）；
在自己的$(SolutionDir)\Debug或$(SolutionDir)\Release中添加arpack_win32.dll，以及MinGW\Bin下面的libgfortran-3.dll和libquadmath-0.dll。
用MinGW build之后再用VS2010的lib.exe打包，好像就会产生这种依赖，这在中也有提及。在用f2c转换的LAPACK的VS2005版本中就没有这个问题，令人非常不爽（）。不管怎么样好用就行了，效率应该也没影响。
但此时人生还不完整，因为没有头文件。
7. 中场休息
啊不，应该叫伤害测试才对。。。虽然还没有官方头文件，我们可以自己捏造个，然后ARPACK应该可以正常使用了。比如这位提供了两个自己写的ARPACK头文件和测试程序，我们可以以这里的特征值问题（,&）为例子作测试。当然啦众所周知哈佛==不靠谱，所以不用下了，这两个文件显然是编不过的。其中的改动有：（1）&fstream.h&要改成&iostream&（2）加（3）pow的参数类型产生二义性改成powf（4）头文件的dseupd_接口参数名重复 - 估计这是不靠谱造成的（尼玛，拿调用当定义真是太可怕了），而才是正确的。废话少说，上修改后的源码：
In this header file, I have defined a simplified function
call to the ARPACK solver routine for a simple, symmetric
eigenvalue problem Av = lv.
The looping procedure and
final extraction of eigenvalues and vectors is handled
automatically.
Most of the parameters to the FORTRAN
functions are hidden from the user, since most of them are
determined from user input anyway.
The remaining parameters to the function calls are as follows:
dsaupd(int n, int nev, double *Evals)
dsaupd(int n, int nev, doubel *Evals, double **Evecs)
n: the order of the square matrix A
nev: the number of eigenvalues to be found, starting at the
Note that the highest eigenvalues, or some
other choices, can be found.
For now, the choice of
the lowest nev eigenvalues is hard-coded.
Evals: a one-dimensional array of length nev to hold the
eigenvalues.
Evecs: a two-dimensional array of size nev by n to hold the
eigenvectors.
If this argument is not provided, the
eigenvectors are not calculated.
Note that the
vectors are stored as columns of Evecs, so that the
elements of vector i are the values Evecs[j][i].
The function is overdefined, so that you can call it in either
fashion and the appropriate version will be run.
To use these function calls, there must be a function
defined in the calling program of the form
av(int n, double *in, double *out)
where the function av finds out = A.in, and n is the order of
the matrix A.
This function must be defined before the
statement that includes this header file, as it needs to know
about this function.
It is used in the looping procedure.
30 August 1999
45 #include &math.h&
47 extern "C" void dsaupd_(int *ido, char *bmat, int *n, char *which,
int *nev, double *tol, double *resid, int *ncv,
double *v, int *ldv, int *iparam, int *ipntr,
double *workd, double *workl, int *lworkl,
int *info);
53 extern "C" void dseupd_(int *rvec, char *All, int *select, double *d,
double *z, int *ldz, double *sigma,
char *bmat, int *n, char *which, int *nev,
double *tol, double *resid, int *ncv, double *v,
int *ldv, int *iparam, int *ipntr, double *workd,
double *workl, int *lworkl, int *ierr);
60 void dsaupd(int n, int nev, double *Evals)
int ido = 0; /* Initialization of the reverse communication
parameter. */
char bmat[2] = "I"; /* Specifies that the right hand side matrix
should be this makes
the problem a standard eigenvalue problem.
Setting bmat = "G" would have us solve the
problem Av = lBv (this would involve using
some other programs from BLAS, however). */
char which[3] = "SM"; /* Ask for the nev eigenvalues of smallest
magnitude.
The possible options are
LM: largest magnitude
SM: smallest magnitude
LA: largest real component
SA: smallest real compoent
LI: largest imaginary component
SI: smallest imaginary component */
double tol = 0.0; /* S tol&=0 specifies
machine precision */
resid = new double[n];
int ncv = 4* /* The largest number of basis vectors that will
be used in the Implicitly Restarted Arnoldi
Work per major iteration is
proportional to N*NCV*NCV. */
if (ncv&n) ncv =
double *v;
v = new double[ldv*ncv];
iparam = new int[11]; /* An array used to pass information to the routines
about their functional modes. */
iparam[0] = 1;
// Specifies the shift strategy (1-&exact)
iparam[2] = 3*n; // Maximum number of iterations
iparam[6] = 1;
/* Sets the mode of dsaupd.
1 is exact shifting,
2 is user-supplied shifts,
3 is shift-invert mode,
4 is buckling mode,
5 is Cayley mode. */
ipntr = new int[11]; /* Indicates the locations in the work array workd
where the input and output vectors in the
callback routine are located. */
workd = new double[3*n];
workl = new double[ncv*(ncv+8)];
int lworkl = ncv*(ncv+8); /* Length of the workl array */
int info = 0; /* Passes convergence information out of the iteration
routine. */
int rvec = 0; /* Specifies that eigenvectors should not be calculated */
int *select;
select = new int[ncv];
double *d;
d = new double[2*ncv]; /* This vector will return the eigenvalues from
the second routine, dseupd. */
/* Here we enter the main loop where the calculations are
performed.
The communication parameter ido tells us when
the desired tolerance is reached, and at that point we exit
and extract the solutions. */
dsaupd_(&ido, bmat, &n, which, &nev, &tol, resid,
&ncv, v, &ldv, iparam, ipntr, workd, workl,
&lworkl, &info);
if ((ido==1)||(ido==-1)) av(n, workd+ipntr[0]-1, workd+ipntr[1]-1);
} while ((ido==1)||(ido==-1));
/* From those results, the eigenvalues and vectors are
extracted. */
if (info&0) {
cout && "Error with dsaupd, info = " && info && "\n";
cout && "Check documentation in dsaupd\n\n";
dseupd_(&rvec, "All", select, d, v, &ldv, &sigma, bmat,
&n, which, &nev, &tol, resid, &ncv, v, &ldv,
iparam, ipntr, workd, workl, &lworkl, &ierr);
if (ierr!=0) {
cout && "Error with dseupd, info = " && ierr && "\n";
cout && "Check the documentation of dseupd.\n\n";
} else if (info==1) {
cout && "Maximum number of iterations reached.\n\n";
} else if (info==3) {
cout && "No shifts could be applied during implicit\n";
cout && "Arnoldi update, try increasing NCV.\n\n";
/* Before exiting, we copy the solution information over to
the arrays of the calling program, then clean up the
memory used by this routine.
For some reason, when I
don't find the eigenvectors I need to reverse the order of
the values. */
for (i=0; i& i++) Evals[i] = d[nev-1-i];
delete select;
190 void dsaupd(int n, int nev, double *Evals, double **Evecs)
int ido = 0;
char bmat[2] = "I";
char which[3] = "SM";
double tol = 0.0;
resid = new double[n];
int ncv = 4*
if (ncv&n) ncv =
double *v;
v = new double[ldv*ncv];
iparam = new int[11];
iparam[0] = 1;
iparam[2] = 3*n;
iparam[6] = 1;
ipntr = new int[11];
workd = new double[3*n];
workl = new double[ncv*(ncv+8)];
int lworkl = ncv*(ncv+8);
int info = 0;
int rvec = 1;
// Changed from above
int *select;
select = new int[ncv];
double *d;
d = new double[2*ncv];
dsaupd_(&ido, bmat, &n, which, &nev, &tol, resid,
&ncv, v, &ldv, iparam, ipntr, workd, workl,
&lworkl, &info);
if ((ido==1)||(ido==-1)) av(n, workd+ipntr[0]-1, workd+ipntr[1]-1);
} while ((ido==1)||(ido==-1));
if (info&0) {
cout && "Error with dsaupd, info = " && info && "\n";
cout && "Check documentation in dsaupd\n\n";
dseupd_(&rvec, "All", select, d, v, &ldv, &sigma, bmat,
&n, which, &nev, &tol, resid, &ncv, v, &ldv,
iparam, ipntr, workd, workl, &lworkl, &ierr);
if (ierr!=0) {
cout && "Error with dseupd, info = " && ierr && "\n";
cout && "Check the documentation of dseupd.\n\n";
} else if (info==1) {
cout && "Maximum number of iterations reached.\n\n";
} else if (info==3) {
cout && "No shifts could be applied during implicit\n";
cout && "Arnoldi update, try increasing NCV.\n\n";
for (i=0; i& i++) Evals[i] = d[i];
for (i=0; i& i++) for (j=0; j&n; j++) Evecs[j][i] = v[i*n+j];
delete select;
dsaupd.cpp
This is an example program using the header file
dsaupd.h, my implementation of the ARPACK sparse matrix
solver routine.
See that header file for more
documentation on its use.
For this example, we will work with the four by four
matrix H[i,j] = i^j + j^i.
The output of this program
should be:
Eigenvalue
0: 0.151995
Eigenvector 0: (0.610828, -0.7065, -0.0527271)
Eigenvalue
1: 1.80498
Eigenvector 1: (0.23, -0..103072)
Eigenvalue
2: 17.6352
Eigenvector 2: (-0.216885, -0.547084, -0..26195)
7 September 1999
// Begin with some standard include files.
#include &math.h&
#include &stdio.h&
#include &iostream&
using namespace
To use the routines from dsaupd.h, we need to define our
own matrix multiplication routine.
Here I show one method
that would exploit sparsity by storing only the non-zero
matrix elements (though for this example we won't actually
have a sparse matrix to work with).
There is a routine for the multiplication, and a global
variable T that will hold the non-zero matrix elements
and their indicies.
The global integer tlen will hold the
number of non-zero elements.
The matrix multiplication function
needs to have the form below, and must be defined before
the file dsaupd.h is included.
void av(int n, double *in, double *out);
#include "dsaupd.h"
double **T;
void main()
int n, nev, i,
double **Evecs, *E
n = 4; // The order of the matrix
Here we generate the matrix T for the multiplication.
It helps if we know the number of non-zero matrix
elements before we find them, so that we can save on
storage space.
If not, generate enough storage space
for T to cover an upper limit.
tlen = n*n;
T = new double*[tlen];
for (i=0; i& i++) T[i] = new double[3];
for (i=0; i&n; i++) for (j=0; j&n; j++) {
T[tlen][0] =
T[tlen][1] =
T[tlen][2] = powf(i+1, j+1) + powf(j+1, i+1);
We will calculate both the eigenvalues and the
eigenvectors in this example.
To not calculate the
vectors, just leave that argument out of the call to
We will find only three of the eigenvalues, to make
it more clear how the eigenvectors come out of the
Besides that, the ARPACK routine won't
calculate all of the eigenvalues.
I'm not sure if
that is intentional, or if I am missing something.
nev = 3; // The number of values to calculate
Evals = new double[nev];
Evecs = new double*[n];
for (i=0; i&n; i++) Evecs[i] = new double[nev];
dsaupd(n, nev, Evals, Evecs);
// Now we output the results.
for (i=0; i& i++) {
cout && "Eigenvalue
" && i && ": " && Evals[i] && "\n";
cout && "Eigenvector " && i && ": ("
&& Evecs[0][i] && ", "
&& Evecs[1][i] && ", "
&& Evecs[2][i] && ", "
&& Evecs[3][i] && ")\n";
void av(int n, double *in, double *out)
for (i=0; i&n; i++) out[i] = 0;
for (i=0; i& i++) out[(int)T[i][0]] += in[(int)T[i][1]] * T[i][2];
如果VS2010在步骤6中配置好，那可以正常编译运行，得到的结果和注释中一样。如果不能，那之前&一定是什么地方搞错了&，这时候应该睡一会:D。
对于仅仅用来求对称稀疏矩阵特征值的简单需求而言，安装到此为止即可。
8. SUPERLU
在ARPACK++的README里提到了对SUPERLU的引用，不过里面给的SUPERLU下载地址已经404了。由于ARPACK++2003年已经不更新，使用的SUPERLU可以推断出是之前的版本，也就是2.0版。事实上高版本的SuperLU接口写的不一样，函数参数有变，会出一堆空引用。
首先解压，然后新建一个临时VS 2010工程，把SUPERLU\SRC中的头文件和代码导入。
SUPERLU附带了CBLAS，这个BLAS和CLAPACK的BLAS又不一样，所以我们要奇葩地把这个\CBLAS里的文件同样混到VS2010工程里，重复的选择覆盖。这样SUPERLU编译出来的包就不引用BLAS了。接下来，在临时工程的编译选项中：
项目属性-&配置属性-&常规-&配置类型选成.lib文件；
项目属性-&配置属性-&C/C++-&代码生成-&运行库选成/MTd（或/MT）；
项目属性-&配置属性-&C/C++-&代码生成-&基本运行时检查，对于Debug，要改成/RTCs。因为SuperLU代码中有未初始化变量作为参数传递，会被VS的Run-Time Checking认为是非法的，所以把那个关掉。对于Release不用管，本来就没有RTC。我不是的~
找到superlu_time.c，内容改成如下所示（windows缺乏sys/time.h, sys/times.h，所以忽略之）：
superlu_time.c
Returns the time in seconds used by the process.
* Note: the timer function call is machine dependent. Use conditional
compilation to choose the appropriate function.
#ifdef SUN
It uses the system call gethrtime(3C), which is accurate to
nanoseconds.
#include &sys/time.h&
double SuperLU_timer_() {
return ( (double)gethrtime() / 1e9 );
#elif _WIN32
double SuperLU_timer_() {
return 0.0;
#include &sys/types.h&
#include &sys/times.h&
#include &time.h&
#include &sys/time.h&
#ifndef CLK_TCK
#define CLK_TCK 60
double SuperLU_timer_()
times(&use);
tmp = use.tms_
tmp += use.tms_
return (double)(tmp) / CLK_TCK;
返回0.0是很猥琐的，但我也不是的啦。。。
再然后，由于LAPACK和SuperLU共用了6个同名函数（dlamch, dlamc1, dlamc2, dlamc3, dlamc4, dlamc5, 好像是检测硬件的），共存会打架（没有namespace的后果~）。所以我们还要把SuperLU中这个6个函数名改了。具体做法就是在当前项目中&全部替换&，比如我把这6个函数统统加上_lu的后缀。作为sanity check，替换的次数是：
dlamch：43个，
dlamc1：8个，
dlamc2：8个，
dlamc3：42个，
dlamc4：9个，
dlamc5：7个
最后编译，输出.lib文件，比如我的是SuperLUd_win32.lib。那么在相应的VS2010工程配置中：
项目属性-&配置属性-&链接器-&常规-&附加库目录中添加SUPERLU\你的lib文件夹目录；
项目属性-&配置属性-&链接器-&输入-&附加依赖项中添加SuperLUd_win32.lib；
9. ARPACK++
ARPACK++是RICE那几个人开发的ARPACK的C++包装，依赖于ARPACK，LAPACK，SUPERLU等，文档和例程都比较详细。实际上我们只需要它的头文件（或者看看例程什么的），makefile和readme都无视好了。当然问题仍然是比较古老没人维护，永远的beta版，所以现代c++编译器无论哪个都是编不过的。这时候一位MIT的出手了，搞了个.diff，可以在g++下编译通过。当然这和VS还不一样。
首先解压ARPACK++，我的路径是E:\ARPACK++\。
下载，放在和ARPACK++同级的目录，比如我的E:根目录。
打开MinGW Shell：
$ mingw-get install msys-patch
（联网装一个patch utility，可以apply这个.diff文件，以前装过就免了）
$ patch -p 1 -d arpack++ & arpackpp1.2.patch.diff
ARPACK++\INCLUDE里就是我们要的头文件。记得在VS2010工程中，在
项目属性-&配置属性-&C/C++-&常规-&附加包含目录中添加ARPACK++\INCLUDE目录位置
10. 建立VS2010工程&测试
现在ARPACK++能不能用呢？还有选项要改。
首先看VS2010工程属性：
项目属性-&配置属性-&C/C++-&常规-&附加包含目录下面，应该有2个文件夹路径，分别是LAPACK和ARPACK++
项目属性-&配置属性-&C/C++-&预处理器-&预处理器定义，加入__GNUG__，因为这个版本的SuperLU根本没考虑到什么Win32，要VS假装自己是Gnu..
项目属性-&配置属性-&C/C++-&代码生成-&运行库应该是/MT、/MTd，说了很多次了
项目属性-&配置属性-&链接器-&常规-&附加库目录下面，应该有3个文件夹，分别是LAPACK，ARPACK和SUPERLU
项目属性-&配置属性-&链接器-&输入-&附加依赖项下面，应该有6个lib：LAPACK4个，ARPACK和SUPERLU各一个
最后给一套测试文件（直接从arpack++\exmaples里挖来的），可以正常运行：
ARPACK++ v1.2 2/18/2000
c++ interface to ARPACK code.
MODULE AREig.h.
Example functions AREig. This function returns eigenvalues, EigVal, and
eigenvectors, EigVec, of real symmetric, real nonsymmetric and complex
problems in various modes using ARluSymStdEig, ARluNonSymStdEig and
ARluNonSymGenEig classes.
The SuperLU package is employed to solve the linear systems that appear
when the shift-and-invert spectral transformation is being used.
There are eighteen different versions of AREig, as shown below. The type
and the meaning of each AREig parameter is briefly described in section
II. For a complete description of all parameters, see the Appendix of
"ARPACK++: a c++ implementation of ARPACK eigenvalue package."
I) How to call AREig:
AREig is a c++ overloaded function, which means that there are
several definitions for this single function. Each definition is
related to a different problem and a different kind of data. All
twenty six available AREig functions are listed below. They are
divided according to the problem type (standard or generalized),
the type of the involved matrices (real or complex), the
computational mode being used (regular or shift-and-invert) and
the desired output data (only eigenvalues or eigenvalues and
eigenvectors).
1) Real symmetric standard problems:
a) Regular mode:
To obtain only eigenvalues:
int AREig(EigVal, n, nnz, A, irow, pcol, uplo, nev,
which, ncv, tol, maxit, resid, AutoShift)
ii) To obtain eigenvalues and eigenvectors:
int AREig(EigVal, EigVec, n, nnz, A, irow, pcol, uplo,
nev, which, ncv, tol, maxit, resid, AutoShift)
b) Real shift-and-invert mode:
To obtain only eigenvalues:
int AREig(EigVal, n, nnz, A, irow, pcol, uplo, sigma,
nev, which, ncv, tol, maxit, resid, AutoShift)
ii) To obtain eigenvalues and eigenvectors:
int AREig(EigVal, EigVec, n, nnz, A, irow, pcol, uplo, sigma,
nev, which, ncv, tol, maxit, resid, AutoShift)
2) Real symmetric generalized problems:
a) Regular mode:
To obtain only eigenvalues:
int AREig(EigVal, n, nnzA, A, irowA, pcolA, nnzB, B,
irowB, pcolB, uplo, nev, which, ncv, tol,
maxit, resid, AutoShift)
ii) To obtain eigenvalues and eigenvectors:
int AREig(EigVal, EigVec, n, nnzA, A, irowA, pcolA,
nnzB, B, irowB, pcolB, uplo, nev, which,
ncv, tol, maxit, resid, AutoShift)
b) Shift-and-invert, buckling and Cayley modes:
To obtain only eigenvalues:
int AREig(EigVal, n, nnzA, A, irowA, pcolA, nnzB, B,
irowB, pcolB, uplo, InvertMode, sigma, nev,
which, ncv, tol, maxit, resid, AutoShift)
ii) To obtain eigenvalues and eigenvectors:
int AREig(EigVal, EigVec, n, nnzA, A, irowA, pcolA, nnzB,
B, irowB, pcolB, uplo, InvertMode, sigma, nev,
which, ncv, tol, maxit, resid, AutoShift)
3) Real nonsymmetric standard problems:
a) Regular mode:
To obtain only eigenvalues:
int AREig(EigValR, EigValI, n, nnz, A, irow, pcol, nev,
which, ncv, tol, maxit, resid, AutoShift)
ii) To obtain eigenvalues and eigenvectors:
int AREig(EigValR, EigValI, EigVec, n, nnz, A, irow, pcol,
nev, which, ncv, tol, maxit, resid, AutoShift)
b) Real shift-and-invert mode:
To obtain only eigenvalues:
int AREig(EigValR, EigValI, n, nnz, A, irow, pcol, sigma,
nev, which, ncv, tol, maxit, resid, AutoShift)
ii) To obtain eigenvalues and eigenvectors:
int AREig(EigValR, EigValI, EigVec, n, nnz, A, irow, pcol,
sigma, nev, which, ncv, tol, maxit, resid, AutoShift)
4) Real nonsymmetric generalized problems:
a) Regular mode:
To obtain only eigenvalues:
int AREig(EigValR, EigValI, n, nnzA, A, irowA, pcolA,
nnzB, B, irowB, pcolB, nev, which, ncv, tol,
maxit, resid, AutoShift)
ii) To obtain eigenvalues and eigenvectors:
int AREig(EigValR, EigValI, EigVec, n, nnzA, A, irowA,
pcolA, nnzB, B, irowB, pcolB, nev, which, ncv,
tol, maxit, resid, AutoShift)
b) Real shift-and-invert mode:
To obtain only eigenvalues:
int AREig(EigValR, EigValI, n, nnzA, A, irowA, pcolA,
nnzB, B, irowB, pcolB, sigma, nev, which, ncv,
tol, maxit, resid, AutoShift)
ii) To obtain eigenvalues and eigenvectors:
int AREig(EigValR, EigValI, EigVec, n, nnzA, A, irowA,
pcolA, nnzB, B, irowB, pcolB, sigma, nev, which,
ncv, tol, maxit, resid, AutoShift)
c) Complex shift-and-invert mode:
To obtain only eigenvalues:
int AREig(EigValR, EigValI, n, nnzA, A, irowA, pcolA,
nnzB, B, irowB, pcolB, part, sigmaR, SigmaI,
nev, which, ncv, tol, maxit, resid, AutoShift)
ii) To obtain eigenvalues and eigenvectors:
int AREig(EigValR, EigValI, EigVec, n, nnzA, A, irowA, pcolA,
nnzB, B, irowB, pcolB, part, sigmaR, SigmaI,
nev, which, ncv, tol, maxit, resid, AutoShift)
5) Complex standard problems:
a) Regular mode:
To obtain only eigenvalues:
int AREig(EigVal, n, nnz, A, irow, pcol, nev,
which, ncv, tol, maxit, resid, AutoShift)
ii) To obtain eigenvalues and eigenvectors:
int AREig(EigVal, EigVec, n, nnz, A, irow, pcol, nev,
which, ncv, tol, maxit, resid, AutoShift)
b) Shift-and-invert mode:
To obtain only eigenvalues:
int AREig(EigVal, n, nnz, A, irow, pcol, sigma,
nev, which, ncv, tol, maxit, resid, AutoShift)
ii) To obtain eigenvalues and eigenvectors:
int AREig(EigVal, EigVec, n, nnz, A, irow, pcol, sigma,
nev, which, ncv, tol, maxit, resid, AutoShift)
6) Complex generalized problems:
a) Regular mode:
To obtain only eigenvalues:
int AREig(EigVal, n, nnzA, A, irowA, pcolA, nnzB, B,
irowB, pcolB, nev, which, ncv, tol, maxit,
resid, AutoShift)
ii) To obtain eigenvalues and eigenvectors:
int AREig(EigVal, EigVec, n, nnzA, A, irowA, pcolA,
nnzB, B, irowB, pcolB, nev, which, ncv, tol,
maxit, resid, AutoShift)
b) Shift-and-invert mode:
To obtain only eigenvalues:
int AREig(EigVal, n, nnzA, A, irowA, pcolA, nnzB, B,
irowB, pcolB, sigma, nev, which, ncv, tol,
maxit, resid, AutoShift)
ii) To obtain eigenvalues and eigenvectors:
int AREig(EigVal, EigVec, n, nnzA, A, irowA, pcolA,
nnzB, B, irowB, pcolB, sigma, nev, which,
ncv, tol, maxit, resid, AutoShift)
II) AREig parameters:
Each AREig function contains a variable number of parameters, because
some parameters are not mandatory. Optional parameters should only be
defined by the used if the default values supplied by the function are
not suitable.
In the description of AREig parameters, the following notation is used
to describe eigenvalue problems:
a) standard problems:
A*EigVec = Eigvec*EigVal,
b) generalized problems: A*EigVec = B*EigVec*EigVal.
1) Compulsory input parameters:
Dimension of the problem.
Number of nonzero elements in matrix A.
Same as nnz.
Number of nonzero elements in matrix B.
Array of nonzero elements in matrix A.
TYPE must be one of "float", "double",
"arcomplex&float&" or "arcomplex&double&".
Array of nonzero elements in matrix B.
TYPE must be one of "float", "double",
"arcomplex&float&" or "arcomplex&double&".
int irow[]
Array of row indices of the nonzero elements in A.
int irowA[]
Same as irow.
int irowB[]
Array of row indices of the nonzero elements in B.
int pcol[]
Array of pointers to the beginning of columns in
A and irow. pcol must have n+1 elements and the
last element must be nnz.
int pcolA[]
Same as pcol.
int pcolB[]
Array of pointers to the beginning of columns in
B and irowB. pcol must have n+1 elements and the
last element must be nnzB.
A parameter used only if the problem is symmetric.
uplo indicates whether the lower (uplo = "L") or the
upper triangular (uplo = "U") part of A (and also
B, if the problem is a generalized one) is being
supplied by the user.
char InvertMode
Spectral transformation used when solving symmetric
generalized problems. To use the shift and invert
mode, the user must set InvertMode to 'S'.
and Cayley modes are represented by 'B' and 'C',
respectively.
TYPE sigma
The shift (when shift-and invert mode is being used).
TYPE must be one of "float" or "double" if the
problem is nonsymmetric and sigma is real, or one of
"arcomplex&float&" or "arcomplex&double&" if the
problem is complex.
ARFLOAT sigmaR
Real part of the shift if the problem is real but
the shift is complex. ARFLOAT must be one of "float"
or "double".
ARFLOAT sigmaI
Imaginary part of the shift if the problem is real
but the shift is complex. ARFLOAT must be one of
"float" or "double".
A parameter that characterizes which part (real or
imaginary) of the vector y = OP*x will be used by
ARPACK++ when the problem is real but the shift is
complex. "part" must be set to one of 'R' (real
part) or 'I' (imaginary part).
Number of eigenvalues to be computed.
2) Optional input parameters:
char* which
A parameter thar specifies which of the Ritz values
are to be computed. "which" must be set to one of:
LM: to find eigenvalues wi
SM: to find eigenvalues wit
LR: to find eigenvalues wi
SR: to find eigenvalues wit
LI: to find eigenvalues with la
SI: to find eigenvalues with smallest imaginary part.
Default: LM.
Number of Arnoldi vectors generated at each
iteration of ARPACK. Default: 2*nev+1.
ARFLOAT tol
Stopping criterion (relative accuracy of Ritz
values). ARFLOAT must be one of "float" or "double".
Default: machine precision.
Maximum number of Arnoldi update iterations allowed.
Default: 100*nev.
TYPE* resid
A pointer to an array that contains the initail
vector. Default: a random vector.
TYPE must be one of "float", "double",
"arcomplex&float&" or "arcomplex&double&".
bool AutoShift
A parameter that indicates if exact shifts for
implicit restarting of the Arnoldi method are to
be generated internally by ARPACK++ or shifts are
being supplied by the user. Default: true (exact
shifts are being used).
3) Output parameters:
A vector that contains the "converged" eigenvalues
when the problem is symmetric or complex. TYPE
must be one of float, double, "arcomplex&float&" or
"arcomplex&double&".
ARFLOAT EigValR[] A vector that contains the real part of the
"converged" eigenvalues (when the problem is real).
ARFLOAT must be one of "float" or "double".
ARFLOAT EigValI[] A vector that contains the imaginary part of the
"converged" eigenvalues (when the problem is real).
ARFLOAT must be one of "float" or "double".
A vector that stores all "converged" eigenvectors
consecutively. For real problems, complex eigenvectors
are given as two consecutive vectors. The first
contains the real part of the eigenvector, while
the imaginary part is stored in the second vector.
TYPE must be one of "float", "double",
"arcomplex&float&" or "arcomplex&double&".
ARPACK Authors
Richard Lehoucq
Danny Sorensen
Dept. of Computational & Applied Mathematics
Rice University
Houston, Texas
359 #ifndef AREIG_H
360 #define AREIG_H
363 #include "arcomp.h"
364 #include "arlsmat.h"
365 #include "arlnsmat.h"
366 #include "arlssym.h"
367 #include "arlgsym.h"
368 #include "arlsnsym.h"
369 #include "arlgnsym.h"
370 #include "arlscomp.h"
371 #include "arlgcomp.h"
374 template &class ARFLOAT&
375 int AREig(arcomplex&ARFLOAT& EigVal[], int n, int nnz, arcomplex&ARFLOAT& A[],
int irow[], int pcol[], int nev, char* which = "LM", int ncv = 0,
ARFLOAT tol = 0.0, int maxit = 0, arcomplex&ARFLOAT&* resid = NULL,
bool AutoShift = true)
// Creating a matrix in ARPACK++ format.
ARluNonSymMatrix&arcomplex&ARFLOAT&, ARFLOAT& matrix(n, nnz, A, irow, pcol);
// Defining the eigenvalue problem.
ARluCompStdEig&ARFLOAT& prob(nev, matrix, which, ncv, tol,
maxit, resid, AutoShift);
// Finding eigenvalues.
return prob.Eigenvalues(EigVal);
394 } // complex standard problem, only eigenvalues, regular mode.
397 template &class ARFLOAT&
398 int AREig(arcomplex&ARFLOAT& EigVal[], arcomplex&ARFLOAT& EigVec[], int n,
int nnz, arcomplex&ARFLOAT& A[], int irow[], int pcol[],
int nev, char* which = "LM", int ncv = 0, ARFLOAT tol = 0.0,
int maxit = 0, arcomplex&ARFLOAT&* resid = NULL,
bool AutoShift = true)
// Creating a matrix in ARPACK++ format.
ARluNonSymMatrix&arcomplex&ARFLOAT&, ARFLOAT& matrix(n, nnz, A, irow, pcol);
// Defining the eigenvalue problem.
ARluCompStdEig&ARFLOAT& prob(nev, matrix, which, ncv, tol,
maxit, resid, AutoShift);
// Finding eigenvalues and eigenvectors.
return prob.EigenValVectors(EigVec, EigVal);
418 } // complex standard problem, values and vectors, regular mode.
421 template &class ARFLOAT&
422 int AREig(arcomplex&ARFLOAT& EigVal[], int n, int nnz, arcomplex&ARFLOAT& A[],
int irow[], int pcol[], arcomplex&ARFLOAT& sigma, int nev,
char* which = "LM", int ncv = 0, ARFLOAT tol = 0.0, int maxit = 0,
arcomplex&ARFLOAT&* resid = NULL, bool AutoShift = true)
// Creating a matrix in ARPACK++ format.
ARluNonSymMatrix&arcomplex&ARFLOAT&, ARFLOAT& matrix(n, nnz, A, irow, pcol);
// Defining the eigenvalue problem.
ARluCompStdEig&ARFLOAT& prob(nev, matrix, sigma, which, ncv,
tol, maxit, resid, AutoShift);
// Finding eigenvalues.
return prob.Eigenvalues(EigVal);
441 } // complex standard problem, only eigenvalues, shift-and-invert.
444 template &class ARFLOAT&
445 int AREig(arcomplex&ARFLOAT& EigVal[], arcomplex&ARFLOAT& EigVec[], int n,
int nnz, arcomplex&ARFLOAT& A[], int irow[], int pcol[],
arcomplex&ARFLOAT& sigma, int nev, char* which = "LM",
int ncv = 0, ARFLOAT tol = 0.0, int maxit = 0,
arcomplex&ARFLOAT&* resid = NULL, bool AutoShift = true)
// Creating a matrix in ARPACK++ format.
ARluNonSymMatrix&arcomplex&ARFLOAT&, ARFLOAT& matrix(n, nnz, A, irow, pcol);
// Defining the eigenvalue problem.
ARluCompStdEig&ARFLOAT& prob(nev, matrix, sigma, which, ncv,
tol, maxit, resid, AutoShift);
// Finding eigenvalues and eigenvectors.
return prob.EigenValVectors(EigVec, EigVal);
465 } // complex standard problem, values and vectors, shift-and-invert.
468 template &class ARFLOAT&
469 int AREig(arcomplex&ARFLOAT& EigVal[], int n, int nnzA,
arcomplex&ARFLOAT& A[], int irowA[], int pcolA[], int nnzB,
arcomplex&ARFLOAT& B[], int irowB[], int pcolB[], int nev,
char* which = "LM", int ncv = 0, ARFLOAT tol = 0.0, int maxit = 0,
arcomplex&ARFLOAT&* resid = NULL, bool AutoShift = true)
// Creating two matrices in ARPACK++ format.
ARluNonSymMatrix&arcomplex&ARFLOAT&,ARFLOAT& matrixA(n,nnzA,A,irowA,pcolA);
ARluNonSymMatrix&arcomplex&ARFLOAT&,ARFLOAT& matrixB(n,nnzB,B,irowB,pcolB);
// Defining the eigenvalue problem.
ARluCompGenEig&ARFLOAT& prob(nev, matrixA, matrixB, which,
ncv, tol, maxit, resid, AutoShift);
// Finding eigenvalues.
return prob.Eigenvalues(EigVal);
490 } // complex generalized problem, only eigenvalues, regular mode.
493 template &class ARFLOAT&
494 int AREig(arcomplex&ARFLOAT& EigVal[], arcomplex&ARFLOAT& EigVec[], int n,
int nnzA, arcomplex&ARFLOAT& A[], int irowA[], int pcolA[],
int nnzB, arcomplex&ARFLOAT& B[], int irowB[], int pcolB[],
int nev, char* which = "LM", int ncv = 0, ARFLOAT tol = 0.0,
int maxit = 0, arcomplex&ARFLOAT&* resid = NULL,
bool AutoShift = true)
// Creating two matrices in ARPACK++ format.
ARluNonSymMatrix&arcomplex&ARFLOAT&,ARFLOAT& matrixA(n,nnzA,A,irowA,pcolA);
ARluNonSymMatrix&arcomplex&ARFLOAT&,ARFLOAT& matrixB(n,nnzB,B,irowB,pcolB);
// Defining the eigenvalue problem.
ARluCompGenEig&ARFLOAT& prob(nev, matrixA, matrixB, which,
ncv, tol, maxit, resid, AutoShift);
// Finding eigenvalues and eigenvectors.
return prob.EigenValVectors(EigVec, EigVal);
516 } // complex generalized problem, values and vectors, regular mode.
519 template &class ARFLOAT&
520 int AREig(arcomplex&ARFLOAT& EigVal[], int n, int nnzA, arcomplex&ARFLOAT& A[],
int irowA[], int pcolA[], int nnzB, arcomplex&ARFLOAT& B[],
int irowB[], int pcolB[], arcomplex&ARFLOAT& sigma, int nev,
char* which = "LM", int ncv = 0, ARFLOAT tol = 0.0, int maxit = 0,
arcomplex&ARFLOAT&* resid = NULL, bool AutoShift = true)
// Creating two matrices in ARPACK++ format.
ARluNonSymMatrix&arcomplex&ARFLOAT&,ARFLOAT& matrixA(n,nnzA,A,irowA,pcolA);
ARluNonSymMatrix&arcomplex&ARFLOAT&,ARFLOAT& matrixB(n,nnzB,B,irowB,pcolB);
// Defining the eigenvalue problem.
ARluCompGenEig&ARFLOAT& prob(nev, matrixA, matrixB, sigma, which,
ncv, tol, maxit, resid, AutoShift);
// Finding eigenvalues.
return prob.Eigenvalues(EigVal);
541 } // complex generalized problem, only eigenvalues, shift-and-invert mode.
544 template &class ARFLOAT&
545 int AREig(arcomplex&ARFLOAT& EigVal[], arcomplex&ARFLOAT& EigVec[], int n,
int nnzA, arcomplex&ARFLOAT& A[], int irowA[], int pcolA[],
int nnzB, arcomplex&ARFLOAT& B[], int irowB[], int pcolB[],
arcomplex&ARFLOAT& sigma, int nev, char* which = "LM",
int ncv = 0, ARFLOAT tol = 0.0, int maxit = 0,
arcomplex&ARFLOAT&* resid = NULL, bool AutoShift = true)
// Creating two matrices in ARPACK++ format.
ARluNonSymMatrix&arcomplex&ARFLOAT&,ARFLOAT& matrixA(n,nnzA,A,irowA,pcolA);
ARluNonSymMatrix&arcomplex&ARFLOAT&,ARFLOAT& matrixB(n,nnzB,B,irowB,pcolB);
// Defining the eigenvalue problem.
ARluCompGenEig&ARFLOAT& prob(nev, matrixA, matrixB, sigma, which,
ncv, tol, maxit, resid, AutoShift);
// Finding eigenvalues and eigenvectors.
return prob.EigenValVectors(EigVec, EigVal);
567 } // complex generalized problem, values and vectors, shift-and-invert mode.
570 template &class ARFLOAT&
571 int AREig(double EigValR[], ARFLOAT EigValI[], int n, int nnz,
ARFLOAT A[], int irow[], int pcol[], int nev, char* which = "LM",
int ncv = 0, ARFLOAT tol = 0.0, int maxit = 0,
ARFLOAT* resid = NULL, bool AutoShift = true)
// Creating a matrix in ARPACK++ format.
ARluNonSymMatrix&ARFLOAT, ARFLOAT& matrix(n, nnz, A, irow, pcol);
// Defining the eigenvalue problem.
ARluNonSymStdEig&ARFLOAT& prob(nev, matrix, which, ncv, tol,
maxit, resid, AutoShift);
// Finding eigenvalues.
return prob.Eigenvalues(EigValR, EigValI);
590 } // real nonsymmetric standard problem, only eigenvalues, regular mode.
593 template &class ARFLOAT&
594 int AREig(float EigValR[], ARFLOAT EigValI[], int n, int nnz,
ARFLOAT A[], int irow[], int pcol[], int nev, char* which = "LM",
int ncv = 0, ARFLOAT tol = 0.0, int maxit = 0,
ARFLOAT* resid = NULL, bool AutoShift = true)
// Creating a matrix in ARPACK++ format.
ARluNonSymMatrix&ARFLOAT, ARFLOAT& matrix(n, nnz, A, irow, pcol);
// Defining the eigenvalue problem.
ARluNonSymStdEig&ARFLOAT& prob(nev, matrix, which, ncv, tol,
maxit, resid, AutoShift);
// Finding eigenvalues.
return prob.Eigenvalues(EigValR, EigValI);
613 } // real nonsymmetric standard problem, only eigenvalues, regular mode.
616 template &class ARFLOAT&
617 int AREig(ARFLOAT EigValR[], ARFLOAT EigValI[], ARFLOAT EigVec[], int n,
int nnz, ARFLOAT A[], int irow[], int pcol[], int nev,
char* which = "LM", int ncv = 0, ARFLOAT tol = 0.0,
int maxit = 0, ARFLOAT* resid = NULL, bool AutoShift = true)
// Creating a matrix in ARPACK++ format.
ARluNonSymMatrix&ARFLOAT, ARFLOAT& matrix(n, nnz, A, irow, pcol);
// Defining the eigenvalue problem.
ARluNonSymStdEig&ARFLOAT& prob(nev, matrix, which, ncv, tol,
maxit, resid, AutoShift);
// Finding eigenvalues and eigenvectors.
return prob.EigenValVectors(EigVec, EigValR, EigValI);
636 } // real nonsymmetric standard problem, values and vectors, regular mode.
639 template &class ARFLOAT&
640 int AREig(double EigValR[], ARFLOAT EigValI[], int n, int nnz,
ARFLOAT A[], int irow[], int pcol[], ARFLOAT sigma, int nev,
char* which = "LM", int ncv = 0, ARFLOAT tol = 0.0,
int maxit = 0, ARFLOAT* resid = NULL, bool AutoShift = true)
// Creating a matrix in ARPACK++ format.
ARluNonSymMatrix&ARFLOAT, ARFLOAT& matrix(n, nnz, A, irow, pcol);
// Defining the eigenvalue problem.
ARluNonSymStdEig&ARFLOAT& prob(nev, matrix, sigma, which, ncv,
tol, maxit, resid, AutoShift);
// Finding eigenvalues.
return prob.Eigenvalues(EigValR, EigValI);
659 } // real nonsymmetric standard problem, only eigenvalues, shift-and-invert.
662 template &class ARFLOAT&
663 int AREig(float EigValR[], ARFLOAT EigValI[], int n, int nnz,
ARFLOAT A[], int irow[], int pcol[], ARFLOAT sigma, int nev,
char* which = "LM", int ncv = 0, ARFLOAT tol = 0.0,
int maxit = 0, ARFLOAT* resid = NULL, bool AutoShift = true)
// Creating a matrix in ARPACK++ format.
ARluNonSymMatrix&ARFLOAT, ARFLOAT& matrix(n, nnz, A, irow, pcol);
// Defining the eigenvalue problem.
ARluNonSymStdEig&ARFLOAT& prob(nev, matrix, sigma, which, ncv,
tol, maxit, resid, AutoShift);
// Finding eigenvalues.
return prob.Eigenvalues(EigValR, EigValI);
682 } // real nonsymmetric standard problem, only eigenvalues, shift-and-invert.
685 template &class ARFLOAT&
686 int AREig(ARFLOAT EigValR[], ARFLOAT EigValI[], ARFLOAT EigVec[], int n,
int nnz, ARFLOAT A[], int irow[], int pcol[], ARFLOAT sigma,
int nev, char* which = "LM", int ncv = 0, ARFLOAT tol = 0.0,
int maxit = 0, ARFLOAT* resid = NULL, bool AutoShift = true)
// Creating a matrix in ARPACK++ format.
ARluNonSymMatrix&ARFLOAT, ARFLOAT& matrix(n, nnz, A, irow, pcol);
// Defining the eigenvalue problem.
ARluNonSymStdEig&ARFLOAT& prob(nev, matrix, sigma, which, ncv,
tol, maxit, resid, AutoShift);
// Finding eigenvalues and eigenvectors.
return prob.EigenValVectors(EigVec, EigValR, EigValI);
705 } // real nonsymmetric standard problem, values and vectors, shift-and-invert.
708 template &class ARFLOAT&
709 int AREig(double EigValR[], ARFLOAT EigValI[], int n, int nnzA,
ARFLOAT A[], int irowA[], int pcolA[], int nnzB,
ARFLOAT B[], int irowB[], int pcolB[], int nev,
char* which = "LM", int ncv = 0, ARFLOAT tol = 0.0,
int maxit = 0, ARFLOAT* resid = NULL, bool AutoShift = true)
// Creating two matrices in ARPACK++ format.
ARluNonSymMatrix&ARFLOAT, ARFLOAT& matrixA(n, nnzA, A, irowA, pcolA);
ARluNonSymMatrix&ARFLOAT, ARFLOAT& matrixB(n, nnzB, B, irowB, pcolB);
// Defining the eigenvalue problem.
ARluNonSymGenEig&ARFLOAT& prob(nev, matrixA, matrixB, which,
ncv, tol, maxit, resid, AutoShift);
// Finding eigenvalues.
return prob.Eigenvalues(EigValR, EigValI);
730 } // real nonsymmetric generalized problem, only eigenvalues, regular mode.
733 template &class ARFLOAT&
734 int AREig(float EigValR[], ARFLOAT EigValI[], int n, int nnzA,
ARFLOAT A[], int irowA[], int pcolA[], int nnzB,
ARFLOAT B[], int irowB[], int pcolB[], int nev,
char* which = "LM", int ncv = 0, ARFLOAT tol = 0.0,
int maxit = 0, ARFLOAT* resid = NULL, bool AutoShift = true)
// Creating two matrices in ARPACK++ format.
ARluNonSymMatrix&ARFLOAT, ARFLOAT& matrixA(n, nnzA, A, irowA, pcolA);
ARluNonSymMatrix&ARFLOAT, ARFLOAT& matrixB(n, nnzB, B, irowB, pcolB);
// Defining the eigenvalue problem.
ARluNonSymGenEig&ARFLOAT& prob(nev, matrixA, matrixB, which,
ncv, tol, maxit, resid, AutoShift);
// Finding eigenvalues.
return prob.Eigenvalues(EigValR, EigValI);
755 } // real nonsymmetric generalized problem, only eigenvalues, regular mode.
758 template &class ARFLOAT&
759 int AREig(ARFLOAT EigValR[], ARFLOAT EigValI[], ARFLOAT EigVec[], int n,
int nnzA, ARFLOAT A[], int irowA[], int pcolA[],
int nnzB, ARFLOAT B[], int irowB[], int pcolB[],
int nev, char* which = "LM", int ncv = 0, ARFLOAT tol = 0.0,
int maxit = 0, ARFLOAT* resid = NULL, bool AutoShift = true)
// Creating two matrices in ARPACK++ format.
ARluNonSymMatrix&ARFLOAT, ARFLOAT& matrixA(n, nnzA, A, irowA, pcolA);
ARluNonSymMatrix&ARFLOAT, ARFLOAT& matrixB(n, nnzB, B, irowB, pcolB);
// Defining the eigenvalue problem.
ARluNonSymGenEig&ARFLOAT& prob(nev, matrixA, matrixB, which,
ncv, tol, maxit, resid, AutoShift);
// Finding eigenvalues and eigenvectors.
return prob.EigenValVectors(EigVec, EigValR, EigValI);
780 } // real nonsymmetric generalized problem, values and vectors, regular mode.
783 template &class ARFLOAT&
784 int AREig(double EigValR[], ARFLOAT EigValI[], int n, int nnzA,
ARFLOAT A[], int irowA[], int pcolA[], int nnzB,
ARFLOAT B[], int irowB[], int pcolB[], ARFLOAT sigma,
int nev, char* which = "LM", int ncv = 0, ARFLOAT tol = 0.0,
int maxit = 0, ARFLOAT* resid = NULL, bool AutoShift = true)
// Creating two matrices in ARPACK++ format.
ARluNonSymMatrix&ARFLOAT, ARFLOAT& matrixA(n, nnzA, A, irowA, pcolA);
ARluNonSymMatrix&ARFLOAT, ARFLOAT& matrixB(n, nnzB, B, irowB, pcolB);
// Defining the eigenvalue problem.
ARluNonSymGenEig&ARFLOAT& prob(nev, matrixA, matrixB, sigma, which,
ncv, tol, maxit, resid, AutoShift);
// Finding eigenvalues.
return prob.Eigenvalues(EigValR, EigValI);
805 } // real nonsymmetric generalized problem, only eigenvalues,
// real shift-and-invert mode.
809 template &class ARFLOAT&
810 int AREig(float EigValR[], ARFLOAT EigValI[], int n, int nnzA,
ARFLOAT A[], int irowA[], int pcolA[], int nnzB,
ARFLOAT B[], int irowB[], int pcolB[], ARFLOAT sigma,
int nev, char* which = "LM", int ncv = 0, ARFLOAT tol = 0.0,
int maxit = 0, ARFLOAT* resid = NULL, bool AutoShift = true)
// Creating two matrices in ARPACK++ format.
ARluNonSymMatrix&ARFLOAT, ARFLOAT& matrixA(n, nnzA, A, irowA, pcolA);
ARluNonSymMatrix&ARFLOAT, ARFLOAT& matrixB(n, nnzB, B, irowB, pcolB);
// Defining the eigenvalue problem.
ARluNonSymGenEig&ARFLOAT& prob(nev, matrixA, matrixB, sigma, which,
ncv, tol, maxit, resid, AutoShift);
// Finding eigenvalues.
return prob.Eigenvalues(EigValR, EigValI);
831 } // real nonsymmetric generalized problem, only eigenvalues,
// real shift-and-invert mode.
835 template &class ARFLOAT&
836 int AREig(ARFLOAT EigValR[], ARFLOAT EigValI[], ARFLOAT EigVec[], int n,
int nnzA, ARFLOAT A[], int irowA[], int pcolA[], int nnzB,
ARFLOAT B[], int irowB[], int pcolB[], ARFLOAT sigma, int nev,
char* which = "LM", int ncv = 0, ARFLOAT tol = 0.0,
int maxit = 0, ARFLOAT* resid = NULL, bool AutoShift = true)
// Creating two matrices in ARPACK++ format.
ARluNonSymMatrix&ARFLOAT, ARFLOAT& matrixA(n, nnzA, A, irowA, pcolA);
ARluNonSymMatrix&ARFLOAT, ARFLOAT& matrixB(n, nnzB, B, irowB, pcolB);
// Defining the eigenvalue problem.
ARluNonSymGenEig&ARFLOAT& prob(nev, matrixA, matrixB, sigma, which,
ncv, tol, maxit, resid, AutoShift);
// Finding eigenvalues and eigenvectors.
return prob.EigenValVectors(EigVec, EigValR, EigValI);
857 } // real nonsymmetric generalized problem, values and vectors,
// real shift-and-invert mode.
861 template &class ARFLOAT&
862 int AREig(ARFLOAT EigValR[], ARFLOAT EigValI[], int n, int nnzA, ARFLOAT A[],
int irowA[], int pcolA[], int nnzB, ARFLOAT B[], int irowB[],
int pcolB[], char part, ARFLOAT sigmaR, ARFLOAT sigmaI,
int nev, char* which = "LM", int ncv = 0, ARFLOAT tol = 0.0,
int maxit = 0, ARFLOAT* resid = NULL, bool AutoShift = true)
// Creating two matrices in ARPACK++ format.
ARluNonSymMatrix&ARFLOAT, ARFLOAT& matrixA(n, nnzA, A, irowA, pcolA);
ARluNonSymMatrix&ARFLOAT, ARFLOAT& matrixB(n, nnzB, B, irowB, pcolB);
// Defining the eigenvalue problem.
ARluNonSymGenEig&ARFLOAT& prob(nev, matrixA, matrixB, part,
sigmaR, sigmaI, which, ncv,
tol, maxit, resid, AutoShift);
// Finding eigenvalues.
return prob.Eigenvalues(EigValR, EigValI);
884 } // real nonsymmetric generalized problem, only eigenvalues,
// complex shift-and-invert mode.
888 template &class ARFLOAT&
889 int AREig(ARFLOAT EigValR[], ARFLOAT EigValI[], ARFLOAT EigVec[], int n,
int nnzA, ARFLOAT A[], int irowA[], int pcolA[], int nnzB,
ARFLOAT B[], int irowB[], int pcolB[], char part, ARFLOAT sigmaR,
ARFLOAT sigmaI, int nev, char* which = "LM", int ncv = 0,
ARFLOAT tol = 0.0, int maxit = 0, ARFLOAT* resid = NULL,
bool AutoShift = true)
// Creating two matrices in ARPACK++ format.
ARluNonSymMatrix&ARFLOAT, ARFLOAT& matrixA(n, nnzA, A, irowA, pcolA);
ARluNonSymMatrix&ARFLOAT, ARFLOAT& matrixB(n, nnzB, B, irowB, pcolB);
// Defining the eigenvalue problem.
ARluNonSymGenEig&ARFLOAT& prob(nev, matrixA, matrixB, part,
sigmaR, sigmaI, which, ncv,
tol, maxit, resid, AutoShift);
// Finding eigenvalues and eigenvectors.
return prob.EigenValVectors(EigVec, EigValR, EigValI);
912 } // real nonsymmetric generalized problem, values and vectors,
// complex shift-and-invert mode.
916 template &class ARFLOAT&
917 int AREig(ARFLOAT EigVal[], int n, int nnz, ARFLOAT A[], int irow[],
int pcol[], char uplo, int nev, char* which = "LM",
int ncv = 0, ARFLOAT tol = 0.0, int maxit = 0,
ARFLOAT* resid = NULL, bool AutoShift = true)
// Creating a matrix in ARPACK++ format.
ARluSymMatrix&ARFLOAT& matrix(n, nnz, A, irow, pcol, uplo);
// Defining the eigenvalue problem.
ARluSymStdEig&ARFLOAT& prob(nev, matrix, which, ncv, tol,
maxit, resid, AutoShift);
// Finding eigenvalues.
return prob.Eigenvalues(EigVal);
936 } // real symmetric standard problem, only eigenvalues, regular mode.
939 template &class ARFLOAT&
940 int AREig(ARFLOAT EigVal[], ARFLOAT EigVec[], int n, int nnz, ARFLOAT A[],
int irow[], int pcol[], char uplo, int nev,
char* which = "LM", int ncv = 0, ARFLOAT tol = 0.0,
int maxit = 0, ARFLOAT* resid = NULL, bool AutoShift = true)
// Creating a matrix in ARPACK++ format.
ARluSymMatrix&ARFLOAT& matrix(n, nnz, A, irow, pcol, uplo);
// Defining the eigenvalue problem.
ARluSymStdEig&ARFLOAT& prob(nev, matrix, which, ncv, tol,
maxit, resid, AutoShift);
// Finding eigenvalues and eigenvectors.
return prob.EigenValVectors(EigVec, EigVal);
959 } // real symmetric standard problem, values and vectors, regular mode.
962 template &class ARFLOAT&
963 int AREig(ARFLOAT EigVal[], int n, int nnz, ARFLOAT A[], int irow[],
int pcol[], char uplo, ARFLOAT sigma, int nev,
char* which = "LM", int ncv = 0, ARFLOAT tol = 0.0,
int maxit = 0, ARFLOAT* resid = NULL, bool AutoShift = true)
// Creating a matrix in ARPACK++ format.
ARluSymMatrix&ARFLOAT& matrix(n, nnz, A, irow, pcol, uplo);
// Defining the eigenvalue problem.
ARluSymStdEig&ARFLOAT& prob(nev, matrix, sigma, which, ncv,
tol, maxit, resid, AutoShift);
// Finding eigenvalues.
return prob.Eigenvalues(EigVal);
982 } // real symmetric standard problem, only eigenvalues, shift-and-invert.
985 template &class ARFLOAT&
986 int AREig(ARFLOAT EigVal[], ARFLOAT EigVec[], int n, int nnz, ARFLOAT A[],
int irow[], int pcol[], char uplo, ARFLOAT sigma,
int nev, char* which = "LM", int ncv = 0, ARFLOAT tol = 0.0,
int maxit = 0, ARFLOAT* resid = NULL, bool AutoShift = true)
// Creating a matrix in ARPACK++ format.
ARluSymMatrix&ARFLOAT& matrix(n, nnz, A, irow, pcol, uplo);
// Defining the eigenvalue problem.
ARluSymStdEig&ARFLOAT& prob(nev, matrix, sigma, which, ncv,
tol, maxit, resid, AutoShift);
// Finding eigenvalues and eigenvectors.
return prob.EigenValVectors(EigVec, EigVal);
1005 } // real symmetric standard problem, values and vectors, shift-and-invert.
1008 template &class ARFLOAT&
1009 int AREig(ARFLOAT EigVal[], int n, int nnzA, ARFLOAT A[], int irowA[],
int pcolA[], int nnzB, ARFLOAT B[], int irowB[], int pcolB[],
char uplo, int nev, char* which = "LM", int ncv = 0,
ARFLOAT tol = 0.0, int maxit = 0, ARFLOAT* resid = NULL,
bool AutoShift = true)
// Creating two matrices in ARPACK++ format.
ARluSymMatrix&ARFLOAT& matrixA(n, nnzA, A, irowA, pcolA, uplo);
ARluSymMatrix&ARFLOAT& matrixB(n, nnzB, B, irowB, pcolB, uplo);
// Defining the eigenvalue problem.
ARluSymGenEig&ARFLOAT& prob(nev, matrixA, matrixB, which,
ncv, tol, maxit, resid, AutoShift);
// Finding eigenvalues.
return prob.Eigenvalues(EigVal);
1030 } // real symmetric generalized problem, only eigenvalues, regular mode.
1033 template &class ARFLOAT&
1034 int AREig(ARFLOAT EigVal[], ARFLOAT EigVec[], int n, int nnzA, ARFLOAT A[],
int irowA[], int pcolA[], int nnzB, ARFLOAT B[], int irowB[],
int pcolB[], char uplo, int nev, char* which = "LM",
int ncv = 0, ARFLOAT tol = 0.0, int maxit = 0,
ARFLOAT* resid = NULL, bool AutoShift = true)
// Creating two matrices in ARPACK++ format.
ARluSymMatrix&ARFLOAT& matrixA(n, nnzA, A, irowA, pcolA, uplo);
ARluSymMatrix&ARFLOAT& matrixB(n, nnzB, B, irowB, pcolB, uplo);
// Defining the eigenvalue problem.
ARluSymGenEig&ARFLOAT& prob(nev, matrixA, matrixB, which,
ncv, tol, maxit, resid, AutoShift);
// Finding eigenvalues and eigenvectors.
return prob.EigenValVectors(EigVec, EigVal);
1055 } // real symmetric generalized problem, values and vectors, regular mode.
1058 template &class ARFLOAT&
1059 int AREig(ARFLOAT EigVal[], int n, int nnzA, ARFLOAT A[], int irowA[],
int pcolA[], int nnzB, ARFLOAT B[], int irowB[], int pcolB[],
char uplo, char InvertMode, ARFLOAT sigma, int nev,
char* which = "LM", int ncv = 0, ARFLOAT tol = 0.0,
int maxit = 0, ARFLOAT* resid = NULL, bool AutoShift = true)
// Creating two matrices in ARPACK++ format.
ARluSymMatrix&ARFLOAT& matrixA(n, nnzA, A, irowA, pcolA, uplo);
ARluSymMatrix&ARFLOAT& matrixB(n, nnzB, B, irowB, pcolB, uplo);
// Defining the eigenvalue problem.
ARluSymGenEig&ARFLOAT& prob(InvertMode, nev, matrixA, matrixB, sigma,
which, ncv, tol, maxit, resid, AutoShift);
// Finding eigenvalues.
return prob.Eigenvalues(EigVal);
1080 } // real symmetric generalized problem, only eigenvalues,
// shift-and-invert, buckling and Cayley modes.
1084 template &class ARFLOAT&
1085 int AREig(ARFLOAT EigVal[], ARFLOAT EigVec[], int n, int nnzA, ARFLOAT A[],
int irowA[], int pcolA[], int nnzB, ARFLOAT B[], int irowB[],
int pcolB[], char uplo, char InvertMode, ARFLOAT sigma,
int nev, char* which = "LM", int ncv = 0, ARFLOAT tol = 0.0,
int maxit = 0, ARFLOAT* resid = NULL, bool AutoShift = true)
// Creating two matrices in ARPACK++ format.
ARluSymMatrix&ARFLOAT& matrixA(n, nnzA, A, irowA, pcolA, uplo);
ARluSymMatrix&ARFLOAT& matrixB(n, nnzB, B, irowB, pcolB, uplo);
// Defining the eigenvalue problem.
ARluSymGenEig&ARFLOAT& prob(InvertMode, nev, matrixA, matrixB, sigma,
which, ncv, tol, maxit, resid, AutoShift);
// Finding eigenvalues and eigenvectors.
return prob.EigenValVectors(EigVec, EigVal);
1106 } // real symmetric generalized problem, values and vectors,
// shift-and-invert, buckling and Cayley modes.
1110 #endif // AREIG_H
ARPACK++ v1.2 2/18/2000
c++ interface to ARPACK code.
MODULE ASymSol.cc
Template functions that exemplify how to print information
about symmetric standard and generalized eigenvalue problems.
ARPACK Authors
Richard Lehoucq
Danny Sorensen
Dept. of Computational & Applied Mathematics
Rice University
Houston, Texas
18 #ifndef ASYMSOL_H
19 #define ASYMSOL_H
21 #include &cmath&
22 #include "blas1c.h"
23 #include "lapackc.h"
24 #include "arlsmat.h"
26 template&class ARFLOAT, class ARINT&
27 void Solution(ARINT nconv, ARINT n, ARINT nnz, ARFLOAT A[], ARINT irow[],
ARINT pcol[], char uplo, ARFLOAT EigVal[], ARFLOAT* EigVec = 0)
Prints eigenvalues and eigenvectors of symmetric eigen-problems
on standard "std::cout" stream.
ARluSymMatrix&ARFLOAT& matrix(n, nnz, A, irow, pcol, uplo);
std::cout && std::endl && std::endl && "Testing ARPACK++ function AREig" && std::
std::cout && "Real symmetric eigenvalue problem: A*x - lambda*x \n \n";
std::cout && "Dimension of the system
std::cout && "Number of 'converged' eigenvalues
: " && nconv && std::endl && std::
// Printing eigenvalues.
std::cout && "Eigenvalues:" && std::
for (i=0; i& i++) {
std::cout && "
lambda[" && (i+1) && "]: " && EigVal[i] && std::
std::cout && std::
// Printing eigenvectors.
if (EigVec != 0) {
// Finding the residual norm || A*x - lambda*x ||
// for the nconv accurately computed eigenvectors.
= new ARFLOAT[n];
ResNorm = new ARFLOAT[nconv+1];
for (i=0; i& i++) {
matrix.MultMv(&EigVec[i*n], Ax);
axpy(n, -EigVal[i], &EigVec[i*n], 1, Ax, 1);
ResNorm[i] = nrm2(n, Ax, 1)/fabs(EigVal[i]);
for (i=0; i& i++) {
std::cout && "||A*x(" && (i+1) && ") - lambda(" && (i+1);
std::cout && ")*x(" && (i+1) && ")||: " && ResNorm[i] && std::
std::cout && std::
delete[] Ax;
delete[] ResN
83 } // Solution.
86 template&class ARFLOAT, class ARINT&
87 void Solution(ARINT nconv, ARINT n, ARINT nnzA, ARFLOAT A[], ARINT irowA[],
ARINT pcolA[], ARINT nnzB, ARFLOAT B[], ARINT irowB[],
ARINT pcolB[], char uplo, ARFLOAT EigVal[], ARFLOAT* EigVec = 0)
Prints eigenvalues and eigenvectors of symmetric generalized
eigen-problem on standard "std::cout" stream.
ARluSymMatrix&ARFLOAT& matrixA(n, nnzA, A, irowA, pcolA, uplo);
ARluSymMatrix&ARFLOAT& matrixB(n, nnzB, B, irowB, pcolB, uplo);
std::cout && std::endl &&std::endl && "Testing ARPACK++ function AREig" &&std::
std::cout && "Real symmetric generalized eigenvalue problem: A*x - lambda*B*x";
std::cout && std::endl &&std::
std::cout && "Dimension of the system
std::cout && "Number of 'converged' eigenvalues
: " && nconv && std::endl &&std::
// Printing eigenvalues.
std::cout && "Eigenvalues:" && std::
for (i=0; i& i++) {
std::cout && "
lambda[" && (i+1) && "]: " && EigVal[i] && std::
std::cout && std::
// Printing eigenvectors.
if (EigVec != 0) {
// Printing the residual norm || A*x - lambda*B*x ||
// for the nconv accurately computed eigenvectors.
= new ARFLOAT[n];
= new ARFLOAT[n];
ResNorm = new ARFLOAT[nconv+1];
for (i=0; i& i++) {
matrixA.MultMv(&EigVec[i*n], Ax);
matrixB.MultMv(&EigVec[i*n], Bx);
axpy(n, -EigVal[i], Bx, 1, Ax, 1);
ResNorm[i] = nrm2(n, Ax, 1)/fabs(EigVal[i]);
for (i=0; i& i++) {
std::cout && "||A*x(" && i && ") - lambda(" &&
std::cout && ")*B*x(" && i && ")||: " && ResNorm[i] && std::
std::cout && std::
delete[] Ax;
delete[] Bx;
delete[] ResN
149 } // Solution.
152 #endif // ASYMSOL_H
lsmatrxb.h
ARPACK++ v1.2 2/20/2000
c++ interface to ARPACK code.
MODULE LSMatrxB.h
Function template for the one dimensional discrete Laplacian
on the interval [0, 1], with zero Dirichlet boundary conditions.
ARPACK Authors
Richard Lehoucq
Danny Sorensen
Dept. of Computational & Applied Mathematics
Rice University
Houston, Texas
18 #ifndef LSMATRXB_H
19 #define LSMATRXB_H
21 #include &cmath&
23 template&class ARFLOAT, class ARINT&
24 void SymmetricMatrixB(ARINT n, ARINT& nnz, ARFLOAT* &A,
ARINT* &irow, ARINT* &pcol, char uplo = 'L')
// Defining internal variables.
// Defining constants.
= (ARFLOAT(n+1)*ARFLOAT(n+1));
dd = 2.0*h2;
// Defining the number of nonzero elements in A.
nnz = 2*n-1;
// Creating output vectors.
= new ARFLOAT[nnz];
irow = new ARINT[nnz];
pcol = new ARINT[n+1];
// Defining matrix A.
pcol[0] = 0;
if (uplo == 'U') {
for (j = 0; j & j++) {
irow[i++] = j-1;
irow[i++] =
pcol[j+1] =
for (j = 0; j & j++) {
irow[i++] =
if (n-j-1) {
irow[i++] = j+1;
pcol[j+1] =
77 } // SymmetricMatrixB.
79 #endif // LSMATRXB_H
asymshf.cpp
ARPACK++ v1.2 2/18/2000
c++ interface to ARPACK code.
MODULE ASymShf.cc.
Example program that illustrates how to solve a real symmetric
standard eigenvalue problem in shift and invert mode using the
AREig function.
1) Problem description:
In this example we try to solve A*x = x*lambda in shift and
invert mode, where A is derived from the central difference
discretization of the one-dimensional Laplacian on [0, 1]
with zero Dirichlet boundary conditions.
2) Data structure used to represent matrix A:
{nnz, irow, pcol, A}: upper triangular part of matrix A
stored in CSC format.
3) Library called by this example:
The SuperLU package is called by AREig to solve some linear
systems involving (A-sigma*I). This is needed to implement
the shift and invert strategy.
4) Included header files:
-----------
--------------------------------------------
lsmatrxb.h
SymmetricMatrixB, a function that generates
matrix A in CSC format.
The AREig function definition.
The Solution function.
5) ARPACK Authors:
Richard Lehoucq
Kristyn Maschhoff
Danny Sorensen
Dept. of Computational & Applied Mathematics
Rice University
Houston, Texas
48 #include "lsmatrxb.h"
49 #include "areig.h"
50 #include "asymsol.h"
53 int main()
// Dimension of the problem.
// Number of "converged" eigenvalues.
// Number of nonzero elements in A.
// pointer to an array that stores the row
// indices of the nonzeros in A.
// pointer to an array of pointers to the
// beginning of each column of A in vector A.
double* A;
// pointer to an array that stores the
// nonzero elements of A.
double EigVal[101];
// Eigenvalues.
double EigVec[1001]; // Eigenvectors stored sequentially.
// Variable that indicates whether the upper
// (uplo='U') ot the lower (uplo='L') part of
// A will be stored in A, irow and pcol.
// Creating a 100x100 matrix.
uplo = 'U';
SymmetricMatrixB(n, nnz, A, irow, pcol, uplo);
// Finding the four eigenvalues of A nearest to 1.0 and the
// related eigenvectors.
nconv = AREig(EigVal, EigVec, n, nnz, A, irow, pcol, uplo, 1.0, 4);
// Printing solution.
Solution(nconv, n, nnz, A, irow, pcol, uplo, EigVal, EigVec);
88 } // main
结果非零，误差极小。收工~~
这里有个给了一组各个数值求解库的lib/dll编译版，目测是神器级别的（i.e.&妈妈再也不用担心我Windows下不会XX了&），但没有亲自测试过。
以后有时间，我索性用DBank把dll和lib都上传上来好了。嗯
感觉Opensource的东东都是这个样子，毫不考虑windows平台（应该说是没放在眼里），代码的维护也时有时无，有时候真觉得怒其不争。难怪intel mkl啊matlab啊，内核都是blas/lapack/arpack，但确实做得方便省事，所以人没有竞争对手，能卖上百上千刀一个copy。写这个文档也是希望在各位配置LAPACK/ARPACK时有所帮助，不要在这些蛋疼的问题上浪费时间。}