和是一个求解大规模稠密/稀疏矩阵问题的库，最近在做特征值问题时用到。
：下载左边链接的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：
clapack.vcxproj打开工程项目文件
解压CLAPACK-3.1.1-VisualStudio，它是自带某版本的Visual
Studio工程的。
注意，/LIB/Win32下已经有一些lib了。大家最好把他们都先删除了，以免新旧文件混淆。
然我们根据个人需要，编译成debug模式，或者release模式。
为了能在自己的程序调用中方便的进行debug，我这里就以debug模式为例说明。
首先编译F2CLIBS，用于将fortran转换为c语言。选择libf2c子项目，直接生成之。编译过程中可能会有一些warning，不要理会他们。编译成功后，输出文件是libf2cd.lib。
接着编译tmglib。这是一个矩阵库。这个库在TESTINGMATGEN里面，选择他生成就好了。输出文件还是在LIB里面。文件名是tmglibd.lib。
然后是编译BLAS，选择项目blas, 编译之。输出文件BLASd.lib。
最后是编译CLAPACK，生成clapackd.lib。
使用LAPACK：
项目的属性 --& C/C++选项卡 --&常规
--&&附加包含目录 --& 添加INCLUDE目录。
项目的属性
--& 链接器 --&常规 --&&附加包含目录
--& 添加LIB目录/Win32目录。
项目属性&--&&链接器&--&&输入&--&&附加依赖项&--&&添加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&);
编译结果为：
s[ 0 ] = 34.000000
s[ 1 ] = 17.888544
s[ 2 ] = 4.472136
s[ 3 ] = 0.000000
请按任意键继续. . .
3.安装MigGW:
可参考网页：
安装:执行MinGW安装程序双击安装。Vista/Win7用户请右键用管理员身份安装。请尽量装在默认的
C:\MinGW 下。记得至少得勾上C++，gfortran编译器和MSYS。
设置“环境”变量：
右击 “我的电脑”，选择“属性”，在弹出的对话框中选择“高级”选型卡，在选择“环境变量”，在弹出的对话框中可以设置环境变量，当中分系统变量和用户变量，系统变量就是所有用户可用，用户变量就是所用用户可用，对于个人电脑一般建议设置系统变量，这样不管用哪个用户名登录都可以用。
变量：PATH
&值：C:\MinGW\bin
变量：LIBRARY_PATH
&值：C:\MinGW\C:\MinGW\lib\gcc\mingw32\4.6.2 好像可以不用写
变量：C_INCLUDE_PATH
&值：C:\MinGW\C:\MinGW\lib\gcc\mingw32\4.6.2\include&好像可以不用写
变量:&CPLUS_INCLUDE_PATH
&好像可以不用写
值：C:\MinGW\C:\MinGW\lib\gcc\mingw32\4.6.2\include\c++;D:\MinGW\lib\gcc&
\mingw32\4.6.2\include\c++\C:\MinGW\lib\gcc\mingw32\4.6.2\inclu de\c++\mingw32
注：Mingw Shell.exe在C:\ProgramData\Microsoft\Windows\开始菜单\Programs\Mingw下
默认的home目录在 C:\MinGW\msys\1.0\home 或 C:\MinGW\msys\1.0\home\用户名下
4. 解压并改动ARPACK源码
先把arpack96.tar.gz解压，我的解压后的目录是D:\ProgrammFile\ARPACK。
将arpack96_patch.tar.gz进行解压，将该文件夹中的文件复制到D:\ProgrammFile\ARPACK，同名文件进行覆盖。
然后纠结的来了，这里有三个文件要修改，不改动会导致现代gfortran编译不过。
第一个是ARPACK\UTIL\second.f，里面的ETIME编译时有问题，在网上找的资料都说将整个子过程替换成：
SUBROUTINE
SECOND( T )
REAL & & &T & &&
CPU_TIME(T)& & &
但是在我的电脑上编译的时候却说有问题，经过尝试改成以下形式：
SUBROUTINE SECOND( T )
-- LAPACK auxiliary routine (preliminary version) --
Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
Courant Institute, Argonne National Lab, and Rice University
July 26, 1991
SECOND returns the user time for a process in seconds.
This version gets the time from the system function ETIME.
.. Local Scalars ..
CALL CPU_TIME(T)
.. Local Arrays ..
.. External Functions ..
.. Executable Statements ..
End of SECOND
第二个是ARPACK\BLAS\dnrm2.f，里面的ASSIGN遍不过。不过我们可以在netlib上找到BLAS里最新版本的。下载并直接覆盖之。
第三个是ARPACK\BLAS\snrm2.f，问题同上，找到最新版本的覆盖。
改动ARPACK\ARmake.inc里面需要改动的有3处：
home = D:/ProgrammFile/ARPACK （改成自己的ARPACK目录）
FC = gfortran
FFLAGS = -O2
= SUN4就可以，不用改成Win32
6. 编译&打包&VS2010配置
打开MinGW Shell，输入
$ cd ProgrammFile/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
可能有警告提示：
dllwrap.exe: no export definition file provided.
Creating one, but that may not be what you want
不用管它，关闭当前MSYS shell，打开vs2010 --& 工具 --&
Visual Studio命令提示 --&&弹出一窗口，敲入命令
d:\vs2010\vc\bin&d:
d:\vs2010\vc\bin&cd d:\ProgrammFile\ARPACK
d:\ProgrammFile\ARPACK&lib /machine:i386 /def:arpack_win32.def （注意空格，否则出错）
会发现在ARPACK目录下会新生成一些文件，文件中内容变为：
其中最重要的是新生成的arpack_win32.dll和arpack_win32.lib。
在vs2010中配置ARPACK:
项目属性 --& 配置属性 --& 链接器 --& 常规 --& 附加库目录 --& 添加lib文件夹目录D:\ProgrammFile\ARPACK；
项目属性 --& 配置属性 --& 链接器 --&&输入 --& 附加依赖项 --& 添加arpack_win32.lib；
项目属性 --& 配置属性 --&&C/C++ --& 代码生成 --& 运行库选成/MTd（或/MT）；
在自己新建立的工程中找到Debug文件夹，我的工程建在D盘，起名为ar，如图：
测试程序为：
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
#include &math.h&
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);
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);
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 */
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];
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
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];
}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)
September 1999
// Begin with some standard include files.
#include &math.h&
#include &stdio.h&
#include &iostream&
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];
在所建项目ar下找到Debug文件，在其中放入以下3个dll文件：ARPACK文件夹中的arpack_win32.dll，以及MinGW\Bin文件夹中的libgfortran-3.dll和libquadmath-0.dll。这时候就可以用了，下面给出测试例子：
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
#include &math.h&
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);
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);
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 */
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];
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
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];
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)
September 1999
// Begin with some standard include files.
#include &math.h&
#include &stdio.h&
#include &iostream&
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];
编译结果为：
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)
请按任意键继续. . .
版权声明：本文为博主原创文章，未经博主允许不得转载。
* 以上用户言论只代表其个人观点，不代表CSDN网站的观点或立场
访问：1302次
排名：千里之外